2  激斗?“机人” vs “袋人”!

你们可能要熟悉一下我的这种标题党的风格,毕竟这些题目是我花了很久的时间想的。社会心理学的研究者,最擅长的就应该是起这种看上去很有趣的陷阱题目,等你开始读就会发现上当了。

心理统计的真正起点,其实就在两组数据的比较上,从这里开始,就是真正会在发表出来的文章中大量使用的统计方法了。对于心理学这门学科来说,这一点也不奇怪,毕竟我们是一门自诩为科学的学科。那么最自然的,验证因果关系的研究设计,就是一个实验组和一个控制组,通过操纵自变量的水平,探究其对因变量的影响。下面我们就开始回忆一下对这样数据的分析。

首先,我们得有一个数据。由于我是一个很懒的人,不想另外去找一个数据,还要单独提供下载方式,我们直接用下面的代码生成一个“假”数据。因为这不是一本工具书,所以不可能展开讲每段代码具体做了什么,感兴趣的读者可以复制代码后问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+运动。由于这两个认同名称实在有点长,后文中我们分别简称他们为“机人”和“袋人”。

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

2.1 独立样本t检验

针对上面这样的研究问题,拿去问心理统计课程最终及格的本科同学应该使用什么统计分析,想必他们都会不假思索地回答独立样本t检验(One-Sample T test)。没错,我们可以针对这组数据做独立样本t检验。

//补充抽样分布原理

我们可以直接用下面地代码进行独立样本t检验。如果你还记得心理统计的老师教的流程性知识的话,应该还会记得,在做t检验之前要先保证方差齐性这一假设。我们先使用下面的代码检验方差齐性。

## 方差齐性检验
var.test(df1.1[df1.1$iden=="武装直升机",c("magic_score")],
         df1.1[df1.1$iden=="某品牌塑料袋",c("magic_score")])

    F test to compare two variances

data:  df1.1[df1.1$iden == "武装直升机", c("magic_score")] and df1.1[df1.1$iden == "某品牌塑料袋", c("magic_score")]
F = 0.74977, num df = 49, denom df = 49, p-value = 0.3168
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.4254777 1.3212395
sample estimates:
ratio of variances 
         0.7497719 

可以看到,针对方差齐性的检验显示 F = 0.74977, num df = 49, denom df = 49, p-value = 0.3168。检验结果不显著,说明无法拒绝方差齐性的假设,我们就默认它成立了。3 因此,我们接着进行传统的独立样本t检验。

## 独立样本t检验,假设方差齐性
## 下面的代码的意思大致是,魔法分数(magic_score)作为因变量,
## 自我认同(iden)作为自变量进行t检验。
## df1.1是两个变量所属的数据框。
## var.equal = TRUE 即假设方差齐性满足进行检验。
t.test(magic_score~iden, df1.1, var.equal = TRUE)

    Two Sample t-test

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

可以看到,t检验的结果是 t = 2.8824, df = 98, p-value = 0.00485。作为一个熟练的心理学p值检验机器,想比大家第一时间就已经把目光对准了0.00485这个数,并在脑中自动化生成了“显著”这样一个金灿灿的结果。由于是我们创造的“假数据”,这个结果并不让人意外。如果我的

2.2 点二列相关

\[ r_{点二列}=\frac{\bar{Y_1}−\bar{Y_0}}{S_Y}\sqrt{\frac{n_0n_1}{N(N-1)}} \tag{2.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}} \tag{2.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"

将上面的点二列相关系数公式 Equation 2.1 代入对其显著性进行检验的t检验公式 Equation 2.2 中,即可得到下面的独立t检验公式。

\[ t = \frac{\bar{Y}_1 - \bar{Y}_0}{\sqrt{\frac{(n_1-1)S_1^2 + (n_0-1)S_0^2}{n_0 + n_1 - 2} \left( \frac{1}{n_0} + \frac{1}{n_1} \right)}} \]

这个推导十分简单,感兴趣的读者可以自己找一张草稿纸代入计算,并以此入眠。“可惜此处的空白太小”4,我就不详细推导了。读者只需要注意,N=n0 + n1 即可。

2.3 皮尔逊相关

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)

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

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 

2.4 一元(简单)线性回归

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 \]

思考

可以看出,标准化之后改变了我们原来对于分类变量的赋值。当然,统计学的好的同学都知道,标准化变量并不会改变它参与的统计检验的结果。那么如果我随意给分类变量的两个水平赋值,结果又会如何呢?感兴趣的读者可以自己试着理解并运行下面的代码,变换不同的赋值,并尝试解释结果。

table(scale(df1.1$iden_n))

-0.99498743710662  0.99498743710662 
               50                50 
df1.1$iden_nn <- c(rep(233,50),rep(7788,50))
lm2 <- lm(magic_score~iden_nn, df1.1)
summary(lm2)

是的,在上一章中,我说我们要放弃直线,走一段弯路。这段弯路,就是让我们运用”直线“——以线性关系来统合你之前学过的心理统计。很有意思的转折,不是吗?

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
[1] -2.882442

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

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

  3. 其实这种做法有些争议。如果你还记得假设检验的逻辑的话,我们不能单独从不显著的结果中得出结论。这也是大多数针对统计假设的检验手段存在的问题。有些观点认为,我们应该默认一些统计假设无法满足,而直接使用稳健性估计处理。↩︎

  4. 费马语。↩︎