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

Kevin Wright kw.stat at gmail.com
Sat Mar 4 00:08:40 CET 2017

```Try this example to simulate AR1xAR1 variance structure.

Kevin Wright

library(mvtnorm)
library(asreml)
set.seed(300)
nr <- 10; nc <- 8 # number of rows and columns
colcor <- .9      # correlation for columns, rows
rowcor <- .1
sigAR <- diag(nr)
sigAR <- rowcor^ abs(row(sigAR) - col(sigAR))
sigAC <- diag(nc)
sigAC <- colcor ^ abs(row(sigAC) - col(sigAC))
sig <- 100 * kronecker(sigAR, sigAC) # scale 100 doesn't really matter

yy <- rmvnorm(1, mean=rep(0, nr*nc), sig)
dat <- data.frame(y=as.vector(yy), row=rep(1:nr, each=nc), col=rep(1:nc,
nr))
dat\$row=factor(dat\$row)
dat\$col=factor(dat\$col)
m1 <- asreml(y ~ 1, data=dat, rcov= ~ ar1(row):ar1(col))
summary(m1)\$varcomp

##                 gamma   component   std.error    z.ratio    constraint
## R!variance 1.00000000 95.60247530 35.62104468  2.6838762      Positive
## R!row.cor  0.08306616  0.08306616  0.11421749  0.7272631 Unconstrained
## R!col.cor  0.91306348  0.91306348  0.03451116 26.4570518 Unconstrained

On Mon, Nov 28, 2016 at 12:03 AM, yz <yzhlinscau at 163.com> wrote:

> 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]]
>
>
> ______________________________________________
> R-help at 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.
>

--
Kevin Wright

[[alternative HTML version deleted]]

```