내 혼합 모델에서 lme4 예측 함수를 사용하는 데 문제가 있습니다.

Question

혼합 모델에서 lme4 예측 함수를 사용하려는 데 어려움을 겪고 있습니다. 술어를 만들 때, 내 설명 변수 중 일부를 특정 레벨로 설정할 수 있지만 다른 변수의 평균을 설정할 수 있기를 원합니다.

a <- data.frame(
    TLR4=factor(rep(1:3, each=4, times=4)), 
    repro.state=factor(rep(c("a","j"),each=6,times=8)), 
    month=factor(rep(1:2,each=8,times=6)), 
    sex=factor(rep(1:2, each=4, times=12)), 
    year=factor(rep(1:3, each =32)), 
    mwalkeri=(sample(0:15, 96, replace=TRUE)), 
    AvM=(seq(1:96)) 
) 

AVM 번호가 물 밭쥐 식별 번호는 다음과 같습니다

다음은 내 원래의 데이터 세트의 단순화, 말도 안되는 버전입니다 몇 가지로 구성 데이터를합니다. 응답 변수 ( mwalkeri)는 각 벼룩에있는 벼룩의 개수입니다. 관심있는 주요 변수는 Tlr4로 3 가지 유전자형 (1, 2 및 3 코드)이있는 유전자입니다. 다른 설명 변수로는 생식 상태 (성인 또는 청소년), 달 (1 또는 2), 성별 (1 또는 2) 및 연도 (1, 2 또는 3)가 있습니다. 모두를 차지하면서

install.packages("lme4") 
library(lme4) 
mm <- glmer(mwalkeri~TLR4+repro.state+month+sex+year+(1|AvM), data=a, 
    family=poisson,control=glmerControl(optimizer="bobyqa"))` 
summary(mm) 

내가 각각 다른 TLR4 유전자형에 대한 기생충 부담에 대한 예측을 할 : 내 모델 (물론이 모델은 지금까지 만들어진 데이터를 부적절하지만 그건 중요하지한다)과 같다 다른 공변량. 이를 위해 나는에 설명 변수의 각을 설정하기 위해 원하는 레벨을 지정하기 위해 새로운 데이터 집합을 생성하고, 예측 기능 사용 :

b <- data.frame(
    TLR4=factor(1:3), 
    repro.state=factor(c("a","a","a")), 
    month=factor(rep(1, times=3)), 
    sex=factor(rep(1, times=3)), 
    year=factor(rep(1, times=3)) 
) 
predict(mm, newdata=b, re.form=NA, type="response") 

이 일을 한을하지만 난 정말 평균에 걸쳐 년에 선호하는 대신에 설정하는 것 특정 수준으로 1 년. 그러나 평균 년을 시도 할 때마다 다음과 같은 오류 메시지가 나타납니다.

model.frame.default (delete.response (terms), newdata, na.action = na.action, : factor year has new에 오류가 있습니다. 수준

내가 지정한 수준을 선택하는 대신 몇 년 동안 평균을 낼 수 있습니까? 또한 이러한 예측과 관련된 표준 오류를 얻는 방법을 찾지 못했습니다. (lsmeans 패키지의) 함수를 사용하여 예측에 대한 표준 오류를 얻으려는 경우 :

c <- lsmeans(mm, "TLR4", type="response") 
summary(c, type="response") 

자동으로 표준 오류가 생성됩니다. 그러나 이것은 다른 모든 설명 변수를 평균하여 생성됩니다. 나는 아마도 그것을 바꿀 수 있다고 확신하지만 가능하다면 나는 predict() 함수를 사용하려고합니다. 내 목표는 x 축에 Tlr4 유전자형 그래프를 만들고 y 축에 기생충 부담을 예측하여 다른 모든 중요한 공변량을 고려하면서 각 유전자형에 대한 기생충 부담의 차이를 예측합니다.

Answer 1

반 구체의 데이터 세트를 작성한 다음 새로운 데이터에 대한 예측을 작성하여 변수가 결과에 미치는 실질적인 영향을 조사 할 수있는 두 가지 기능이 포함 된 merTools 패키지에 관심이있을 수 있습니다. 이에 대한 좋은 예가 README 및 패키지 비 네트에서 제공됩니다.

범주와 연속적 예측 자 사이에 상호 작용 용어가있는 모델의 영향을 탐색하고자하는 경우를 생각해 봅시다.

data(VerbAgg) 
fmVA <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 + 
     (1|id) + (1|item), family = binomial, 
     data = VerbAgg) 

이제 우리는 merTools에서 draw 기능을 사용하여 데이터를 수험 : 첫째, 우리는 상호 작용 모델을 맞습니다. 여기서 우리는 모델 프레임으로부터의 평균 관측치를 그립니다.그런 다음 var 매개 변수로 지정된 변수의 값이 다른 동일한 관찰을 반복하도록 데이터 프레임을 확장하여 데이터를 wiggle으로 만듭니다. 여기에서 데이터 세트를 btype, situ 및 Anger의 모든 값으로 확장합니다.

# Select the average case 
newData <- draw(fmVA, type = "average") 
newData <- wiggle(newData, var = "btype", values = unique(VerbAgg$btype)) 
newData <- wiggle(newData, var = "situ", values = unique(VerbAgg$situ)) 
newData <- wiggle(newData, var = "Anger", values = unique(VerbAgg$Anger)) 

head(newData, 10) 

#> r2 Anger Gender btype situ id  item 
#> 1 N 20  F curse other 5 S3WantCurse 
#> 2 N 20  F scold other 5 S3WantCurse 
#> 3 N 20  F shout other 5 S3WantCurse 
#> 4 N 20  F curse self 5 S3WantCurse 
#> 5 N 20  F scold self 5 S3WantCurse 
#> 6 N 20  F shout self 5 S3WantCurse 
#> 7 N 11  F curse other 5 S3WantCurse 
#> 8 N 11  F scold other 5 S3WantCurse 
#> 9 N 11  F shout other 5 S3WantCurse 
#> 10 N 11  F curse self 5 S3WantCurse 

이제 우리는 단순히 이러한 counterfactuals에 대한 예측을 생성하기 위해 predictInterval이 새 데이터 세트를 전달합니다. 그런 다음 두 개의 범주 형 변수 situ 및 btype에 대한 연속 변수 Anger 및 패싯 및 그룹에 대해 예상 값을 각각 그립니다.

plotdf <- predictInterval(fmVA, newdata = newData, type = "probability", 
     stat = "median", n.sims = 1000) 
plotdf <- cbind(plotdf, newData) 
ggplot(plotdf, aes(y = fit, x = Anger, color = btype, group = btype)) + 
    geom_point() + geom_smooth(aes(color = btype), method = "lm") + 
    facet_wrap(~situ) + theme_bw() + 
    labs(y = "Predicted Probability")