I want to run MC simulation of AR1(auto-regression) matrix in R, which would be
as residuals in linear mixed model.
AR1 matrix with the following character:
��r,��c is the auto-correlation parameter in the row and column direction.
And I wrote one function in R ( under following). But I run the function, it
seems only work for the first Row auto-corr. When setting different a set of
Row auto-corr values, the simulated dataset would change with the same value.
But it did not work for the second column auto-corr parameter, even if setting
different col atuo-corr, the simulated dataset seemd no changed in col
auto-corr value that nearly is zero all the time. Would someone please help me
to find the questions that the R function codes somewhere got wrong? Thanks a
lots.
####### simulation codes for AR1 model
multi_norm <- function(data_num,Pr,Pc) {
require(MASS)
# data_num for row/col number; Pr for row auto-corr; Pc for colum auto-corr.
V <- matrix(data=NA, nrow=data_num, ncol=data_num)
R.mat=diag(data_num)
C.mat=diag(data_num)
set.seed(2016)
means <- runif(1, min=0, max=1)
means1=rep(means,data_num*data_num)
# variance
set.seed(2016)
var <- runif(1, min=0, max=1)
for (i in 1:data_num) {
# a two-level nested loop to generate AR matrix
for (j in 1:data_num) {
if (i == j) {
# covariances on the diagonal
V[i,j] <- 1 #varsmodule[i]
} else if(i<j){
# covariances
R.mat[i,j]<- V[i,i]*(Pr^(j-i))
C.mat[i,j]<- V[i,i]*(Pc^(j-i))
}else {R.mat[i,j]=R.mat[j,i];C.mat[i,j]=C.mat[j,i]}
}
}
V=var*kronecker(C.mat,R.mat)
# simulate multivariate normal distribution
# given means and covariance matrix
X <- t(mvrnorm(n = data_num, means1, V))
aam=X[1:data_num,]
aad=data.frame()
for(i in 1:data_num){
for(j in 1:data_num){
aad[j+data_num*(i-1),1]=i
aad[j+data_num*(i-1),2]=j
aad[j+data_num*(i-1),3]=aam[i,j]
}
}
names(aad)=c('Row','Col','y')
for(i in 1:2) aad[,i]=factor(aad[,i])
return(aad)
}
The simulation results as following:
> aam=multi_norm(30,0.6,0.01)
> mm2=asreml(y~1,rcov=~ar1(Row):ar1(Col),data=aam,trace=F,maxit=30)
> summary(mm2)$varcomp
gamma component std.error z.ratio constraint
R!variance 1.00000000 0.16267855 0.01038210 15.6691353 Positive
R!Row.cor 0.55811722 0.55811722 0.02734902 20.4072085 Unconstrained
R!Col.cor 0.01735573 0.01735573 0.03368048 0.5153055 Unconstrained
> aam=multi_norm(30,0.6,0.3)
> mm2=asreml(y~1,rcov=~ar1(Row):ar1(Col),data=aam,trace=F,maxit=30)
> summary(mm2)$varcomp
gamma component std.error z.ratio constraint
R!variance 1.00000000 0.17491494 0.01199393 14.583624 Positive
R!Row.cor 0.62097328 0.62097328 0.02534858 24.497358 Unconstrained
R!Col.cor -0.03744104 -0.03744104 0.03380648 -1.107511 Unconstrained
> aam=multi_norm(30,0.6,0.6)
> mm2=asreml(y~1,rcov=~ar1(Row):ar1(Col),data=aam,trace=F,maxit=30)
> summary(mm2)$varcomp
gamma component std.error z.ratio constraint
R!variance 1.000000000 0.180804271 0.01227539 14.7289994 Positive
R!Row.cor 0.581797580 0.581797580 0.02861663 20.3307541 Unconstrained
R!Col.cor 0.007598536 0.007598536 0.03448510 0.2203426 Unconstrained
> aam=multi_norm(30,0.3,0.6)
> mm2=asreml(y~1,rcov=~ar1(Row):ar1(Col),data=aam,trace=F,maxit=30)
> summary(mm2)$varcomp
gamma component std.error z.ratio constraint
R!variance 1.000000000 0.177888691 0.008979462 19.8106171 Positive
R!Row.cor 0.269572147 0.269572147 0.031823892 8.4707474 Unconstrained
R!Col.cor -0.004159379 -0.004159379 0.035830577 -0.1160846 Unconstrained
> aam=multi_norm(30,0.9,0.6)
> mm2=asreml(y~1,rcov=~ar1(Row):ar1(Col),data=aam,trace=F,maxit=30)
> summary(mm2)$varcomp
gamma component std.error z.ratio constraint
R!variance 1.00000000 0.19194479 0.02674158 7.1777654 Positive
R!Row.cor 0.91213667 0.91213667 0.01247011 73.1458677 Unconstrained
R!Col.cor 0.01203907 0.01203907 0.03474589 0.3464891 Unconstrained
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
