关键词：
【中文标题】glmer logit - 概率尺度上的交互效应（用 `predict` 复制`effects`）【英文标题】：glmer logit - interaction effects on probability scale (replicating `effects` with `predict`) 【发布时间】：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 的初始语法来绘制反映预测概率的图？
【参考方案1】：
要获得与您的问题中提供的effect 函数类似的结果，您只需使用您提供的转换将预测值和置信区间的边界从 logit 标度反向转换为原始标度： exp(x)/(1+exp(x))。
这个转换可以用plogis函数在base R中完成：
> a <- 1:5
> plogis(a)
[1] 0.7310586 0.8807971 0.9525741 0.9820138 0.9933071
> exp(a)/(1+exp(a))
[1] 0.7310586 0.8807971 0.9525741 0.9820138 0.9933071
所以使用来自@eipi10 的建议，使用丝带代替虚线（我也发现这个演示文稿更具可读性）：
ggplot(newdat, aes(LengthofStay, fill=factor(CancerStage), color=factor(CancerStage))) +
geom_ribbon(aes(ymin = plogis(plo), ymax = plogis(phi)), colour=NA, alpha=0.2) +
geom_line(aes(y = plogis(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()
结果是一样的（effects_3.1-2和lme4_1.1-13）：
> compare <- merge(newdat, eff.m)
> compare[, c("remission", "plo", "phi")] <-
+ sapply(compare[, c("remission", "plo", "phi")], plogis)
> head(compare)
CancerStage LengthofStay remission Experience plo phi fit se lower upper
1 1 10 0.20657613 17.64129 0.12473504 0.3223392 0.20657613 0.3074726 0.12473625 0.3223368
2 1 2 0.35920425 17.64129 0.27570456 0.4522040 0.35920425 0.1974744 0.27570598 0.4522022
3 1 4 0.31636299 17.64129 0.26572506 0.3717650 0.31636299 0.1254513 0.26572595 0.3717639
4 1 6 0.27642711 17.64129 0.22800277 0.3307300 0.27642711 0.1313108 0.22800360 0.3307290
5 1 8 0.23976445 17.64129 0.17324422 0.3218821 0.23976445 0.2085896 0.17324530 0.3218805
6 2 10 0.09957493 17.64129 0.06218598 0.1557113 0.09957493 0.2609519 0.06218653 0.1557101
> compare$remission-compare$fit
[1] 8.604228e-16 1.221245e-15 1.165734e-15 1.054712e-15 9.714451e-16 4.718448e-16 1.221245e-15 1.054712e-15 8.326673e-16
[10] 6.383782e-16 4.163336e-16 7.494005e-16 6.383782e-16 5.689893e-16 4.857226e-16 2.567391e-16 1.075529e-16 1.318390e-16
[19] 1.665335e-16 2.081668e-16
置信边界之间的差异更高但仍然很小：
> compare$plo-compare$lower
[1] -1.208997e-06 -1.420235e-06 -8.815678e-07 -8.324261e-07 -1.076016e-06 -5.481007e-07 -1.429258e-06 -8.133438e-07 -5.648821e-07
[10] -5.806940e-07 -5.364281e-07 -1.004792e-06 -6.314904e-07 -4.007381e-07 -4.847205e-07 -3.474783e-07 -1.398476e-07 -1.679746e-07
[19] -1.476577e-07 -2.332091e-07
但是，如果我使用正态分布的实分位数cmult <- qnorm(0.975) 而不是cmult <- 1.96，我也会在这些边界上获得非常小的差异：
> compare$plo-compare$lower
[1] 5.828671e-16 9.992007e-16 9.992007e-16 9.436896e-16 7.771561e-16 3.053113e-16 9.992007e-16 8.604228e-16 6.938894e-16
[10] 5.134781e-16 2.289835e-16 4.718448e-16 4.857226e-16 4.440892e-16 3.469447e-16 1.006140e-16 3.382711e-17 6.765422e-17
[19] 1.214306e-16 1.283695e-16
【讨论】：
谢谢！这很有帮助！不幸的是，尽管这两个图之间仍然存在细微差别，但我将它们带到了相同的比例，因此它在曲线中可见（我添加了xlim 和ylim）。您还可以看到差异，例如compare <- merge(newdat, eff.m) head(compare) compare$remission-compare$fit确实，在这个例子中差异非常小，但我想了解偏差来自哪里，所以我可以在我的研究中消除它。 PS：我编辑了这些图并添加了plyr 包。感谢您的回答！
查看编辑后的回复。我无法复制任何显着差异。也许软件包版本有所不同？注意，您还应该在代码中添加library(effects) 并删除您的第一个图的ylim（此图在 logit 标度上，因此 0,0.5 限制超出了图的范围）移动通信笔记-小尺度衰落和多径效应
移动无线电传播：小尺度衰落和多径效应移动无线电传播：小尺度衰落和多径效应 查看详情
如何绘制新数据的预测与 R 中的 gee、lme、glmer 和 gamm4 相匹配？
...行比较。我使用geepack拟合GEE模型，使用lme(nlme)在log(count)上的线性混合效应模型，使用g 查看详情
如何使用超过 5000 万个观测值的样本计算具有固定效应的 logit 模型的边际效应
】如何使用超过5000万个观测值的样本计算具有固定效应的logit模型的边际效应【英文标题】：Howtocalculatemarginaleffectsoflogitmodelwithfixedeffectsbyusingasampleofmorethan50millionobservations【发布时间】：2022-01-1116:02:27【问题描述】：我有超过500... 查看详情
固定效应 logit：R 中调整的 r square-bife 包
】固定效应logit：R中调整的rsquare-bife包【英文标题】：Fixedeffectlogit:adjustedrsquare-bifepackageinR【发布时间】：2020-09-2413:22:38【问题描述】：我正在使用R中的bife包研究我的固定效应logit模型。但问题是我需要调整后的r平方，目前我... 查看详情
物体检测的尺度效应实验
...下：如果把图切成两块方形的：可以看出检测效果和图像尺度还是有关系的。因为不管图像大小网络首先将图像缩放到320*320，第一张是长方形，所以缩放后变形大，识别的精度也就低。把SUV识别为truck，右上角白车的遮挡导致识... 查看详情
如何将 SVM 类概率转换为 logits？
】如何将SVM类概率转换为logits？【英文标题】：HowcanconvertSVMclassprobabilitiestologits?【发布时间】：2019-02-2014:17:51【问题描述】：我想将SVM输出的概率类转换为logits。为了得到每个类的概率model=svm.SVC(probability=True)model.fit(X,Y)results=mod... 查看详情
[培训-无线通信基础-3]：窄带无线信道（大小尺度衰落多普勒效应）
...cle/details/118719351目录引言：第1部分无线信道概述与大尺度衰落1.1 无线信道概述1.2大尺度衰落的成因与数学模型（正态分布）第 查看详情
为啥在 sotfmax_cross_entropy_with_logits 中将 logit 解释为“未缩放的对数概率”？
...sotfmax_cross_entropy_with_logits中将logit解释为“未缩放的对数概率”？【英文标题】：whyexplainlogitas\'unscaledlogprobabililty\'insotfmax_cross_entropy_with_logits?为什么在sotfmax_cross_entropy_with_logits中将logit解释为“未缩放的对数概率”？【发布时间... 查看详情
Keras - 如何获得非标准化的logits而不是概率
】Keras-如何获得非标准化的logits而不是概率【英文标题】：Keras-howtogetunnormalizedlogitsinsteadofprobabilities【发布时间】：2018-04-1216:21:34【问题描述】：我正在Keras中创建一个模型，并想计算我自己的指标（困惑度）。这需要使用非标... 查看详情
glmer - 使用二项式数据进行预测（cbind 计数数据）
...在我的二项式数据上运行的glmer模型随时间变化的值（x轴上的天数）。TotalAlive和TotalDead是计数数据。这是我的模型，下面是相应的步骤。full.model.dredge<-g 查看详情
spss广义线性混合效应模型中的随机效应怎么交互
...效应来描述数据的个体差异或重复度量数据，并通过添加交互项来探讨随机效应之间的关系。一般来说，在模型中包含多个随机效应时，应谨慎考虑它们之间的交互作用，并进行适当的解释，以便评估其对数据的影响。 查看详情
如果存在主效应和交互效应，解释整体效应？
】如果存在主效应和交互效应，解释整体效应？【英文标题】：Interpretingoveralleffectifmainandinteractioneffectsarepresent?【发布时间】：2021-10-3021:05:49【问题描述】：假设，我有三个独立的分类变量e、f和g，并且想估计因变量y。经过一... 查看详情
无线信道传输特性
...性是瞬息万变的，通常认为无线信道衰落可以分为大尺度衰落和小尺度衰落，又称为快衰落和慢衰落。小尺度衰落传统上从两个角度来描述其信道特性:即时域上对应的多径效应;频率域上对应的多普勒频移效应引发的频率... 查看详情
获取 glmer 模型的标准化系数？
】获取glmer模型的标准化系数？【英文标题】：Gettingstandardizedcoefficientsforaglmermodel?【发布时间】：2021-01-1211:49:21【问题描述】：有人要求我为glmer模型提供标准化系数，但我不确定如何获得它们。不幸的是，beta函数不适用于glmer... 查看详情
混杂中介效应调节效应
...度和车祸发生，即使相同速度，喝酒比不喝酒发生车祸的概率要大。调节效应多用分层来解决。 查看详情
glmer 的 MUMIn pdredge 错误
】glmer的MUMInpdredge错误【英文标题】：MUMInpdredgeerrorwithglmer【发布时间】：2015-07-2711:52:00【问题描述】：我正在尝试使用dredge的并行计算版本（包MUMIn）来选择我的完整glmer模型的模型：modmer.pom.full<-glmer(cbind(TEST,CONTROL)~G+MS+l+MS*l+... 查看详情
图像特征提取方法
...征2、对边缘光滑的图像难以准确提取特征点原理：1、在尺度空间（例如高斯金字塔）上搜寻keypoints兴趣点（对于尺度和旋转不变）2、筛选上一步获得的兴趣点（1）对空间中的极值点进行精确定位（2）用Hessian矩阵消除边缘效应... 查看详情
tf-logtis模型
...发生与该事件不发生的比值的对数。假设一个事件发生的概率为p，那么该事件的logits为.现在来看一下这个式子和softmax有啥关系。在dl中，softmax层会对输入进行归一化处理以得到概率分布：如下面式子所述就是tensorflow所称的logit... 查看详情