《Discovering Statistics Using R》

2020-04-04 本文已影响0人新云旧雨

笔记说明

读《Discovering Statistics Using R》第六章 Correlation中的6.3-6.5.4.2节做的笔记。另外涉及第四章散点图的内容。主要是介绍Pearson相关系数。

示例数据

有个心理学家对考试焦虑对考试成绩的影响比较感兴趣。她设计了一个量表评估考试焦虑程度。考试前用量表测量学生的焦虑程度（变量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" ...

散点图

进行两个定量变量的相关性分析的第一步应当是做两个变量的散点图：

#散点图
library(ggplot2)
scatter <- ggplot(examData, aes(Anxiety, Exam)) + geom_point()

从散点图可以看出一些趋势。低水平的焦虑与好成绩相关联，而高水平的焦虑对应成绩的方差很大。

协方差与Pearson相关系数

协方差(covariance)可以衡量成对出现的两变量间线性关系的强度。
$cov(x,y)=\frac{\Sigma(x_i-\overline{x}) (y_i-\overline{y}) }{N-1}$
计算两个变量的协方差可以评估两个变量间的线性相关关系。
协方差为正表示当x偏离 $\overline{x}$ 时，（平均来说）y以同样的方向偏离 $\overline{y}$ 。
协方差为负表示当x偏离 $\overline{x}$ 时，（平均来说）y以相反的方向偏离 $\overline{y}$ 。
可以用cov()计算协方差：

> cov(examData$Anxiety,examData$Exam)
[1] -196.554

协方差受到量纲影响且不容易直接解读，对协方差进行标准化后得到线性相关系数:
$r=cor(X,Y)=\frac{cov(X,Y) }{sd(X)sd(Y)}$
线性相关系数消除了量纲的影响。取值范围为[-1, 1]，相关系数绝对值越接近1表示线性相关的程度越强，相关系数为0表示无线性相关关系。
由上式定义的r称为皮尔森相关系数（pearson correlation coefficient)。有时也用R来表示。一般来说在回归的语境下使用大写的R代表复相关系数（multiple correlation coefficient）；另外，当我们使用 $r^2$ 时也经常使用大写的R。样本person相关系数用r表示，总体person相关系数用ρ表示。

Pearson相关系数的假设检验与区间估计

对Pearson相关系数的假设检验的无效假设为总体相关系数ρ=0，备择假设为总体相关系数ρ≠0。有两种方法进行假设检验。实际应用的时候通常使用t检验来进行假设检验，利用Fisher-Z变换来进行区间估计：
1、Z检验。
如果抽样分布符合正态分布，可以利用Z值进行假设检验。应用到pearson相关系数时有一个问题：r的抽样分布不符合正态分布。根据Fisher大神，对r进行Fisher-Z变换后得到统计量 $z_r$ 近似服从正态分布：
$z_r=\frac{1}{2} ln(\frac{1+r}{1-r})$ 或 $z_r=tanh^{-1}r$
其中 $tanh^{-1}$ 为反双曲正切函数
变换后的 $z_r$ 近似服从均值为 $\frac{1}{2} ln(\frac{1+ρ}{1-ρ})$ ，标准差为 $\frac{1} {\sqrt{N-3}}$ 的正态分布。
进而可以对 $z_r$ 进行Z检验：
$z=\frac{z_r}{SE_{z_r}}=z_r\sqrt{N-3}$
利用Fisher-Z变换得到的 $z_r$ 我们可以对相关系数的进行区间估计。
$z_r$ 的95%置信区间的上下限可以由 $z_r±1.96SE_{z_r}$ 求得。
$r=tanhz_r$
双曲正切函数是单调的，因此r的95%置信区间的上下限可以由 $z_r$ 的95%置信区间的上下限通过双曲正切函数求出。

2、t检验
实际上Pearson相关系数的假设检验通常不使用上述方法，而是用t检验完成。
检验统计量 $t_r$ 服从自由度为N-2的t分布：
$t_r=\frac{r-0}{S_r}$
其中 $S_r=\sqrt{\frac {1-r^2}{N-2}}$

Pearson相关系数的假设

