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

Viechtbauer, Wolfgang (SP) wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Wed Mar 3 13:43:02 CET 2021


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