4 描述关系
4.1 何为变量关系？
在大多数研究情境下，学者关注的焦点并非单一变量的分布特征1，而是两个或多个变量之间呈现的关联模式。
两个变量存在关联关系，本质上意味着：通过其中一个变量的信息可推测另一个变量的特征。
以儿童身高与年龄为例：通常年龄越大，身高越高。因此，若已知一名儿童13岁而另一名6岁，即可较准确地推测两者的相对身高。
我们可以将身高与年龄的关系定义为正相关——即一个变量的取值越高，另一个变量的预期取值也越高（年龄增长伴随身高增加）。与之相对的负相关则表现为一个变量取值越高，另一变量取值越低（年龄增长伴随哭闹减少）。此外还存在零相关，即变量间无任何关联（儿童是否居住在法国与年龄无关）。
其他复杂关系模式包括：
时变相关：某些时段呈正相关，其他时段呈负相关
强度变化：初期高度正相关，后期相关性减弱
分类变量相关：当涉及分类变量时（如交通工具选择），关系表现为组间差异而非数值高低（年长儿童更倾向骑自行车出行）
本章的核心目标在于建立描述双变量关系的方法论体系，从而精准呈现研究问题中的数据关联特征——需再次强调，研究问题的本质通常正蕴含于变量关系之中。掌握变量关系的描述方法后，我们便能在后续章节中验证所发现的关联模式是否真正解答了研究问题。
本章将采用\(Emily \space Oster\)基于美国国家健康与营养调查\(NHANES\)的研究数据作为范例2。其核心研究问题是：药物推荐的健康效益是否存在测量偏误——即健康状况更佳的人群本身就更可能遵循用药建议，从而导致效益被高估？
为验证这一假设，\(Oster\)选取了维生素\(E\)补充剂作为研究对象——该营养素仅在特定时段被列为推荐用药。她通过分析以下变量间的动态关联来解答研究问题：(1)服用维生素E的行为；(2)其他健康关注指标（如戒烟）；(3)死亡率等健康结果；(4)这些关联在维生素\(E\)推荐期前、中、后的变化模式3。
我们可以从一个非常直观的方法开始，展示两个连续变量之间的关系——这就是散点图，如图4.1所示。散点图简单地展示了所有的数据点。这种图表非常便于我们全面观察数据，并从中直观地看出两个变量之间存在什么样的关系。数据整体是向上倾斜的吗？还是大幅向下倾斜？或者像图4.1中那样只是轻微向下倾斜？又或者是上下波动的？
散点图是展示两个连续变量之间关系的常用方法，正如第3章中密度图用于展示单一连续变量一样4。它通常是我们描述变量关系的最佳起点。
Figure 4.1: Age and Heart Health, 150 Observations
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(vtable)## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(purrr)
library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(Cairo)
library(extrafont)## Registering fonts with Rlibrary(haven)
library(mvtnorm)
library(ggpubr)##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:cowplot':
##
## get_legendlibrary(modelsummary)oster <- read_dta('DescribingRelationships/nhanes_summary_cleaned.dta') %>%
dplyr::mutate(supplement_vite_single = case_when(
!supplement_vite_single ~ 'No Vitamin E',
TRUE ~ 'Took Vitamin E'
))
ggplot2::ggplot(oster %>% slice(150:300), aes(x=age, y=heart_health)) +
ggplot2::geom_point(na.rm=TRUE) +
ggplot2::labs(x='Age', y='Heart Health Score') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=12, family="sans"),
axis.title.x = element_text(size=12, family="sans"),
axis.title.y = element_text(size=12, family= "sans"))
ggplot2::ggsave('DescribingRelationships/scatterplot.pdf',width=7, height=6, units='in', device=cairo_pdf)set.seed(20200702)
s_x <- 10000
s_y <- 1.5
rho <- -.5
sigma <- base::matrix(c(s_x^2, s_x*s_y*rho, s_x*s_y*rho, s_y^2), 2)
dat <- dplyr::as_tibble(rmvnorm(n = 500, mean = c(50000, 15), sigma = sigma)) %>%
dplyr::rename(Income = V1, Depression = V2)## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ggplot2::ggplot(dat, aes(x=Income, y=Depression)) +
ggplot2::geom_jitter(alpha = .5) +
ggplot2::scale_x_continuous(labels = scales::dollar) +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"))
ggplot2::ggsave('DescribingRelationships/depression_income.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)散点图除了直观展示数据外，还隐含两层含义：一个缺点，一个优点。缺点是人们很容易从散点图中看到变量关系就误认为\(X\)轴导致了\(Y\)轴的变化，即使我们知道这种因果关系并不成立，这种诱惑依然强烈。优点则是散点图能促使我们思考其他描述变量关系的方法，从而以更易理解的方式获取所需信息。本章后续内容将围绕这一点展开。
4.2 条件分布
第三章全面探讨了变量的分布特征，但所述内容均属于无条件分布范畴5。
条件分布是指在给定另一个变量取值的情况下，某个变量的分布 。
我们可以从一个更基础的概念开始讨论——条件概率。一个人是女性的概率大约是50%。但对于名叫Sarah的人来说，是女性的概率则远高于50%。也可以这样表述：“在所有叫Sarah的人中，女性占多大比例？”这种情况我们称之为”在名字是Sarah的条件下，某人是女性的概率”。
得知某人名叫Sarah会改变我们对其性别为女性的概率判断。条件分布的工作原理与此类似，不同之处在于，这次发生变化的不是一个单一概率，而是整个概率分布。
以图4.2为例。这张图表展示了服用维生素\(E\)的人群中，其服用剂量的分布情况。我们进一步根据调查对象过去一个月是否进行过剧烈运动，对该分布进行了分组比较。
Figure 4.2: Distribution of Amount of vitamin E Taken by Exercise Level
ggplot2::ggplot(subset(oster, vite > 0), aes(x=vite, linetype = factor(vigorous_exercise_month))) +
# kernel density
ggplot2::geom_density(linewidth=1.5, na.rm=TRUE) +
ggplot2::scale_x_log10() +
ggplot2::labs(x='Vitamin E Taken (Log Scale)',y='Density') +
ggplot2::annotate(geom='label', x=3.5, y=.75, label='No Vigorous Exercise\nLast Month', hjust=1, family='sans', size=10/.pt) +
ggplot2::annotate(geom='label', x=15, y=.75, label='Vigorous Exercise\nLast Month', hjust=0, family='sans', size=10/.pt) +
ggplot2::guides(linetype=FALSE) +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family="sans"),
panel.grid = element_blank())## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot2::ggsave('DescribingRelationships/exercise_vite.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)我们可以观察到运动人群与非运动人群的分布存在轻微差异6。具体而言，经常运动的人群在服用维生素\(E\)时倾向于选择更大剂量。运动者与非运动者之间的分布差异表明，在这组数据中维生素\(E\)的服用与运动行为存在相关性。
上述示例针对的是连续变量，但同样适用于分类变量。现在我们不再关注服用剂量大小，而是考察是否服用维生素\(E\)这一行为本身！\(Oster\)的研究假设认为：服用维生素\(E\)的人群更可能从事其他健康行为（如运动），因为这两种行为都受到健康意识水平的驱动。
图4.3展示了这种情况的一个实例。图中分别呈现了吸烟人群与非吸烟人群服用维生素\(E\)的分布情况。两组分布存在明显差异：非吸烟人群中服用维生素\(E\)的比例更高——这完全符合\(Oster\)的预期。
Figure 4.3: Distribution of Whether vitamin E is Taken by Whether you Smoke
ggplot2::ggplot(oster %>%
dplyr::select(smoke,supplement_vite_single) %>%
stats::na.omit() %>%
dplyr::mutate(smoke=as.factor(case_when(smoke==0 ~ 'No Smoking',TRUE~'Smoking'))) %>%
dplyr::group_by(smoke) %>%
dplyr::mutate(BigN=n()) %>%
dplyr::group_by(smoke, supplement_vite_single) %>%
dplyr::summarize(N = n()/dplyr::first(BigN), .groups = "drop") %>%
dplyr::group_by(smoke, supplement_vite_single) %>%
dplyr::mutate(sup_label = case_when(
row_number() == 1 ~ supplement_vite_single,
TRUE ~ NA_character_
)),
aes(x = smoke, y = N, group = supplement_vite_single)) +
ggplot2::geom_col(aes(fill = supplement_vite_single), position = 'dodge', fill = 'white', color = 'blue') +
ggplot2::geom_text(aes(label = sup_label, y = N + .05),position = position_dodge(0.9), family = 'sans', size = 10/.pt) +
ggplot2::guides(fill = FALSE) +
ggplot2::labs(x=NULL, y='Proportion') +
ggplot2::scale_fill_manual(values = c('gray','red'))+
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family = 'sans'),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank())
ggplot2::ggsave('DescribingRelationships/smoking_vite.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)# 不同年级学生获得财政援助情况的频率表
freq_table <-
dplyr::tribble(
~Classification, ~Aid, ~freq,
"Freshman", "Y", 508,
"Freshman", "N", 371,
"Sophomore", "Y", 349,
"Sophomore", "N", 337,
"Junior", "Y", 425,
"Junior", "N", 384,
"Senior", "Y", 288,
"Senior", "N", 338,
)
freq_table %>%
dplyr::rename(`Financial Aid`=Aid) %>%
tidyr::pivot_wider(names_from=Classification, values_from=freq) %>%
vtable::dftoLaTeX(title = 'Financial Aid by Student Standing',
anchor = 'tab:describingrelationships-finaid',
file = 'DescribingRelationships/financial_aid_table.tex')## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{Financial Aid by Student Standing}\\label{tab:describingrelationships-finaid}\n\\begin{tabular}{lllll}\n\\hline\n\\hline\nFinancial Aid & Freshman & Sophomore & Junior & Senior \\\\ \n\\hline\nY & 508 & 349 & 425 & 288 \\\\ \nN & 371 & 337 & 384 & 338\\\\ \n\\hline\n\\hline\n\\end{tabular}\n\\end{table}\n"4.3 条件均值
在掌握条件分布的概念后，我们自然可以计算基于另一变量取值的分布特征：维生素\(E\)服用量的总体第95百分位数与吸烟人群的第95百分位数分别是多少？中位数是多少？对于服用量处于第90百分位和第10百分位的人群，其死亡率的标准差又有何差异？ 当这些可能性尚待探讨时，我们将聚焦于条件均值的研究：在给定\(X\)的特定取值时，\(Y\)的期望均值是多少7？
一旦获得条件均值，我们就能较好地描述两个变量之间的关系。如果\(Y\)的均值在\(X\)取值较高时也较高，那么\(Y\)与\(X\)就呈正相关。进一步地，我们可以绘制\(Y\)在每个\(X\)值上的所有条件均值，从而完整展现一个变量的均值如何与另一个变量的取值相关联。
在某些情况下，这种计算非常简单。如果条件变量是离散型（或分类）变量，只需计算该取值下所有观测值的均值即可。以图4.4为例，该图显示了在维生素\(E\)被推荐之前、推荐期间和推荐之后三个时期中服用维生素\(E\)的比例8。具体操作是：首先提取推荐前数据中的所有观测值，计算服用维生素\(E\)的比例；然后对推荐期间和推荐后的数据重复相同计算。
Figure 4.4: Proportion Taking Vitamin E Before It Was Recommended, During, and After
oster %>%
dplyr::filter(year < 2015 & !is.na(year)) %>%
dplyr::arrange(year) %>%
dplyr::mutate(time=fct_inorder(case_when(
year <= 1995 ~ 'Before Recommendation',
year <= 2004 ~ 'During Recommendation',
TRUE ~ 'After Recommendation'
))) %>%
dplyr::group_by(time) %>%
dplyr::summarize(vite=mean(supplement_vite_single=='Took Vitamin E')) %>%
ggplot2::ggplot(aes(x=time, y=vite)) +
ggplot2::geom_col(fill='white', color='black') +
ggplot2::geom_text(aes(label=time, y=vite+.005), family='sans', size=10/.pt) +
ggplot2::labs(x='Year', y='Proportion Taking Vitamin E') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank())
ggplot2::ggsave('DescribingRelationships/rec_vite.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)图4.4展示了维生素\(E\)服用行为与推荐时点之间的关系。我们可以观察到：维生素\(E\)服用与推荐政策实施呈正相关（推荐期间服用比例更高）。同时还能发现，维生素\(E\)服用与时间的关系最初为正相关（随着推荐政策生效而上升），随后转为负相关（随着推荐取消而下降）。
当条件变量为连续变量时，情况会变得稍复杂些。毕竟，我们无法准确计算年收入84,325美元人群的维生素\(E\)服用比例——因为这个精确数值对应的个体很可能不超过一人。事实上，对于多数具体数值，我们根本没有任何数据。
我们有两种处理方案：其一，对条件变量采用取值范围而非单一数值；其二，通过某种形状或线条拟合无观测数据的区间。
首先我们关注取值区间法。表4.1展示了以体重指数\(BMI\)为条件的维生素\(E\)服用比例。由于\(BMI\)是连续变量，我将其划分为十个等距区间（分组），并计算每个区间内的服用比例。在实际研究中，这种通过数据分组获取条件均值的方法并不常用，但它能直观展现后续其他方法的核心思想。
Table 4.1: Proportion Taking Vitamin E by Range of Body Mass Index Values
oster %>%
dplyr::filter(bmi < 100) %>%
# 将 bmi 等分
dplyr::mutate(bmi_cut = cut(bmi, 8)) %>%
dplyr::group_by(bmi_cut) %>%
dplyr::summarize(vite = mean(supplement_vite_single == 'Took Vitamin E')) %>%
dplyr::mutate(vite = scales::number(vite, accuracy = .001)) %>%
dplyr::rename(`BMI Bin` = bmi_cut,`Proportion Taking Vitamin E` = vite) %>%
# save to the latex
vtable::dftoLaTeX(
anchor = 'tab:describingrelationships-cut',
file = 'DescribingRelationships/vite_by_cut_bmi.tex',
title = 'Proportion Taking Vitamin E by Range of Body Mass Index Values'
)## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{Proportion Taking Vitamin E by Range of Body Mass Index Values}\\label{tab:describingrelationships-cut}\n\\begin{tabular}{ll}\n\\hline\n\\hline\nBMI Bin & Proportion Taking Vitamin E \\\\ \n\\hline\n(11.6,20.6] & 0.133 \\\\ \n(20.6,29.5] & 0.159 \\\\ \n(29.5,38.4] & 0.171 \\\\ \n(38.4,47.3] & 0.178 \\\\ \n(47.3,56.2] & 0.203 \\\\ \n(56.2,65.1] & 0.243 \\\\ \n(65.1,74] & 0.067 \\\\ \n(74,83] & 0.143\\\\ \n\\hline\n\\hline\n\\end{tabular}\n\\end{table}\n"这些区间同样可通过图形呈现（如图4.5）。图中的水平线段表明：我们对\(BMI\)同一区间内的所有观测值赋予相同均值，即该\(BMI\)分组对应的条件均值。由此可知，\(BMI\)与维生素\(E\)服用量呈正相关关系，但当\(BMI\)超过70后，条件均值开始下降。
Figure 4.5: Proportion Taking Vitamin E by Range of Body Mass Index Values
oster %>%
dplyr::filter(bmi < 100) %>%
dplyr::mutate(bmi_cut = cut(bmi, 8)) %>%
dplyr::group_by(bmi_cut) %>%
dplyr::summarize(vite = mean(supplement_vite_single=='Took Vitamin E'), bmi=min(bmi)) %>%
ggplot2::ggplot(aes(x = bmi, y = vite)) +
ggplot2::geom_step(color='black') +
ggplot2::geom_text(aes(label = bmi_cut, y = vite + .01), family = 'sans', size = 10/.pt) +
ggplot2::labs(x = 'Body Mass Index', y = 'Proportion Taking Vitamin E') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank())
ggplot2::ggsave('DescribingRelationships/vite_by_cut_bmi.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)当然，尽管这种分组方法简单直观，却也相当随意：我凭空选择了10个分组（而非9个、11个或其他数量）；使用等距区间同样具有随意性——并无实质性依据；更关键的是，这种划分方式过于生硬。难道真能认为某个区间上限的个体，与其区间下限的个体相似度，反而高于相邻区间的下限个体吗？
我们可以转而采用局部均值法 ，通过\(X\)值的邻近区间来计算\(Y\)的条件均值。具体而言，要计算\(X=2.5\)时的\(Y\)条件均值，只需取\(X\)值接近2.5的所有观测点的\(Y\)均值。这种方法需明确两个关键参数：邻近范围究竟多近才算”接近”？对不同程度接近的观测值是否赋予同等权重？
实现此类分析的常用方法是LOESS曲线（亦称LOWESS）910。该技术通过以下方式提供局部预测：针对\(X\)轴上的每个取值分别拟合曲线形态，且在拟合过程中对邻近程度不同的观测值赋予差异化权重——距离越近的观测值权重越高。最终得到的拟合曲线具有优良的平滑特性 。
图4.6展示了维生素E服用比例与BMI关系的LOESS拟合曲线。
Figure 4.6: Proportion Taking Vitamin E by BMI with a LOESS Curve
oster %>%
dplyr::filter(bmi < 100) %>%
dplyr::mutate(vite=(supplement_vite_single=='Took Vitamin E')*1) %>%
ggplot2::ggplot(aes(x=bmi, y=vite)) +
ggplot2::geom_smooth(color='black', se=FALSE, method='gam') +
ggplot2::labs(x='Body Mass Index', y='Proportion Taking Vitamin E') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank())## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
ggsave('DescribingRelations/vite_by_bmi_loess.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
从图4.6可以观察到明显的相关性：\(BMI\)值越高，服用维生素\(E\)的人群比例越大。这种关联起初非常显著，随后趋于平缓，但始终保持正相关趋势11。该曲线通过以下方式生成：仅计算特定\(BMI\)区间内人群的维生素\(E\)服用比例，并令这个”特定区间”逐步右移，从而构建出一系列条件均值。
4.4 线性拟合
通过\(X\)的局部值展示\(Y\)的均值虽具价值，能细致呈现\(X\)与\(Y\)的关联特征，但该方法仍存在局限：其一，数据缺口可能难以填补；其二，观察到的关联模式往往难以简明表述12。
由此引入线性拟合概念（亦称回归分析）13。
不同于局部思维下的条件均值估计，我们可以假设\(Y\)与\(X\)的本质关系可通过某种数学形态表征。基础回归分析中，该形态表现为直线。例如直线方程：
\[Y=3 + 4 \times X \space
\space (公式4.1) \]
表明：当\(X=5\)时，Y的条件均值为\(3+4\times 5=23\)；同时揭示\(X\)每增加1个单位，\(Y\)的条件均值将相应增加4个单位。
在图4.7中，我们重现了先前维生素\(E\)与\(BMI\)的关系数据，但此次采用直线拟合。该拟合直线的斜率为0.002，意味着\(BMI\)每增加1个单位，服用维生素\(E\)补充剂的概率将比低1个单位的人群高出0.2个百分点。
这种方法具有显著优势：其一，它能提供任意\(X\)值对应的\(Y\)条件均值——即使该特定值并无实际观测数据14；其二，可清晰描述\(Y\)与\(X\)的关系。若\(X\)的斜率系数（如维生素\(E\)/\(BMI\)回归中的0.002）为正，则两变量呈正相关；若为负，则呈负相关。
采用拟合直线具有多重实用优势，其统计层面的益处更为显著：由于直线估计基于全体数据而非局部数据，结果精度更高；此外，该方法可轻松扩展至多元变量分析（下节将详细探讨）。
Figure 4.7: Proportion Taking Vitamin E by BMI with a Fitted Straight Line
oster %>%
dplyr::filter(bmi < 100) %>%
dplyr::mutate(vite = (supplement_vite_single == 'Took Vitamin E')*1) %>%
ggplot2::ggplot(aes(x = bmi, y = vite)) +
ggplot2::geom_smooth(method = 'lm', color = 'black', se = FALSE) +
ggplot2::labs(x = 'Body Mass Index', y = 'Proportion Taking Vitamin E') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank())## `geom_smooth()` using formula = 'y ~ x'
ggplot2::ggsave('DescribingRelationships/vite_by_bmi_lm.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)## `geom_smooth()` using formula = 'y ~ x'当然，这种方法也存在缺陷。最显著的局限在于：直线拟合本质上强制要求关系呈线性——我们必须预先设定关联形态（直线？曲线？还是波动线？）。拟合程序虽能优化给定形态的参数，但若初始形态设定错误，条件均值估计将完全失真。试想用直线描述图4.6中的非线性关系！
问题的关键不在于直线拟合必然有误（回归分析同样允许使用曲线拟合），而在于必须预先准确判断何时该用曲线，并确定具体的曲线类型。
当然，这些缺陷需与其显著优势权衡——正因优势突出，拟合方法已成为所有应用统计领域极其普遍的做法。那么，具体该如何操作呢？
普通最小二乘法(OLS)是直线拟合中最著名的应用方法 。该算法通过选择使残差平方和最小的直线进行拟合。残差即观测值实际数值与拟合直线所赋条件均值之间的差异15
以前文所述\(Y=3+4X\)直线为例：当\(X=5\)时，我们计算出\(Y\)的条件均值为\(3+4 \times 5=23\)。但若观测到\(X=5\)而\(Y=25\)的个体，其残差即为\(25-23=2\)。\(OLS\)算法将此残差平方得4，再对数据集中所有预测值执行相同运算并求和。最终，该算法会确定直线\(Y=\beta_0+\beta_1X\)中的系数\(\beta_0\)和\(\beta_1\)，使残差平方和最小化（如图4.8所示）。
Figure 4.8: Fitting an OLS Line to Four Points
tb2 <- dplyr::tibble(X=c(1,3,5,6.5), Y=c(2,6,3.6, 8)) %>%
dplyr::mutate(line1=1.5+X, line2=2+.6*X, line3=1.7324+.8175*X)
p1 <- ggplot2::ggplot(tb2, aes(x=X, y=Y)) +
ggplot2::geom_point(size=2) +
ggplot2::expand_limits(x = c(0,7), y = c(1,9))+
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank(),
plot.title.position = 'plot') +
ggplot2::coord_fixed() +
labs(title = 'Let\'s fit a line to four points')
p2 <- ggplot2::ggplot(tb2, aes(x=X, y=Y)) +
ggplot2::geom_point(size=2) +
ggplot2::expand_limits(x=c(0,7), y=c(1,9))+
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank(),
plot.title.position = 'plot') +
ggplot2::coord_fixed() +
ggplot2::geom_line(aes(x=X, y=line3), size=1.0) +
ggplot2::labs(title = 'Add the OLS line')## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p3 <- ggplot2::ggplot(tb2, aes(x=X, y=Y)) +
ggplot2::geom_point(size=2) +
ggplot2::expand_limits(x=c(0,7), y=c(1,9))+
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank(),
plot.title.position = 'plot') +
ggplot2::coord_fixed() +
ggplot2::geom_segment(aes(x=X, xend=X, y=Y, yend=line3), size=1.0, color='blue') +
ggplot2::geom_line(aes(x=X, y=line3), linewidth=1.0) +
ggplot2::labs(title = 'Residuals are from point to line')
p4 <- ggplot2::ggplot(tb2, aes(x=X, y=Y)) +
ggplot2::geom_point(size=2) +
ggplot2::expand_limits(x=c(0,7), y=c(1,9))+
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.x = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family= "sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank(),
plot.title.position = 'plot') +
ggplot2::geom_segment(aes(x=X, xend=X, y=Y, yend=line3), size=1.0) +
ggplot2::geom_rect(aes(xmin=X, xmax=X+abs(line3-Y), ymin=Y, ymax=line3), alpha=.5, fill='red') +
ggplot2::geom_line(aes(x=X, y=line3), linewidth=1.0) +
ggplot2::labs(title = 'Goal: minimize squared residuals')+
coord_fixed()
plot_grid(p1,p2,p3,p4, nrow = 2)
ggsave('DescribingRelationships/ols_fit.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)其实现原理在于利用两变量协同变动的信息——这些信息编码在协方差之中16
若回顾第3章所学的方差概念，您应记得计算变量X方差的步骤：(1)用X减去其均值；(2)将差值平方；(3)对所有观测值求和；(4)除以样本量减一。最终所得的方差值反映了该变量的实际变异程度。
协方差的计算原理与方差完全相同，仅存在两处关键差异：(1)步骤(a)需对两个变量分别去均值；(2)步骤(b)需将两变量的差值相乘。最终所得的协方差度量两变量的协同变动：若两变量同时高于或低于均值，多数观测值的乘积结果为正，协方差增大；若两变量互不相关，则乘积结果正负各半，在步骤(c)相互抵消，使协方差趋近于0。
\(OLS\)如何利用协方差推导\(Y\)与\(X\)的关系？其核心公式仅为：\(\frac{cov(X,Y)}{var(X)}\)17。该计算本质上是在回答”在\(X\)的所有变异中，有多少与\(Y\)协同变动？“这一命题18。确定斜率后，算法通过调整截距使残差均值（非平方残差）归零，从而确保条件均值至少在整体上准确无误。
\(OLS\)最终生成带截距和斜率的直线方程（如\(Y=3+4X\)）。通过代入\(X\)值即可求得\(Y\)的条件均值。关键在于，该斜率能清晰描述变量关系：由于方程含\(4X\)项，可知\(X\)每增加1个单位，\(Y\)平均增加4个单位。
有时我们需要对\(OLS\)结果进行标准化处理，这就引入了相关系数的概念。具体而言，皮尔逊相关系数通过对\(OLS\)斜率进行标准化计算获得：将斜率乘以\(X\)的标准差，再除以\(Y\)的标准差。该计算等价于用\(X\)与\(Y\)的协方差同时除以两变量的标准差。
相关系数同样基于直线拟合原理19，只是呈现方式不同：虽然损失了斜率在\(\frac{X}{Y}\)原单位的解释性，但能更直观判断关联强度。相关系数取值限定于[-1,1]区间，且不受原变量单位影响——越接近-1表示变量反向变动（负斜率）越显著；越接近1则同向变动（正斜率）越显著。
对于维生素\(E\)与\(BMI\)的关系，\(OLS\)估计得出以下线性方程：
\[维生素E = \beta_0 +
\beta_1BMI \space\space\space (公式4.3)\]
经参数优化后得到最佳拟合方程：
\[ 维生素E = 0.110 + 0.002BMI \space\space\space (公式4.3) \]\(BMI\)每增加1个单位，维生素\(E\)的条件均值预计增加0.002 。鉴于维生素\(E\)是二元变量，该条件均值的0.002增量可理解为：服用维生素\(E\)的人群比例将上升0.2个百分点。
由于服用维生素\(E\)的标准差为0.369，\(BMI\)的标准差为6.543，两者的皮尔逊相关系数为：0.002 × (6.543/0.369) ≈ 0.0355。
线性假设有时并不适用。虽然\(OLS\)强制拟合直线，但许多变量间本质上是非线性关联——如图4.6所示，维生素\(E\)与\(BMI\)的关系正是如此。此时该如何处理？
此时有两大方法可破解此局20：
首个方法竟由”反派”\(OLS\)摇身变为救星——事实上，\(OLS\)并非必须拟合直线。其关键要求是”系数线性”，即斜率系数仅需与变量相乘，无需更复杂的运算关系。
\(OLS\)完全可以估计\(Y=\beta_0+\beta_1X\)这类方程，同样适用于\(Y=\beta_0+\beta_1+\beta X^2\)（非线性！）或\(Y=\beta_0+\beta_1ln(X)\)（亦非线性！）等模型。而诸如\(Y=\beta_0+X^{\beta_1}\)或\(Y=\frac{\beta_0}{1+\beta_1X}\)等方程则不符合”系数线性”要求。
因此，图4.6中看似复杂的曲线并非难题——只要预先分析数据形态（需二次项构成抛物线？需对数项实现快速上升后趋平？），我们就能通过\(OLS\)模拟此类非线性关系。针对图4.6，采用\(Y=\beta_0+\beta_1+\beta X^2\)方程即可兼顾\(LOESS\)的灵活性与\(OLS\)的拟合优势。
set.seed(20200703)
dat <- dplyr::tibble(x=rnorm(100)) %>%
dplyr::mutate(y=1+2*x+2*x^2+rnorm(100))
a <- ggplot2::ggplot(dat, aes(x=x, y=y)) +
ggplot2::geom_jitter(alpha=.5) +
ggplot2::stat_smooth(formula=y~x, method="lm", se=F, color='black') +
ggplot2::labs(x="Math Exam Scores", y="Intelligence Scores") +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"))
b <- ggplot2::ggplot(dat, aes(x=x, y=y)) +
ggplot2::geom_jitter(alpha=.5) +
ggplot2::stat_smooth(formula=y~x, method="loess", se=F, color='black') +
ggplot2::labs(x="Math Exam Scores", y="Intelligence Scores") +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"))
cowplot::plot_grid(a, b, labels = "AUTO")
ggplot2::ggsave('DescribingRelationships/parabolic.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)第二种方法是非线性回归——其形态极其多样，通常表现为\(Y=F(\beta_0+\beta_1X)\)的形式，其中函数\(F(.)\)的具体形式取决于建模需求。
非线性回归尤其适用于因变量取值受限的场景。以维生素\(E\)（二元变量：服用/不服用）与\(BMI\)的关系为例：\(OLS\)直线拟合无法真实反映二元跳变特性——即便采用\(维生素E=\beta_0+\beta_1BMI+\beta_2BMI^2\)的曲线拟合，当将因变量解释为服用概率时，仍会导致荒谬预测（极高\(BMI\)者概率>100%，极低\(BMI\)者概率<0%）。
解决之道在于采用概率单元函数\(probit\)或逻辑函数\(logit\)等\(F()\)——这类函数能确保输出值始终介于0%-100%之间（详见第13章）。
实际上，针对不同类型因变量的取值特性，还存在众多其他函数形式。本书虽不深入探讨，但读者应意识到：这些方法为非线性拟合提供了重要工具。
set.seed(20200715)
s_x <- 1
s_y <- 1
rho <- .6
sigma <- base::matrix(c(s_x^2, s_x*s_y*rho, s_x*s_y*rho, s_y^2),2)
dat <- dplyr::as_tibble(rmvnorm(n = 5, mean = c(5, 5), sigma = sigma)) %>%
dplyr::rename(get_along = V1, sat = V2)
mod <- stats::lm(sat~get_along, data=dat)
dat <- dat %>%
dplyr::mutate(pred_sat=predict(mod)) %>%
dplyr::mutate_all(round, 2) %>%
dplyr::mutate(id = 1:5) %>%
dplyr::select(id, everything())
dat %>%
dplyr::select(-id) %>%
rename(GetAlong = get_along,Satisfaction = sat,Prediction = pred_sat) %>%
dftoLaTeX(title = 'Satisfaction and Getting Along with Coworkers',
anchor = 'fig:describingrelationships-getalong',
file = 'DescribingRelationships/get_along.tex')## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{Satisfaction and Getting Along with Coworkers}\\label{fig:describingrelationships-getalong}\n\\begin{tabular}{lll}\n\\hline\n\\hline\nGetAlong & Satisfaction & Prediction \\\\ \n\\hline\n4.7 & 5.07 & 4.72 \\\\ \n4.21 & 4.05 & 4.28 \\\\ \n5.42 & 5.33 & 5.38 \\\\ \n4.14 & 4.02 & 4.22 \\\\ \n3.3 & 3.59 & 3.45\\\\ \n\\hline\n\\hline\n\\end{tabular}\n\\end{table}\n"4.5 双重条件均值—-变量控制
让我们进入统计建模的”未解释之地”——这里特指残差领域。
无论采用何种方法计算\(Y\)基于\(X\)的条件均值，本质上都将观测值分解为两部分：\(X\)可解释部分（条件均值）与\(X\)未解释部分（残差）。例如，当\(X=5\)时\(Y\)的条件均值为10，而实际观测到\(X=5\)且\(Y=13\)的数据点，则预测值为10，残差为13-10。图4.9直观展示了条件均值与残差的区分方式21。
Figure 4.9: An OLS Line and its Residuals
set.seed(1234)
tib <- base::data.frame(x=runif(20), eps=rnorm(20)) %>% dplyr::mutate(y = .5*x + eps)
tib <- tib %>% dplyr::mutate(pred=predict(lm(y~x, data=tib)))
ggplot2::ggplot(tib, aes(x=x, y=y)) +
ggplot2::geom_point(size=2) +
ggplot2::geom_line(aes(x=x, y=pred), size=1, linetype='solid') +
ggplot2::geom_segment(aes(x=x, xend=x, y=y, yend=pred), linetype='dashed') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family="sans"),
panel.grid = element_blank(),
axis.ticks = element_blank()) +
ggplot2::labs(x="X", y="Y") +
ggplot2::geom_rect(aes(xmin=.46, xmax=.85, ymin=-2.42, ymax = -1.5), fill = 'white', color = 'black') +
ggplot2::annotate('segment', x = .48, xend = .58, y= -2.2, yend = -2.2, linetype='dashed') +
ggplot2::annotate('segment', x = .48, xend = .58, y= -1.8, yend = -1.8, linetype='solid') +
ggplot2::annotate(geom = 'text', x = .6, y = -1.8,
label='Conditional Mean (from OLS)',
hjust=0,
size=10/.pt,
lineheight=.7,
family='sans') +
ggplot2::annotate(geom = 'text', x = .6, y = -2.2, label='Residual', hjust = 0, size = 10/.pt, family = 'sans')
ggplot2::ggsave('DescribingRelationships/ols_resid.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)表面看来，残差似乎只是无法预测的干扰项或模型缺陷。但事实上，它们蕴含着重要信息——残差本质上反映了\(Y\)中与\(X\)无关的部分。例如：当条件均值为10而实际值为13时，\(X\)仅能解释10这个部分，超出的3必定源于数据生成过程中的其他因素。
残差分析的实际价值何在？其应用场景多样：例如，我们可能不仅关注维生素\(E\)的整体变异，更想量化其中不受\(BMI\)影响的部分。通过分析图4.7的残差，即可精确解答这一问题。
当我们同时考察两个变量的残差时，分析将变得极具启发性。若从两个不同变量中剔除已解释部分会如何？现引入第三个变量\(Z\)扩展分析框架（保持\(X/Y/Z\)的简洁结构），具体操作方法如下22：
计算\(Y\)基于\(Z\)的条件均值
求\(Y\)与该条件均值之差，获得\(Y\)的残差\(Y^R\)
计算\(X\)基于\(Z\)的条件均值
求\(X\)与该条件均值之差，获得\(X\)的残差\(X^R\)
分析\(Y^R\)与\(X^R\)之间的关联关系
由于\(Y^R\)和\(X^R\)已剔除\(Z\)可解释的部分，两者间的关系即反映\(Y\)与\(X\)关系中不受\(Z\)影响的成分。换言之，我们获得的是基于\(Z\)条件下的\(Y\)对\(X\)的条件均值——该方法能有效过滤掉\(X/Y\)关系中由\(Z\)解释的部分。
通过此操作，我们剔除了与\(Z\)相关的所有变异，实质上固定了\(Z\)的变动。正因如此，该过程被称为“控制\(Z\)”（尽管”调整\(Z\)“的表述或许更为精确） 。
以冰淇淋消费量与短裤穿着量为例：二者在观测中呈现正相关，但气温变量同时影响着这两个指标。
若需确证冰淇淋消费是否直接影响短裤穿着，关键在于量化二者关系中不受气温解释的部分。具体步骤：
计算冰淇淋消费量基于气温的条件均值，提取残差（获得仅与气温无关的冰淇淋消费变异）
计算短裤穿着量基于气温的条件均值，提取残差（获得仅与气温无关的短裤穿着变异）
分析短裤穿着残差对冰淇淋消费残差的条件均值
若短裤穿着均值随冰淇淋残差变化不大，则原关联可完全归因于气温；若仍存在强相关，则可能揭示真实效应。
实现双重条件均值的最简方法是通过多元回归。只需在方程中加入控制变量\(Z\)，即可自动完成对\(Z\)的统计控制。原单变量模型：
\[Y=\beta_0+\beta_1X \space\space\space\space
(公式4.4)\]
升级为多元模型：
\[Y=\beta_0+\beta_1X+\beta_2Z \space\space\space\space
(公式4.5)\]
此时，\(OLS\)对\(\beta_1\)的估计将自动执行”剔除\(Z\)的影响后分析\(Y\)与\(X\)关系”的全过程。
更强大的是，该方法可扩展至多重变量控制。例如添加变量\(W\)构建方程：
\[Y=\beta_0+\beta_1X+\beta_2Z+\beta_3W \space\space\space\space
(公式4.6)\]
此时\(OLS\)估计的\(\beta_1\)即表示同时控制\(Z\)和\(W\)时\(Y\)与\(X\)的关系。
以维生素\(E\)与\(BMI\)关系为例：性别和年龄可能同时影响这两个变量。将这些变量纳入回归后，我们将获得更精确的关联估计。
stats::lm((supplement_vite_single=='Took Vitamin E')~bmi, data = oster %>% dplyr::filter(bmi < 100))##
## Call:
## stats::lm(formula = (supplement_vite_single == "Took Vitamin E") ~
## bmi, data = oster %>% dplyr::filter(bmi < 100))
##
## Coefficients:
## (Intercept) bmi
## 0.109926 0.001843stats::sd(oster$supplement_vite_single=='Took Vitamin E', na.rm = TRUE)## [1] 0.3687025stats::sd(oster$bmi, na.rm = TRUE)## [1] 6.543197# HH income control
m1 <- oster %>%
dplyr::filter(bmi < 100) %>%
dplyr::mutate(vite=(supplement_vite_single=='Took Vitamin E')*1) %>%
stats::lm(vite~bmi, data=.)
m2 <- oster %>%
dplyr::filter(bmi < 100) %>%
dplyr::mutate(vite=(supplement_vite_single=='Took Vitamin E')*1) %>%
stats::lm(vite~bmi+age+I(gender==2), data=.)
modelsummary::msummary(list(m1,m2), stars=TRUE)| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 0.110*** | -0.006 |
| (0.007) | (0.008) | |
| bmi | 0.002*** | 0.001*** |
| (0.000) | (0.000) | |
| age | 0.002*** | |
| (0.000) | ||
| I(gender == 2)TRUE | 0.016*** | |
| (0.003) | ||
| Num.Obs. | 59018 | 59018 |
| R2 | 0.001 | 0.016 |
| R2 Adj. | 0.001 | 0.016 |
| AIC | 49676.0 | 48800.0 |
| BIC | 49703.0 | 48844.9 |
| Log.Lik. | -24834.998 | -24394.985 |
| RMSE | 0.37 | 0.37 |
此前仅含\(BMI\)的模型估计结果为：
\[维生素E = 0.110 +
0.002BMI \space\space\space\space (公式4.7)\]
而加入性别和年龄变量后，新模型为：
\[ 维生素E = -0.006 + 0.0018BMI + 0.002Age + 0.016Female \space\space\space\space (公式4.8) \]\(BMI\)的系数从0.002降至0.001，表明原\(BMI\)与维生素\(E\)的部分关联实际由年龄和/或性别因素解释。新模型还显示：
年龄效应：年龄每增加1岁，维生素\(E\)服用比例预计上升0.2个百分点
性别效应：女性比男性的服用比例高1.6个百分点（“Female”变量从0→1的变化）
需注意的是，\(OLS\)本身并不识别处理变量——年龄效应的估计已控制\(BMI\)和性别的影响，性别效应的估计亦控制\(BMI\)和年龄。但需避免过度解读控制变量的系数，因为这些估计通常未经过严格的因果识别检验（详见第5章）。
回归分析如何实现这一功能？让我们深入其数学本质：
线性代数视角（了解者可阅，否则请跳至下段）：
多元\(OLS\)估计公式为\((A^\top A)^{-1}A^\top Y\) ，其中矩阵\(A\)包含除\(Y\)外的所有变量（包括目标变量\(X\)）。该公式通过协方差矩阵运算，自动消除所有非\(X\)变量对\(X/Y\)关系的影响。
几何视角：若将双变量\(OLS\)直线\(Y=\beta_0+\beta_1X\)视为二维空间中的直线，则三变量\(OLS\)模型可理解为三维空间中的平面（变量继续增加则对应更高维超平面）。我们可通过三维图像的三个侧面逐一观察这种关系。
图4.10左上角展示\(X-Y\)轴关系，右侧为\(Z-Y\)轴，下方是\(Z-X\)轴。需注意\(Z-X\)轴的坐标方向进行了翻转——虽然此处计算的是基于\(Z\)的\(X\)条件均值，但为保持与\(X-Y\)图的一致性，仍将\(X\)置于横轴。\(Z-Y\)与\(Z-X\)轴上的上升斜率表明：\(Z\)能同时解释\(X\)和\(Y\)的部分变异，剔除这些影响后即可聚焦于残差分析。
而在图4.11中，这些解释性影响被”熨平”——上升斜率经调整后消失，连带移动\(X\)和\(Y\)的数据点。如图所示，通过数学剔除\(Z\)解释的部分后，\(X/Y\)关系中已完全不包含\(Z\)的影响。\(Z\)在双向维度上被完全中和，不再对\(X/Y\)散点图产生任何”抬升”效应。此时\(X/Y\)平面上呈现的，已是剔除\(Z\)后的纯净关联。
Figure 4.10: A ThreeVariable Regression from All Three Dimensions
Figure 4.11: A ThreeVariable Regression from All Three Dimensions After Removing the Variation Explained by Z
# 3-D
set.seed(1015114)
df <- dplyr::tibble(eps=rnorm(100), nu=rnorm(100), Z=rnorm(100)) %>%
dplyr::mutate(X=Z+nu) %>%
dplyr::mutate(Y=X+2*Z+eps)
y_model <- stats::lm(Y~Z, data=df)
x_model <- stats::lm(X~Z, data=df)
df <- df %>%
dplyr::mutate(X_res = stats::residuals(x_model),Y_res = stats::residuals(y_model))
pmain <- ggplot2::ggplot(df, aes(x=X, y=Y)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(formula = 'y ~ x',se=FALSE, color='black', size=1, method='lm') +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size = 10, family="sans"),
axis.line = element_line(linewidth=1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.ticks = element_blank())
pywithline <- ggplot2::ggplot(df, aes(x=Z, y=Y)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(formula = 'y ~ x', se=FALSE, method='lm', color='black', size=1) +
ggplot2::scale_x_continuous(limits=c(min(df$Y), max(df$Y))) +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.line = element_line(linewidth=1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.ticks = element_blank())
pxwithline <- ggplot2::ggplot(df, aes(x=Z, y=X)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(formula = 'y ~ x', se=FALSE, method='lm', color='black', size=1) +
ggplot2::scale_y_continuous(limits=c(min(df$X), max(df$X))) +
ggpubr::theme_pubr() +
ggplot2::coord_flip() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.line = element_line(linewidth=1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family="sans"),
axis.title.y = element_text(size=10, family="sans"),
panel.grid = element_blank(),
axis.ticks = element_blank())
cowplot::plot_grid(pmain, pywithline, pxwithline, ncol=2)
ggplot2::ggsave('DescribingRelationships/control_anim_1.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)pmain <- ggplot2::ggplot(df, aes(x=X_res, y=Y_res)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(formula = 'y ~ x',se=FALSE, color='black', size=2, method='lm') +
ggpubr::theme_pubr() +
ggplot2::scale_y_continuous(limits=c(min(df$Y), max(df$Y))) +
ggplot2::scale_x_continuous(limits=c(min(df$X), max(df$X))) +
ggplot2::labs(x='X Residual',y='Y Residual')+
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.line = element_line(size = 1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.ticks = element_blank())## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.pywithline <- ggplot2::ggplot(df, aes(x=Z, y=Y_res)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(formula = 'y ~ x', se=FALSE, method='lm', color='black', size=1) +
ggplot2::scale_y_continuous(limits=c(min(df$Y), max(df$Y))) +
ggpubr::theme_pubr() +
ggplot2::theme(text = element_text(size=10, family="sans"),
axis.line = element_line(size=1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.ticks = element_blank()) +
ggplot2::labs(y='Y Residual')
pxwithline <- ggplot(df, aes(x=Z, y=X_res)) +
geom_point() +
geom_smooth(formula = 'y ~ x', se=FALSE, method='lm', color='black', size=2) +
scale_y_continuous(limits=c(min(df$X), max(df$X))) +
theme_pubr() +
coord_flip() +
theme(text = element_text(size=10, family="sans"),
axis.line = element_line(size = 1),
axis.text.y = element_blank(),
axis.text.x = element_blank(),
axis.title.x = element_text(size=10, family= "sans"),
axis.title.y = element_text(size=10, family= "sans"),
panel.grid = element_blank(),
axis.ticks = element_blank())+
labs(y="X Residual")
plot_grid(pmain, pywithline, pxwithline, ncol=2)
ggsave('DescribingRelationships/control_anim_2.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)4.6 未涵盖内容
前一章关于变量描述的内容已基本涵盖核心知识点，而本章在变量关系描述方面存在更多未涉及的方法。 这种取舍主要基于研究重点的考量。本书聚焦于研究设计——只有确立研究设计后，才需处理后续统计问题。诸如具体概率分布（正态/对数正态/\(t\)分布/泊松分布等）、函数形式（\(OLS/probit/logit\)等）或标准误与假设检验等细节，反而可能干扰研究核心问题的解决思路。
有一处省略并非基于聚焦考量，而是为后续更系统的探讨：请注意前文虽以\(Oster\)论文为例说明维生素\(E\)与健康指标的关系随时间变化，却未展示\(BMI\)关系的时序变化23。实际上，研究关联性在不同情境下的演变本身构成重要研究设计——这部分内容将在本书第二部深入展开。
必须明确，这些知识至关重要。本书第13章及第二部其他章节将深入探讨，读者亦可参考《真实计量经济学》或《计量经济学导论》等经典教材。但就观测数据而言，上述内容通常应在确定研究设计并准备实施后考量。
现在，我希望您思考如何处理数据——您的研究设计需要怎样的变量描述、需要揭示何种变量关系、涉及哪些条件均值与双重条件均值分析。先规划数据的分析路径，正如旅行前确定行程，出发时再收拾行装。
4.7 软件中的变量关系分析
本节将演示如何在三种编程环境中分析变量关系：R语言/Stata/Python
以下代码可能依赖需安装的扩展包：
R语言：
library(X)或X::表示调用包X，未安装时需执行install.packages('X')Python：
import X或from X import需通过pip/conda等包管理器安装Stata：扩展包无需重复加载，但示例中将标注需安装的包
书所有案例数据集均收录于causaldata扩展包（支持三种语言环境）：
R：执行
install.packages('causaldata')Stata：执行
ssc install causaldataPython：执行
pip install causaldata（如使用pip）
现在开始代码演示！虽然\(Oster\)数据可免费下载但需特殊授权，故改用\(Mroz\)(1987)的研究数据24——该数据集包含1975年女性劳动力参与率与收入信息。
以下为各语言统一分析流程：
数据加载：导入\(Mroz\)数据集
散点图绘制：在职女性的对数收入与家庭其他收入对数关系25
条件均值分析：按是否接受高等教育分组计算女性收入均值
分组条件均值：按家庭其他收入区间分组计算女性收入均值
曲线拟合：
\(LOESS\)平滑曲线
线性回归直线（均基于女性收入对数对家庭其他收入对数的条件均值）
回归建模：
基准模型：女性收入对数 ~ 家庭其他收入对数
扩展模型：增加高等教育参与率和5岁以下子女数量控制变量
library(tidyverse)
library(modelsummary)
df <- causaldata::Mroz %>%
# Keep just working women
dplyr::filter(lfp == TRUE & inc > 0) %>%
# Get unlogged earnings
dplyr::mutate(earn = exp(lwg))
# 1. Draw a scatterplot
ggplot2::ggplot(df, aes(x = inc, y = earn)) +
ggplot2::geom_point() +
# Use a log scale for both axes
# We'll get warnings as it drops the 0s, that's ok
ggplot2::scale_x_log10() +
ggplot2::scale_y_log10()
# 2. Get the conditional mean by college attendance
df %>%
# wc is the college variable
dplyr::group_by(wc) %>%
# Functions besides mean could be used here to get other conditionals
dplyr::summarize(earn = mean(earn))## # A tibble: 2 × 2
## wc earn
## <lgl> <dbl>
## 1 FALSE 3.58
## 2 TRUE 5.34# 3. Get the conditional mean by bins
df %>%
# use cut() to cut the variable into 10 bins
dplyr::mutate(inc_cut = cut(inc, 10)) %>%
dplyr::group_by(inc_cut) %>%
dplyr::summarize(earn = mean(earn))## # A tibble: 10 × 2
## inc_cut earn
## <fct> <dbl>
## 1 (1.11,10.2] 3.03
## 2 (10.2,19.2] 3.95
## 3 (19.2,28.1] 4.99
## 4 (28.1,37.1] 4.63
## 5 (37.1,46.1] 4.28
## 6 (46.1,55.1] 6.91
## 7 (55.1,64.1] 3.83
## 8 (64.1,73] 5.23
## 9 (73,82] 7.02
## 10 (82,91.1] 3.67# 4. Draw the LOESS and linear regression curves
ggplot2::ggplot(df, aes(x = inc, y = earn)) +
ggplot2::geom_point() +
## geom_smooth by default draws a LOESS; we don't want standard errors
ggplot2::geom_smooth(se = FALSE) +
ggplot2::scale_x_log10() +
ggplot2::scale_y_log10()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# Linear regression needs a 'lm' method
ggplot2::ggplot(df, aes(x = inc, y = earn)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(method = 'lm', se = FALSE) +
ggplot2::scale_x_log10() +
ggplot2::scale_y_log10()## `geom_smooth()` using formula = 'y ~ x'
# 5. Run a linear regression, by itself and including controls
model1 <- stats::lm(lwg ~ log(inc), data = df)
# k5 is number of kids under 5 in the house
model2 <- stats::lm(lwg ~ log(inc) + wc + k5, data = df)
# And make a nice table
modelsummary::msummary(list(model1, model2))| (1) | (2) | |
|---|---|---|
| (Intercept) | 0.569 | 0.707 |
| (0.185) | (0.184) | |
| log(inc) | 0.221 | 0.135 |
| (0.065) | (0.066) | |
| wcTRUE | 0.335 | |
| (0.075) | ||
| k5 | -0.067 | |
| (0.087) | ||
| Num.Obs. | 427 | 427 |
| R2 | 0.027 | 0.070 |
| R2 Adj. | 0.024 | 0.064 |
| AIC | 928.7 | 913.1 |
| BIC | 940.9 | 933.4 |
| Log.Lik. | -461.340 | -451.536 |
| F | 11.651 | 10.673 |
| RMSE | 0.71 | 0.70 |
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import seaborn as sns
import seaborn.objects as so
from causaldata import Mroz
# Read in data
dt = Mroz.load_pandas().data
# Keep just working women
dt = dt.query('lfp == True')
# Create unlogged earnings
dt = dt.assign(earn = lambda x: np.exp(x['lwg']))
# 1. Draw a scatterplot
(so.Plot(dt, x = 'inc', y = 'earn').add(so.Dot()).scale(x = 'log', y = 'log'))## <seaborn._core.plot.Plot object at 0x0000013C8CD86850># 2. Get the conditional mean by college attendance
# wc is the college variable
dt.groupby('wc').agg({'earn': 'mean'})## earn
## wc
## False 3.583539
## True 5.349448# 3. Get the conditional mean by bins
# Use cut to get 10 bins
dt = dt.assign(inc_bin = lambda x: pd.cut(x['inc'],10))
dt.groupby('inc_bin').agg({'earn': 'mean'})## <string>:1: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
## earn
## inc_bin
## (-0.12, 9.074] 3.091594
## (9.074, 18.177] 3.826430
## (18.177, 27.28] 4.954454
## (27.28, 36.383] 4.563332
## (36.383, 45.485] 4.241790
## (45.485, 54.588] 6.139739
## (54.588, 63.691] 3.825000
## (63.691, 72.794] 5.232558
## (72.794, 81.897] 7.015306
## (81.897, 91.0] 3.666667# 4. Draw the LOESS and linear regression curves
# Do log beforehand for these axes
dt = dt.assign(linc = lambda x: np.log(x['inc']))## C:\Users\capen\AppData\Roaming\Python\Python313\site-packages\pandas\core\arraylike.py:399: RuntimeWarning: invalid value encountered in log
## result = getattr(ufunc, method)(*inputs, **kwargs)sns.regplot(x = 'linc', y = 'lwg',data = dt, lowess = True, ci = None)
# 5. Run a linear regression, by itself and including controls
m1 = smf.ols(formula = 'lwg ~ linc', data = dt).fit()
print(m1.summary())## OLS Regression Results
## ==============================================================================
## Dep. Variable: lwg R-squared: 0.027
## Model: OLS Adj. R-squared: 0.024
## Method: Least Squares F-statistic: 11.65
## Date: 周日, 17 8月 2025 Prob (F-statistic): 0.000703
## Time: 22:38:48 Log-Likelihood: -461.34
## No. Observations: 427 AIC: 926.7
## Df Residuals: 425 BIC: 934.8
## Df Model: 1
## Covariance Type: nonrobust
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 0.5690 0.185 3.079 0.002 0.206 0.932
## linc 0.2207 0.065 3.413 0.001 0.094 0.348
## ==============================================================================
## Omnibus: 65.243 Durbin-Watson: 1.947
## Prob(Omnibus): 0.000 Jarque-Bera (JB): 163.890
## Skew: -0.761 Prob(JB): 2.58e-36
## Kurtosis: 5.626 Cond. No. 17.1
## ==============================================================================
##
## Notes:
## [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.# k5 is number of kids under 5 in the house
m2 = smf.ols(formula = 'lwg ~ linc + wc + k5', data = dt).fit()
print(m2.summary())## OLS Regression Results
## ==============================================================================
## Dep. Variable: lwg R-squared: 0.070
## Model: OLS Adj. R-squared: 0.064
## Method: Least Squares F-statistic: 10.67
## Date: 周日, 17 8月 2025 Prob (F-statistic): 8.92e-07
## Time: 22:38:48 Log-Likelihood: -451.54
## No. Observations: 427 AIC: 911.1
## Df Residuals: 423 BIC: 927.3
## Df Model: 3
## Covariance Type: nonrobust
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 0.7065 0.184 3.833 0.000 0.344 1.069
## wc[T.True] 0.3353 0.075 4.450 0.000 0.187 0.483
## linc 0.1351 0.066 2.041 0.042 0.005 0.265
## k5 -0.0675 0.087 -0.775 0.439 -0.239 0.104
## ==============================================================================
## Omnibus: 82.950 Durbin-Watson: 2.001
## Prob(Omnibus): 0.000 Jarque-Bera (JB): 253.726
## Skew: -0.889 Prob(JB): 8.02e-56
## Kurtosis: 6.332 Cond. No. 17.6
## ==============================================================================
##
## Notes:
## [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
* Don't forget to install causaldata with ssc install causaldata
* if you haven't yet.
causaldata Mroz.dta, use clear download
* Keep just working women
keep if lfp == 1
* Get unlogged earnings
g earn = exp(lwg)
* Drop negative other earnings
drop if inc < 0
* 1. Draw a scatterplot
twoway scatter inc earn, yscale(log) xscale(log)
* 2. Get the conditional mean college attendance
* Note the syntax is different in older versions of Stata, see help table
table wc, stat(mean earn)
* 3. Get the conditional mean by bins
* Create the cut variable with ten groupings
egen inc_cut = cut(inc), group(10) label
table inc_cut, stat(mean earn)
* 4. Draw the LOESS and linear regression curves
* Create the logs manually for the fitted lines
g loginc = log(inc)
twoway scatter loginc lwg || lowess loginc lwg
twoway scatter loginc lwg || lfit loginc lwg
* 5. Run a linear regression, by itself and including controls
reg lwg loginc
reg lwg loginc wc k5第三章请退场吧，此处不再需要你的分析。↩︎
Emily Oster. Health recommendations and selection in health behaviors. American Economic Review: Insights, 2(2):143–60, 2020b.
↩︎本章将在\(Oster\)原研究的基础上，补充若干同质性的分析案例——这些新增内容虽非原始研究组成部分，但能更有效地阐释变量关系的描述方法。说来有趣，原作者竟将研究目的置于服务教科书范例之上，这般”本末倒置”着实令人莞尔。↩︎
不过与密度图不同，当数据量过大时，散点图的可读性会显著降低。正因如此，我在绘制该图时仅使用了150个观测值，而非全部数据。↩︎
这些分布亦可称为”边际分布”——尽管我个人强烈反对这一术语，因其字面含义与实际概念恰恰相反。↩︎
差异看似不大，但这正是社会科学中许多显著差异的典型表现。这种向右的偏移可能比看起来更具实质性！↩︎
为何选择均值？首先，均值在小样本中表现更稳定——当我们按X的特定值分组分析时，样本量往往会变小。其次，均值有助于加权预测误差，从而找到最小化误差的方法。简言之，它确实实用。↩︎
“比例”就是二元变量的均值。↩︎
关于LOESS与LOWESS的界定，学界存在两种观点：一种认为二者完全等同，另一种则指出其在局部预测的估计方法上存在细微差异——但实践中这两个名称常被混用于指代各类局部预测变体。↩︎
Locally Estimated Scatterplot Smoothing
↩︎为何该曲线未如图4.5那样在末端下降？这是因为极高\(BMI\)区间的观测值极少。\(LOESS\)算法不会让这少量观测值过度影响整体趋势，因而在某种程度上忽略了它们——这与图4.5的分组均值法形成鲜明对比。↩︎
更不必说，虽然并非完全不可行，但要用这些方法实现”条件中的条件均值”分析（参见”条件均值”章节）仍颇具难度。↩︎
严格而言，这两个概念并非完全等同，但在大多数应用场景中已足够近似。此外，尽管本节反复讨论条件均值，线性拟合同样存在可输出条件中位数、百分位数等其他统计量的变体方法。↩︎
但若该取值远离现有数据范围，则尝试求取其条件均值显然不妥。↩︎
.换言之，即实际观测值与模型预测值之间的差异。↩︎
微积分是方法之一，但除此之外↩︎
当前讨论仅针对双变量版本，更复杂的模型将在后续展开。↩︎
这种计算方式极强的直观性，或许能解释为何我们选择最小化残差平方和（而非四次方残差、乘积残差或绝对残差和）。尽管\(OLS\)因其限制性假设在统计界备受争议，但其无处不在的应用场景与广泛的理论关联——恕我直言——堪称多元统计方法中的”圆周率\(\pi\)“。其数学机理之精妙，甚至值得专章论述。看啊，我竟为一个比值如此着迷。↩︎
相关系数：衡量两个变量线性关联程度的标准化指标，其取值介于-1到1之间。直线斜率表示”Y单位变化/X单位变化”，其中”每(per)“即除法含义。乘以X的标准差（X单位）后，与分母的”每X单位”相抵消，仅余Y单位；再除以Y的标准差（Y单位），最终实现完全无量纲化。↩︎
本书对这两种方法的探讨将止步于满足研究设计需求的程度，更多细节可参阅第13章，或专项回归教材如《真实计量经济学》↩︎
残差：观测值的实际值与预测值之差。↩︎
这套特定计算方法在线性回归中被称为Frisch-Waugh-Lovell定理，但并不完全适用于参数非线性的回归模型（如前述\(logit/probit\)模型）。不过，其核心思想在这些模型中仍然适用。↩︎
控制时间变量仅能剔除时间解释的部分关联，却无法展现关系随时间的变化模式。↩︎
Thomas A Mroz. The sensitivity of an empirical model of married women’s hours of work to economic and statistical assumptions. Econometrica, 55(4):765–799, 1987.
↩︎为何在多数步骤中使用收入对数？请结合第3章对数运算的知识深入思考。↩︎