[R-sig-ME] lmer code for multiple random slopes
Peter R Law
pr|db @end|ng |rom protonm@||@com
Mon Mar 1 23:34:11 CET 2021
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
Sent with ProtonMail Secure Email.
‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Thursday, February 25, 2021 10:26 AM, Phillip Alday <me using phillipalday.com> wrote:
> On 25/2/21 5:16 am, Peter R Law via R-sig-mixed-models wrote:
>
> > Alas I am still puzzled. I have extracted some data in trial.txt (no categorical variables) and attached some code in trial.R. I am only using this data to test code, which I hope to apply to a larger data set obtained from multiple populations, so with more structure. So the results themselves are not important, only whether the code does what I think it should. It occurred to me to test each model with lme and lmer.
> > For a model with only a random intercept plus fixed effects, lme and lmer returned the same results, except that lme returns estimates with more decimal places (so here lmer returned zero variance for the random intercpt and lme a very small number).
>
> This is to be expected: lme can't handle singular models (e.g. models
> with an exactly zero variance component), but lmer can. So lme just gets
> as close as it can.
>
> > Adding a random slope, the main difference in the results is that lme and lmer returned very different estimates of the slope variance. Is that surprising?
>
> It depends on what's going wrong. Generally these should return
> identical fits for converged models (with the fine print that there are
> certain models that you can specify in one but not the other without a
> lot of effort).
>
> > For a model with two random slopes, lme returned results as expected but lmer still claims there are 78 variance components, which is the number one would get from a 12 x 12 covariance matrix. Where the number 12 comes from is a mystery to me, especially as lme does what I expected (I ran this model with both raw data and normalized predictor values to see if something fishy was happening there but no):
>
> The list strips most attachment types, so your script didn't make it
> through. But just your data was enough for me to find a few problems
>
> First, you're fitting and displaying two different models here:
>
> > > M5n <- lme(Response~normP+normA, random=~normP+normA|Group,data=Trial, method="ML")
> > > summary(M2n)
>
> Second, when I try to fit the model that's given by the summary here, I
> get a numerical error in lme:
>
> > lme(Response~P+A, random=~P+A|Group,data=Trial, method="ML")
>
> Error in lme.formula(Response ~ P + A, random = ~P + A | Group, data =
> Trial, :
> nlminb problem, convergence error code = 1
> message = false convergence (8)
>
> This makes me think that something is wrong wrong in both bits of
> software, but maybe the versions on your local machine are catching it
>
> > 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
> >
> > > 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.
>
> This screams that your normP and normA variables aren't being handled as
> continuous but rather categorical/factors. This means that you have a
> lot more slopes being estimated, which the model survives, but not the
> extra correlation parameters (hence the success when you suppress them
> with the two syntax variants you mention below).
>
> > It didn't matter which pair of the three predictors I used as far as the error message lmer returned but lme only returned results for P+A, not any pair involving PR, which apparently created computational issues, distinct to the interpretation of the code.
> > In lmer, replacing the code (x+y|Group) by either (x+y||Group) or (X|Group) + (Y|Group) returned the expected results (in that lme4 interpreted the code as expected). Is there lme code for either of these models?
> > Many thanks for any help.
> > Peter
> > P. S. I just noticed that responding to Ben's email sends my email to Ben, not to r-sig. I hope that sending my email to r-sig is the right thing to do and doesn't break the chain to my previous email.
>
> Yep! Frequently, we have the reverse problem: people forget to keep the
> list in CC. So thanks! The matching in the threading is handled mostly
> by the subject line.
>
> > Sent with ProtonMail Secure Email.
> > ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
> > On Wednesday, February 17, 2021 10:40 PM, Peter R Law prldb using protonmail.com wrote:
> >
> > > Thanks for your quick response. I take it that the code should mean what I thought it would and that somehow lmer is not interpreting what I wrote in my actual example as intended. None of the variables are categorical but I'll give it some more thought and see if I can figure it out. If not, I'll provide more details.
> > > Peter
> > > Sent with ProtonMail Secure Email.
> > > ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
> > > On Tuesday, February 16, 2021 10:04 AM, Ben Bolker bbolker using gmail.com wrote:
> > >
> > > > I second Phillip's point. The example below works as expected (gets
> > > > a singular fit, but there are 6 covariance parameters as expected).
> > > > Based on what you've told us so far, the most plausible explanation is
> > > > that one or both of your covariates (x and/or z) are factors
> > > > (categorical) rather than numeric.
> > > > Ben Bolker
> > > >
> > > > ==========================================================================================================================================================================================================================================================================================================================
> > > >
> > > > set.seed(101)
> > > > dd <- data.frame(x=rnorm(500),z=rnorm(500),
> > > > g=factor(sample(1:6,size=500,replace=TRUE)))
> > > > form <- y ~ x + z + (x+z|g)
> > > > dd$y <- simulate(form[-2],
> > > > newdata=dd,
> > > > newparams=list(beta=rep(0,3),
> > > > theta=rep(1,6),
> > > > sigma=1))[[1]]
> > > > library(lme4)
> > > > m1 <- lmer(form, data=dd)
> > > > VarCorr(m1)
> > > > On 2/16/21 8:18 AM, Phillip Alday wrote:
> > > >
> > > > > I suspect we'll need to know a bit more about your data to answer this
> > > > > question. Can you share it in any form (e.g. variables renamed and
> > > > > levels of factors changed to something opaque) ?
> > > > > Best,
> > > > > Phillip
> > > > > On 16/2/21 4:02 am, Peter R Law via R-sig-mixed-models wrote:
> > > > >
> > > > > > I am trying to fit a model with two covariates, x and z say, for response y, with a random factor g and want each of x and y to have a random slope. I expected
> > > > > > lmer(y ~ x + z + (x+z|g),...)
> > > > > > to fit a model with 6 random variance components, the intercept, two slopes and three correlations. But I got an error message saying there were 74 random variance components and my data was insufficient to fit the model. Yet
> > > > > > lmer(y ~ x + z + (x+z||g),...)
> > > > > > returned what I expected, a model with the random intercept and two slopes but no correlations. How is lmer interpreting the first line of code above and how I would code for what I want. I have not been able to find any examples in the literature or online that help me but I may have easily missed something so if anyone knows of a useful link that'd be great. The only examples of multiple random slopes I've seen take the form
> > > > > > lmer(y~x + z +(x|g) + (z|g),...)
> > > > > > specifically excluding correlations between the random slopes and intercept of the two predictors. Even if the latter is a more sensible approach I'd like to understand the coding issue.
> > > > > > Thanks.
> > > > > > Peter
> > > > > > Sent with ProtonMail Secure Email.
> > > > > > [[alternative HTML version deleted]]
> > > > > > R-sig-mixed-models using r-project.org mailing list
> > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > > > >
> > > > > R-sig-mixed-models using r-project.org mailing list
> > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > > >
> > > > R-sig-mixed-models using r-project.org mailing list
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> > R-sig-mixed-models using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Trial.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20210301/5e7bf7c6/attachment.txt>
More information about the R-sig-mixed-models
mailing list