第 55 章 统计推断

Statistical Inference: A Tidy Approach

55.1 案例1:你会给爱情片还是动作片高分?

这是一个关于电影评分的数据集8

library(tidyverse)
d <- ggplot2movies::movies
d
## # A tibble: 58,788 × 24
##   title     year length budget rating votes    r1    r2
##   <chr>    <int>  <int>  <int>  <dbl> <int> <dbl> <dbl>
## 1 $         1971    121     NA    6.4   348   4.5   4.5
## 2 $1000 a…  1939     71     NA    6      20   0    14.5
## 3 $21 a D…  1941      7     NA    8.2     5   0     0  
## 4 $40,000   1996     70     NA    8.2     6  14.5   0  
## 5 $50,000…  1975     71     NA    3.4    17  24.5   4.5
## 6 $pent     2000     91     NA    4.3    45   4.5   4.5
## # … with 58,782 more rows, and 16 more variables:
## #   r3 <dbl>, r4 <dbl>, r5 <dbl>, r6 <dbl>, r7 <dbl>,
## #   r8 <dbl>, r9 <dbl>, r10 <dbl>, mpaa <chr>,
## #   Action <int>, Animation <int>, Comedy <int>,
## #   Drama <int>, Documentary <int>, Romance <int>,
## #   Short <int>

数据集包含58788 行 和 24 变量

variable description
title 电影名
year 发行年份
budget 预算金额
length 电影时长
rating 平均得分
votes 投票人数
r1-10 各分段投票人占比
mpaa MPAA 分级
action 动作片
animation 动画片
comedy 喜剧片
drama 戏剧
documentary 纪录片
romance 爱情片
short 短片

我们想看下爱情片与动作片(不是爱情动作片)的平均得分是否显著不同。

  • 首先我们简单的整理下数据,主要是剔除既是爱情片又是动作片的电影
movies_genre_sample <- d %>%
  select(title, year, rating, Action, Romance) %>%
  filter(!(Action == 1 & Romance == 1)) %>% 
  mutate(genre = case_when(
    Action == 1 ~ "Action",
    Romance == 1 ~ "Romance",
    TRUE ~ "Neither"
  )) %>%
  filter(genre != "Neither") %>%
  select(-Action, -Romance) %>%
  group_by(genre) %>%
  #slice_sample(n = 34) %>%  # sample size = 34
  slice_head(n = 34) %>%  
  ungroup()

movies_genre_sample
## # A tibble: 68 × 4
##   title                  year rating genre 
##   <chr>                 <int>  <dbl> <chr> 
## 1 $windle                2002    5.3 Action
## 2 'A' gai waak           1983    7.1 Action
## 3 'A' gai waak juk jaap  1987    7.2 Action
## 4 'Crocodile' Dundee II  1988    5   Action
## 5 'Gator Bait            1974    3.5 Action
## 6 'Sheba, Baby'          1975    5.5 Action
## # … with 62 more rows
  • 先看下图形
movies_genre_sample %>%
  ggplot(aes(x = genre, y = rating)) +
  geom_boxplot() +
  geom_jitter()
  • 看下两种题材电影评分的分布
movies_genre_sample %>%
  ggplot(mapping = aes(x = rating)) +
  geom_histogram(binwidth = 1, color = "white") +
  facet_grid(vars(genre))
  • 统计两种题材电影评分的均值
summary_ratings <- movies_genre_sample %>%
  group_by(genre) %>%
  summarize(
    mean = mean(rating),
    std_dev = sd(rating),
    n = n()
  )
summary_ratings
## # A tibble: 2 × 4
##   genre    mean std_dev     n
##   <chr>   <dbl>   <dbl> <int>
## 1 Action   4.78   1.56     34
## 2 Romance  5.91   0.994    34

55.1.1 传统的基于频率方法的t检验

假设:

  • 零假设:

    • \(H_0: \mu_{1} - \mu_{2} = 0\)
  • 备选假设:

    • \(H_A: \mu_{1} - \mu_{2} \neq 0\)

两种可能的结论:

  • 拒绝 \(H_0\)
  • 不能拒绝 \(H_0\)
t_test_eq <- t.test(rating ~ genre,
  data = movies_genre_sample,
  var.equal = TRUE
) %>%
  broom::tidy()
t_test_eq
## # A tibble: 1 × 10
##   estimate estimate1 estimate2 statistic  p.value
##      <dbl>     <dbl>     <dbl>     <dbl>    <dbl>
## 1    -1.12      4.78      5.91     -3.54 0.000736
## # … with 5 more variables: parameter <dbl>,
## #   conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>
t_test_uneq <- t.test(rating ~ genre,
  data = movies_genre_sample,
  var.equal = FALSE
) %>%
  broom::tidy()
