如何將glmer（）調用轉換爲lme（）;幷包括用於隨機效應的列表（）

Question

我以前在軟件包lme4中使用glmer（）運行混合模型分析。我現在想在包nlme中使用lme（）運行相同的分析。這是因爲隨後使用的函數需要輸出或調用lme（）混合模型。

隨後使用的函數嘗試使用函數segmented.lme（）在數據中查找斷點。這個函數的代碼可以在這裏找到：https://www.researchgate.net/publication/292986444_segmented_mixed_models_in_R_code_and_data

以前，我所使用的函數：

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

對於重複的例子，請參見下文。

請注意：隨機效應是：（1 |塊/橫切），即爲了解決塊與塊之間的橫切之間的交互作用。

現在，我不知道如何重寫lme（）中的隨機效果部分以完全匹配glmer（）中使用的代碼，特別是因爲segmented.lme（）似乎需要「列表」。我曾嘗試以下方法：

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

請注意：我只對predictor1的數據中的潛在斷點感興趣。

需要的軟件包：lme4，NLME

參考工作文件，請訪問：https://www.researchgate.net/publication/292629179_Segmented_mixed_models_with_random_changepoints_in_R

這是數據的一個子集：

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") 

提前任何建議非常感謝。

Answer 1

我對segmented.lme不熟悉，但是如果它的功能與nlme相同（您的問題的開頭似乎表明了這一點），那麼您可以指定隨機效果，如下所示。

我用我自己的一些數據作爲例子，因爲你的數據集沒有足夠的信息來估計一個模型。你應該能夠推導出你自己的數據集所需的模型。

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) 

您的模型指示塊和橫斷面上的橫斷面嵌套在塊中的隨機截距。我的數據具有相同的結構，但LeerlingID嵌套在SCHOOL中。我用lmer代替glmer（因爲警告信息會告訴你：calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly）。但是lmer和glmer的想法是一樣的。輸出如下：

> 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 

你可以看到，隨機和固定效應的估計是一樣的，「在收斂REML標準」 = -2 * logLik。總之，您可以指定隨機結構爲random= list(block= ~1, transect= ~ 1)以獲得相同的模型。

編輯：pdDiag是標準pdMat類的一部分，用於指定隨機效應的方差 - 協方差矩陣。您的原始模型只能在兩個級別上指定一個隨機攔截，所以pdDiag不會執行任何操作。如果指定隨機斜率和隨機截距，則pdDiag將斜率截距相關設置爲0.有關詳細信息，請參閱Pinheiro（2000）的Bates &。