审稿人：你的模型结果汇报的太含糊！请给出预测变量效应图，谢谢！

这篇文章介绍如何画出广义线性模型的预测变量效应图（predictor effect plot；关于这个概念的详细解释，可以查看文献1）[1]，此方法可以清晰的展示自变量与因变量的关系，同时将其他变量考虑在内。

将介绍两种模型的应用，包含Logistic回归以及线性回归模型。

如果在实践中被问到，这篇文章的帮助将会是巨大！

首先安装并且载入R包：

install.packages("effects")
library(effects)
R包get！

将使用这个R包中的自带数据集Arrests，查看下数据集概况：

summary(Arrests)
str(Arrests)


其中的released为一个二分类变量，将作为后续Logistic回归模型中的因变量。关于数据集其他变量，可以使用帮助文档进行详细了解。

下一步，建立一个Logistic回归模型作为例子，released作为因变量，自变量包含checks, employed, 以及两个交互项colour*year以及colour*age，代码如下：

mydata <- Arrests
mymodel <- glm(released ~ checks + employed + colour*year + colour*age, 
                family = binomial, data = mydata)

summary(mymodel)

模型建立好了，那如何查看某一个自变量对因变量的效应呢？

这里以checks变量为例，它是一个连续变量，代码如下：

checks_effects <- predictorEffect("checks", mymodel)
checks_effects

从代码可知，需要在函数中输入感兴趣的自变量名字，加引号，紧跟着这个模型的名字。这是最简洁的输入代码的方式。

默认情况下， 将连续变量切割成50份（checks的值范围为0-6），数值下方分别展示相应的预测值。

还可以分别展示更加详细的信息，比如95%置信区间，代码如下：

summary(checks_effects)

与之前的结果相比，多了两段内容，分别展示了置信区间的两侧数值。

表格或者数值的展示有时并不清晰和直观（如上面的结果），可以选择将上面的数值做成图片，代码如下：

plot(checks_effects)

一行代码出图！

现在回到summary()输出结果，在之前的checks_effects中，默认情况下被分割成了50份。为了让输出结果更加简洁，可以只分割成5份，代码如下：

checks5_effects <- predictorEffect("checks", mymodel,
                                   focal.levels = 5)
summary(checks5_effects)

同理，也可以分别展示其他自变量的效应图，代码如下：

# 查看employed的效应
employed_effects <- predictorEffect("employed", mymodel)
plot(employed_effects)

# 查看age的效应
age_effects <- predictorEffect("age", mymodel)
plot(age_effects)

# 查看year的效应
year_effects <- predictorEffect("year", mymodel)
plot(year_effects)
上述代码输出的图片分别如下： 




继续以age为例子，介绍几个可以修改图片的参数，比如修改x轴标签，题目标签，删除rug plot等，代码如下：

plot(age_effects, 
     main = "",
     xlab = "Age", 
     ylab = "Probability (released)",
     rug = FALSE)

第二个例子为线性回归模型。 

将使用这个R包自带的数据集，名为Prestige，查看下数据集的概况： 

mydata1 <- Prestige
summary(mydata1)
str(mydata1)


其中，将使用prestige作为因变量，它是一个连续变量。再选取几个变量作为自变量纳入方程。因为income的分布偏态非常严重， 因此将其进行log转化。代码如下：

mymodel1 <- lm(prestige ~ education + women + log(income)*type, data = mydata1)
summary(mymodel1)

与之前的例子不同的是，这次将使用predictorEffects()函数计算预测变量效应值，而不是preditorEffect()，请注意，两者的差别就是一个字母s。

predictorEffects()函数多了一个s，意味着这个函数可以一下子计算出所有预测变量效应值。

下一步，使用这个函数将所有效应图一起画出，代码如下： 

eall_model1 <- predictorEffects(mymodel1)
plot(eall_model1)

这对于希望快速的了解所有预测变量效应图将会非常有帮助！

但是，如果只想展示某一预测变量，或者希望对单个的图进行修饰，可以这么做。还是使用predictorEffects()函数，代码如下：

plot(predictorEffects(mymodel1, ~ education))
plot(predictorEffects(mymodel1, ~ women))
plot(predictorEffects(mymodel1, ~ income))
plot(predictorEffects(mymodel1, ~ type))





以income为例，如果不希望将三种type分开画，可以将三张图片合并在一起，以三条线条的形式展示，代码如下：

plot(predictorEffects(mymodel1, ~ income), 
     lines = list(multiline = TRUE))

好啦，今天的内容就到这里。如果有帮助，记得分享给需要的人！

参考文献

[1] Visualizing Fit and Lack of Fit in Complex Regression Models with Predictor Effect Plots and Partial Residuals


公众号的线上课程
1. 《R语言和统计新手课程》
2. 《回归：从入门到进阶》
【通过公众号菜单栏--线上课程--新手课程/回归课程】

统计咨询
《服务介绍和经典合作案例》

公众号核心成员的成果发表
《SCI医学1区影响因子9分论文》

公众号核心成员担任SCI杂志Associate Editor!
《JAD杂志Associate editor》
《Frontiers in Neuroscience, Frontiers in Neurology and Frontiers in Psychiatry杂志的神经退行性病变板块》


▌本文由R语言和统计首发
▌课程相关咨询可添加R师妹微信: kefu_rstats
▌编辑：June
▌邮箱：contact@rstats.cn
▌网站：www.rstats.cn
▌我们致力于让R语言和统计变得简单！


2024-04-15 14:43
首页  ꄲ  公众号  ꄲ  审稿人：你的模型结果汇报的太含糊！请给出预测变量效应图，谢谢！