有序回归

logistic回归

• 二元logistic回归：Y为定类且为2个，比如是否购买(1购买；0不购买)
• 多分类logistic回归：Y为定类且选项大于2个，比如总统候选人偏好(特朗普、希拉里、卢比奥)
• 有序logistic回归：Y为定类且有序，幸福感(不幸福、比较幸福和非常幸福)

生活中的有序logistic回归

• 人们在肯德基里点餐，一般都会买可乐，可乐有四种型号(small, medium, large or extra large)，选择何种型号的可乐会与哪些因素有关呢？是否购买了汉堡、是否购买了薯条，消费者的年龄等。我们这里考察的被解释变量，可乐的大小就是一个有序的值。
• 问卷调查。问大三的学生是否申请读研究生，有三个选项：1不愿意，2有点愿意，3非常愿意。那么这里的被解释变量是三个有序的类别，影响读研意愿的因素可能与父母的教育水平、本科阶段学习成绩、经济压力等有关。

案例

教育代际传递。通俗点说子女的教育程度是否受到父母教育程度的影响。我这个案例思路参考了南京大学池彪的《教育人力资本的代际传递研究》硕士论文，这篇文章思路很清晰，建议大家可以看看。根据文中提供的数据来源，我们下载2016年的中国家庭追踪调查数据CFPS，并整理了部分数据。

tb <- readr::read_rds("./demo_data/cfps.rds")
head(tb)

tb %>% count(edu)
tb %>% count(edu_f)
tb %>% count(edu_m)

为了方便处理，减少分类，我们将大专以及大专以上的都归为一类

df <- tb %>%
  dplyr::mutate(
    across(
      starts_with("edu"),
      ~ case_when(
        . %in% c(5, 6, 7, 8) ~ 5,
        TRUE ~ .
      )
    )
  )
df

看起很复杂？那我写简单点

tb %>%
  dplyr::mutate(
    across(
      starts_with("edu"),
      ~ if_else(. %in% c(5, 6, 7, 8), 5, .)
    )
  )

df %>% count(edu)
df %>% count(edu_f)
df %>% count(edu_m)

问题的提出

问题的提出：

• 学历上父母是否门当户对？
• 父母的受教育程度对子女的受教育水平是正向影响？
• 父亲和母亲谁的影响大？
• 对男孩影响大？还是对女孩影响大？
• 以上情况城乡有无差异？

父母门当户对？

数据还是比较有意思的，我们来看看父母是否门当户对

多大比例选择门当户对?

df

df %>%
  dplyr::summarise(
    eq_n = sum(edu_m == edu_f),
    n = n()
  ) %>%
  dplyr::mutate(prop = eq_n / n)

df %>%
  dplyr::count(edu_m, edu_f) %>%
  dplyr::group_by(edu_m) %>%
  dplyr::mutate(prop = n / sum(n)) %>%
  dplyr::ungroup()

df %>%
  dplyr::count(edu_m, edu_f) %>%
  dplyr::group_by(edu_m) %>%
  dplyr::mutate(percent = n / sum(n)) %>%
  dplyr::select(-n) %>%
  tidyr::pivot_wider(
    names_from = edu_f,
    values_from = percent
  )

母亲的教育程度对子女的影响

library(ggridges)
df %>%
  dplyr::mutate(
    across(edu_m, as.factor)
  ) %>%
  ggplot(aes(x = edu, y = edu_m)) +
  geom_density_ridges() +
  scale_x_continuous(limits = c(0, 6), breaks = c(1:5)) +
  labs(
    title = "The influence of mother's education on children in the family",
    subtitle = "The greater the number, the higher the level of education",
    x = "Education level of children",
    y = "Education level of Mother"
  )

父亲和母亲谁的影响大

这里需要用到有序logistic回归。 为了理解模型的输出，我们需要先简单介绍下模型的含义。假定被解释变量$Y$有$J$类且有序，那么$Y$ 小于等于某个具体类别$j$的累积概率，可以写为$P(Y \le j)$，这里$j = 1, \cdots, J-1$.

从而，小于等于某个具体类别$j$的比率就可以定义为

$$\frac{P(Y \le j)}{P(Y>j)}$$ 对这个比率取对数，就是我们熟知的logit

$$log \frac{P(Y \le j)}{P(Y>j)} = logit (P(Y \le j)).$$

，有序logistic回归的数学模型就是

