在lmer迴歸中預測（），但我只需要它2個類別

Question

我正在嘗試估計多級模型。我的代碼是：

fullModel2 <- lmer(pharmexp_2001 ~ gdp_1000_gm + health_exp_per_cap_1000_gm + life_exp + 
        labour_cost_1000_gm + (year_gm|lowerID), data=adat, REML=F) 

導致以下模型：

Linear mixed model fit by maximum likelihood ['lmerMod'] 
Formula: pharmexp_2001 ~ gdp_1000_gm + health_exp_per_cap_1000_gm + life_exp +  
     labour_cost_1000_gm + (year_gm | lowerID) 
    Data: adat 

    AIC  BIC logLik deviance df.resid 
    1830.2 1859.9 -906.1 1812.2  191 

Scaled residuals: 
    Min  1Q Median  3Q  Max 
-2.5360 -0.6853 -0.0842 0.4923 4.0051 

Random effects: 
Groups Name  Variance Std.Dev. Corr 
lowerID (Intercept) 134.6851 11.6054  
      year_gm  0.4214 0.6492 -1.00 
Residual    487.5324 22.0801  
Number of obs: 200, groups: lowerID, 2 

Fixed effects: 
          Estimate Std. Error t value 
(Intercept)    -563.7924 75.4125 -7.476 
gdp_1000_gm     -0.9050  0.2051 -4.413 
health_exp_per_cap_1000_gm 37.5394  6.3943 5.871 
life_exp      8.8571  0.9498 9.326 
labour_cost_1000_gm   -1.3573  0.4684 -2.898 

Correlation of Fixed Effects: 
      (Intr) g_1000 h____1 lif_xp 
gdp_1000_gm -0.068      
hl____1000_ 0.374 -0.254    
life_exp -0.996 0.072 -0.393  
lbr_c_1000_ -0.133 -0.139 -0.802 0.142 

我知道，這是一個問題，即相關係數爲-1通過隨機效應，但我有一個更大的問題。我必須繪製我的結果，但只有我需要2行：當lowerID=0和lowerID=1。所以我想繪製y軸上的pharmaexp_2001與x軸上的year，但我只需要2行（通過lowerID）。我知道我必須使用predict.merMod，但如何繪製這些結果，僅繪製這兩行？目前我的情節有21行（因爲我分析了21個國家的藥品支出）。

Answer 1

歡迎來到本站，@EszterTakács！

您只需要在newdata中指定兩個ID。這裏是一個基於sleepstudy數據的例子，在R。我假設你想在y軸上繪製預測值。只需將代碼替換爲您的數據和變量，即可獲得lowerID==0和lowerID==1的預測值。然後，您可以使用您的代碼爲兩個ID繪製兩條線。

> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=F)) 
Linear mixed model fit by maximum likelihood ['lmerMod'] 
Formula: Reaction ~ Days + (Days | Subject) 
    Data: sleepstudy 
     AIC  BIC logLik deviance 
1763.9393 1783.0971 -875.9697 1751.9393 
Random effects: 
Groups Name  Std.Dev. Corr 
Subject (Intercept) 23.781  
      Days   5.717 0.08 
Residual    25.592  
Number of obs: 180, groups: Subject, 18 
Fixed Effects: 
(Intercept)   Days 
    251.41  10.47 

> newdata = sleepstudy[sleepstudy$Subject==308 | sleepstudy$Subject==333,] 
> str(p <- predict(fm1,newdata)) # new data, all RE 
Named num [1:20] 254 274 293 313 332 ... 
- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...