预测类型 "response" 的 lme4(二元选择)——置信区间
predict lme4 for type "response" (binary choices) -- Confidence intervals
I 运行 与 lme4 (type="response") 的混合效应 logistig 回归。现在我使用了预测功能并且还想确定置信区间。
我发现这个代码(http://glmm.wikidot.com/faq)用于预测,它有效,但 CI 不适合二进制响应(我的预测在 0 和 1 之间,突然 CI 在 -3 和 3 之间.有谁知道在哪里调整这个?
library(lme4)
library(ggplot2) # Plotting
fm1<- glmer(choice~rating + indi + rating*indi + (1|ID),data=z,family="binomial")
newdat<-data.frame(indi=factor(c(1,1,1,1,1,1,2,2,2,2,2,2)),
rating=factor(1:6), ID=factor(rep(c(1:30), each=12)), choice=0)
newdat$prob<-predict(fm1,newdata=newdat, re.form=NULL, type="response")
mm <- model.matrix(terms(fm1),newdat)
newdat$choice<- predict(fm1,newdat)
## or newdat$choice<- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$ID[1] ## must be adapted
tvar1 <-
newdat <- data.frame(
newdat
, plo = newdat$choice-2*sqrt(pvar1)
, phi = newdat$choice+2*sqrt(pvar1)
, tlo = newdat$choice-2*sqrt(tvar1)
, thi = newdat$choice+2*sqrt(tvar1)
)
#plot confidence
g0 <- ggplot(newdat, aes(x=rating, y=choice, colour=indi))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
opts(title="CI based on FE uncertainty + RE variance")
非常感谢!
您需要使用反函数-link,在这种情况下plogis()
:
newdat <- transform(newdat,
plo = plogis(plo),
phi = plogis(phi),
tlo = plogis(tlo),
thi = plogis(thi))
更一般地说,如果你已经拟合了一个模型gm1
,反link函数将存储在gm1@resp$family$linkinv
中(尽管像这样处理对象的内部结构不是保证与未来版本保持兼容)。
I 运行 与 lme4 (type="response") 的混合效应 logistig 回归。现在我使用了预测功能并且还想确定置信区间。
我发现这个代码(http://glmm.wikidot.com/faq)用于预测,它有效,但 CI 不适合二进制响应(我的预测在 0 和 1 之间,突然 CI 在 -3 和 3 之间.有谁知道在哪里调整这个?
library(lme4)
library(ggplot2) # Plotting
fm1<- glmer(choice~rating + indi + rating*indi + (1|ID),data=z,family="binomial")
newdat<-data.frame(indi=factor(c(1,1,1,1,1,1,2,2,2,2,2,2)),
rating=factor(1:6), ID=factor(rep(c(1:30), each=12)), choice=0)
newdat$prob<-predict(fm1,newdata=newdat, re.form=NULL, type="response")
mm <- model.matrix(terms(fm1),newdat)
newdat$choice<- predict(fm1,newdat)
## or newdat$choice<- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$ID[1] ## must be adapted
tvar1 <-
newdat <- data.frame(
newdat
, plo = newdat$choice-2*sqrt(pvar1)
, phi = newdat$choice+2*sqrt(pvar1)
, tlo = newdat$choice-2*sqrt(tvar1)
, thi = newdat$choice+2*sqrt(tvar1)
)
#plot confidence
g0 <- ggplot(newdat, aes(x=rating, y=choice, colour=indi))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
opts(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
opts(title="CI based on FE uncertainty + RE variance")
非常感谢!
您需要使用反函数-link,在这种情况下plogis()
:
newdat <- transform(newdat,
plo = plogis(plo),
phi = plogis(phi),
tlo = plogis(tlo),
thi = plogis(thi))
更一般地说,如果你已经拟合了一个模型gm1
,反link函数将存储在gm1@resp$family$linkinv
中(尽管像这样处理对象的内部结构不是保证与未来版本保持兼容)。