R：如何从 lmer() 模型的回归结果中提取协变量 p 值列表？

Question

我有一个带有 8 个预测变量的混合 lmer 模型示例，我想提取协变量的名称、它们的系数、它们的标准误差和它们的 p 值，并将它们放入一个矩阵中，这样我就可以将它们写到一个 . .csv。

我已将前 3 个提取到列中，但我不知道如何提取 p 值。你怎么做到这一点？它是 vcov 还是 getME() 的变体？

模型和摘要如下所示：

mod <- lmer(outcome ~ predictor1 + etc...
summary(mod)

Generalized linear mixed model fit by the Laplace approximation 
Formula: Freq ~ pm.lag0 + pm.lag1 + pm.lag2 + pm.lag3 + pm.lag4 + pm.lag5 
+ temp13 + temp013 + rh13 + rh013 + (1 | county) 
   Data: dt 
  AIC  BIC logLik deviance
 3574 3636  -1775     3550
Random effects:
 Groups Name        Variance Std.Dev.
 county (Intercept) 1.6131   1.2701  
Number of obs: 1260, groups: county, 28

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.9356504  0.2614892  11.227  < 2e-16 ***
pm.lag0      0.0012996  0.0005469   2.376 0.017494 *  
pm.lag1      0.0005021  0.0005631   0.892 0.372568    
pm.lag2      0.0009126  0.0005596   1.631 0.102893    
pm.lag3     -0.0007073  0.0005678  -1.246 0.212896    
pm.lag4      0.0031566  0.0005316   5.939 2.88e-09 ***
pm.lag5      0.0019598  0.0005359   3.657 0.000255 ***
temp13      -0.0028040  0.0007315  -3.833 0.000126 ***
temp013     -0.0023532  0.0009683  -2.430 0.015087 *  
rh13         0.0058769  0.0009909   5.931 3.01e-09 ***
rh013       -0.0028568  0.0006070  -4.706 2.52e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
        (Intr) pm.lg0 pm.lg1 pm.lg2 pm.lg3 pm.lg4 pm.lg5 temp13 tmp013 rh13  
pm.lag0 -0.025                                                               
pm.lag1 -0.032 -0.154                                                        
pm.lag2 -0.021  0.044 -0.179                                                 
pm.lag3  0.002  0.003  0.033 -0.176                                          
pm.lag4  0.016  0.102 -0.016  0.041 -0.176                                   
pm.lag5  0.008  0.027  0.090 -0.002  0.040 -0.186                            
temp13  -0.316  0.026  0.027  0.004 -0.019 -0.055 -0.035                     
temp013  0.030 -0.015  0.051  0.015 -0.015  0.002 -0.069 -0.205              
rh13    -0.350  0.043  0.078  0.056 -0.012 -0.042 -0.030  0.430  0.055       
rh013    0.193 -0.008 -0.021  0.011  0.030  0.101 -0.028 -0.278  0.025 -0.524

我已经在这里为 p-value 列留了一个空格并为它输入了一个 colname，所以这个代码示例无法操作：

mixed.results <- mod
cbind(names(fixef(mod)),as.numeric(fixef(mod)),sqrt(diag(vcov(mod))),  ????  )
mixed.results
colnames(mixed.results) <- c("Pred", "Coef", "St. Error", "Pr(>|z|)")
mixed.results
write.csv(mixed.results, file="mixedmod1.csv")

谢谢！

Answer 1

这只是coef(summary(model))，我相信：

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                   data = cbpp, family = binomial)
cc <- coef(summary(gm1))
str(cc)
# num [1:4, 1:4] -1.376 -1.058 -1.196 -1.638 0.205 ...
# - attr(*, "dimnames")=List of 2
#  ..$ : chr [1:4] "(Intercept)" "period2" "period3" "period4"
#  ..$ : chr [1:4] "Estimate" "Std. Error" "z value" "Pr(>|z|)"
cc[,4] ## or cc[,"Pr(>|z)"] to be more explicit
# (Intercept)      period2      period3      period4 
#1.907080e-11 1.996120e-41 4.634385e-43 4.657952e-47 

我使用了lme4的开发版，但我认为这已经有一段时间了。