1  放弃直线?

还在学生时代的时候,学习心理统计是一件痛苦的事情。那个时候完全不理解心理统计为什么要这么上,我不想知道抽样分布是什么,统计的假设和前提是什么,为什么要手算一遍那些公式,方差分析对变异的分解对我以及我的很多同学来说都如听天书。这个困惑我也不会归因于所谓国内的统计教育上。因为同样的困惑到我留学之后还依然存在,存在于我做助教的心理统计课堂的学生上,存在在我自己的博士同学上。

这个困惑的高潮到来于我们一群上统计课的研究生组了一个课后学习小组,一个女同学哭着喊:“我不明白我学这些干什么,告诉我怎么点不就行了。”(大意;原话大致是I don’t understand why I have to spend the XXXX time to study this. Just tell me where to click.)

那时我突然想起我在本科时也是这么跟同学抱怨的。我们作为应用统计的学科,为什么还要去学习这些统计检验的原理?统计学家们为什么不能把统计的流程标准化,我们只需要把数据整理好,假设有没有满足,结果有没有显著,要报告什么,生成表格全都交给软件完成。这就是我那时候学习统计时的困惑与愿望。

直到我自己开始做一个大学教师,开始上一点跟统计和方法相关的课程。有一天突然看到学生在分享一个思维导图,或者说是一个决策树,根据自变量、因变量数据类型决定应该使用什么数据分析方法与具体步骤。那一刻,上面那些经历就重新涌进我的脑海。一方面,我确实敬佩每一代学生的努力与新的资料的进步,似乎真的在向我上面描述的那个愿望努力。另一方面,也感叹果然每一代学生学习统计的需求都是类似的。

想要那个图的读者只能说不好意思了,我并没有要来那张图保存。对我来说,这样一张图对我来说没有什么帮助。我仔细思考了很多要怎么否定这种做法的说辞,但我觉得这种否定其实并不必要也太过居高临下。事实上,很多领域的研究,确实可以使用一个相对简洁的统计模型里纯靠图形化界面完成统计,比如一个实验组一个对照组的实验研究。而且,即使你完全不懂统计,你也完全可以把这部分工作“外包”给会的人做,比如你的合作者、你的男朋友或女朋友、你的导师、你的同学,甚至在闲鱼1上找一个出卖自己统计技能的“流水工人”。对于统计原理的理解程度完全掌握在学习者自己的手里。就如同现在在看这本书的读者,我能用来诱惑你的理由就是,如果你并不喜欢知其然而不知其所以然,做统计的时候只能按图索骥的感觉。如果你更想创造一种知道自己在做什么的错觉,这就是我写这本书想要带给你的东西。

比起直接学习决策树那样的操作指南,我想传达的学习方式完全相反——我们应该从另一种视角来看待统计。至少在心理学中常见的大部分统计方法,它们并不是待在终点的一个个独立的方框,一组数据也并不是只有一种统计分析方法。我想要带你看到的,是它们内在存在的共通性,让你即使在学习新的陌生的统计方法时,能够利用这种共通性,假装自己能更快地掌握这种新的方法。而这个过程,可能意味着相比直接打开图形化的统计软件开始学习怎么点击,你要多花一些时间,走一段弯路。

那么接下来,如果你还保持着兴趣,我们的弯路会从比较两组之间的差异开始。

2 理解“直线”

2.1 比较两组差异

2.1.1 独立样本t检验

首先,我们得有一个数据。由于我是一个很懒的人,不想另外去找一个数据,还要单独提供下载方式,我们直接用下面的代码生成一个“假”数据。因为这不是一本工具书,所以不可能展开讲每段代码具体做了什么,感兴趣的读者可以复制代码后问AI这段代码在做什么。

## 生成随机种子
set.seed(233)
## 被试ID
ID <- c(1:100)
## 分类变量
iden <- as.factor(c(rep("武装直升机",50),rep("某品牌塑料袋",50)))
## 连续变量
magic_score <- c(rnorm(n = 50, mean = 80, sd = 8),
                 rnorm(n = 50, mean = 85, sd = 8))
## 生成数据框
df1.1 <- data.frame(ID,iden,magic_score)
## 看看数据框长啥样
head (df1.1)
  ID       iden magic_score
1  1 武装直升机    87.16906
2  2 武装直升机    85.85956
3  3 武装直升机    77.53476
4  4 武装直升机    72.38858
5  5 武装直升机    77.15004
6  6 武装直升机    90.52906
tail (df1.1)
     ID         iden magic_score
95   95 某品牌塑料袋    88.65119
96   96 某品牌塑料袋    82.54543
97   97 某品牌塑料袋    91.21250
98   98 某品牌塑料袋    82.36802
99   99 某品牌塑料袋    97.20007
100 100 某品牌塑料袋    80.15866

这就是一个我们在学习心理统计时最开始会看到的数据形式了。为了方便大家理解,一个两水平分类变量当然使用性别最为直观。但是我又不想使用男女这两个现实里纷争不断的性别,避免无形中再次强化读者们的性别意识。所以使用了“武装直升机”和“XXX塑料袋”这两个网上常见的梗,由于不想在书里为任何品牌宣传,所以又把品牌名字隐去,变成了“某品牌塑料袋”。由此可见我是一个性格上相当谨慎的人了。出于这份谨慎,我还想强调一下使用这两个性别认同只是善意的开个多元性别的玩笑,我本人并不反对lgbtq+运动。由于这两个认同名称实在有点长,后文中我们分别简称他们为“机人”和“袋人”。