$$logit (P(Y \le j)) = \alpha_{j} + \beta_{1}x_1 + \beta_{2}x_2 $$ $\alpha$ 是截距 $\beta$ 是回归系数，注意到有序分类 logistic 回归模型中就有 $J-1$ 个 logit 模型。对于每个模型，系数是相同的，截距不同。

在R语言中，我们可以使用MASS::polr函数，但需要注意的是，使用这个函数，对应的模型表达式为^[感谢@huanfachen指出我之前的错误]，斜率符号写为负号

$$logit (P(Y \le j)) = \alpha_{j} - \beta_{1}x_1 - \beta_{2}x_2 $$

下面我们通过代码来演示

library(MASS)
df1 <- df %>%
  dplyr::mutate(
    across(c(edu, sex, urban), as.factor),
    across(edu, ~ fct_inseq(., ordered = TRUE))
  )

mod_mass <- polr(edu ~ edu_f + edu_m + sex + num_siblings + urban,
  data = df1,
  method = c("logistic")
)

summary(mod_mass)

输出结果得到有序分类 logistic 回归模型中截距和回归系数的最大似然估计值，确定出回归方程为：

$$ \begin{aligned} \text{logit}(\hat{P}(Y \le 1))&= \text{logit}\left(p_{1}\right) = \ln \left(\frac{p_{1}}{1 - p_{1}}\right) = -0.8385 - 0.46 \times \text{edu_f} - 0.51 \times\text{edu_m} -(-0.46)\times\text{sex1} -(-0.15)\times \text{num_siblings} -0.96 \times\text{urban1} \\ \text{logit}(\hat{P}(Y \le 2))&= \text{logit}\left(p_{1} + p_{2}\right) = \ln \left(\frac{p_{1} + p_{2}}{1 - p_{1} - p_{2}}\right) = 0.6742 - 0.46 \times \text{edu_f} - 0.51 \times\text{edu_m} -(-0.46)\times\text{sex1} -(-0.15)\times \text{num_siblings} -0.96 \times\text{urban1} \\ \text{logit}(\hat{P}(Y \le 3))&= \text{logit}\left(p_{1} + p_{2} + p_{3}\right) = \ln \left(\frac{p_{1} + p_{2} + p_{3}}{1 - p_{1} - p_{2} - p_{3}}\right) = 2.5093 - 0.46 \times \text{edu_f} - 0.51 \times\text{edu_m} -(-0.46)\times\text{sex1} -(-0.15)\times \text{num_siblings} -0.96 \times\text{urban1}\\ \text{logit}(\hat{P}(Y \le 4))&= \text{logit}\left(p_{1} + p_{2} + p_{3} + p_{4}\right) = \ln \left(\frac{p_{1} + p_{2} + p_{3} + p_{4}}{1 - p_{1} - p_{2} - p_{3} - p_{4}}\right) = 3.5454 - 0.46 \times \text{edu_f} - 0.51 \times\text{edu_m} -(-0.46)\times\text{sex1} -(-0.15)\times \text{num_siblings} -0.96 \times\text{urban1}\\ \end{aligned} $$

系数的解释

关于系数的解释，推荐您阅读这里。

先将系数转换成odds ratios(OR)

coef(mod_mass) %>% exp()

• 在其它因素不变的情况下，父亲教育程度每增加一个等级（比如，大专到本科）， 会增加子女教育程度向上提高一个级别的概率1.58倍，也就是增加了58%。
• 在其它因素不变的情况下，母亲教育程度每提高一个等级，会增加提升子女教育水平的概率1.66倍.
• 从子女的性别差异来看, 在其它因素不变的情况下，女性的受教育水平向上提高一个级别的概率更大，是男性的(1/0.630)倍，或者说，男性受教育水平向上提高一个级别的概率比女性减少37%(1 - 0.63).
• 从城乡差异来看，城镇子女提升教育水平的概率是农村的2.6倍

边际效应

library(margins)
# me_mass <- marginal_effects(mod_mass, variables = "sex")
me_mass <- marginal_effects(mod_mass, variables = "edu_m")
me_mass %>% 
  head()

从边际效应图可以看到，随着父母教育程度的增加，子女低学历的的概率减少，高学历的概率增加

其他宏包

ordinal 包

library(ordinal)
mod_ordinal <- clm(edu ~ edu_f + edu_m + sex + num_siblings + urban,
  data = df1,
  link = "logit",
  thresholds = "flexible"
)

broom::tidy(mod_ordinal)

# remove the objects
# rm(list=ls())
rm(df, df1, me_mass, mod_mass, mod_ordinal, tb)

pacman::p_unload(pacman::p_loaded(), character.only = TRUE)