dimanche 19 février 2017

How to rewrite glmer() function using lme() instead and including list() for random effects

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: http://ift.tt/2m2YCwR

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: http://ift.tt/2lBDX5C

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.




Aucun commentaire:

Enregistrer un commentaire