평균 모델에서 불연속 변수에 대한 플롯 모델 적합

Question

선형 혼합 모델 집합이 있고 평균 모델을 만들었습니다. 평균 모델에 포함 된 두 가지 수준의 요인에 대한 모델 적합성을 플로팅하고 싶습니다. 간단한 예 :

library(lme4) 
library(MuMIn) 

mtcars2 <- mtcars 
mtcars2$vs <- factor(mtcars2$vs) 

gl <- lmer(mpg ~ am + disp + hp + qsec + (1 | cyl), mtcars2, 
      REML = FALSE, na.action = 'na.fail') 
d <- dredge(gl) 

av <- model.avg(d, subset = cumsum(weight) <= 0.95) 
summary(av) 

Call: 
model.avg(object = d, subset = cumsum(weight) <= 0.95) 

Component model call: 
lme4::lmer(formula = mpg ~ <7 unique rhs>, data = mtcars2, REML = FALSE, na.action = na.fail) 

Component models: 
    df logLik AICc delta weight 
13 5 -77.81 167.92 0.00 0.37 
123 6 -76.34 168.05 0.13 0.35 
134 6 -77.54 170.43 2.51 0.11 
1234 7 -76.25 171.16 3.24 0.07 
23 5 -79.85 172.00 4.08 0.05 
2  4 -81.63 172.75 4.83 0.03 
124 6 -78.99 173.34 5.42 0.02 

Term codes: 
    am disp hp qsec 
    1 2 3 4 

Model-averaged coefficients: 
(full average) 
      Estimate Std. Error Adjusted SE z value Pr(>|z|)  
(Intercept) 25.457505 6.467643 6.648016 3.829 0.000129 *** 
am   4.103425 1.861593 1.898182 2.162 0.030636 * 
hp   -0.043829 0.017926 0.018265 2.400 0.016415 * 
disp  -0.009419 0.011834 0.011983 0.786 0.431821  
qsec   0.081973 0.284147 0.292015 0.281 0.778929  

(conditional average) 
      Estimate Std. Error Adjusted SE z value Pr(>|z|)  
(Intercept) 25.45751 6.46764  6.64802 3.829 0.000129 *** 
am   4.46519 1.46823  1.51835 2.941 0.003273 ** 
hp   -0.04651 0.01471  0.01515 3.070 0.002140 ** 
disp  -0.01793 0.01068  0.01099 1.632 0.102634  
qsec   0.40421 0.51757  0.53873 0.750 0.453075  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Relative variable importance: 
        hp am disp qsec 
Importance:   0.94 0.92 0.53 0.20 
N containing models: 5 5 5 3 

나는 전체 평균 모델에 의해 추정으로 am의 효과를 플롯합니다.

일반적으로 lsmeans::lsmeans(gl, ~am) 또는 lmerTest::lsmeansLT(gl, 'am')을 사용하고 두 그룹에 대한 최소 제곱 평균과 신뢰 구간을 그립니다.

평균 모델의 경우 어떻게합니까?

Answer 1

(일부 토론 및 추가 연구 결과는. 나는 emmeans 패키지 저자있어 참고 한 후이, 수정 답변입니다.) 여기

작업이 나타납니다 무언가이다. 일이 없기 때문에

library(emmeans) 

terms.averaging = function(x, ...) 
    terms(x$formula) 

recover_data.averaging = emmeans:::recover_data.lm 
    ### NOTE: still have to provide 'data' argument 

emm_basis.averaging = function(object, trms, xlev, grid, ...) { 
    bhat = coef(object, full = TRUE) 
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev) 
    X = model.matrix(trms, m, contrasts.arg = object$contrasts) 
    V = vcov(object, full = TRUE) 
    dffun = function(k, dfargs) NA 
    dfargs = list() 
    list(X=X, bhat=bhat, nbasis=estimability::all.estble, V=V, 
     dffun=dffun, dfargs=dfargs, misc=list()) 
} 

terms 방법이 필요하다 :

첫째, emmeans 패키지에 필요한 방법을 정의한다. 다른 두 개는 lm 개체의 기존 메서드에서 수정되었습니다. 이제 한 가지 잡을 수 있습니다 : vcov() 호출에는 객체에 NULL"modelList"이 아닌 속성이 있어야합니다. 그리고 av 개체가 실패합니다. 그러나 model.avg에 대한 도움말 페이지 하단의 예는 무엇을 보여줍니다

cs95 = get.models(d, cumsum(weight) <= .95) 
AV = model.avg(cs95) 

지금, AV가 필요한 속성이 있습니다. 이제 우리는 얻을 : 대조 결과가 모델 요약 테이블에 av에 대한 추정 및 보정되지 않은 SE 일치

em = emmeans(AV, ~ am, at = list(am = c("0", "1")), data = mtcars) 
em 
## am emmean  SE df asymp.LCL asymp.UCL 
## 0 15.42665 2.985460 NA 9.575257 21.27805 
## 1 19.53008 1.986149 NA 15.637297 23.42286 

pairs(em) 
## contrast estimate  SE df z.ratio p.value 
## 0 - 1 -4.103425 1.861593 NA -2.204 0.0275 

하는 것으로.

참고coef(..., full = FALSE)를 사용 vcov(... full = FALSE) 부정적인 분산 결과가 아닌 포지티브 명확한 공분산 행렬을 수득는 EMM을위한 추정한다.

그리고 이것은 계산적으로 작동하는 것처럼 보이지만 그 대답이 옳다는 것을 의미하지는 않습니다.