lmer code for multiple random slopes

Wed, Mar 3, 2021 4:43 AM

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 at 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 at 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

lmer code for multiple random slopes

Thread (5 messages)