这篇文章介绍如何画出广义线性模型的预测变量效应图（predictor effect plot；关于这个概念的详细解释，可以查看文献1）[1]，此方法可以清晰的展示自变量与因变量的关系，同时将其他变量考虑在内。将介绍两种模型的应用，包含Logistic回归以及线性回归模型。
install.packages("effects")
library(effects)
将使用这个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中，默认情况下被分割成了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
【通过公众号菜单栏--线上课程--新手课程/回归课程】公众号核心成员担任SCI杂志Associate Editor!▌课程相关咨询可添加R师妹微信: kefu_rstats
