python 多元回归 python多元回归分类变量

转载

mob64ca13fdd43c 2023-08-24 12:26:37

文章标签 python 多元回归 r library car r语言回归分析分类变量分类变量哑变量矩阵指标矩阵 文章分类 Python 后端开发

在前面几讲，我们介绍了线性回归及R的实现。今天的课程将继续带大家学习多元线性回归。当我们提到“线性”回归时，特指的是因变量(结果变量)为连续性变量，与自变量(预测变量)有线性关系，而对自变量(预测变量)并没有要求一定要是连续性变量。前面我们已经提到，当自变量是连续变量时，线性回归可以写成一个线性方程式y = b0 + b1*x1 + b2*x2 + …

那么，当自变量是分类变量时，回归分析时如何处理的呢？我们能不能把各个分类的类别像血压、血糖数值一样，对应为响应的数值大小来处理呢？

这显然是不可以的，因为分类变量各类别间，可能并没有一个高低之分。今天我们将带你找到答案。

1. 分类变量

分类变量(也称为因子或定性变量)是将观察指标分类的变量。它们具有数量有限的不同值，称为级别。例如，性别是可以分为两个级别的分类变量：男性或女性。

回归分析时需要用数值变量。因此，当需要将分类变量用在回归模型中时，需要对分类变量进行补充处理，以使结果可解释。

通常，我们需要将分类变量进行重新编码，使成为一系列二进制的变量，被称为对比矩阵。这个新的编码，被称为“哑变量”。

2. 加载所需的R包

#tidyverse 便于数据操作和可视化

library(tidyverse)

2.1 数据集示例

我们将使用在car软件包中的Salaries数据集，其中包含2008-09年度某学校助理教授，副教授和教授的9个月工资情况。收集这些数据是该学校行政部门为监测男女教职员工薪资差异而收集的数据。

获取数据

data("Salaries", package = "car")
data("Salaries", package = "car")

检查数据，随机抽取3个样本显示。

sample_n(Salaries, 3)
sample_n(Salaries, 3)

输出结果：

rank discipline yrs.since.phd yrs.service    sex salary115 Prof          A            12           0 Female 105000313 Prof          A            29          19   Male  94350162 Prof          B            26          19   Male 176500

2.2 具有两个级别的分类变量～二分类变量

基于预测变量(x)预测结果变量(y)的回归方程可以写为

y = b0 + b1*x。b0和b1是回归beta系数，分别代表截距和斜率。

假设，我们希望调查男性和女性之间的工资差异。基于性别变量，我们可以创建一个新的虚拟变量，其值如下：

·1 是男性

·0 是女性

并将此变量用作回归方程式的预测变量，从而得出以下模型：

·b0 + b1 如果是男性 (y = b0 + b1 * 1)

·b0 如果是女性(x = b0 + b1 * 0)

系数可以解释如下：1)

b0 是女性的平均工资。2) b0+ b1 是男性的平均工资。3) 而 b1在男性和女性之间的工资的平均差异。 我们以Salaries数据集中男女工资差异为例，进行解说：建立一个简单的线性回归模型来解释男性和女性之间的工资差异。

虽然数据中sex是Male和Female，但是R会自动创建哑变量，于是模型结果如下：

model ~ sex, data = Salaries)summary(model)$coef

输出结果

Estimate Std. Error t value Pr(>|t|)(Intercept)   101002       4809   21.00 2.68e-66sexMale        14088       5065    2.78 5.67e-03

从上面的输出中，女性的平均工资估计为101002，而男性的平均工资估计为101002 + 14088 = 115090。哑变量sexMale的p值 = 5.670e-3，非常显着，表明有统计证据表明性别之间的平均工资存在差异。

使用contrasts()函数可以返回R创建的哑变量编码：

contrasts(Salaries$sex)

输出结果

MaleFemale    0Male      1

R创建了一个sexMale虚拟变量，如果性别为Male，则其值为1，否则为0。将男性编码为1，女性编码为0(基线)的决定是任意的，并且对回归计算没有影响，但会改变系数的解释方向。

您可以使用函数relevel()将基线类别设置为男性，如下所示：

