[R-sig-ME] 3 level multivariate multilevel model

[Thierry Onkelinx]

Dear Aydin,

The lme4 package has no infrastructure to fit variance and correlation
structures. You need nlme for those.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2017-01-02 5:43 GMT+01:00 Aydin,Burak <burakaydin op ufl.edu>:

> Hello all,
>
> This is my first post.
>
> I am trying to extract covariance matrices.
>
> I can do it with lme (the model is given by Snijders), please see below;
>
> However, the way I do it is just ugly and inefficient. BTW It is a lot
> easier with mplus.
>
>
>
> Can you help me to get the same covariance matrices using lmer?
>
> Best
>
>
> # 3 level multivariate empty model
> ## Attempts to extract the covariance matrices
>
> ## load Rdata from an online repository
> urlfile2='https://github.com/burakaydin/materyaller/blob/
> gh-pages/EGEintRo/tempdata.Rdata?raw=true'
> load(url(urlfile2))
> rm(urlfile2)
>
> head(tempdat)
>
> ##reshape data
> library(tidyr)
> longeb = gather(tempdat, which, response, x:x2,factor_key=TRUE)
> head(longeb)
> longeb <-
>   longeb[order(longeb$schid,
>                longeb$cid2,
>                longeb$id),]
>
> ## run multivariate multilevel model from Snijder's website
> library(nlme)
> fit <- lme(response ~ - 1 + which, random = ~ -1 + which|schid/cid2,
>            weights=varIdent(form=~1|which),
>            corr=corSymm(form=~as.numeric(which)|schid/cid2/id),
>            data=longeb, method="REML",
>            control = list(maxIter=500, msMaxIter=500, tolerance=1e-8,
>                           niterEM=250))
>
> ##extract covariance matrices
>
> ###level3
>
> l3x1var=as.numeric(VarCorr(fit)[[2,1]])
> l3x2var=as.numeric(VarCorr(fit)[[3,1]])
> l3cor=as.numeric(VarCorr(fit)[[3,3]])
> l3cov=l3cor*sqrt(l3x1var)*sqrt(l3x2var)
>
> ####level 3 cov matrix
> sigma33=matrix(c(l3x1var,l3cov,l3cov,l3x2var),ncol=2)
> sigma33
>
> ###level 2
> l2x1var=as.numeric(VarCorr(fit)[[5,1]])
> l2x2var=as.numeric(VarCorr(fit)[[6,1]])
> l2cor=as.numeric(VarCorr(fit)[[6,3]])
> l2cov=l2cor*sqrt(l2x1var)*sqrt(l2x2var)
>
> #### level 2 matrix
> sigma22=matrix(c(l2x1var,l2cov,l2cov,l2x2var),ncol=2)
> sigma22
>
> ### Level-1 (shame)
>
> l1x1var=as.numeric(VarCorr(fit)[[7,1]])
> l1x2var=l1x1var
> capture.output( summary(fit) , file='fit.txt')
> mystring <- read.table("fit.txt",sep="\n")
> l1cor=factor(mystring[22,1])
> l1cor=levels(l1cor)
> l1cor=strsplit(l1cor," ")
> l1cor=as.numeric(l1cor[[1]][2])
> l1cov=l1cor*sqrt(l1x1var)*sqrt(l1x2var)
>
> #### level 1 covariance matrix
> sigma11=matrix(c(l1x1var,l1cov,l1cov,l1x2var),ncol=2)
> sigma11
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]