[R] How to mimic pdMat of lme under lmer?
Douglas Bates
dmbates at gmail.com
Tue Sep 20 01:27:51 CEST 2005
On 9/19/05, Joris De Wolf <joris.dewolf at cropdesign.com> wrote:
> Dear members,
>
> I would like to switch from nlme to lme4 and try to translate some of my
> models that worked fine with lme.
> I have problems with the pdMat classes.
>
> Below a toy dataset with a fixed effect F and a random effect R. I gave
> also 2 similar lme models.
> The one containing pdLogChol (lme1) is easy to translate (as it is an
> explicit notation of the default model)
> The more parsimonious model with pdDiag replacing pdLogChol I cannot
> reproduce with lmer. The obvious choice for me would be my model lmer2,
> but this is yielding different result.
>
> Somebody any idea?
> Thanks,
> Joris
>
> I am using R version 2.1.0 for Linux
> and the most recent downloads of Matrix and nlme
>
> #dataset from McLean, Sanders and Stroup
> F <- factor(c(1,1,1,1,1,1,2,2,2,2,2,2))
> R <- factor(c(1,1,2,2,3,3,1,1,2,2,3,3))
> s <-
> c(51.43,51.28,50.93,50.75,50.47,50.83,51.91,52.43,52.26,52.33,51.58,51.23)
> DS <- data.frame(F,R,s)
> DS$F <- as.factor(DS$F)
> DS$R <- as.factor(DS$R)
>
> library(nlme)
> lme1 <- lme(data = DS,s ~ F,random = list(R = pdLogChol(~F)))
> lme2 <- lme(data = DS,s ~ F,random = list(R = pdDiag(~F)))
> summary(lme1)
> summary(lme2)
>
> library(lme4)
> lmer1 <- lmer(data = DS,s ~ F + (F|R))
> lmer2 <- lmer(data = DS,s ~ F + (1|R) + (1|F))
> summary(lmer1)
> summary(lmer2)
You need to generate the two columns of the indicator matrix for F as
separate model matrices. This will involve some awkward construction
like
> (fm1 <- lmer(s ~ F +(F|R), DS))
Linear mixed-effects model fit by REML
Formula: s ~ F + (F | R)
Data: DS
AIC BIC logLik MLdeviance REMLdeviance
20.69259 23.60203 -4.346297 6.03915 8.692593
Random effects:
Groups Name Variance Std.Dev. Corr
R (Intercept) 0.108796 0.32984
F2 0.102008 0.31939 -0.014
Residual 0.048525 0.22028
# of obs: 12, groups: R, 3
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 50.9483 0.2106 10 241.9188 < 2.2e-16
F2 1.0083 0.2240 10 4.5014 0.001141
> (fm2 <- lmer(s ~ F + (0+as.numeric(F==1)|R) + (0+as.numeric(F==2)|R), DS))
Linear mixed-effects model fit by REML
Formula: s ~ F + (0 + as.numeric(F == 1) | R) + (0 + as.numeric(F ==
2) | R)
Data: DS
AIC BIC logLik MLdeviance REMLdeviance
19.62621 22.05075 -4.813107 7.439584 9.626215
Random effects:
Groups Name Variance Std.Dev.
R as.numeric(F == 1) 0.108796 0.32984
R as.numeric(F == 2) 0.207896 0.45596
Residual 0.048525 0.22028
# of obs: 12, groups: R, 3; R, 3
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 50.94833 0.21060 10 241.9188 < 2e-16
F2 1.00833 0.34891 10 2.8899 0.01611
There are probably more elegant ways of doing this but I don't think
there is any really clean way. If you want to assume that all the
variances are equal then you can estimate the model using an
interaction.
> (fm3 <- lmer(s ~ F + (1|F:R), DS))
Linear mixed-effects model fit by REML
Formula: s ~ F + (1 | F:R)
Data: DS
AIC BIC logLik MLdeviance REMLdeviance
17.77917 19.7188 -4.889587 7.669022 9.779175
Random effects:
Groups Name Variance Std.Dev.
F:R (Intercept) 0.158346 0.39793
Residual 0.048525 0.22028
# of obs: 12, groups: F:R, 6
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 50.94833 0.24672 10 206.5049 < 2e-16
F2 1.00833 0.34891 10 2.8899 0.01611
>
>
>
> confidentiality notice:
> The information contained in this e-mail is confidential and...{{dropped}}
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
More information about the R-help
mailing list