《Discovering Statistics Using R》
笔记说明
读《Discovering Statistics Using R》第六章 Correlation中的6.7节 Comparing correlations做的笔记。
示例数据
有个心理学家对考试焦虑对考试成绩的影响比较感兴趣。她设计了一个量表评估考试焦虑程度。考试前用量表测量学生的焦虑程度（变量Anxiety），用成绩百分位数反映考试表现（变量Exam）。数据在这里：Exam Anxiety
Revise变量表示修改所花的小时数。
#数据导入
library(rio)
examData <- import("data/Exam Anxiety.dat")
str(examData)
## 'data.frame': 103 obs. of 5 variables:
## $ Code : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Revise : int 4 11 27 53 4 22 16 21 25 18 ...
## $ Exam : int 40 65 80 80 40 70 20 55 50 40 ...
## $ Anxiety: num 86.3 88.7 70.2 61.3 89.5 ...
## $ Gender : chr "Male" "Female" "Male" "Male" ...
Gender变量表示学生性别。为后续处理方便将其变为factor类型：
examData$Gender <- factor(examData$Gender)
str(examData)
## 'data.frame': 103 obs. of 5 variables:
## $ Code : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Revise : int 4 11 27 53 4 22 16 21 25 18 ...
## $ Exam : int 40 65 80 80 40 70 20 55 50 40 ...
## $ Anxiety: num 86.3 88.7 70.2 61.3 89.5 ...
## $ Gender : Factor w/ 2 levels "Female","Male": 2 1 2 2 2 1 1 1 1 1 ...
看一下男女各有多少人：
table(examData$Gender)
## Female Male
## 51 52
两组独立样本Pearson相关系数的差异性检验
假设我们想考察焦虑程度和考试表现之间的相关性关系（用Pearson相关系数表示）在男性女性之间是否有差异。
先在男性、女性人群中分别画出散点图并计算样本Pearson相关系数：
library(ggplot2)
scatter <- ggplot(examData, aes(Anxiety, Exam)) + geom_point() + facet_grid(.~ Gender)
scatter
by(examData[,c("Exam","Anxiety")], examData$Gender, cor)
## examData$Gender: Female
## Exam Anxiety
## Exam 1.0000000 -0.3813845
## Anxiety -0.3813845 1.0000000
## ------------------------------------------------------------------------------
## examData$Gender: Male
## Exam Anxiety
## Exam 1.0000000 -0.5056874
## Anxiety -0.5056874 1.0000000
,
对两独立样本相关系数进行差异性检验，我们利用《Discovering Statistics Using R》笔记7-Pearson相关系数中提到的Fisher-Z变换：
或
(其中为反双曲正切函数)
变换后得到的统计量近似服从均值为
，标准差为
的正态分布。
将两独立样本的Pearson相关系数都进行Z变换得到和
。这两个统计量都服从正态分布，且相互独立，则二者差值统计量
也服从正态分布。差值统计量的均值为
。差值统计量的方差为
、
的方差之和，即
。对差值统计量进行标准化得到：
在无效假设下，
服从标准正态分布。由此可计算P值。
R中没有现成的函数进行两独立样本Pearson相关系数的差异性检验，但我们可以根据上述原理自己写一个实现此功能的函数：
zdifference<-function(r1, r2, n1, n2){
zd <- (atanh(r1) - atanh(r2)) / sqrt(1/(n1 - 3) + 1/(n2 - 3))
p <- 1 - pnorm(abs(zd))
print(paste("Z Difference: ", zd))
print(paste("One-Tailed P-Value: ", p))
}
zdifference(-0.506, -0.381, 52, 51)
## [1] "Z Difference: -0.768709306290097"
## [1] "One-Tailed P-Value: 0.221032949510287"
Z检验下，双侧检验的P值即为单侧检验P值的两倍，本例中即为0.442.此检验结果表示不同性别人群间，考试焦虑程度和考试成绩间的Pearson相关系数的差异无统计学意义。
同组样本的不同Pearson相关系数的差异性检验
具体的检验场景是这样的：同一组样本有三个定量变量x,y,z，比较x,y间的Pearson相关系数和z,y间的Pearson相关系数是否有差异。
示例数据中，假设我们想考察考试焦虑和考试成绩之间的Pearson相关系数与修改时间和考试成绩之间的Pearson相关系数是否有差异。
首先看一下散点图和各自相关系数：
pairs(examData[,c("Exam","Anxiety","Revise")])
cor(examData[,c("Exam","Anxiety","Revise")])
## Exam Anxiety Revise
## Exam 1.0000000 -0.4409934 0.3967207
## Anxiety -0.4409934 1.0000000 -0.7092493
## Revise 0.3967207 -0.7092493 1.0000000
我们可以用t检验来检验同组样本的两不同pearson相关系数是否不同：
对应自由度为N-3
式子看起来比较复杂，实际上只需要以下数字：xyz三个变量两两间的相关系数和样本量N。
R中没有现成的函数进行同组样本不同Pearson相关系数的差异性检验，但我们可以根据上述原理自己写一个实现此功能的函数：(注意比较的是rxy和rzy)
tdifference <- function(rxy, rxz, rzy, n){
df <- n - 3
td <- (rxy - rzy) * sqrt((df * (1 + rxz))/(2*(1 -rxy^2- rxz^2- rzy^2 + (2*rxy*rxz*rzy))))
p <-pt(td, df)
print(paste("t Difference: ", td))
print(paste("One-Tailed P-Value: ", p))
}
tdifference(-0.441, -0.709, 0.397, 103)
## [1] "t Difference: -5.09576822523987"
## [1] "One-Tailed P-Value: 8.21913727738007e-07"
P值很小，考试焦虑和考试成绩之间的Pearson相关系数与修改时间和考试成绩之间的Pearson相关系数的差异有统计学意义。
