# [R] About MC simulation of AR1 model in R

yz yzhlinscau at 163.com
Mon Nov 28 07:03:33 CET 2016

```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,]

for(i in 1:data_num){
for(j in 1:data_num){
}
}

}

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]]

```