t_test_uneq
## # A tibble: 1 × 10
##   estimate estimate1 estimate2 statistic  p.value
##      <dbl>     <dbl>     <dbl>     <dbl>    <dbl>
## 1    -1.12      4.78      5.91     -3.54 0.000810
## # … with 5 more variables: parameter <dbl>,
## #   conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>

55.1.2 infer:基于模拟的检验

所有的假设检验都符合这个框架9:

Hypothesis Testing Framework

图 55.1: Hypothesis Testing Framework

  • 实际观察的差别
library(infer)

obs_diff <- movies_genre_sample %>%
  specify(formula = rating ~ genre) %>%
  calculate(
    stat = "diff in means",
    order = c("Romance", "Action")
  )
obs_diff
## Response: rating (numeric)
## Explanatory: genre (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1  1.12
  • 模拟
null_dist <- movies_genre_sample %>%
  specify(formula = rating ~ genre) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 5000, type = "permute") %>% 
  calculate(
    stat = "diff in means",
    order = c("Romance", "Action")
  )
head(null_dist)
## Response: rating (numeric)
## Explanatory: genre (factor)
## Null Hypothesis: independence
## # A tibble: 6 × 2
##   replicate    stat
##       <int>   <dbl>
## 1         1  0.0529
## 2         2 -0.447 
## 3         3  0.341 
## 4         4 -0.271 
## 5         5 -0.0941
## 6         6 -0.0706
  • 可视化
null_dist %>%
  visualize()
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "both")
# shade_p_value(bins = 100, obs_stat = obs_diff, direction = "both")
  • 计算p值
pvalue <- null_dist %>%
  get_pvalue(obs_stat = obs_diff, direction = "two_sided")

pvalue
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1  0.0004
  • 结论

在构建的虚拟(\(\Delta = 0\))的平行世界里,出现实际观察值(1.1235)的概率为(4^{-4})。 如果以(p < 0.05)为标准,我们看到p_value < 0.05, 那我们有足够的证据证明,H0不成立,即爱情电影和动作电影的评分均值存在显著差异,具体来说,动作电影的平均评分要比爱情电影低些。

55.2 案例2: 航天事业的预算有党派门户之见?

美国国家航空航天局的预算是否存在党派门户之见?

gss <- read_rds("./demo_data/gss.rds")

gss %>%
  select(NASA, party) %>%
  count(NASA, party) %>%
  head(8)
## # A tibble: 8 × 3
##   NASA        party     n
##   <fct>       <fct> <int>
## 1 TOO LITTLE  Dem       8
## 2 TOO LITTLE  Ind      13
## 3 TOO LITTLE  Rep       9
## 4 ABOUT RIGHT Dem      22
## 5 ABOUT RIGHT Ind      37
## 6 ABOUT RIGHT Rep      17
## # … with 2 more rows
gss %>%
  ggplot(aes(x = party, fill = NASA)) +
  geom_bar()

假设:

  • 零假设 \(H_0\):

    • 不同党派对预算的态度的构成比(TOO LITTLE, ABOUT RIGHT, TOO MUCH) 没有区别
  • 备选假设 \(H_a\):

    • 不同党派对预算的态度的构成比(TOO LITTLE, ABOUT RIGHT, TOO MUCH) 存在区别

两种可能的结论:

  • 拒绝 \(H_0\)
  • 不能拒绝 \(H_0\)

55.2.1 传统的方法

chisq.test(gss$party, gss$NASA)
## 
##  Pearson's Chi-squared test
## 
## data:  gss$party and gss$NASA
## X-squared = 1.3, df = 4, p-value = 0.9

或者

gss %>%
  chisq_test(NASA ~ party) %>%
  dplyr::select(p_value) %>%
  dplyr::pull()
## [1] 0.8569

55.2.2 infer:Simulation-based tests

obs_stat <- gss %>%
  specify(NASA ~ party) %>%
  calculate(stat = "Chisq")
obs_stat
## Response: NASA (factor)
## Explanatory: party (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1  1.33
null_dist <- gss %>%
  specify(NASA ~ party) %>%                     # (1)
  hypothesize(null = "independence") %>%        # (2)
    generate(reps = 5000, type = "permute") %>% # (3)
  calculate(stat = "Chisq")                     # (4)
