---
title: "使用 Tidyverse 清洗 CFPS 数据（四）"
author: "Rui"
date: "2022-09-27"
categories: [R, 数据处理, 一些尝试, tidyverse, plm]
image: "yellow_ice.jpg"
format:
html:
code-fold: true
code-tools: true
---
```{r setup, include = FALSE}
# 设置默认参数
knitr::opts_chunk$set(
echo = TRUE,
fig.align = "center",
message = FALSE,
warning = FALSE,
collapse = TRUE
)
```
## 导入数据
```{r}
#| column: page
## 导入数据
library(tidyverse)
file_path <- "data/operated_data/clean_data.sav"
clean_data <- haven::read_sav(file_path) %>% as_tibble()
clean_data %>% head(10) %>% knitr::kable()
```
## 数据转换
将变量转换为论文需要的变量，其中 +1 是为了防止变量为 0 的情况从而导致无法转换。
```{r}
trans_data <- clean_data %>%
mutate(lnexp16 = log(expense16 + 1),
lninc16 = log(fincome16 + 1),
lnasset16 = log(family_asset16 + 1),
lev16 = family_debts16/(family_asset16 + 1),
lnexp18 = log(expense18 + 1),
lninc18 = log(fincome18 + 1),
lnasset18 = log(family_asset18 + 1),
lev18 = family_debts18/(family_asset18 + 1),
year16 = 2016,
year18 = 2018) %>%
select(-c(expense16, fincome16, family_asset16, family_debts16,
expense18, fincome18, family_asset18, family_debts18))
```
## 剔除 lev 变量的极端值
```{r}
library(magrittr)
trans_data %<>% filter(lev16 < 50 & lev18 < 50)
```
## 拆分并对列重命名
上一节中将两个年份的数据合并是为了方便同时对两年的数据进行处理，在这里需要再将数据按照年份拆分开。
```{r}
temp16 <- trans_data %>%
select(ends_with("16")) %>%
cbind(trans_data$id, .) %>%
as_tibble() %>%
rename(year = year16,
id = "trans_data$id",
provcd = provcd16,
urban = urban16,
gender = gender16,
marriage = marriage16,
age = age16,
age2 = age2_16,
health = health16,
fml = fml16,
lnexp = lnexp16,
lev = lev16,
lninc = lninc16,
lnasset = lnasset16)
temp18 <- trans_data %>%
select(ends_with("18")) %>%
cbind(trans_data$id, .) %>%
as_tibble() %>%
rename(year = year18,
id = "trans_data$id",
provcd = provcd18,
urban = urban18,
gender = gender18,
marriage = marriage18,
age = age18,
age2 = age2_18,
health = health18,
fml = fml18,
lnexp = lnexp18,
lev = lev18,
lninc = lninc18,
lnasset = lnasset18)
```
查看数据：
```{r}
temp16 %>% head(10) %>% knitr::kable()
temp18 %>% head(10) %>% knitr::kable()
```
## 论文中的数据
```{r}
## 论文中的数据
delta_lnexp <- temp18$lnexp - temp16$lnexp
delta_lninc <- temp18$lninc - temp16$lninc
delta_lnasset <- temp18$lnasset - temp16$lnasset
delta_fml <- temp18$fml - temp16$fml
data <- data.frame(delta_lnexp = delta_lnexp,
delta_lninc = delta_lninc,
delta_lnasset = delta_lnasset,
delta_fml) %>%
cbind(temp16[c(13, 2:3, 5:12)]) %>%
as_tibble()
```
## 数据类型规整
```{r}
data$provcd %<>% as.factor()
data$urban %<>% as.factor()
data$gender %<>% as.factor()
data$marriage %<>% as.factor()
data$health %<>% as.factor()
data$delta_fml %<>% as.integer()
data$age %<>% as.integer()
data$age2 %<>% as.integer()
```
## 建立模型
```{r}
## 多元线性回归
model1 <- lm(delta_lnexp ~
lev + delta_lninc + lnasset +
provcd + urban + age + age2 +
gender + marriage + health + delta_fml,
data = data)
library(plm)
library(lmtest)
coeftest(model1, vcov = vcovHC, type = "HC1")
```
## 问题
可能存在的问题：
- 数据中 lev 变量大量为 0，原数据中用于计算 lev 的 fameliy_debts 也大量为 0。可能是漏了计算家庭负债变量的相关变量。
- 论文中使用的是 2010 和 2012 年的 CFPS 数据。
- 论文中控制了其他与家庭相关的协变量，如家庭中未婚男性个数等。而本文没有纳入模型中加以控制。
- 论文中将 health 变量粗分为 3 级，本文未作处理。
## 构建面板数据
```{r}
## 面板数据
panel_data <- rbind(temp16, temp18)
## 数据类型规整
panel_data$provcd %<>% as.factor()
panel_data$urban %<>% as.factor()
panel_data$gender %<>% as.factor()
panel_data$marriage %<>% as.factor()
panel_data$health %<>% as.factor()
## 完成
panel_data %<>% select(id, year, everything()) %>% arrange(id, year)
# 先按id排序，再按year排序
panel_data %>% head(10) %>% knitr::kable()
```
## 保存面板数据
```{r}
haven::write_sav(panel_data, "data/operated_data/两年面板数据.sav")
write_excel_csv(panel_data, "data/operated_data/两年面板数据.csv")
```
## 其他尝试
论文中使用的模型相当于一个截面数据模型（具体原因参见论文）。而我们现在拿到了面板数据，可以尝试使用面板数据建立双向固定效应模型。
```{r}
coplot(lnexp ~ year|provcd, type="l", data=panel_data) # Lines
```
```{r}
library(car)
scatterplot(lnexp ~ year|provcd, boxplots=FALSE, smooth=TRUE, reg.line=FALSE, panel_data)
```
```{r}
library(gplots)
plotmeans(lnexp ~ provcd, main="Heterogeineity across provinces", data=panel_data)
plotmeans(lnexp ~ year, main="Heterogeineity across years", data=panel_data)
```
### 普通最小二乘回归
不考虑个体的异质性，仅使用普通最小二乘回归
```{r}
ols <-lm(lnexp ~ lev + lninc + lnasset, data=panel_data)
summary(ols)
```
```{r}
yhat <- ols$fitted # 获取拟合值
plot(panel_data$lev, panel_data$lnexp, pch=19, xlab="lev", ylab="lnexp")
abline(lm(panel_data$lnexp ~ panel_data$lev), lwd=3, col="red")
```
### LSDV
先不考虑个人，将省份作为个体。
用 LSDV 模型控制省份固定效应（其实就是将因子型变量 provcd 转换为多个 0-1 虚拟变量）：
```{r}
fixed.dum <-lm(lnexp ~ lev + lninc + lnasset + factor(provcd) - 1, data=panel_data)
summary(fixed.dum)
```
```{r}
yhat <- fixed.dum$fitted # 获取拟合值
scatterplot(yhat ~ panel_data$lev|panel_data$provcd, boxplots=FALSE, xlab="lev", ylab="yhat",smooth=FALSE)
abline(lm(panel_data$lnexp ~ panel_data$lev),lwd=3, col="red") # 普通最小二乘回归作为对照
```
### 使用 FE
首先先看单向的个体固定效应，即考虑个体差异忽视时间差异。将每个人视为一个个体，控制个体固定效应。个体非常多，无法像 LSDV 模型那样设置 0-1 虚拟变量。
```{r}
library(plm)
fixed1 <- plm(lnexp ~ lev + lninc + lnasset, data=panel_data, index=c("id", "year"), effect = "individual", model="within")
summary(fixed1)
```
使用双向固定效应来同时控制个体与时间：
```{r}
fixed2 <- plm(lnexp ~ lev + lninc + lnasset, data=panel_data, index=c("id", "year"), effect = "twoways", model="within")
summary(fixed2)
```
加入控制变量：
```{r}
fixed3 <- plm(
lnexp ~ lev + lninc + lnasset + provcd + urban + fml + age + age2 + gender + marriage + health,
data = panel_data,
index=c("id", "year"),
effect = "twoways",
model="within"
)
summary(fixed3)
```
```{r}
library(lmtest)
coeftest(fixed3, vcov = vcovHC, type = "HC1")
```
使用 `stargazer` 来输出回归结果
```{r}
library(stargazer)
output <- stargazer(fixed2, fixed3, type="latex", title = "双向固定效应")
```
更多关于 R 语言创建面板模型的教程参考：<https://blog.csdn.net/weixin_43718786/article/details/112854803?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-1-112854803-blog-109094100.235%5Ev31%5Epc_relevant_default_base3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-1-112854803-blog-109094100.235%5Ev31%5Epc_relevant_default_base3&utm_relevant_index=2>
