R语言和统计：如何画出广义线性模型的预测变量效应图?

kexue

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

如果在实践中被问到，这篇文章的帮助将会是巨大！
首先安装并且载入R包：

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

复制代码

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

1. summary(Arrests)
   
2. str(Arrests)

复制代码

其中的released为一个二分类变量，将作为后续Logistic回归模型中的因变量。关于数据集其他变量，可以使用帮助文档进行详细了解。
下一步，建立一个Logistic回归模型作为例子，released作为因变量，自变量包含checks, employed, 以及两个交互项colour*year以及colour*age，代码如下：

1. mydata <- Arrests
   
2. mymodel <- glm(released ~ checks + employed + colour*year + colour*age,
   
3.                 family = binomial, data = mydata)
   
4. 
   
5. summary(mymodel)

复制代码

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

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

1. checks_effects <- predictorEffect("checks", mymodel)
   
2. checks_effects

复制代码

从代码可知，需要在函数中输入感兴趣的自变量名字，加引号，紧跟着这个模型的名字。这是最简洁的输入代码的方式。
默认情况下， 将连续变量切割成50份（checks的值范围为0-6），数值下方分别展示相应的预测值。

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

1. summary(checks_effects)

复制代码

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

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

1. plot(checks_effects)

复制代码

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

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

复制代码

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

1. # 查看employed的效应
   
2. employed_effects <- predictorEffect("employed", mymodel)
   
3. plot(employed_effects)
   
4. 
   
5. # 查看age的效应
   
6. age_effects <- predictorEffect("age", mymodel)
   
7. plot(age_effects)
   
8. 
   
9. # 查看year的效应
   
10. year_effects <- predictorEffect("year", mymodel)
    
11. plot(year_effects)上述代码输出的图片分别如下：

复制代码

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

1. plot(age_effects,
   
2.      main = "",
   
3.      xlab = "Age",
   
4.      ylab = "Probability (released)",
   
5.      rug = FALSE)

复制代码

第二个例子为线性回归模型。
将使用这个R包自带的数据集，名为Prestige，查看下数据集的概况：

1. mydata1 <- Prestige
   
2. summary(mydata1)
   
3. str(mydata1)

复制代码

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

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

复制代码

与之前的例子不同的是，这次将使用predictorEffects()函数计算预测变量效应值，而不是preditorEffect()，请注意，两者的差别就是一个字母s。
predictorEffects()函数多了一个s，意味着这个函数可以一下子计算出所有预测变量效应值。
下一步，使用这个函数将所有效应图一起画出，代码如下：

1. eall_model1 <- predictorEffects(mymodel1)
   
2. plot(eall_model1)

复制代码

这对于希望快速的了解所有预测变量效应图将会非常有帮助！
但是，如果只想展示某一预测变量，或者希望对单个的图进行修饰，可以这么做。还是使用predictorEffects()函数，代码如下：

1. plot(predictorEffects(mymodel1, ~ education))
   
2. plot(predictorEffects(mymodel1, ~ women))
   
3. plot(predictorEffects(mymodel1, ~ income))
   
4. plot(predictorEffects(mymodel1, ~ type))

复制代码

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

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

复制代码

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

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

▌本文由R语言和统计首发▌