如果只要Pearson相关系数能够准确衡量两个变量之间的线性相关关系，那么只要求数据满足定距变量即可。（定距变量的数据可以进行分类、排序、加减运算，不要求乘除运算。定距变量值之间的差值有实际意义）
如果要说明Pearson相关系数有统计学意义。则需要满足更多的假设：需要变量满足定距变量，并且服从正态分布（我学过的教材中讲相关系数检验假设两个变量满足二元正态分布。二元正态分布的数据在散点图上显示出椭圆形的轮廓形状）。
如果不符合这两条，则应该考虑使用其他相关系数或者使用bootstrapping.

用R计算相关系数

R中常用的可以计算相关系数的函数有：cor()， cor.test()， rcorr()。
前两个函数是R自带的，rcorr()是Hmisc包中的函数，使用前需先加载包。
这三个函数能够实现的相关功能见下图：

表中的Spearman,Kendall是另外两种相关系数，会在之后介绍。
本次笔记省略rcorr()的介绍。
cor()用法：
cor(x, y, use = "string", method = "correlation type")

x 为数值型变量或一个由数值型变量构成的dataframe
y 为数值型变量（若x为dataframe则y可不指定）
use:指定缺失数据如何处理：
（1）"everything" ：如果分析涉及的变量有缺失值，则函数返回NA。默认如此处理。
（2）"all.obs"：使用所有样本，若有缺失值则会报错。
（3）"complete.obs"：使用所有变量均未缺失的样本进行相关系数计算。
（4）"pairwise.complete.obs"：计算某2个变量相关系数时使用这2个变量均未缺失的样本计算。
method:指定计算哪种相关系数 "pearson" "spearman" "kendall"。默认为pearson。若想计算多种可以用c()如：c("pearson","spearman")
例:使用cor()计算各数值型变量之间的Pearson相关系数。
由于变量Gender不是数值型变量，需舍弃。

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

可以使用R自带的pairs()生成散点图矩阵：

pairs(examData[,c("Exam","Anxiety","Revise")])

cor.test()用法：
cor.test(x, y, alternative = "string", method = "correlation type", conf.level = 0.95)

x为数值型变量
y为数值型变量 x y长度需相同。
alternative：指定备择假设类型：
（1）"two.sided"指定双侧检验，默认为此
（2）"less"指定单侧检验，预计相关系数为负
（3）”greater“指定单侧检验，预计相关系数为正
method:指定计算哪种相关系数 "pearson" "spearman" "kendall"。默认为pearson。若想计算多种可以用c()如：c("pearson","spearman")
conf.level：指定相关系数置信区间长度默认为0.95 只计算pearson相关系数的置信区间且需要至少4对数据。
例：

cor.test(examData$Anxiety,examData$Exam)

##  Pearson's product-moment correlation
## 
## data:  examData$Anxiety and examData$Exam
## t = -4.938, df = 101, p-value = 3.128e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5846244 -0.2705591
## sample estimates:
##        cor 
## -0.4409934

cor.test()对pearson相关系数的假设检验使用的方法为上面介绍的t检验，计算置信区间使用的方法是上面介绍的Fisher-Z变换。

警惕对相关系数进行因果关系解读

解读相关系数时需要注意相关并不意味着因果关系（causality）。原因有二：
1、混杂变量问题。在相关分析中，可能存在其他变量影响结果。
2、因果关系的方向问题。相关分析中无法分辨两个变量哪个是因哪个是果。

决定系数

Pearson相关系数的平方 $R^2$ 称为决定系数(coefficient of determination).决定系数表示一个变量的总变异中与另一个变量共享的部分占比（the amount of variability in one variable that is shared by the other）。
举例来说：在刚才考试焦虑程度和考试表现的分析中，考试表现在不同人之间存在变异。可以用离均差平方和来代表考试表现的总变异程度。 $R^2$ 表示了考试表现得总变异中与考试焦虑共享的比例。在本例中 $R^2$ =0.194.
当解释 $R^2$ 时很多人会解读为因果关系，“考试焦虑能解释考试表现19.4%的变异”（the variance in y accounted for by x或the variation in one variable explained by the other）。虽然决定系数是效应实际重要性的很有用的测量，但注意它并不能用来推导因果关系。
计算决定系数：

cor(examData[,c("Exam","Anxiety","Revise")])^2

##             Exam   Anxiety    Revise
## Exam    1.0000000 0.1944752 0.1573873
## Anxiety 0.1944752 1.0000000 0.5030345
## Revise  0.1573873 0.5030345 1.0000000