【发布时间】：2017-08-19 10:08:03
【问题描述】：
我正在使用 lme4 包运行 glmer logit 模型。我对各种二向和三向交互效应及其解释感兴趣。为简化起见，我只关心固定效应系数。
我设法想出了一个代码来计算并在 logit 尺度上绘制这些影响，但我无法将它们转换为预测的概率尺度。最终我想复制effects 包的输出。
该示例依赖于UCLA's data on cancer patients。
library(lme4)
library(ggplot2)
library(plyr)
getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
facmin <- function(n) {
min(as.numeric(levels(n)))
}
facmax <- function(x) {
max(as.numeric(levels(x)))
}
hdp <- read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv")
head(hdp)
hdp <- hdp[complete.cases(hdp),]
hdp <- within(hdp, {
Married <- factor(Married, levels = 0:1, labels = c("no", "yes"))
DID <- factor(DID)
HID <- factor(HID)
CancerStage <- revalue(hdp$CancerStage, c("I"="1", "II"="2", "III"="3", "IV"="4"))
})
到这里为止，就是我需要的所有数据管理、功能和包。
m <- glmer(remission ~ CancerStage*LengthofStay + Experience +
(1 | DID), data = hdp, family = binomial(link="logit"))
summary(m)
这是模型。这需要一分钟，然后会出现以下警告：
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0417259 (tol = 0.001, component 1)
尽管我不太确定是否应该担心警告，但我还是使用估计来绘制感兴趣的交互作用的平均边际效应。首先，我准备要输入 predict 函数的数据集，然后使用固定效应参数计算边际效应和置信区间。
newdat <- expand.grid(
remission = getmode(hdp$remission),
CancerStage = as.factor(seq(facmin(hdp$CancerStage), facmax(hdp$CancerStage),1)),
LengthofStay = seq(min(hdp$LengthofStay, na.rm=T),max(hdp$LengthofStay, na.rm=T),1),
Experience = mean(hdp$Experience, na.rm=T))
mm <- model.matrix(terms(m), newdat)
newdat$remission <- predict(m, newdat, re.form = NA)
pvar1 <- diag(mm %*% tcrossprod(vcov(m), mm))
cmult <- 1.96
## lower and upper CI
newdat <- data.frame(
newdat, plo = newdat$remission - cmult*sqrt(pvar1),
phi = newdat$remission + cmult*sqrt(pvar1))
我相当有信心这些是 logit 量表上的正确估计，但也许我错了。总之，剧情是这样的：
plot_remission <- ggplot(newdat, aes(LengthofStay,
fill=factor(CancerStage), color=factor(CancerStage))) +
geom_ribbon(aes(ymin = plo, ymax = phi), colour=NA, alpha=0.2) +
geom_line(aes(y = remission), size=1.2) +
xlab("Length of Stay") + xlim(c(2, 10)) +
ylab("Probability of Remission") + ylim(c(0.0, 0.5)) +
labs(colour="Cancer Stage", fill="Cancer Stage") +
theme_minimal()
plot_remission
我认为现在 OY 量表是在 logit 量表上测量的，但为了理解它，我想将其转换为预测概率。基于wikipedia，类似exp(value)/(exp(value)+1) 的东西应该可以达到预测概率。虽然我可以做到 newdat$remission <- exp(newdat$remission)/(exp(newdat$remission)+1) 我不确定我应该如何为置信区间做到这一点？
最终我想得到effects 包生成的相同情节。那就是：
eff.m <- effect("CancerStage*LengthofStay", m, KR=T)
eff.m <- as.data.frame(eff.m)
plot_remission2 <- ggplot(eff.m, aes(LengthofStay,
fill=factor(CancerStage), color=factor(CancerStage))) +
geom_ribbon(aes(ymin = lower, ymax = upper), colour=NA, alpha=0.2) +
geom_line(aes(y = fit), size=1.2) +
xlab("Length of Stay") + xlim(c(2, 10)) +
ylab("Probability of Remission") + ylim(c(0.0, 0.5)) +
labs(colour="Cancer Stage", fill="Cancer Stage") +
theme_minimal()
plot_remission2
尽管我可以只使用 effects 包，但遗憾的是它无法与我必须为自己的工作运行的许多模型一起编译：
Error in model.matrix(mod2) %*% mod2$coefficients :
non-conformable arguments
In addition: Warning message:
In vcov.merMod(mod) :
variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX
解决这个问题需要调整估算程序，目前我想避免这种情况。另外，我也很好奇 effects 在这里实际做了什么。
如果有任何关于如何调整我的初始语法以达到预测概率的建议，我将不胜感激！
【问题讨论】：
-
我认为如果你这样做，你的情节会更容易阅读：
ggplot(newdat, aes(LengthofStay, fill=factor(CancerStage), color=factor(CancerStage))) + geom_ribbon(aes(ymin=plo, ymax=phi), colour=NA, alpha=0.2) + geom_line(aes(y = remission), size=1.2) + xlab("Length of Stay") + ylab("Probability of Remission") + labs(colour="Cancer Stage", fill="Cancer Stage") + theme_minimal() -
你绝对应该担心收敛警告。
-
我真的不明白为什么这是一个不可能回答的问题......我要求的内容是否不清楚？
-
我同意@JacobSocolar 的观点。我认为您的模型不收敛的事实将导致虚假的模型估计。所以要小心。
-
当然，谢谢！但这是一个相当侧面的观点。如何使用基于
predict的初始语法来绘制反映预测概率的图？
标签： r graphics predict lme4 marginal-effects