How to translate glmer() call to lme(); and including list() for random effects
我之前使用 lme4 包中的 glmer() 运行了混合模型分析。我现在想使用 nlme 包中的 lme() 来运行相同的分析。这是因为随后使用的函数需要输出或调用 lme() 混合模型。
随后使用的函数尝试使用函数segmented.lme() 在数据中查找断点。这个函数的代码可以在这里找到：https://www.researchgate.net/publication/292986444_segmented_mixed_models_in_R_code_and_data
之前，我使用了函数：
|
1
|
global.model <- glmer(response ~ predictor1*predictor2*predictor3*predictor4 + covariate1 + covariate2 + covariate3 + (1|block/transect), data=dat, family=”gaussian”, na.action=”na.fail”)
|
有关可重现的示例，请参见下文。
请注意：随机效应是： (1|block/transect) ，即考虑块之间的交互效应和块内的样带。
现在，我不确定如何重写 lme() 的随机效果部分以完全匹配 glmer() 中使用的代码，特别是因为 segmented.lme() 似乎需要一个”列表”。我尝试了以下方法：
|
1
|
random = list(block = pdDiag(~ 1 + predictor1))
|
请注意：我只对 predictor1 数据中的潜在断点感兴趣。
需要的包：lme4, nlme
参考工作文件可在此处获得：https://www.researchgate.net/publication/292629179_Segmented_mixed_models_with_random_changepoints_in_R
这是数据的一个子集：
|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 |
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”) |
非常感谢您的任何建议。
我对 segmented.lme 不熟悉，但如果它的功能与 nlme 相同(您的问题的开头似乎暗示了这一点)，那么您可以按如下方式指定随机效果。
我以我自己的一些数据为例，因为您的数据集没有包含足够的信息来估计模型。您应该能够为您自己的数据集推断出所需的模型。
|
1
2 3 4 5 6 7 |
library(lme4)
global.model <- lmer(Schaalscore ~ Leeftijd + (1|SCHOOL/LeerlingID),data = Data_RW5, na.action =”na.exclude”) summary(global.model) library(nlme) |
您的模型指示了块和样带上的随机截距，其中样带嵌套在块内。我的数据具有相同的结构，但 LeerlingID 嵌套在 SCHOOL 中。我使用 lmer 而不是 glmer(因为警告消息会显示：calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly)。但是 lmer 和 glmer 的想法是一样的。输出如下：
|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 |
> 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: Random effects: Fixed effects: Correlation of Fixed Effects:
> summary(global.model2) Random effects: Formula: ~1 | LeerlingID %in% SCHOOL Fixed effects: Schaalscore ~ Leeftijd Standardized Within-Group Residuals: Number of Observations: 7798 |
您可以看到随机效应和固定效应估计值相同，并且”REML 收敛标准”等于 -2 * logLik。综上，可以指定随机结构为random= list(block= ~1, transect= ~ 1)，得到相同的模型。
edit：pdDiag 是标准 pdMat 类的一部分，用于指定随机效应的方差-协方差矩阵。您的原始模型仅在两个级别上指定随机截距，因此 pdDiag 不执行任何操作。如果您指定随机斜率和随机截距 pdDiag 将斜率-截距相关性设置为 0。请参阅 Bates
- 我也尝试过：random = list(block = pdDiag(~ 1 predictor1), transect = pdDiag(~ 1 predictor1))。 segmented.lme() 的结果与@Niek 在此处描述的结果相同。这支持关于 pdDiag 的评论。我非常感谢这方面的建议。
来源：https://www.codenong.com/42338979/