null_dist
## Response: NASA (factor)
## Explanatory: party (factor)
## Null Hypothesis: independence
## # A tibble: 5,000 × 2
##   replicate  stat
##       <int> <dbl>
## 1         1  1.65
## 2         2  3.52
## 3         3  3.43
## 4         4 11.0 
## 5         5  2.44
## 6         6  8.44
## # … with 4,994 more rows
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = obs_stat, method = "both", direction = "right")
null_dist %>%
  get_pvalue(obs_stat = obs_stat, direction = "greater")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.870

看到 p_value > 0.05,不能拒绝 \(H_0\),我们没有足够的证据证明党派之间有显著差异

55.2.3 using ggstatsplot::ggbarstats()

library(ggstatsplot)

gss %>% 
  ggbarstats(
    x = party,
    y = NASA
)

55.3 案例3:原住民中的女学生多?

案例 quine 数据集有 146 行 5 列,包含学生的生源、文化、性别和学习成效,具体说明如下

  • Eth: 民族背景:原住民与否 (是”A”; 否 “N”)
  • Sex: 性别
  • Age: 年龄组 (“F0”, “F1,” “F2” or “F3”)
  • Lrn: 学习者状态(平均水平 “AL”, 学习缓慢 “SL”)
  • Days:一年中缺勤天数
td <- MASS::quine %>%
  as_tibble() %>%
  mutate(
    across(c(Sex, Eth), as_factor)
  )
td
## # A tibble: 146 × 5
##   Eth   Sex   Age   Lrn    Days
##   <fct> <fct> <fct> <fct> <int>
## 1 A     M     F0    SL        2
## 2 A     M     F0    SL       11
## 3 A     M     F0    SL       14
## 4 A     M     F0    AL        5
## 5 A     M     F0    AL        5
## 6 A     M     F0    AL       13
## # … with 140 more rows

从民族背景有两组(A, N)来看,性别为 F 的占比 是否有区别?

td %>% count(Eth, Sex)
## # A tibble: 4 × 3
##   Eth   Sex       n
##   <fct> <fct> <int>
## 1 A     F        38
## 2 A     M        31
## 3 N     F        42
## 4 N     M        35

55.3.1 传统方法

prop.test(table(td$Eth, td$Sex), correct = FALSE)
## 
##  2-sample test for equality of proportions
##  without continuity correction
## 
## data:  table(td$Eth, td$Sex)
## X-squared = 0.0041, df = 1, p-value = 0.9
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1564  0.1670
## sample estimates:
## prop 1 prop 2 
## 0.5507 0.5455

55.3.2 基于模拟的方法

被解释变量 sex 中F的占比,解释变量中两组(A,N)

obs_diff <- td %>%
  specify(Sex ~ Eth, success = "F") %>% 
  calculate(
    stat = "diff in props",
    order = c("A", "N")                  
  )

obs_diff
## Response: Sex (factor)
## Explanatory: Eth (factor)
## # A tibble: 1 × 1
##      stat
##     <dbl>
## 1 0.00527
null_distribution <- td %>%
  specify(Sex ~ Eth, success = "F") %>%
  hypothesize(null = "independence") %>%
  generate(reps = 5000, type = "permute") %>%
  calculate(stat = "diff in props", order = c("A", "N"))
null_distribution %>%
  visualize()
pvalue <- null_distribution %>%
  get_pvalue(obs_stat = obs_diff, direction = "both")

pvalue
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       1
null_distribution %>%
  get_ci(level = 0.95, type = "percentile")
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   -0.160    0.170

55.4 宏包infer

我比较喜欢infer宏包的设计思想,它把统计推断分成了四个步骤

下图可以更好的帮助我们理解infer的工作流程
  • specify() 指定解释变量和被解释变量 (y ~ x)

  • hypothesize() 指定零假设 (比如, independence= yx 彼此独立)

  • generate() 从基于零假设的平行世界中抽样:

    • 指定每次重抽样的类型,通俗点讲就是数据洗牌,重抽样type = "bootstrap" (有放回的),对应的零假设往往是null = “point” ; 重抽样type = "permuting" (无放回的),对应的零假设往往是null = “independence”, 指的是y和x之间彼此独立的,因此抽样后会重新排列,也就说原先 value1-group1 可能变成了value1-group2,(因为我们假定他们是独立的啊,这种操作,也不会影响y和x的关系)
    • 指定多少组 (reps = 1000)
  • calculate() 计算每组(reps)的统计值 (stat = "diff in props")

  • visualize() 可视化,对比零假设的分布与实际观察值.

下面是我自己对重抽样的理解