Skip to content

3 level multivariate multilevel model

2 messages · Aydin,Burak, Thierry Onkelinx

Original
Aydin,Burak
# 
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
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 at ufl.edu>: