[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,]
>
> 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 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]]
More information about the R-help
mailing list