How to translate glmer() call to lme(); and including list() for random effects

Question

I have previously run mixed model analyses using glmer() in package lme4. I would now like to run the very same analyses using lme() in package nlme instead. This is because a subsequently used function requires the output or call of a lme() mixed model.

The subsequently used function attempts to find a breakpoint in the data using the function segmented.lme(). The code for this function can be found here: https://www.researchgate.net/publication/292986444_segmented_mixed_models_in_R_code_and_data

Previously, I used the function:

global.model <- glmer(response ~ predictor1*predictor2*predictor3*predictor4 + covariate1 + covariate2 + covariate3 + (1|block/transect), data=dat, family="gaussian", na.action="na.fail")

For a reproducible example please see below.

Please note: the random effect is: (1|block/transect) , i.e. to account for interaction effects between blocks and transects within blocks.

Now, I am not sure how to rewrite the random effects part for lme() to match exactly the code used in glmer(), and in particular because segmented.lme() seems to require a 'list'. I have tried the following:

random = list(block = pdDiag(~ 1 + predictor1))

Please note: I am interested in a potential breakpoint in the data of predictor1 only.

Required packages: lme4, nlme

A reference working paper is available here: https://www.researchgate.net/publication/292629179_Segmented_mixed_models_with_random_changepoints_in_R

This is a subset of the data:

structure(list(block = structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = c("B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8"), class = "factor"), transect = structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = c("B1L", 
"B1M", "B1S", "B2L", "B2M", "B2S", "B3L", "B3M", "B3S", "B4L", 
"B4M", "B4S", "B5L", "B5M", "B5S", "B6L", "B6M", "B6S", "B7L", 
"B7M", "B7S", "B8L", "B8M", "B8S"), class = "factor"), predictor1 = c(28.63734661, 
31.70995133, 27.40407982, 25.48842992, 21.81094637, 24.02032756
), predictor2 = c(5.002945364, 6.85567854, 0, 22.470422, 
0, 0), predictor3 = c(3.72, 3.55, 3.66, 3.65, 3.53, 3.66), 
predictor4 = c(504.8, 547.6, 499.7, 497.8, 473.8, 467.5), 
covariate1 = c(391L, 394L, 351L, 336L, 304L, 335L), covariate2 = c(0.96671086, 
2.81939707, 0.899512367, 1.024730094, 1.641161861, 1.419433714
), covariate3 = c(0.787505444, 0.641693911, 0.115804751, 
-0.041146951, 1.983567486, -0.451039179), response = c(0.81257636, 
0.622662116, 0.490330786, 0.709929461, -0.156398286, -1.185175095
)), .Names = c("block", "transect", "predictor1", "predictor2", "predictor3", "predictor4", "covariate1", "covariate2", "covariate3", "response"), row.names = c(NA, 6L), class = "data.frame")

Many thanks in advance for any advice.

Niek Niek · Accepted Answer · 2017-03-02T15:16:04

I'm not familiar with segmented.lme but if it functions in the same manner as nlme (which the beginning of your question seems to suggest) then you can specify the random effects as follows.

I used some of my own data as an example because your dataset does not contain sufficient information to estimate a model. You should be able to deduce the desired model for your own dataset.

library(lme4)
    global.model <- lmer(Schaalscore ~ Leeftijd + (1|SCHOOL/LeerlingID),data = Data_RW5, na.action = "na.exclude")
    summary(global.model)

library(nlme)
global.model2 <- lme(Schaalscore ~ Leeftijd, random= list(SCHOOL = ~1, LeerlingID = ~ 1) ,data = Data_RW5, na.action = "na.exclude")
summary(global.model2)

Your model indicates a random intercept on block and transect where transect is nested within blocks. My data has the same structure but then LeerlingID is nested within SCHOOL. I used lmer instead of glmer (as the warning message will show you: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly). But the idea for lmer and glmer is the same. The output is as follows:

> summary(global.model)
Linear mixed model fit by REML ['lmerMod']
Formula: Schaalscore ~ Leeftijd + (1 | SCHOOL/LeerlingID)
   Data: Data_RW5

REML criterion at convergence: 58562.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2088 -0.5855 -0.0420  0.5380  4.6893 

Random effects:
 Groups            Name        Variance Std.Dev.
 LeerlingID:SCHOOL (Intercept) 213.46   14.610  
 SCHOOL            (Intercept)  28.39    5.328  
 Residual                       62.35    7.896  
Number of obs: 7798, groups:  LeerlingID:SCHOOL, 1384; SCHOOL, 59

Fixed effects:
            Estimate Std. Error t value
(Intercept) -89.0261     1.2116  -73.48
Leeftijd     18.3646     0.1081  169.86

Correlation of Fixed Effects:
         (Intr)
Leeftijd -0.725




> summary(global.model2)
Linear mixed-effects model fit by REML
 Data: Data_RW5 
       AIC      BIC    logLik
  58572.08 58606.89 -29281.04

Random effects:
 Formula: ~1 | SCHOOL
        (Intercept)
StdDev:    5.327848

 Formula: ~1 | LeerlingID %in% SCHOOL
        (Intercept) Residual
StdDev:    14.61033  7.89634

Fixed effects: Schaalscore ~ Leeftijd 
                Value Std.Error   DF   t-value p-value
(Intercept) -89.02613 1.2116148 6413 -73.47726       0
Leeftijd     18.36460 0.1081172 6413 169.85827       0
 Correlation: 
         (Intr)
Leeftijd -0.725

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.2087839 -0.5855190 -0.0420062  0.5379625  4.6892515 

Number of Observations: 7798
Number of Groups: 
                SCHOOL LeerlingID %in% SCHOOL 
                    59                   1384

You can see that the random and fixed effect estimates are the same and 'REML criterion at convergence' is equal to -2 * logLik. In summary, you can specify the random structure as random= list(block= ~1, transect= ~ 1) to get the same model.

edit: pdDiag is part of the standard pdMat classes used to specify the variance-covariance matrix for the random effects. Your original model only specifies a random intercept on two levels, so pdDiag does not do anything. If you specify a random slope and random intercept pdDiag sets the slope-intercept correlation to 0. See Bates & Pinheiro (2000) for details.

How to translate glmer() call to lme(); and including list() for random effects

1 Answers