Salaries %mutate(sex = relevel(sex, ref = "Male"))

回归拟合的输出变为：

model ~ sex, data = Salaries)summary(model)$coef

输出结果

Estimate Std. Error t value  Pr(>|t|)(Intercept)   115090       1587   72.50 2.46e-230sexFemale     -14088       5065   -2.78  5.67e-03

sexFemale回归输出中的系数为负。这一事实表明，女性与工资的下降(相对于男性)相关。

现在对b0和b1的估计分别为115090和-14088。解释为：

男性的平均工资为11509，女性的平均预测为115090-14088 = 101002。

此外，除了0/1编码外，我们还可以创建一虚拟变量-1(male)/ 1(female)。结果模型：

·b0 - b1 如果是男性

·b0 + b1 如果是女性

2.3 具有两个以上级别的分类变量 – 多分类变量

通常，将具有n个级别的分类变量转换为n-1个哑变量。这n-1个哑变量包含的信息与单个变量相同。重新编码将创建矩阵，称为对比矩阵。

例如，Salaries数据集中的rank具有三个级别：“ AsstProf”，“ AssocProf”和“ Prof”。可以将该变量编码为两个哑变量，一个称为AssocProf，另一个称为Prof：

·如果rank = AssocProf，则列AssocProf的编码为1，Prof的编码为0。

·如果rank = Prof，则将AssocProf列编码为0，将Prof列编码为1。

·如果rank = AsstProf，则“ AssocProf”和“ Prof”两列都将编码为0。

该变量的编码由R自动执行。您可以使用函数model.matrix()为分类变量创建对比矩阵，进行学习：

res head(res[, -1])

输出结果

rankAssocProf rankProf1             0        12             0        13             0        04             0        15             0        16             1        0

建立线性模型时，有多种方法来编码分类变量，称为对比编码系统。R中的默认选项是使用分类变量的第一个级别作为参考，并用以解释相对于该级别的其余级别，比如，结果解释为：相对于该参考级别，某级别是参考级别的几倍。

请注意，ANOVA(方差分析)是线性模型的一种特殊情况，其中预测变量是分类变量。并且，由于R理解ANOVA和回归都是线性模型的示例，因此您可以使用anova()函数或Anova()函数[car包中] ，从回归模型中提取经典ANOVA表。我们通常建议使用此Anova()功能，因为它会自动处理不平衡设计。

下面显示了使用多元回归预测薪水的结果。

library(car)model2 Anova(model2)

输出结果

Anova Table (Type II tests)Response: salary              Sum Sq  Df F value  Pr(>F)   yrs.service 3.24e+08   1    0.63    0.43   rank        1.03e+11   2  100.26 < 2e-16 ***discipline  1.74e+10   1   33.86 1.2e-08 ***sex         7.77e+08   1    1.51    0.22   Residuals   2.01e+11 391        ---Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

考虑到其他变量(服务年数，职级和学科)，可以看出分类变量性别不再与个人之间的薪资差异显着相关。重要变量是教职等级和学科。

如果要查看详细内容，请输入以下内容：

summary(model2)

输出结果

Call:lm(formula = salary ~ yrs.service + rank + discipline + sex,  data = Salaries)Residuals:Min     1Q Median     3Q    Max-64202 -14255  -1533  10571  99163Coefficients:              Estimate Std. Error t value Pr(>|t|)   (Intercept)    73122.9     3245.3   22.53  < 2e-16 ***yrs.service      -88.8      111.6   -0.80  0.42696   rankAssocProf  14560.4     4098.3    3.55  0.00043 ***rankProf       49159.6     3834.5   12.82  < 2e-16 ***disciplineB    13473.4     2315.5    5.82  1.2e-08 ***sexFemale      -4771.2     3878.0   -1.23  0.21931 ---Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 22700 on 391 degrees of freedomMultiple R-squared:  0.448,  Adjusted R-squared:  0.441F-statistic: 63.4 on 5 and 391 DF,  p-value:<2e-16< span="">

例如，可以看出，与学科A(理论系，基线)相比，来自学科B(应用系)的工资平均增加了13473.38。