[R] About MC simulation of AR1 model in R

[yz]

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