好!所以现在我们有了我们的自变量性别认同。2 那么因变量是什么呢,因为我不想用各种成绩或者智商之类的变量,所以把因变量称呼为魔法天赋分数。你可以理解成魔法天赋分数越高的人,越有可能有猫头鹰带着录取信敲你窗户。显而易见的,自我认同为武装直升机的机人有着更强的机械亲和,而应该在魔法天赋上逊于自我认同为某品牌塑料袋的袋人。这就是我们的备择假设,机人与袋人的魔法天赋水平有差异。3

接下来的独立样本t检验相信熟练掌握统计知识的大家都不陌生了。呈现在下面:

t.test(magic_score~iden, df1.1)

    Welch Two Sample t-test

data:  magic_score by iden
t = 2.8824, df = 96.036, p-value = 0.004869
alternative hypothesis: true difference in means between group 某品牌塑料袋 and group 武装直升机 is not equal to 0
95 percent confidence interval:
 1.321362 7.166416
sample estimates:
mean in group 某品牌塑料袋   mean in group 武装直升机 
                  85.23551                   80.99163 

作为一个熟练的心理学p值检验机器,相比大家第一时间就已经把目光对准了0.004869这个数,并在脑中自动化生成了“显著”这样一个金灿灿的结果。

\[r_{点二列}=\frac{\bar{Y_1}−\bar{Y_0}}{S_Y}\sqrt{\frac{n_0n_1}{N(N-1)}}\]

y1_m <- mean(df1.1$magic_score[df1.1$iden=="武装直升机"])
y2_m <- mean(df1.1$magic_score[df1.1$iden=="某品牌塑料袋"])
y_s <- sd(df1.1$magic_score)
N <- 100
r_pbis <- (y1_m-y2_m)/y_s*sqrt((N/2)*(N/2)/(N*(N-1)))
r_pbis
[1] -0.279561

\[t_{点二列}=r_b\sqrt{\frac{N-2}{1-r_b^2}}\]

t_pbis <- r_pbis*sqrt((N-2)/(1-r_pbis^2))
t_pbis
[1] -2.882442
## 乘以2计算双尾
p_pbis <- 2*(1- pt(abs(t_pbis), df=N-2))
p_pbis
[1] 0.004849607
print(paste("点二列相关的双尾检验: r=", round(r_pbis, 4),",t=",round(t_pbis, 4),",p=",round(p_pbis, 4),",df=",N-2))
[1] "点二列相关的双尾检验: r= -0.2796 ,t= -2.8824 ,p= 0.0048 ,df= 98"

当然,你也可以使用as.numeric()命令直接改变iden变量的属性。但因为我被外星人劫持了,如果我不这么写代码的话他们就会把全世界的薯片都变成酸甜的番茄口味。为了全世界人类的福祉,我决定向外星人屈服。

df1.1$iden_n <- df1.1$iden
levels(df1.1$iden_n) <- c(1,0)
df1.1$iden_n <- as.numeric(df1.1$iden_n)
cor.test(df1.1$iden_n,df1.1$magic_score)

    Pearson's product-moment correlation

data:  df1.1$iden_n and df1.1$magic_score
t = -2.8824, df = 98, p-value = 0.00485
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.45120317 -0.08797359
sample estimates:
      cor 
-0.279561 
aov1 <- aov(magic_score~iden, df1.1)
summary(aov1)
            Df Sum Sq Mean Sq F value  Pr(>F)   
iden         1    450   450.3   8.308 0.00485 **
Residuals   98   5311    54.2                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t_pbis^2
[1] 8.30847
lm1 <- lm(magic_score~iden_n, df1.1)
summary(lm1)

Call:
lm(formula = magic_score ~ iden_n, data = df1.1)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.7313  -5.0949  -0.9581   4.9255  20.1798 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   89.479      2.328  38.437  < 2e-16 ***
iden_n        -4.244      1.472  -2.882  0.00485 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.362 on 98 degrees of freedom
Multiple R-squared:  0.07815,   Adjusted R-squared:  0.06875 
F-statistic: 8.308 on 1 and 98 DF,  p-value: 0.00485
lm1 <- lm(scale(magic_score)~scale(iden_n), df1.1)
summary(lm1)

Call:
lm(formula = scale(magic_score) ~ scale(iden_n), data = df1.1)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9311 -0.6679 -0.1256  0.6457  2.6453 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
(Intercept)    7.352e-16  9.650e-02   0.000  1.00000   
scale(iden_n) -2.796e-01  9.699e-02  -2.882  0.00485 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.965 on 98 degrees of freedom
Multiple R-squared:  0.07815,   Adjusted R-squared:  0.06875 
F-statistic: 8.308 on 1 and 98 DF,  p-value: 0.00485
r_pbis
[1] -0.279561

\[ y \]


  1. 一个阿里系下的二手交易平台。↩︎

  2. 准确一点来说,应该叫预测变量,因为性别认同并不能改变;但本书为了方便起见,有时会把预测变量和结果变量叫成自变量和因变量。↩︎

  3. 严格意义上来说,我们的备择假设可以是机人的魔法天赋水平低于袋人。也就是一个单尾的检验。这里我们为了方便后续的呈现,使用双尾检验。事实上,实践中也有些专家并不喜欢单尾检验,认为双尾检验更加“保守”,毕竟数值上p值翻了倍,会要求将单尾检验改为双尾。↩︎