本章将完全关于如何描述一个变量。这似乎是一个奇怪的目标1。本书的开篇都在于提出研究问题以及实证研究如何帮助我们理解世界。然后我们就直接从那里跳到了描述变量?这是为什么呢?
事实证明,实证研究问题实际上完全归结为描述统计变量的密度分布。嗯,这实际上就是定量实证研究的全部内容了。抱歉。
也许这对我来说是个错误的方法。或许我应该这么说:你所听说过的所有有趣的实证研究成果——无论是在物理学、社会学、生物学、医学、经济学、政治学等等领域——都可以用一根线串起来。这根线被巧妙地铺展在大量的概率之上。它铺展时所呈现的形态就是一个统计变量的密度,它将整个宇宙中、永远存在的所有实证知识串联在了一起。
这样是不是更好了?我上《纽约时报》非虚构类畅销书排行榜榜首了吗?
听着,为了让我们所拥有的数据变得有意义,我们必须知道如何对一些观测结果进行分析并描述它们。我们做到这一点的方法是描述我们所拥有的变量类型以及这些变量所呈现的分布情况。这种描述的一部分内容将以描述不同变量之间如何相互作用的形式呈现,这部分内容会在第四章讲解。在本章中,我们将单独描述各个变量。这可能没有第四章的内容有趣,但很遗憾,它更为重要。
在实证研究的语境中,一个变量指的是对同一测量指标的一系列观测值——比如433名南非人的月收入、1984年至2014年期间法国每年发生的企业并购数量、对744名儿童进行访谈得到的心理“神经质”得分、532朵花的颜色、《华盛顿邮报》连续2348天的头版头条内容。成功地描述一个变量意味着能够依据这些观测值,清晰地解释所观测到的内容,而无需让别人自己去查看那744个神经质得分。这可比听起来要难。
弄清楚如何描述一个变量的第一步是弄清楚它是哪种类型的变量。
虽然总会有例外情况,但一般来说,你会遇到的最常见的变量类型有:
连续变量。 连续变量是可以取任何值的变量(可能在某个范围内)。例如,一位南非人的月收入就是一个连续变量。它可以是20,000 ZAR,或者它可以是34,123.32 ZAR,或者介于两者之间的任何值,或者从 0 到无穷大。不存在“下一个最高值”,因为变量会连续变化。20,000 ZAR 之后不是 20,001 ZAR2,因为 20,000.5 ZAR 介于两者之间。在到达那里之前,你必须经过 20,000.25 ZAR 和 20,000.10 ZAR,等等。
计数变量。 计数变量是那些,嗯,计算某些东西的变量。也许是某事发生了多少次或有多少东西在那里。法国在一年内的企业合并数量就是一个计数变量的例子。计数变量不能为负数,而且它们当然不能取分数。与连续变量相比,它们可能更难处理。有时,如果一个计数变量取许多不同的值,那么它的行为很像一个连续变量,因此研究人员经常将它们视为连续变量。
有序变量。 有序变量是指某些值“更多”,而另一些值“更少”的变量,但并不一定有规则说明“更多”到底多多少。一个具有“低水平神经质”、“中等水平神经质”和“高水平神经质”选项的“神经质”得分就是一个有序变量的例子。高于低,但高多少呢?这并不清楚。我们甚至不知道“低”和“中”之间的差异是否与“中”和“高”之间的差异相同。另一个可以更清楚地说明这一点的有序变量的例子是“最终完成的学业水平”,选项包括“小学”、“初中”、“高中”和“大学”。当然,完成高中学业意味着你比完成初中学业的人接受了更多的教育。但多多少呢? 是……多了两年学业吗? 实际上并不是这样运作的。它只是“更多”。 这就是一个有序变量。
分类变量。 分类变量是记录观察结果属于哪个类别的变量 - 非常简单!花的颜色就是一个分类变量的例子。这朵花是白色、橙色还是红色? 这些选项中没有一个比其他选项“更多”;它们只是不同。分类变量在社会科学研究中非常常见,在社会科学研究中,我们感兴趣的许多事物,例如宗教信仰、种族或地理位置,最好被描述为类别而不是数字。
分类变量的一种特殊形式是二元变量,它是只取两个值的分类变量。通常,这些值是“是”和“否”。例如,“你曾经服过兵役吗?” 是或否。“这只动物被给予药物了吗?” 是或否。二元变量很实用,因为它们比分类变量更容易处理,因为它们在询问治疗效果方面很有用(你接受治疗了吗?是或否),而且因为分类变量可以转化为一系列二元变量。我们的宗教信仰变量原本是包含 “基督教”“犹太教”“伊斯兰教” 等选项的分类变量,现在我们可以有一系列二元变量——“你是基督教徒吗?” 是或否。“你是犹太教徒吗?” 是或否。我们为什么要这样做呢?在阅读本书的过程中你会发现,这样做只是碰巧比较方便。此外,这样还能允许数据中的某个人既是基督教徒又是犹太教徒。
定性变量 定性变量是一种针对其他所有事物的包罗万象的类别。它们本质上不是数值型的,也不是分类的。《华盛顿邮报》的标题就是一个定性变量的例子。这些变量处理和描述起来可能非常棘手,因为这类变量往往包含大量难以精简和概括的细节。通常,为了概括这些变量,它们首先会被转换成上面提到的其他变量类型之一。例如,与其试图将《华盛顿邮报》的标题作为一个整体来描述,不如问“这个标题中提到总统多少次?”——这是一个计数变量——然后概括它。或者,既然我们现在生活在未来,我们可以让大型语言模型(AI)为我们阅读文本,并标记给定主题是否存在,或者根据某些特征对其进行评分。
一旦我们对所处理的是何种变量有了概念,下一步就是观察该变量的分布情况 。
一个变量的分布是对不同取值出现频率的描述。就是这么简单!例如,抛硬币结果的分布是正面朝上的概率为 50%,反面朝上的概率也为 50%。再比如,“一个人拥有的肢体数量”的分布是,大多数情况下这个数量是 4,而 0、1、2、3 以及 5 个及以上肢体的情况出现的频率相对较低。
当涉及到分类变量或顺序变量时,只需给出每个类别或取值中观测值所占的百分比,就可以描述该变量的分布情况。完整的分布可以用频数表或条形图来展示,而频数表或条形图展示的正是样本或总体中每个取值所占的百分比。
Table 3.1: Distribution of Kinds of Degrees US Colleges Award
library(conflicted)
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
library(vtable)
## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
library(purrr)
library(cowplot)
library(Cairo)
library(extrafont)
## Registering fonts with R
library(haven)
library(fGarch)
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")
library(ggpubr)
# import the windows' fonts
library(showtext)
## Loading required package: sysfonts
## Loading required package: showtextdb
##
## Attaching package: 'showtextdb'
##
## The following object is masked from 'package:extrafont':
##
## font_install
font_add('TNR','C:/Windows/Fonts/times.ttf')
showtext_auto()
## 加载数据
utils::data(Scorecard, package='pmdplyr')
# Categorical variable
catvar <- Scorecard %>%
ggpubr::group_by(unitid) %>%
dplyr::summarize(
# get the first value as all values are same
pred_degree_awarded_ipeds = dplyr::first(pred_degree_awarded_ipeds)
) %>%
# add column
dplyr::mutate(
`Primary Degree Type Awarded` = base::factor(
case_when(
pred_degree_awarded_ipeds == 1 ~ 'Less than Two-Year Degree',
pred_degree_awarded_ipeds == 2 ~ 'Two-Year Degree',
pred_degree_awarded_ipeds == 3 ~ 'Four Year or More'
),
levels = c('Less than Two-Year Degree', 'Two-Year Degree', 'Four Year or More')
)
) %>%
dplyr::select(`Primary Degree Type Awarded`)
vtable::st(catvar,
title = 'Distribution of Kinds of Degrees US Colleges Award',
# Character variable to be used to set an anchor link in HTML tables, or a label tag in LaTeX.
anchor = 'tab:describingvariables-frequencytable',
# comments
note = 'Data from College Scorecard',
# Determines where the completed table is sent
out='viewer',
# save file
file = 'DescribingVariables/frequency_table.html'
)
## NULL
这些表格告诉你所有你需要知道的信息。从表3.1我们可以看到,在我们的数据中的7424所大学中3,有3495所(47.1%)主要授予需要不到两年时间完成的学位,有1647所(22.2%)主要授予需要两年时间完成的学位,以及有2282所(30.7%)主要授予需要四年或更长时间完成的学位或高等学位。图3.1以图形格式显示了完全相同的信息 。
Figure 3.1: Distribution of Kinds of Degrees US Colleges Award
# 绘图
ggplot2::ggplot(catvar %>%
ggpubr::group_by(`Primary Degree Type Awarded`) %>%
dplyr::summarize(N=n()) %>%
dplyr::mutate(`Primary Degree Type Awarded`= base::factor(`Primary Degree Type Awarded`,
levels=c('Less than Two-Year Degree','Two-Year Degree','Four Year or More'))),
aes(x=`Primary Degree Type Awarded`, y=N)
) +
ggplot2::geom_col(fill='white', color='black') +
ggplot2::geom_text(aes(label=`Primary Degree Type Awarded`, x=`Primary Degree Type Awarded`, y=N+70),
vjust=0,
family='sans') +
ggpubr::theme_pubr() +
# label value
ggplot2::labs(x='Primary Degree Type Awarded',y='Number of Colleges') +
ggplot2::theme(text=element_text(size = 10, family="sans"),
axis.text.x = element_blank(),
# the size and family
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans"),
panel.grid = element_blank(),
axis.ticks.x = element_blank())
ggsave('DescribingVariables/college_type_distribution.pdf', width = 6, height = 5, device = cairo_pdf)
需要考虑的可能性只有这么多,而表格和图表分别展示了每种可能性出现的频率。一旦完成这一步,你就完整地描述了该变量的分布情况——这个变量能提供的信息真的就只有这些了!如果我们想展示更多细节(比如每所大学倾向于专攻哪些具体专业),就需要另一个包含不同变量的数据来源。
连续变量的处理要复杂一些。由于不同观测值几乎不会完全重合,我们无法直接为连续变量制作频数表。例如,虽然24,201南非兰特的收入与24,202兰特极为接近,但两者并不相等,因此在频数表或条形图中仍会对应不同的位置。
对于连续变量,分布的描述并非基于变量取某个确定值 的概率,而是基于变量落在该值附近 的概率。
表达连续变量分布的一种常见方法是使用直方图 。直方图将数据的可能范围划分为若干个区间(\(bin\)),并显示落入每个区间的观测值的比例。这与我们用于分类变量的频数表 或频数图 完全相同,唯一的区别在于,直方图的类别是变量的取值范围,而不是它可能取的全部具体值。
例如,图3.2展示了大学毕业生毕业几年后收入的分布情况,每个数据点代表一个院校的每届毕业生(“队列”)。从图中可以看出,有超过2万个院校队列的毕业生年收入在2万至4万美元之间。收入在1万至2万美元之间的队列较少,大约有4000个。极少数院校队列的平均收入在8万至16万美元之间。而收入介于16万至32万美元的队列数量极少,在图中几乎无法观察到。
Figure 3.2: Distribution of Average Earnings across US College Cohorts
# Histogram/density
ratiovec <- function(n,top) {
rate <- top^(1/(n-1))
v <- c(1)
for (i in 2:n) {v[i] <- v[i-1]*rate}
return(v)
}
ratiovec(6,32)
## [1] 1 2 4 8 16 32
# x轴值
ggplot2::ggplot(Scorecard, aes(x = earnings_med)) +
# 直方图
ggplot2::geom_histogram(color='Blue',fill='white',breaks=ratiovec(6,32)*10000, na.rm=TRUE) +
# x轴刻度标签
ggplot2::scale_x_log10(labels = scales::label_dollar(),breaks=ratiovec(6,32)*10000) +
# x轴和y轴对应标签
ggplot2::labs(x='Average Earnings of This College & Cohort\'s Graduates (Log Scale)',y='Number of College Cohorts') +
ggpubr::theme_pubr() +
ggplot2::theme(text=element_text(size = 10, family="sans", color='Blue'),
axis.title.x = element_text(size = 10, family="sans", color='Blue'),
axis.title.y = element_text(size = 10, family= "sans", color='red'),
panel.grid = element_blank())
ggsave('DescribingVariables/college_earnings_histogram.pdf', width=6, height=5, device=cairo_pdf)
对于一个连续变量,我们可以更进一步,从直方图一直到密度图4。密度图展示了如果直方图的组距变得越来越窄会发生什么,正如你在图3.3中看到的那样5。
Figure 3.3: Distribution of Average Earnings across US College Cohorts
p <- c(6,15,30) %>%
map(function(x)
ggplot2::ggplot(Scorecard, aes(x=earnings_med)) +
# set the breaks
ggplot2::geom_histogram(color='black',fill='white',breaks=ratiovec(x,32)*10000, na.rm=TRUE) +
ggplot2::scale_x_log10(labels=scales::label_dollar(),limits = c(10000,200000)) +
ggplot2::labs(x = 'Average Earnings (Log Scale)',y = 'Number of College Cohorts') +
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()))
p[[4]] <- ggplot2::ggplot(Scorecard, aes(x = earnings_med)) +
ggplot2::geom_density(na.rm=TRUE) +
ggplot2::scale_x_log10(labels = scales::label_dollar(),limits = c(10000,200000)) +
ggplot2::labs(x = 'Average Earnings of This College \n& Cohort\'s Graduates (Log Scale)',y = 'Density') +
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())
plot_grid(p[[1]],p[[2]],p[[3]],p[[4]], nrow=2)
ggsave('histogram_to_density.pdf',width = 7, height = 6, units='in', device=cairo_pdf)
当我们有一个密度图时,我们可以通过查看分布曲线下的面积大小来描述变量处于给定范围内的概率。例如,图3.4展示了收入的分布情况,并对4万美元至5万美元之间的区域进行了阴影标注。相对于整个分布曲线下的面积大小而言,这个阴影区域的面积就是收入处于4万美元至5万美元之间的概率。这个特定的阴影区域约占曲线下总面积的16%,因此,所有群体中有16%的群体平均收入在4万美元至5万美元之间6。
Figure 3.4: Shaded Distribution of Earnings across US College Cohorts
densdata <- base::with(stats::density(log(Scorecard$earnings_med),na.rm=TRUE),
base::data.frame(x=exp(x),y))
ggplot2::ggplot() +
# generate the gray ribbon
ggplot2::geom_ribbon(data=densdata %>% dplyr::filter(x > 40000 & x < 50000),
mapping=aes(ymax=y,ymin=0,x=x),
alpha=.5,
color='gray') +
# draw a line
ggplot2::geom_line(data=densdata, mapping=aes(x=x, y=y), na.rm=TRUE) +
ggplot2::scale_x_log10(labels = c('$40k','$50k'),limits = c(10000,200000),breaks=c(40000,50000)) +
ggplot2::labs(x = 'Average Earnings of This College & Cohort\'s Graduates (Log Scale)',y='Density') +
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"),
axis.text.x = element_text(
# angle = c(0,0,20,0,0,20),
size=10,
# hjust = c(.5,.5,.7,.5,.5,.7)
family='sans'
),
panel.grid = element_blank())
ggplot2::ggsave('DescribingVariables/shaded_density.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
就这样吧。一旦你掌握了变量的分布情况,实际上这就是你能说的全部了7。毕竟,我们拥有什么呢?对于变量可能取的每一个可能值,我们知道该结果发生的可能性有多大,或者至少是类似结果发生的可能性有多大。关于一个变量,你还能说什么呢?
当然,在很多情况下,这些分布细节过于丰富,难以完整呈现。诚然,对于只有几个类别的分类变量,我们可以轻松展示完整的频数表。但对于任何类型的连续变量,即便我们向某人展示了密度图,他们也很难全盘吸收其中的所有信息。
那么,对于这种分布,我们可以怎么做才能让它更易于理解呢?我们选取它的几个关键特征并告知你这些特征。换句话说,我们对它进行总结。
# That particular shaded area makes up about 16% of the area underneath the curve,
# and so 16% of all cohorts have average earnings between $40,000 and $50,000
mean(Scorecard$earnings_med >= 40000 & Scorecard$earnings_med <= 50000, na.rm =TRUE)
## [1] 0.1603287
一旦我们掌握了变量的分布情况,就可以将注意力转向对该变量进行总结。整个分布所包含的信息可能过多,让我们难以有效利用,尤其是对于连续变量而言。因此,我们的目标是选取合适的方法,对整个分布进行处理,从而得出几个能够很好地描述该分布的数值。
可能最广为人知的一个试图描述整个分布的单一数字的例子是平均数。平均数就是你将所有观察值加起来,然后除以观察值的数量所得出的。因此,如果你有2、5、5和6,那么平均数就是(2+5+5+6)/4 = 18/4 = 4.5。
更正式一点来说,均值的计算方式是,它会取你可能得到的每个数值,用你得到该数值的可能性对其进行加权,然后将所有结果相加,依此类推!对于我们包含 2、5、5 和 6 的数据集,其分布情况是怎样的呢?我们的频数表如表 3.2 所示。
Table 3.2: Distribution of a Variable
a <- base::data.frame(`Observed Values`=factor(c(2,5,5,6)))
vtable::st(a, title='Distribution of a Variable',
anchor='tab:describingvariables-minifrequency',
out='viewer',
file='DescribingVariables/mini_frequency_table.html')
## NULL
Table 3.2 给出了我们变量的分布。现在我们可以计算均值。再次, 这根据我们获得每个值的可能性来缩放每个值。在 2、5、5 和 6 中,我们得到 5 的一半,所以我们计算 5 的一半。我们得到 2 的四分之一的,所以我们计算 2 的四分之一。
好的,所以,由于 2 只出现 25%(1/4),我们只计算 2 的 1/4,得到 .5。 接下来,5 出现 50%(1/2),所以我们计算 5 的一半,得到 2.5。我们看到 6 出现 也是 25% (1/4) ,所以我们将 6 缩放 1/4,得到 1.5。将它们全部加起来得到我们的 均值是 .5 + 2.5 + 1.5 = 4.5。
所以,均值 实际上所做的 是观察变量的 分布情况 ,并对其进行总结,将其归结为一个单一的数字。那个数字是什么呢?均值应该代表数据的集中趋势 —— 它处于你可能得到的值的中间位置。更具体地说,它试图得出一个有代表性的 数值。如果变量是“这台老虎机的赔付金额(美元)”,且均值为 4.50 美元,而玩一次老虎机的费用是 4.50 美元,那么如果你玩这台老虎机很多次,你将恰好不赚不赔。
有时候,直截了当地表达更有好处。我们当然可以用均值来描述一个分布,而且在这本书里,我们会多次这么做。但如果目标是向别人描述这个分布,那我们为什么还要费劲去计算均值,而不是直接把分布的情况告诉别人呢?
一个变量分布的第 X 百分位数是这样一个值,在该分布中,有 X%的观测值小于这个值。例如,如果你让 100 个人按身高排队,如果排在队列中前面有 5 个人的那个人身高是 5 英尺 4 英寸,那么 5 英尺 4 英寸就是第 5 百分位数。
我们可以在分布图上看到百分位数。图 3.5 展示了我们之前的大学同届毕业生收入分布情况。我们从分布曲线的左侧开始进行阴影填充,然后持续填充,直到曲线下方 5%的区域都被阴影覆盖。我们停止填充时在 x 轴上对应的点——16400 美元,就是第 5 百分位数。
Figure 3.5: Distribution of Average Earnings across US College Cohorts
ggplot2::ggplot() +
ggplot2::geom_ribbon(data=densdata %>% dplyr::filter(x < quantile(Scorecard$earnings_med,.05, na.rm=TRUE)),
mapping=aes(ymax=y,ymin=0,x=x),
alpha=.5,
color='blue') +
ggplot2::geom_line(data=densdata, mapping=aes(x=x, y=y), na.rm=TRUE, color='red') +
ggplot2::scale_x_log10(
labels=function(x) scales::dollar(x, accuracy=1),
limits=c(10000,200000),
breaks=c(quantile(Scorecard$earnings_med,.05, na.rm=TRUE))
) +
ggplot2::annotate(geom='text',x=15000,y=.1,label='5%',family='sans', size = 10/.pt) +
ggplot2::annotate(geom='text',x=30000,y=.1,label='95%',family='sans', size = 10/.pt) +
ggplot2::labs(x = 'Average Earnings of This College & Cohort\'s Graduates (Log Scale)',y = 'Density') +
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())
ggplot2::ggsave('DescribingVariables/fifth_percentile.pdf',width=7, height=6, units='in', device=cairo_pdf)
var(Scorecard$earnings_med,na.rm=TRUE)
## [1] 153287962
sd(Scorecard$earnings_med,na.rm=TRUE)
## [1] 12380.95
mean(Scorecard$earnings_med,na.rm=TRUE)
## [1] 33347.62
(38000-33348.62)/12380.95
## [1] 0.3756885
我们实际上可以用这种方式完美地描述整个分布。什么是第一个百分位数?好的,现在什么是第二个?等等8。很快,我们将通过每次只多一点阴影来绘制出整个分布。所以百分位数是描述分布的一种相当直接的方式。
有几个百分位数值得特别提及。
第一个是 中位数,即第 50 百分位数。处于这个位置的人正好处于中间——样本中有一半的人比他们高,另一半比他们矮。和均值一样,中位数也是在衡量数据的集中趋势。均值试图得出一个具有代表性的 数值,而中位数则给出一个具有代表性的 观测值。
例如,假设你正在研究一万个人的财富状况,其中有一位是亚马逊创始人杰夫·贝索斯。均值会“这么认为”:“嗯……没错,大多数人没什么财富,但偶尔会出现像杰夫·贝索斯这样的人,这就把整体水平拉高了。所以均值会非常高。”但中位数则表示:“杰夫·贝索斯并不能很好地代表其他所有人。他和其他每个人在排序中所占的比重是一样的。所以中位数相对较低9。”
因此,当您想要描述典型的观测值是什么样子,或者当您有一个高度偏斜的变量,其中包含一些非常大的观测值,例如财富和杰夫·贝佐斯时,通常会使用中位数而不是均值。那个房间的平均财富可能是 15,000,000 美元,但几乎都是杰夫·贝佐斯,并不能真正代表其他人。一旦他走出房间,均值就会像石头一样下降。均值对杰夫非常敏感!但中位数可能更接近 90,000 美元,这对于一个美国家庭来说是一个相当典型的净资产10,而且如果杰夫离开房间,它几乎会保持完全相同。中位数非常适合这样的事情11!
另外两个需要关注的百分位数是 最小值,即第 0 百分位数,以及 最大值,即第 100 百分位数。它们分别是该变量所取的最低值和最高值。它们很有用,因为它们能让你了解该变量会产生哪些类型的值。例如,一大群人的身高最小值和最大值能让你了解人类可能的身高范围,即最矮和最高能达到多少。
最小值和最大值的另一个优点是,我们可以通过计算它们之间的差值来得到 极差。极差是衡量一个变量 变化程度 的一种方式。如果最大值和最小值相距很远,就像房间里有杰夫·贝索斯(世界顶级富豪)时财富的情况那样,你就知道这个变量可以取到非常广泛的值。如果最大值和最小值很接近,就像“你拥有的眼睛数量”这种情况,你就知道这个变量所能取到的值的范围相当小。
顺便说一句,从技术上讲,在百分位数方面存在一些歧义。如果100人排成一列,那么第5个百分位数是队列中的第5个人(因为5%的人处于或低于他们的水平,即“包含式”百分位数),还是队列中的第6个人(因为5%的人低于他们的水平,即“排除式”百分位数)?如果选择的百分位数介于两个人之间呢?在一个由10人组成的群体中,没有人处于第25个百分位数。那么,取队列中的第2个人?第3个?通常你会以某种方式混合这两个,但以什么方式呢(Hyndman和Fan(1996)讨论了九个(!)选项)。在大的样本中,这些都不重要,但在较小的样本中,这些决定可能很重要,你应该考虑什么对你的情况有意义,以及你的软件包默认使用哪个选项。
names <- c('Very Low Variation','A Little More Variation','Even More Variation','Lots of Variation')
set.seed(1000)
# Variance is a measure of variation that focuses on values
dfs <- c(.2,.25,.5,.75) %>% map(function(x) data.frame(x=rnorm(100,sd=x)))
p <- 1:4 %>%
purrr::map(function(x)
ggplot2::ggplot(dfs[[x]],aes(x=x))+
ggplot2::geom_density()+
ggplot2::scale_x_continuous(limits=c(min(dfs[[4]]$x),max(dfs[[4]]$x)))+
ggplot2::scale_y_continuous(limits=c(0,max(density(dfs[[1]]$x)$y)))+
ggplot2::labs(x='Observed Value',y='Density',title=names[x])+
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()))
cowplot::plot_grid(p[[1]],p[[2]],p[[3]],p[[4]], nrow = 2)
ggplot2::ggsave('DescribingVariables/amount_of_variation.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
有些变量的变化幅度小,有些则变化幅度大。例如,以“一个人拥有的孩子数量”为例。对很多人来说,这个数量是零。对于三十多岁的人而言,这个数量的平均值可能在 2 左右。不过,有些人会有很多很多孩子。少数罕见的女性会有十个或更多孩子。世界上有些男性有几十个孩子。一个人拥有的孩子数量可能会有相当大的差异。
与”人的眼睛数量”形成鲜明对比:少数个体可能拥有0只、1只甚至3只眼睛,但总体分布高度集中——约99.7%的人口具有双眼[数据来源:WHO]。该变量的变异系数不足0.5%,堪称生物学中最稳定的形态特征之一。
“子女人数”与”眼睛数量”这两个变量虽拥有相近的均值和中位数,但变异特性截然不同。这说明仅靠集中趋势或百分位数不足以完整描述变量特征,我们还需量化变异程度的指标。
分布图的形态直观展现变异程度:若呈现”高瘦型”,则观测值紧密聚集在均值附近,变异度低;若呈”扁平宽阔型”,则两侧”肥尾”包含大量远离均值的观测值,变异度高(参见图3.6)。具体而言,中部”堆积”面积越大的分布,其变异度越小——远离中心的数据点极少;反之,中部堆积少而尾部厚的分布,则包含大量偏离中心的观测值。
Figure 3.6: Four Variables with Different Levels of Variation
描述变异的方法多样:均值等指标侧重数值特征,而中位数、百分位数等则聚焦观测值的序位特性。
方差作为衡量变异度的指标,以均值为基准聚焦数值特征,其计算步骤如下:’
计算均值:数据组2, 5, 5, 6的均值为 (2+5+5+6)/4 = 4.5
求离均差:各数值减均值得-2.5, 0.5, 0.5, 1.5(即各数据点与均值的偏差)
平方处理12:得6.25, 0.25, 0.25, 2.25
求和:6.25+0.25+0.25+2.25=9
除以n-113:样本方差=9/(4-1)=3
方差值越大,变量变异程度越高。其原理在于:步骤4-5实质是对”离均差平方”求均值——远离均值的观测值经平方运算后权重倍增,从而反映数据整体的平均离散程度。
方差的局限在于其”平方单位”难以直观解释——例如大学毕业生收入变量的方差为153,287,962(美元平方?)。为此我们常通过开平方运算将方差转化为标准差(√153,287,962=12,380.95美元),使其单位与原始变量一致。
已知均值收入为33,348.62美元时,某毕业生群体平均收入38,000美元意味着:(38,000-33,348.62)/12,380.95=0.376,即高于均值0.376个标准差。这不仅显示该群体收入绝对值(高4,651.38美元),更揭示其相对位置——超出典型变异范围的37.6%。
Figure 3.7: Distribution of Average Earnings across US College Cohorts
logd <- Scorecard %>% dplyr::mutate(lx=log(earnings_med)) %>% dplyr::filter(!is.na(lx))
ggplot2::ggplot(logd, aes(x=lx)) +
ggplot2::geom_density(na.rm=TRUE) +
# (mean-1sd, mead+1sd)
ggplot2::geom_vline(aes(xintercept=mean(lx))) +
ggplot2::geom_vline(aes(xintercept=mean(lx)-sd(lx)),linetype='dashed') +
ggplot2::geom_vline(aes(xintercept=mean(lx)+sd(lx)),linetype='dashed') +
ggplot2::scale_x_continuous(
labels = scales::dollar(exp(c(mean(logd$lx),mean(logd$lx)-sd(logd$lx),mean(logd$lx)+sd(logd$lx))), accuracy=1),
limits = c(log(10000),log(200000)),
breaks=c(mean(logd$lx),mean(logd$lx)-sd(logd$lx),mean(logd$lx)+sd(logd$lx))
) +
ggplot2::annotate(geom='label',x=mean(logd$lx)+.01,y=1.25,label='Mean',family='sans',hjust=.5) +
ggplot2::annotate(geom='label',x=mean(logd$lx)-sd(logd$lx)+.01,y=1.25,label='Mean - 1 SD',family='sans',hjust=.5) +
ggplot2::annotate(geom='label',x=mean(logd$lx)+sd(logd$lx)+.01,y=1.25,label='Mean + 1 SD',family='sans',hjust=.5) +
ggplot2::labs(x='Average Earnings of This College & Cohort\'s Graduates (Log Scale)',y='Density') +
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())
ggplot2::ggsave('DescribingVariables/standard_deviation.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
理解”一个标准差”的变异幅度需要实践与直觉,但图示可辅助理解。如图3.7所示,“均值±1个标准差”范围覆盖了样本的68.2%(32.7%在均值左侧,35.5%在右侧),意味着超过60%的样本点集中在距均值不足一个标准差的范围内。
mean(logd$lx < mean(logd$lx)) - mean(logd$lx < mean(logd$lx)-sd(logd$lx))
## [1] 0.3276826
mean(logd$lx < mean(logd$lx)+sd(logd$lx)) - mean(logd$lx < mean(logd$lx))
## [1] 0.3551727
与均值偏离一个标准差究竟有多特殊?简言之,约三分之一样本分布在你与均值之间——其意义可自行解读。
通过对比百分位数也能衡量变量变异度,该方法操作简便:只需选取中位数上方和下方的对称百分位数(如第25与第75百分位数),计算其差值即可。
极差(最大值与最小值之差)虽能反映变异幅度,但对极端值异常敏感——以财富分布为例,杰夫·贝索斯是否在场会导致极差天差地别,故其并非稳健的变异度指标。
最常用的百分位变异度指标是四分位距(IQR)14即第75与第25百分位数之差。其优势在于:1) 涵盖样本最接近中位数的50%数据,能准确反映核心分布密度;2) 不像方差那样易受尾部极值影响。这使其成为刻画观测值序位特征的稳健指标。
图3.8展示了四分位距(IQR)在分布中的位置。本例中,第25百分位数为21,761美元,第75百分位数为44,879美元,故IQR=44,879-21,761=23,118美元。这意味着最接近中位数的50%毕业生群体,其平均收入波动范围为23,118美元。
Figure 3.8: Distribution of Average Earnings across US College Cohorts
logd <- Scorecard %>% dplyr::mutate(lx=log(earnings_med)) %>% dplyr::filter(!is.na(lx))
ggplot2::ggplot(logd, aes(x=lx)) +
ggplot2::geom_density(na.rm=TRUE) +
ggplot2::geom_vline(aes(xintercept = median(lx))) +
ggplot2::geom_vline(aes(xintercept = quantile(lx,.25)),linetype = 'dashed') +
ggplot2::geom_vline(aes(xintercept = quantile(lx,.75)),linetype = 'dashed') +
ggplot2::scale_x_continuous(
labels = scales::dollar(exp(c(mean(logd$lx),mean(logd$lx)-sd(logd$lx),mean(logd$lx)+sd(logd$lx))),accuracy=1),
limits = c(log(10000),log(200000)),
breaks=c(median(logd$lx),quantile(logd$lx,c(.25,.75)))
) +
# add annotate
ggplot2::annotate(geom='label',x=median(logd$lx),y=1.25,label='Median',family='sans',hjust=0.5) +
ggplot2::annotate(geom='label',x=quantile(logd$lx,.25),y=1.15,label='25th Percentile',family='sans',hjust=.7) +
ggplot2::annotate(geom='label',x=quantile(logd$lx,.75),y=1.15,label='75th Percentile',family='sans',hjust=.3) +
# plot the arrow
ggplot2::annotate('segment',
x=quantile(logd$lx,.25)+.02,
xend=quantile(logd$lx,.75)-.02,
y=.4,
yend=.4,
arrow=arrow(ends='both',type='closed',length=unit(.1,'inches')))+
ggplot2::annotate(geom='label',x=median(logd$lx),y=.4,label='IQR',family='sans',hjust=.5)+
ggplot2::labs(x = 'Average Earnings of This College & Cohort\'s Graduates (Log Scale)',y='Density') +
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"),
axis.text.x = element_text(angle=20, hjust = .7),
panel.grid = element_blank())
ggsave('DescribingVariables/iqr.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
除变异度外,分布特征尚有诸多可描述维度。本书将重点探讨其中一个关键特征——偏态。
偏态用于描述分布向一侧倾斜的程度。以年收入为例:多数人的收入集中在相对狭窄区间(0至15万美元),但仍有相当数量人群的收入远超此范围——甚至达到数百万、数千万乃至数亿美元。
就年收入而言,其右尾(即分布曲线右侧部分)显著延伸。如图3.9所示,尽管大部分数据点集中在0值附近,但部分人群的百万美元级收入使得右尾极度延长。而左尾则不存在对称情况——至少在当前数据集中,未出现负收入观测值。
Figure 3.9: Distribution of Personal Income in 2018 American Community Survey
## Income distribution
df <- haven::read_stata(gzfile('DescribingVariables/usa_00016.dta.gz'))
## get the value (>1)
df_noz <- df %>% dplyr::filter(inctot > 1)
## use the perwt as the weight
ggplot2::ggplot(df_noz,aes(x=inctot,weight=perwt)) +
ggplot2::geom_density() +
ggplot2::labs(x='Personal Income',y='Density') +
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"),
axis.text.y = element_blank(),
panel.grid = element_blank(),
axis.ticks.y = element_blank())
ggplot2::ggsave('DescribingVariables/income_distribution.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
我们称此类分布具有”右偏态”——因其存在显著右尾(大量极端大值)但无明显左尾。同理,若左尾存在大量异常观测值则称为左偏态。当两侧尾部对称时,则称该分布为对称分布。
右偏态变量在社会科学研究中极为常见。任何存在不均衡分布的变量(如收入)通常呈现此特征:多数人拥有量较少,而少数人持有量极大——且这些少数个体的数值往往异常突出。
偏态是描述分布特征的重要维度,但若涉及均值与方差计算则可能引发问题——因极端值会显著影响基于数值的统计量。
处理数据偏态的常用方法之一是数据变换——通过对数据施加某种函数转换,削弱极端大值对统计量的影响,从而使均值与方差更具代表性。自然对数变换便是典型方案:取数据的自然对数替代原值,可显著改善数据分布特性15。
如图3.10所示,收入数据经对数变换后仍存在轻微尾部(此次出现在左侧!),但整体已呈现近似对称分布——此时均值作为统计量将更具代表性16。
Figure 3.10: Distribution of Logged Personal Income in 2018 American Community Survey
ggplot2::ggplot(df_noz,aes(x=log(inctot),weight=perwt)) +
ggplot2::geom_density() +
ggplot2::labs(x='Log Personal Income',y='Density') +
ggplot2::scale_x_continuous(labels=scales::number)+
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"),
axis.text.y = element_blank(),
panel.grid = element_blank(),
axis.ticks.y = element_blank())
ggplot2::ggsave('DescribingVariables/log_income_distribution.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
自然对数变换广受青睐的原因之一在于其直观的可解释性:对数变量的单位增量(如0.01)对应原始变量约0.01×100=1%的增幅。例如,当对数收入从10.00增至10.02时,(10.02-10.00)=0.02≈2%即表示实际收入增长约2%。该近似换算在0.01或0.02等微小增量时效果良好,但当增量超过0.1(约10%)时误差显著,此时应避免使用此近似方法17。
若数据经对数变换后仍呈现偏态,则可能隐含严重问题——这常见于存在”厚尾”特征的数据集(即远离均值的极端值频繁出现)。当单侧厚尾现象显著时,数据分析将变得异常棘手。本书第22章将对此进行专题探讨。
现实与真理的区别在于:现实始终如影随形,而真理纵使跋涉千里仍在地平线外。因此,人们往往选择眯眼耸肩,嘟囔着”差不多得了”,然后打道回府——这显然省事得多。
统计学严格区分”真理”与”观测数据”——尽管数据本身确实是实际采集的结果,但其终极意义在于揭示背后更广义的真理。
假设你想探究儿童分享玩具的平均年龄。通过对1000名家长进行访谈,统计得出样本中儿童开始主动分享的平均年龄约为4.2岁。
实际上,你获得的仅是样本中1000名儿童的平均分享年龄(4.2岁)。但研究初衷绝非仅了解这1000个孩子——而是要推断全体儿童的普遍规律。因此,总体真实均值与样本计算均值本质上是两个不同概念。
这正是统计学的核心价值——我们无法普查所有儿童的分享行为发展时点,但基于现有样本数据,能否推断总体真实参数?
解决这一问题需要思考数据在不同真实参数下的表现形态:若儿童平均分享年龄的真实值为3.8岁,会生成何种数据分布?若真实值为5岁,又将产生怎样的数据模式?为此,我们需要将本章讨论的观测分布,与不同真实参数下的理论分布进行配对分析。
在深入探讨前,需明确关键统计符号:\(\beta\)、\(\mu\)、\(\hat\mu\)、\(\bar x\)等常见符号虽需记忆,但其命名体系实则存在内在逻辑规律。这些符号究竟代表什么?
英文字母(拉丁字母)表示实际观测数据:如变量\(x\)可代表1000份家长问卷中的儿童分享年龄数据。
字母修饰符号表示基于真实数据的计算结果:常用上划线”¯“表示均值,故\(\bar x\)即样本\(x\)的均值——如问卷数据计算得出的4.2岁。
希腊字母表征理论真值1819—虽其具体数值未知,但可设定假设。特定希腊字母有固定含义:\(\mu\)通常表示总体均值20,\(\sigma\)代表标准差,\(\rho\)指相关系数,\(\beta\)为回归系数,\(\epsilon\)表示误差项(后续详述)。其核心要义在于:希腊字母始终代表理论层面的真实参数。
希腊字母加修饰符表示对真值的估计——尽管真实参数未知,我们仍可进行最优推测(无论其准确性如何)。最常见的估计表示法是在希腊字母上方加”^“符号:如\(\hat\mu\)表示”对\(\mu\)的估计值”。若采用\(x\)的均值估计\(\mu\),则可表述为\(\hat\mu\)= \(\bar x\)。
理论分布是生成观测数据的底层机制——这一认知视角颇具启发性:理论分布囊括了所有潜在数据的分布特征,包括未实际采集(甚至理论上永远无法采集)的数据。若能获得无限量观测值,其分布形态即为理论分布21。
这一事实向我们揭示了几点要义。首先,它解释了为何我们最初就对理论分布如此关注——因为那正是我们数据的源头!若想探究儿童平均分享年龄的规律,理论分布才是这一数据的真正出处。因此,当我们需要透过既有数据窥见更深层的真相时,就必须凭借这些数据逆向追溯至理论分布。唯有如此,我们才能触及真正有趣的发现!
需要牢记的是:我们真正关心的并非样本均值\(\bar x\),而是总体真实均值\(\mu\)!之所以费心收集数据,正是为了通过估计量\(\hat\mu\)来推断其背后的理论分布特征。
这个”无限观测”事实给我们的第二个启示是:观测数据越多,就越能匹配理论分布。单次观测价值有限,而无限次观测将完美呈现理论分布。现实中我们只能取其中间值——随着观测次数增加,匹配精度将不断提升。
这一点在图3.11中清晰可见。无论观测数据量多少,实线表示的理论分布自然始终保持不变。虽然仅凭10次观测数据难以准确描述该分布,但到100次时拟合效果已显著改善,而1,000次时则已相当吻合。当然,这并非意味着1,000次就”足以完全匹配理论分布”,但在本例中确实效果良好。
Figure 3.11: Trying to Match the Theoretical Distribution
library(tidyverse)
library(ggpubr)
## Normal distribution matching with bigger sample
nobs <- c(10,50,100,1000)
set.seed(500)
p <- nobs %>%
purrr::map(~ {
# 模拟数据并计算密度
x_vals <- stats::rnorm(.x)
dens <- stats::density(x_vals)
df <- base::data.frame(x=dens$x, y=dens$y)
# 添加理论密度值(标准正态分布)
df <- df %>% dplyr::mutate(yn=dnorm(x))
# 绘图
ggplot2::ggplot(df, aes(x=x, y=y)) +
ggplot2::geom_line(linetype='dashed', color="blue", linewidth=1) +
ggplot2::geom_line(aes(y=yn), linetype='solid', color="red", linewidth=1) +
ggplot2::labs(x='Value of Observation',y='Density',title = paste(.x, 'Observations')
) +
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()
)
})
plot_grid(p[[1]],p[[2]],p[[3]],p[[4]], nrow=2)
ggsave('DescribingVariables/approach_limit.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)
这意味着随着观测数据的不断增加,我们得到的经验分布将越来越接近抽样的理论分布。由于后者才是我们真正关心的目标分布,这无疑是件好事!关键在于确保获取足够数量的观测数据22。
理论分布有无限多种,但在实际应用中某些类型会反复出现。一些经典分布被不断重复使用,若数据符合这些分布特征,我们就能借助现成的理论分布完成大量分析工作——这无疑是研究者的福音!
在应用社会科学研究中,我将重点介绍两种至关重要的分布(均如图3.12所示)。当然还有众多未被涵盖的分布类型:均匀分布、泊松分布、二项分布、伽马分布、贝塔分布等等。若您感兴趣,不妨参考更专业的统计学著作——比如在亚马逊搜索”统计学导论”时弹出的那八百万本教材。我确信它们几乎都和本书一样出色。几乎。
首先要介绍的是正态分布。这种对称分布(即左右尾端等长,不存在偏斜)常用于描述存在实际物理限制的特征,比如身高或智力水平。
正态分布在汇总统计量分析中也频繁出现23。虽然原始收入数据可能存在严重右偏,但当我们反复从不同样本中计算收入均值时,这些样本均值的分布将呈现正态特征。具体而言:在第一个样本中计算观测值的收入均值,在第二个样本中计算另一组观测值的收入均值,如此不断重复,最终不同样本均值的分布会趋近正态分布。
从技术上讲,正态分布具有无限延伸的特性——这意味着任何数值在理论上都可能出现,即便概率极低。因此,对于存在取值限制的变量(如高度不可能为负值),严格来说并不符合正态分布。但只要近似程度足够理想,我们通常选择忽略这一差异。
这种宽容态度部分源于正态分布的”薄尾”特性——远离均值的观测值出现概率极低。如图3.12所示,曲线两端会迅速趋近于零。因此,当声称身高服从正态分布时,虽然从技术层面暗示了负值存在的可能性,但实际概率可能低至0.00000001%,这种误差范围在实践层面完全可以接受。
Figure 3.12: Normal and Log-Normal Distributions
# Normal and log-normal
set.seed(500)
df <- base::data.frame(x=seq(-4,4,length.out=100000)) %>% dplyr::mutate(y=dnorm(x),expy=exp(x))
p1 <- ggplot2::ggplot(df, aes(x=x,y=y)) +
ggplot2::geom_line() +
ggplot2::labs(x='Value',y='Density', title='Normal Distribution') +
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"),
axis.text.y = element_blank(),
panel.grid = element_blank(),
axis.ticks.y = element_blank())
p2 <- ggplot2::ggplot(df, aes(x=expy)) +
ggplot2::geom_density() +
ggplot2::labs(x='Value',y='Density', title='Log-Normal Distribution') +
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"),
axis.text.y = element_blank(),
panel.grid = element_blank(),
axis.ticks.y = element_blank())
plot_grid(p1,p2, nrow = 1)
ggsave('DescribingVariables/theoretical.pdf',width = 7, height = 3, units = 'in', device = cairo_pdf)
第二个要介绍的是取巧的”对数正态分布”。这种分布呈现显著右偏特征,但对其取对数后即转化为正态分布——多么便捷的性质!
对数正态分布作为偏态分布中极具实用价值的类型,仅需对数变换即可消除偏斜特性。现实世界中存在大量重偏数据:无论是无上限的不均衡分布(如收入、财富),还是呈现”赢家通吃”特征的现象(如歌曲下载量、游戏时长、企业规模等),都普遍存在偏斜。观察其分布形态,你会发现绝大多数数据密集于左侧,只为极少数右侧的极端值预留空间——这就是偏态的典型特征!
当我们遇到偏态分布时,出于分析便利的考虑,往往希望其服从对数正态分布。但现实中还存在诸多其他类型的偏态分布——尤其是那些具有”厚尾”特征的分布,往往需要特殊方法进行处理。因此,对偏态变量进行对数变换后,务必检验其分布形态是否真正接近正态24。若变换未达预期效果,则可能面临复杂的分析困境!
那么,如何利用经验数据来推断理论分布的特征呢?
请始终牢记:我们观察和描述变量的根本目的,在于更准确地把握理论分布的特征。真正值得关注的并非样本数据本身,而是如何通过样本推断变量的总体行为规律。
实际操作中,我们可以通过样本分布(包括均值、标准差等统计量)来推测理论分布的形态——这无疑是我们所能获得的最佳估计25。但必须清醒认识到,这种推定存在固有局限:难道仅凭手头数据获得的分布,就真能代表潜在的理论分布吗?
相对而言,我们更容易实现的是排除某些理论分布的可能性——虽然无法精确确定理论分布的具体形态,但至少可以否定部分不符合的分布假设。
要判断某个理论分布的可能性程度,可遵循以下步骤:
选定理论分布的描述指标——均值、中位数、标准差等。以均值为例展开说明。
根据理论分布特性及样本容量,确定该统计量的抽样分布——均值通常服从正态分布,且样本量越大,其标准误越小。
计算观测数据的对应统计量——至此我们既掌握理论均值的分布特征,也获得实际观测均值。
利用该统计量的理论分布评估观测结果的极端性——若理论均值分布为\(\mu\)=1、\(\sigma\)=2的正态分布,而观测均值为1.5,则需要计算”从该正态分布中获得≥1.5均值的概率”。
若概率极低,则可初步排除原假设分布。若进行统计显著性检验,可判定观测均值与理论分布均值存在”统计学显著差异”。
让我们通过一个案例逐步解析:假设我们研究篮球运动员的单场得分情况。收集100名球员的数据后,发现观测数据既不符合常规分布形态,也与已知的理论分布不符。但计算得出其均值为102分,标准差为30分。
根据步骤1,我们需要验证”该数据是否可能来自\(\mu\)=90的总体分布”——请注意,此处尚未假定分布形态(正态、对数正态或其他)。当前目标仅是尝试排除均值为90的分布可能。同时留意我们使用希腊字母\(\mu\)表征总体均值(即真实参数),核心问题在于能否排除\(\mu\)=90这一真实情况。
进入步骤2,我们需要确定该统计量的抽样分布。如本章前文所述,样本均值通常服从以理论均值\(\mu\)为中心的正态分布,其标准差计算公式为\(\frac{\sigma}{\sqrt{N}}\)(总体标准差除以样本量N的平方根),这一指标在统计学中称为”标准误”。
这一现象的本质是:若每次调查N名篮球球员并计算均值,通过反复抽样(N名→N名→N名…),每次获得的样本均值都会存在差异。将各次抽样的均值视为随机变量时,其分布将呈现以下特征:(1) 服从正态分布;(2) 中心位置等于总体真实均值\(\mu\);(3) 离散程度由\(\frac{\sigma}{\sqrt{N}}\)决定。其中\(\sqrt{N}\)的倒数关系表明:单次抽样的样本量越大,样本均值逼近真实均值的概率越高。例如:10名球员的样本均值可能严重偏离总体均值,而1,000名球员的样本均值则有极大可能接近真实值——这正是标准误随样本量增大而减小的直观体现。
\(\sigma\)作为希腊字母,代表未知的总体真实标准差。但我们可通过样本标准差(前文计算的30分)对其进行估计。因此,样本均值的抽样分布可表示为:服从均值为\(\mu\)=90、标准误为\(\frac{\sigma}{\sqrt{N}}\)=\(\frac{30}{\sqrt{100}}\)=3的正态分布(注:\(\hat \sigma\)表示对真实\(\sigma\)的估计值,此处即观测数据的标准差)。
执行步骤3时,我们对观测数据进行相同计算(如前所述,样本均值为102分)。
进入步骤4,需要求解”从\(\mu\)=90、\(\sigma\)=3的正态分布中,获得≥102分极端值的概率26“。更准确地说,我们通常关注的是获得偏离\(\mu\)值达到或超过当前程度的概率(即双侧检验)27——102分与90分偏离12分,因此实际计算的是获得≥102分或≤78分(90±12)的联合概率。
通过分析均值为90、标准差为3的正态分布百分位数,可以发现:78分甚至未达到第1百分位(实际约0.004%分位)。这意味着当真实均值\(\mu\)=90、总体标准差\(\sigma\)=30时,在100名球员样本中出现≤78分的均值概率仅为0.004%。同理,102分对应99.996%分位,出现≥102分均值的概率同样仅为0.004%。
在步骤5中,我们将双侧概率相加得出:当真实均值μ=90、真实标准差σ=30时,获得与90分偏离达到或超过102分程度的概率仅为0.008%。这一极端小概率事件表明,现有数据极不可能来源于均值为90的总体分布。
## Chapter Questions
salary_dat <- tibble(Salary=c(85000, 90000, 100000, 120000, 125000, 130000),Frequency=c(5, 4, 1, 2, 3, 2))
dftoLaTeX(salary_dat, title='Economics Professor Salaries', anchor='tab:describingvariables-profsalary', file='professor_salary.tex')
## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{Economics Professor Salaries}\\label{tab:describingvariables-profsalary}\n\\begin{tabular}{ll}\n\\hline\n\\hline\nSalary & Frequency \\\\ \n\\hline\n85000 & 5 \\\\ \n90000 & 4 \\\\ \n1e+05 & 1 \\\\ \n120000 & 2 \\\\ \n125000 & 3 \\\\ \n130000 & 2\\\\ \n\\hline\n\\hline\n\\end{tabular}\n\\end{table}\n"
set.seed(1000)
dat <- tibble(a = rnorm(1000, 10, 3),
b = runif(1000, 0, 20),
c = rnorm(1000, 10, 5),
d = rnorm(1000, 10, 1))
dat %>%
tidyr::pivot_longer(cols=a:d,names_to="var",values_to="val") %>%
ggplot2::ggplot(aes(x=val)) +
ggplot2::geom_density() +
ggplot2::facet_wrap(~ var) +
ggplot2::labs(x="Value", y="Density") +
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"),
axis.text.y = element_blank(),
panel.grid = element_blank())
ggplot2::ggsave('DescribingVariables/question_variability.pdf', width = 5, height = 4, units = 'in', device = cairo_pdf)
exam_scores <- tibble(scores=rsnorm(1000, mean=80, sd=3, xi=1.5))
ggplot2::ggplot(exam_scores, aes(x=scores)) +
ggplot2::geom_density() +
ggplot2::labs(y="Density", x="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"),
axis.text.y = element_blank(),
panel.grid = element_blank(),
axis.ticks.y = element_blank())
ggplot2::ggsave('DescribingVariables/question_students.pdf', width = 4, height = 3, units = 'in', device = cairo_pdf)
dat <- tibble(Classification = c("Freshman", "Sophomore", "Junior", "Senior"),Count = c(1000, 1200, 900, 1500))
dftoLaTeX(
dat,
title = 'Student Standing',
anchor = 'tab:describingvariables-studstanding',
file = 'DescribingVariables/student_standing.tex'
)
## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{Student Standing}\\label{tab:describingvariables-studstanding}\n\\begin{tabular}{ll}\n\\hline\n\\hline\nClassification & Count \\\\ \n\\hline\nFreshman & 1000 \\\\ \nSophomore & 1200 \\\\ \nJunior & 900 \\\\ \nSenior & 1500\\\\ \n\\hline\n\\hline\n\\end{tabular}\n\\end{table}\n"
我们亦可采用假设检验框架重新表述上述过程,具体步骤如下:
步骤1:建立原假设(\(H_0\):\(\mu\)=90)与备择假设(\(H_1\):\(\mu\)≠90)。
步骤2:选择已知分布的检验统计量。由于样本均值服从正态分布,可采用\(Z\)统计量(用于描述正态分布中的点位)。
步骤3:计算观测数据对应的检验统计量值。
步骤4:基于检验统计量的理论分布,计算获得当前统计量值(或更极端值)的概率(即\(p\)值)。
步骤5:确定是否拒绝原假设。这取决于显著性水平\(\alpha\)的设定——若步骤4所得的\(p\)值小于\(\alpha\)阈值(通常设为5%),则拒绝\(\mu\)=90的原假设28。
本节揭示了假设检验的本质与用途。需要明确的是,统计显著性既非结论正确性的保证,亦非效应重要性的度量,其核心价值仅在于能否拒绝预设的理论分布。这一过程本身已极具学术意义——至少,我们得以缩小数据生成机制的潜在范围。而这正是假设检验的根本价值所在!
这又到底是什么意思呢?我们接下来会探讨这个问题。↩︎
ZAR指的是南非兰特,它是南非的官方货币。↩︎
在统计学领域及本书中,“N”统一表示样本容量(即观测值数量)。↩︎
在图形化表示时称为“密度分布”,而在数学描述时称为“概率密度”。↩︎
这之所以可行,仅仅是因为我们有足够多的观测值来填充每个区间。否则,它看起来会有很多起伏。从形式上来说,我们需要让区间变窄,并且观测值的数量增多。↩︎
懂微积分的人会意识到,所有这些操作其实就是对密度曲线求积分。你知道吗,统计学里居然藏着微积分呢!有些事你的教授可不会告诉你,等你发现的时候,都来不及退课啦。不过,在这本书里我们不会涉及微积分相关内容。↩︎
直到你开始将它与其他变量的关系纳入考虑,正如我们将在第四章中讨论的那样。↩︎
好吧,行,实际上你不可能做到 完全 精确地描述分布,除非你还算出第 1.00001 百分位数、第 1.00002 百分位数,以及这两者之间无穷多个百分位数,以此类推。但你懂我的意思。↩︎
那么为什么还要使用均值呢?它也有很多不错的特性!由于一些技术细节过于复杂,这里就不展开说明了,通常在数学上处理均值要容易得多,而且它会出现在各种统计计算中,而不仅仅用于描述一个变量。此外,有时候你确实想要的是代表数值,而不是代表性的观测值。凡事皆有其用武之地。↩︎
是的,真的。孩子们,开始存钱吧。↩︎
有时,如果先对有偏态的数据进行转换,均值也能很好地适用。一种常见的方法是对数据取对数,这样可以控制那些非常大的观测值。我们将在“理论分布”部分更深入地讨论这个问题。↩︎
为何取平方而非其他运算?当然可以使用其他方式!此处计算的方差属于分布的”矩”——均值是一阶矩(线性),方差是二阶矩(平方)反映变异度;三阶矩(立方)体现分布偏斜方向;四阶矩(四次方)度量尾部厚度。经标准化处理后,三阶矩称为偏度,四阶矩称为峰度。↩︎
为何要减1?因均值本身也是估计值,不同样本的均值计算会存在与方差相关的微小差异,导致直接除以n会带来(n-1)/n的偏差。通过除以(n-1)而非n,相当于乘以n/(n-1)的校正因子,即可消除这一偏差。↩︎
另一常见指标是”90/10比率”——第90百分位数与第10百分位数的比值,该指标通过衡量分布顶端与底端的差距,常用于不平等研究领域。↩︎
“自然对数”以常数e为底数。在统计学中,若无特别说明,“log”默认指以e为底的自然对数——这一约定已成学科共识。↩︎
对数变换存在明显局限:无法处理负值或零值(因log(0)无定义)。当前数据集中实际存在大量零收入样本(未图示显示),若需纳入这些样本则必须放弃对数变换。虽然负值处理更为复杂,但对于常见的零值问题,可采用反双曲正弦变换(asinh)——该变换对大数值的拟合效果类似对数变换,且能保持零值不变。↩︎
对于较大增幅,对数增量p实际对应原始变量的增长比例为(e^p - 1)×100%;反之,原始变量增长X%则等价于对数增加log(1+X)。该近似方法成立的原因在于:当p值较小时,e^p - 1与p的数值极为接近。↩︎
由此形成了一套符号体系:英文字母表征”尘世现实”,希腊字母象征”完美真理”。这般精妙设计,究竟源自何人?↩︎
值得指出的是,在一个糟糕的模型中,希腊字母代表的参数未必是”真实的”。但其设计初衷仍是为了指向真相,即便实际上未能达成。↩︎
从技术上讲,它更多代表的是”期望值”而非均值,不过本书不会就此展开深入探讨。↩︎
在统计学中,这被称为”极限分布”——因为在极限情况下(即当观测数据量趋近于无穷大时),得到的就是这个分布。(悄悄说:这又是微积分的概念。微积分正向你逼近呢!)↩︎
但需要明确的是,我们最终得到的理论分布仅反映数据来源的群体特征——这可能并非我们真正关心的目标群体!例如,若儿童分享行为数据仅来自日托中心儿童,那么随着观测数据增加,我们获得的只是日托儿童的理论分布。若研究目标是全体儿童,那么无论采集多少日托儿童数据,都无法推及总体情况。↩︎
这一现象源于所谓的”中心极限定理”。说真的,不妨去读读我之前提到的那些汗牛充栋的统计学导论教材!其中奥妙无穷。↩︎
在偏态分布族中,“幂律分布”尤为棘手——当数据中存在大量极端值时(如股市中频繁出现的暴涨暴跌行情),这类分布就会显现。处理这类分布颇具挑战性,许多理论上的幂律分布甚至不具备明确的均值或方差定义。↩︎
若采用贝叶斯统计方法,我们可将观测分布与先验分布(即数据收集前对理论分布的最佳预估)进行整合分析。↩︎
这里的关键区别在于:3是样本均值抽样分布的标准差(即标准误),而非原始变量自身的标准差30。需要明确的是,\(\frac{\sigma}{\sqrt{N}}\)=30/10=3反映的是样本均值的变异程度,而非个体观测值的离散程度。↩︎
这属于”双侧检验”范畴——若仅评估获得102分及以上值的概率,则称为”单侧检验”。↩︎
将5%作为显著性阈值虽纯属人为设定,却已成为学界普遍标准——犹如QWERTY键盘布局的沿袭,其合理性更多源于历史惯例而非数学必然。↩︎