[R-sig-ME] lmer code for multiple random slopes

Peter R Law pr|db @end|ng |rom protonm@||@com
Thu Mar 4 00:08:27 CET 2021


Many thanks Wolfgang! That does clear up my confusion. That interpretation of 'random effects' didn't occur to me.

Thanks also for pointing out that the differences in the estimates provided by lme and lmer occur because they handle singular models differently and won't occur for non-singular examples.

When I have collected enough data, the model I am actually interested in will have two nested random factors, G1 within G2, and I only want random slopes (say for a covariate x) with respect to G2, but random intercepts for both factors. Am I right in thinking the code will take the form

response ~ x + (1|G2/G1) + (x|G2)?

Much appreciated,

Peter


Peter R Law
Research Associate
Center for African Conservation Ecology
Nelson Mandela University
South Africa

Sent with ProtonMail Secure Email.

‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Wednesday, March 3, 2021 7:43 AM, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

> Dear Peter, Dear All,
>
> The error doesn't say 78 variance components, it says 78 random effects. There are 26 groups in this dataset and you are estimating random intercepts and two sets of random slopes. Hence, there are 3*26 = 78 random effects that are in this model. But the dataset has only 74 rows. This is a check built into lmer() to avoid fitting overly complex models for a given dataset. You can disable this check with (also, use REML=FALSE, not "False"):
>
> M2 <- lmer(Response ~ P + A + (P + A | Group), REML=FALSE, data=Trial, control=lmerControl(check.nobs.vs.nRE="ignore"))
>
> Then it will run and this is the same model as fitted with:
>
> M2n <- lme( Response ~ P + A, random = ~ P + A | Group, data=Trial, method="ML")
>
> The results differ somewhat though, because both packages handle such ill-defined problems in different ways.
>
> Best,
> Wolfgang
>
> > -----Original Message-----
> > From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org] On
> > Behalf Of Peter R Law via R-sig-mixed-models
> > Sent: Monday, 01 March, 2021 23:34
> > To: Phillip Alday; r-sig-mixed-models using r-project.org
> > Subject: Re: [R-sig-ME] lmer code for multiple random slopes
> > Thanks Phillip. Sorry about the goof in the summary code. It doesn't actually
> > affect the behaviour that I described, however. Below is a more complete output
> > from R. You will see the str output for the data file, which doesn't seem to
> > indicate that R is treating any of the variables other than Group as a factor (by
> > the way, I recently updated both nlme and lme4 packages but I haven't noticed any
> > changes to their behaviour as a result. But I just noticed that when lme4 loads
> > there is a warning that it was built under R 3.5.3 Is that a possible problem?).
> > For Model 5n (random slopes for both normP and normA) lme returns an analysis but
> > lmer complains about the 78 variance components. Same if one uses A and P instead
> > (Model2ns and 2). If I use variable PR, however, I get the error message you got
> > with lme (Model 7n for normA and normPR) or a different error message (Model 6n
> > for normPR and normA).
> > Yet, lmer seems to behave properly for ~P+A+(P|Group)+(A|Group) (Model 4) and
> > P+PR+(P+PR||Group) (model M3). While these models remove the correlations between
> > the random factors, if lmer is seeing variables as factors rather than numeric,
> > wouldn't lmer still see the wrong number of random factors in these two models?
> > Finally, Models 1 and 1n compare lme and lmer on a model with a single random
> > slope, showing the differences in parameter estimates but not log-likelihood or
> > AIC.
> > I have attached the trial.txt and trial.r files again. I assume you will receive
> > them directly if you want to open them.
> > R output:
> >
> > > library(lme4)
> > > Loading required package: Matrix
> > > Warning message:
> > > package ‘lme4’ was built under R version 3.5.3
> > > library(nlme)
> >
> > Attaching package: ‘nlme’
> > The following object is masked from ‘package:lme4’:
> > lmList
> >
> > > Trial <- read.table(file="C:/PRL/R/RGFRRIBI/Trial.txt", header = TRUE)
> > > names(Trial)
> > > [1] "Response" "P" "A" "PR" "Group"
> > > str(Trial)
> > > 'data.frame': 74 obs. of 5 variables:
> > > $ Response: int 24 22 24 24 23 24 24 24 35 31 ...
> > > $ P : int 15 82 95 71 88 77 93 91 13 82 ...
> > > $ A : num 44.2 39.6 41.4 34.1 44.6 ...
> > > $ PR : num 121 115 250 200 252 ...
> > > $ Group : Factor w/ 26 levels "G1","G10","G11",..: 1 12 12 20 20 21 21 22 23 24
> > > ...
> > > Trial$normP <- as.double(scale(Trial$P))
> > > Trial$normA <- as.double(scale(Trial$A))
> > > Trial$normPR <- as.double(scale(Trial$PR))
> > > names(Trial)
> > > [1] "Response" "P" "A" "PR" "Group" "normP" "normA"
> > > "normPR"
> > > str(Trial)
> > > 'data.frame': 74 obs. of 8 variables:
> > > $ Response: int 24 22 24 24 23 24 24 24 35 31 ...
> > > $ P : int 15 82 95 71 88 77 93 91 13 82 ...
> > > $ A : num 44.2 39.6 41.4 34.1 44.6 ...
> > > $ PR : num 121 115 250 200 252 ...
> > > $ Group : Factor w/ 26 levels "G1","G10","G11",..: 1 12 12 20 20 21 21 22 23 24
> > > ...
> > > $ normP : num -1.783 0.719 1.205 0.308 0.943 ...
> > > $ normA : num 1.576 0.221 0.754 -1.38 1.708 ...
> > > $ normPR : num -1.415 -1.486 0.251 -0.385 0.281 ...
> >
> > > M5n <- lme(Response~normP+normA, random=~normP+normA|Group,data=Trial,
> > > method="ML")
> > > summary(M5n)
> > > Linear mixed-effects model fit by maximum likelihood
> > > Data: Trial
> > > AIC BIC logLik
> > > 536.4562 559.4969 -258.2281
> >
> > Random effects:
> > Formula: ~normP + normA | Group
> > Structure: General positive-definite, Log-Cholesky parametrization
> > StdDev Corr
> > (Intercept) 0.0002514711 (Intr) normP
> > normP 0.0002146005 0
> > normA 0.0001485695 0 0
> > Residual 7.9298211705
> > Fixed effects: Response ~ normP + normA
> > Value Std.Error DF t-value p-value
> > (Intercept) 29.081081 0.9410966 46 30.901270 0.0000
> > normP -0.784566 0.9620920 46 -0.815479 0.4190
> > normA 0.480397 0.9620920 46 0.499326 0.6199
> > Correlation:
> > (Intr) normP
> > normP 0.000
> > normA 0.000 -0.173
> > Standardized Within-Group Residuals:
> > Min Q1 Med Q3 Max
> > -1.8447131 -0.5563605 -0.2613520 0.2755547 5.5643551
> > Number of Observations: 74
> > Number of Groups: 26
> >
> > > M5 <- lmer(Response~normP+normA+(normP+normA|Group), REML="False", data=Trial)
> > > Error: number of observations (=74) <= number of random effects (=78) for term
> > > (normP + normA | Group); the random-effects parameters and the residual variance
> > > (or scale parameter) are probably unidentifiable
> >
> > > M2n <- lme(Response~P+A, random=~P+A|Group,data=Trial, method="ML")
> > > summary(M2n)
> > > Linear mixed-effects model fit by maximum likelihood
> > > Data: Trial
> > > AIC BIC logLik
> > > 536.4562 559.4969 -258.2281
> >
> > Random effects:
> > Formula: ~P + A | Group
> > Structure: General positive-definite, Log-Cholesky parametrization
> > StdDev Corr
> > (Intercept) 3.061144e-04 (Intr) P
> > P 1.535752e-09 0
> > A 8.517774e-13 0 0
> > Residual 7.929821e+00
> > Fixed effects: Response ~ P + A
> > Value Std.Error DF t-value p-value
> > (Intercept) 25.453398 10.828841 46 2.3505192 0.0231
> > P -0.029307 0.035939 46 -0.8154787 0.4190
> > A 0.140819 0.282018 46 0.4993255 0.6199
> > Correlation:
> > (Intr) P
> > P -0.033
> > A -0.975 -0.173
> > Standardized Within-Group Residuals:
> > Min Q1 Med Q3 Max
> > -1.8447131 -0.5563605 -0.2613520 0.2755547 5.5643551
> > Number of Observations: 74
> > Number of Groups: 26
> >
> > > M2 <- lmer(Response~P+A+(P+A|Group), REML="False", data=Trial)
> > > Error: number of observations (=74) <= number of random effects (=78) for term (P
> >
> > -   A | Group); the random-effects parameters and the residual variance (or scale
> >     parameter) are probably unidentifiable
> >
> >
> > > M6n <- lme(Response~A+PR, random=~A+PR|Group,data=Trial, method="ML")
> > > Error in logLik.lmeStructInt(lmeSt, lmePars) :
> > > NA/NaN/Inf in foreign function call (arg 3)
> > > In addition: There were 50 or more warnings (use warnings() to see the first 50)
> > > summary(M6n)
> > > Error in summary(M6n) : object 'M6n' not found
> >
> > > M7n <- lme(Response~normA+normPR, random=~normA+normPR|Group,data=Trial,
> > > method="ML")
> > > Error in lme.formula(Response ~ normA + normPR, random = ~normA + normPR | :
> > > nlminb problem, convergence error code = 1
> > > message = iteration limit reached without convergence (10)
> > > summary(M7n)
> > > Error in summary(M7n) : object 'M7n' not found
> >
> > > M4 <- lmer(Response~P+A+(P|Group)+(A|Group), REML="False", data=Trial)
> > > boundary (singular) fit: see ?isSingular
> > > summary(M4)
> > > Linear mixed model fit by maximum likelihood ['lmerMod']
> > > Formula: Response ~ P + A + (P | Group) + (A | Group)
> > > Data: Trial
> >
> >     AIC      BIC   logLik deviance df.resid
> >
> >
> > 536.5 559.5 -258.2 516.5 64
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -1.8447 -0.5564 -0.2614 0.2756 5.5644
> > Random effects:
> > Groups Name Variance Std.Dev. Corr
> > Group (Intercept) 1.035e-06 1.017e-03
> > P 5.031e-09 7.093e-05 -1.00
> > Group.1 (Intercept) 5.262e-07 7.254e-04
> > A 1.941e-11 4.406e-06 -1.00
> > Residual 6.288e+01 7.930e+00
> > Number of obs: 74, groups: Group, 26
> > Fixed effects:
> > Estimate Std. Error t value
> > (Intercept) 25.45340 10.60707 2.400
> > P -0.02931 0.03520 -0.833
> > A 0.14082 0.27624 0.510
> > Correlation of Fixed Effects:
> > (Intr) P
> > P -0.033
> > A -0.975 -0.173
> > convergence code: 0
> > boundary (singular) fit: see ?isSingular
> >
> > > M3 <- lmer(Response~P+PR+(P+PR||Group), REML="False", data=Trial)
> > > boundary (singular) fit: see ?isSingular
> > > summary(M3)
> > > Linear mixed model fit by maximum likelihood ['lmerMod']
> > > Formula: Response ~ P + PR + ((1 | Group) + (0 + P | Group) + (0 + PR |
> > > Group))
> > > Data: Trial
> >
> >     AIC      BIC   logLik deviance df.resid
> >
> >
> > 530.5 546.7 -258.3 516.5 67
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -1.8106 -0.5567 -0.2981 0.2984 5.0950
> > Random effects:
> > Groups Name Variance Std.Dev.
> > Group (Intercept) 5.422e-04 0.023285
> > Group.1 P 0.000e+00 0.000000
> > Group.2 PR 7.239e-05 0.008508
> > Residual 5.912e+01 7.688728
> > Number of obs: 74, groups: Group, 26
> > Fixed effects:
> > Estimate Std. Error t value
> > (Intercept) 29.265396 3.747575 7.809
> > P -0.023208 0.035100 -0.661
> > PR 0.006151 0.012288 0.501
> > Correlation of Fixed Effects:
> > (Intr) P
> > P -0.633
> > PR -0.754 0.043
> > convergence code: 0
> > boundary (singular) fit: see ?isSingular
> >
> > > M1 <- lmer(Response~P+A+(A|Group), REML="False", data=Trial)
> > > boundary (singular) fit: see ?isSingular
> > > summary(M1)
> > > Linear mixed model fit by maximum likelihood ['lmerMod']
> > > Formula: Response ~ P + A + (A | Group)
> > > Data: Trial
> >
> >     AIC      BIC   logLik deviance df.resid
> >
> >
> > 530.5 546.6 -258.2 516.5 67
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -1.8447 -0.5564 -0.2614 0.2756 5.5644
> > Random effects:
> > Groups Name Variance Std.Dev. Corr
> > Group (Intercept) 0.000e+00 0.000e+00
> > A 1.657e-10 1.287e-05 NaN
> > Residual 6.288e+01 7.930e+00
> > Number of obs: 74, groups: Group, 26
> > Fixed effects:
> > Estimate Std. Error t value
> > (Intercept) 25.45340 10.60707 2.400
> > P -0.02931 0.03520 -0.833
> > A 0.14082 0.27624 0.510
> > Correlation of Fixed Effects:
> > (Intr) P
> > P -0.033
> > A -0.975 -0.173
> > convergence code: 0
> > boundary (singular) fit: see ?isSingular
> >
> > > M1n <- lme(Response~1+P+A, random=~1+A|Group,data=Trial, method="ML")
> > > summary(M1n)
> > > Linear mixed-effects model fit by maximum likelihood
> > > Data: Trial
> > > AIC BIC logLik
> > > 530.4562 546.5847 -258.2281
> >
> > Random effects:
> > Formula: ~1 + A | Group
> > Structure: General positive-definite, Log-Cholesky parametrization
> > StdDev Corr
> > (Intercept) 5.277473e-04 (Intr)
> > A 1.487869e-18 0
> > Residual 7.929821e+00
> > Fixed effects: Response ~ 1 + P + A
> > Value Std.Error DF t-value p-value
> > (Intercept) 25.453398 10.828841 46 2.3505192 0.0231
> > P -0.029307 0.035939 46 -0.8154787 0.4190
> > A 0.140819 0.282018 46 0.4993255 0.6199
> > Correlation:
> > (Intr) P
> > P -0.033
> > A -0.975 -0.173
> > Standardized Within-Group Residuals:
> > Min Q1 Med Q3 Max
> > -1.8447131 -0.5563605 -0.2613520 0.2755547 5.5643551
> > Number of Observations: 74
> > Number of Groups: 26
> > Peter



More information about the R-sig-mixed-models mailing list