2012-04-04 153 views
2

我有一個示例混合lmer模型與8個預測變量,我想提取協變量的名稱,它們的係數,它們的標準錯誤和它們的p值,並將它們放入一個矩陣中,這樣我就可以將它們寫入一個.csv 。R:如何從lmer()模型的迴歸結果中提取協變量p值列表?

我已經提取了第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值列中的空間,進入了colname的吧,這樣的代碼ISN的這個樣本't操作:

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

謝謝!

+0

我首先看一下彙總函數的源代碼。作爲摘要的一部分,計算它的可能性很大。下面是'lm()'模型的一個相關答案,我想可以很容易地採用它:http://stackoverflow.com/questions/5587676/pull-out-p-values-and-r-squared-from-a-線性迴歸/ 5587781#5587781 – Chase 2012-04-04 22:16:50

+0

嗨大通,是的,我看到那篇文章,但命令似乎並沒有轉移到lmer()。你說得對,p值已經在總結中計算了 - 他們是第4列。有沒有辦法像summary [,4]那樣做? – mEvans 2012-04-04 22:18:50

回答

4

這只是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開發版本,但我認爲這已經有一段工作。

+0

好吧!所以,在此之後,我可以使用as.data.frame(cc),然後調出我想要的特定列!或者......只是使用cc作爲自己......哈哈。謝謝! – mEvans 2012-04-04 22:23:14

相關問題