R语言

Lecture07 : R统计分析

罗益民

2024-09-20

本节内容一览

• Part01 统计学快速回顾
• Part02 假设检验
  • 参数假设检验
  • 非参数假设检验
• Part03 回归分析
• Part04 主成分分析
• Part05 聚类分析

本章依赖的包

运行本章代码之前，请确保已载入以下包：

library(tidyverse)
library(ggplot2)
library(Hmisc)
library(pastecs)
library(psych)
library(MASS)

#我们将使用state.x77数据，确保已将其转为tibble数据框
states <- as_tibble(state.x77)

Part01 统计学快速回顾

统计学讲了什么？

描述性统计量

1. 整体描述性统计

可用内置的summary()函数，以及Hmisc、pastecs、psych等第三方包提供的函数

summary(rnorm(100))

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-2.888068 -0.876392  0.153925 -0.004244  0.782835  2.847946 

summary(mpg)

 manufacturer          model               displ            year     
 Length:234         Length:234         Min.   :1.600   Min.   :1999  
 Class :character   Class :character   1st Qu.:2.400   1st Qu.:1999  
 Mode  :character   Mode  :character   Median :3.300   Median :2004  
                                       Mean   :3.472   Mean   :2004  
                                       3rd Qu.:4.600   3rd Qu.:2008  
                                       Max.   :7.000   Max.   :2008  
      cyl           trans               drv                 cty       
 Min.   :4.000   Length:234         Length:234         Min.   : 9.00  
 1st Qu.:4.000   Class :character   Class :character   1st Qu.:14.00  
 Median :6.000   Mode  :character   Mode  :character   Median :17.00  
 Mean   :5.889                                         Mean   :16.86  
 3rd Qu.:8.000                                         3rd Qu.:19.00  
 Max.   :8.000                                         Max.   :35.00  
      hwy             fl               class          
 Min.   :12.00   Length:234         Length:234        
 1st Qu.:18.00   Class :character   Class :character  
 Median :24.00   Mode  :character   Mode  :character  
 Mean   :23.44                                        
 3rd Qu.:27.00                                        
 Max.   :44.00                                        

sapply(mpg,mean)

manufacturer        model        displ         year          cyl        trans 
          NA           NA     3.471795  2003.500000     5.888889           NA 
         drv          cty          hwy           fl        class 
          NA    16.858974    23.440171           NA           NA 

# 使用Hmisc的describe函数
library(Hmisc)
describe(mpg)

              vars   n    mean    sd median trimmed   mad    min  max range
manufacturer*    1 234    7.76  5.13    6.0    7.68  5.93    1.0   15  14.0
model*           2 234   19.09 11.15   18.5   18.98 14.08    1.0   38  37.0
displ            3 234    3.47  1.29    3.3    3.39  1.33    1.6    7   5.4
year             4 234 2003.50  4.51 2003.5 2003.50  6.67 1999.0 2008   9.0
cyl              5 234    5.89  1.61    6.0    5.86  2.97    4.0    8   4.0
trans*           6 234    5.65  2.88    4.0    5.53  1.48    1.0   10   9.0
drv*             7 234    1.67  0.66    2.0    1.59  1.48    1.0    3   2.0
cty              8 234   16.86  4.26   17.0   16.61  4.45    9.0   35  26.0
hwy              9 234   23.44  5.95   24.0   23.23  7.41   12.0   44  32.0
fl*             10 234    4.63  0.70    5.0    4.77  0.00    1.0    5   4.0
class*          11 234    4.59  1.99    5.0    4.64  2.97    1.0    7   6.0
               skew kurtosis   se
manufacturer*  0.21    -1.63 0.34
model*         0.11    -1.23 0.73
displ          0.44    -0.91 0.08
year           0.00    -2.01 0.29
cyl            0.11    -1.46 0.11
trans*         0.29    -1.65 0.19
drv*           0.48    -0.76 0.04
cty            0.79     1.43 0.28
hwy            0.36     0.14 0.39
fl*           -2.25     5.76 0.05
class*        -0.14    -1.52 0.13

# 使用pastecs的stat.desc()函数
library(pastecs)
stat.desc(mpg)

         manufacturer model      displ         year          cyl trans drv
nbr.val            NA    NA 234.000000 2.340000e+02  234.0000000    NA  NA
nbr.null           NA    NA   0.000000 0.000000e+00    0.0000000    NA  NA
nbr.na             NA    NA   0.000000 0.000000e+00    0.0000000    NA  NA
min                NA    NA   1.600000 1.999000e+03    4.0000000    NA  NA
max                NA    NA   7.000000 2.008000e+03    8.0000000    NA  NA
range              NA    NA   5.400000 9.000000e+00    4.0000000    NA  NA
sum                NA    NA 812.400000 4.688190e+05 1378.0000000    NA  NA
median             NA    NA   3.300000 2.003500e+03    6.0000000    NA  NA
mean               NA    NA   3.471795 2.003500e+03    5.8888889    NA  NA
SE.mean            NA    NA   0.084458 2.948048e-01    0.1053493    NA  NA
CI.mean            NA    NA   0.166399 5.808237e-01    0.2075589    NA  NA
var                NA    NA   1.669158 2.033691e+01    2.5970434    NA  NA
std.dev            NA    NA   1.291959 4.509646e+00    1.6115345    NA  NA
coef.var           NA    NA   0.372130 2.250884e-03    0.2736568    NA  NA
                  cty          hwy fl class
nbr.val   234.0000000  234.0000000 NA    NA
nbr.null    0.0000000    0.0000000 NA    NA
nbr.na      0.0000000    0.0000000 NA    NA
min         9.0000000   12.0000000 NA    NA
max        35.0000000   44.0000000 NA    NA
range      26.0000000   32.0000000 NA    NA
sum      3945.0000000 5485.0000000 NA    NA
median     17.0000000   24.0000000 NA    NA
mean       16.8589744   23.4401709 NA    NA
SE.mean     0.2782199    0.3892672 NA    NA
CI.mean     0.5481481    0.7669333 NA    NA
var        18.1130736   35.4577785 NA    NA
std.dev     4.2559457    5.9546434 NA    NA
coef.var    0.2524439    0.2540358 NA    NA

# 使用psych包的describe()函数
# 注意库名::方法的写法
psych::describe(mpg)

              vars   n    mean    sd median trimmed   mad    min  max range
manufacturer*    1 234    7.76  5.13    6.0    7.68  5.93    1.0   15  14.0
model*           2 234   19.09 11.15   18.5   18.98 14.08    1.0   38  37.0
displ            3 234    3.47  1.29    3.3    3.39  1.33    1.6    7   5.4
year             4 234 2003.50  4.51 2003.5 2003.50  6.67 1999.0 2008   9.0
cyl              5 234    5.89  1.61    6.0    5.86  2.97    4.0    8   4.0
trans*           6 234    5.65  2.88    4.0    5.53  1.48    1.0   10   9.0
drv*             7 234    1.67  0.66    2.0    1.59  1.48    1.0    3   2.0
cty              8 234   16.86  4.26   17.0   16.61  4.45    9.0   35  26.0
hwy              9 234   23.44  5.95   24.0   23.23  7.41   12.0   44  32.0
fl*             10 234    4.63  0.70    5.0    4.77  0.00    1.0    5   4.0
class*          11 234    4.59  1.99    5.0    4.64  2.97    1.0    7   6.0
               skew kurtosis   se
manufacturer*  0.21    -1.63 0.34
model*         0.11    -1.23 0.73
displ          0.44    -0.91 0.08
year           0.00    -2.01 0.29
cyl            0.11    -1.46 0.11
trans*         0.29    -1.65 0.19
drv*           0.48    -0.76 0.04
cty            0.79     1.43 0.28
hwy            0.36     0.14 0.39
fl*           -2.25     5.76 0.05
class*        -0.14    -1.52 0.13

2. 分组描述性统计

1) 使用内置方法

• by(data, INDICES, FUN, ...)是一个通用分组函数，可以对数据集进行分组，并在每个分组上应用任何函数。
• aggregate(x, data, FUN, ...)是将一个或多个变量的值聚合到由另一个变量定义的组中。它通常用于计算每个组的描述性统计量，如平均值、最大值、最小值等。

df <- data.frame(group = rep(c("A", "B"), each = 5), value = 1:10)
# 使用by()计算每个组的平均值
by(df$value, df$group, mean)

df$group: A
[1] 3
------------------------------------------------------------ 
df$group: B
[1] 8

# 使用aggregate()计算每个组的平均值
aggregate(value ~ group, data = df, FUN = mean)

  group value
1     A     3
2     B     8

2)使用dplry

# 注意where和summarise_all的用法
mpg |>  
  #select(where(is.numeric)) |> ## 奇怪奇怪
  group_by(cyl) |> 
  summarise_all(mean)

# A tibble: 4 × 11
    cyl manufacturer model displ  year trans   drv   cty   hwy    fl class
  <int>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1     4           NA    NA  2.15 2003     NA    NA  21.0  28.8    NA    NA
2     5           NA    NA  2.5  2008     NA    NA  20.5  28.8    NA    NA
3     6           NA    NA  3.41 2003.    NA    NA  16.2  22.8    NA    NA
4     8           NA    NA  5.13 2005.    NA    NA  12.6  17.6    NA    NA

频数表和列联表

频数表通常用于显示每个类别在数据集中出现的次数（即频数），又称为一维列联表

评分	频数
差	10
中	20
好	30

列联表又称为交叉表，是一种特殊的频数表，用于展示两个或多个分类变量之间的关系。它通常用于探索变量之间的关联性，通常为二维

	喜欢运动	不喜欢运动	总计
男	20	10	30
女	15	25	40
总计	35	35	70

列联表相关方法：

• table(var1,var2,...,varN)使用N个分类变量创建一个N维列联表
• xtabs(formula, data) 使用一个公式和一个矩阵/数据框创建一个N维列联表

创建列联表

1. 创建一维列联表（频数表）

t_cyl <- table(mtcars$cyl)
t_cyl

 4  6  8 
11  7 14 

# 用prop.table()转化为比例值
prop.table(t_cyl)

      4       6       8 
0.34375 0.21875 0.43750 

2. 创建二维列联表

#使用`table(A,B)`创建二维列联表
with(mpg,table(cyl,drv))

   drv
cyl  4  f  r
  4 23 58  0
  5  0  4  0
  6 32 43  4
  8 48  1 21

#使用`xtabs()`创建列联表
xtabs(~cyl+drv, data=mpg)

   drv
cyl  4  f  r
  4 23 58  0
  5  0  4  0
  6 32 43  4
  8 48  1 21

3. 创建多维列联表

with(mtcars,table(vs,am,gear))

, , gear = 3

   am
vs   0  1
  0 12  0
  1  3  0

, , gear = 4

   am
vs   0  1
  0  0  2
  1  4  6

, , gear = 5

   am
vs   0  1
  0  0  4
  1  0  1

xtabs(~ vs + am + gear,data=mtcars)

, , gear = 3

   am
vs   0  1
  0 12  0
  1  3  0

, , gear = 4

   am
vs   0  1
  0  0  2
  1  4  6

, , gear = 5

   am
vs   0  1
  0  0  4
  1  0  1

Part02 假设检验

假设检验的基本概念

假设检验是统计推断的核心，它使我们能够从样本数据中得出关于总体的结论，并量化这些结论的不确定性。

1. 零假设（H₀）：通常表示“无效应”或“无差异”，是研究者希望通过数据来否定的假设。

2. 备择假设（H₁）：与零假设相对，表示存在某种效应或差异。

Ⅰ类错误：指原假设（H0）本身是正确的，但我们拒绝了它所犯的错误（假阳性） Ⅱ类错误：指原假设（H0）本来是错误的，但我们接受了它所犯的错误（假阴性）

假设检验的步骤

1. 设定假设：明确零假设和备择假设。
2. 选择显著性水平（α）：通常设定为0.05或0.01，表示拒绝零假设的概率阈值。
3. 选择适当的检验统计量：根据数据类型和假设选择合适的统计检验方法，如t检验或z检验。
4. 计算检验统计量：使用样本数据计算检验统计量，并确定其在零假设成立时的分布。
5. 做出决策：根据计算出的检验统计量和显著性水平，决定是否拒绝零假设。如果检验统计量落在拒绝域内，则拒绝零假设；否则，不拒绝零假设。

常见类型的假设检验

1. 参数假设检验：当总体分布已知时，检验总体参数的假设。
2. 非参数假设检验：当对总体分布的了解较少时，使用不依赖于特定分布的检验方法。

参数假设检验

参数假设检验类型

• 两组
  • 独立样本t检验
  • 非独立样本t检验
• 多组
  • 方差分析(ANOVA)

进行参数检验时通常需要满足一些条件：

1. 正态性：尤其是t检验和z检验，通常要求数据服从正态分布。对于大样本，中心极限定理允许一定程度的偏离正态分布，但对于小样本，正态性的要求更为严格。
2. 独立性：观测值之间应该是相互独立的，没有相互影响。这适用于大多数参数检验，如t检验、方差分析（ANOVA）和回归分析。
3. 方差齐性：在方差分析和t检验中，不同组的方差应该相等，这称为方差齐性或方差同质性。

正态性检验

使用Shapiro-Wilk正态性检验来检查数据是否近似正态分布。该检验的零假设是总体呈正态分布。因此，如果p值小于所选的alpha 水平(通常是0.05)，则拒绝零假设，并且有证据表明测试的数据不是正态分布的。

shapiro.test(mpg$hwy)

    Shapiro-Wilk normality test

data:  mpg$hwy
W = 0.95885, p-value = 2.999e-06

也可以通过绘制直方图观察正态性（钟形曲线）

hist(mpg$hwy, freq = FALSE, col = "white")
lines(density(mpg$hwy), lwd = 2, col = "red")

也可以通过绘制QQ图来直观判断数据的正态性。如果数据点在QQ图上近似落在一条直线上，则认为数据是正态分布的

qqnorm(mpg$hwy)
qqline(mpg$hwy)

方差齐性检验

使用Bartlett检验来检查两组数据的方差是否相等，适用于正态分布的数据。零假设是所有组的方差相等，备择假设是至少有一个组的方差与其他组不同。如果检验的p值小于显著性水平（通常是0.05），则拒绝零假设，认为方差不齐。

set.seed(123)
# 创建数据框
df <- data.frame(
  group = rep(c('A', 'B', 'C'), each = 10),
  score = c(85, 86, 88, 75, 78, 94, 98, 79, 71, 80,
           91, 92, 93, 85, 87, 84, 82, 88, 95, 96,
           79, 78, 88, 94, 92, 85, 83, 85, 82, 81)
)

# 执行Bartlett检验
bartlett.test(score ~ group, data = df) 

    Bartlett test of homogeneity of variances

data:  score by group
Bartlett's K-squared = 3.3024, df = 2, p-value = 0.1918

独立样本t检验

任务：将state.x77数据集按照结霜天气的天数（Frost）分为“高、低”两组，比较两组之间的犯罪率（Murder）有无显著性差异。

R语言中进行t检验有两种方式：¹。

1. t.test(x ~ y, data)： \(x\)是一个数值型变量， \(y\) 是个二分变量。适用于数据框
2. t.test(x,y)：\(x\)和\(y\)均为数值型变量。适用于独立数据

state_grouped= as_tibble(state.x77) |> 
  mutate(Frost_Group = ifelse(Frost <= median(Frost, na.rm = TRUE), "Low", "High"))

t.test(Murder ~ Frost_Group, data=state_grouped)

    Welch Two Sample t-test

data:  Murder by Frost_Group
t = -3.9964, df = 47.984, p-value = 0.0002206
alternative hypothesis: true difference in means between group High and group Low is not equal to 0
95 percent confidence interval:
 -5.489359 -1.814641
sample estimates:
mean in group High  mean in group Low 
             5.552              9.204 

# 也可以指定方差齐性
# t.test(Murder ~ Frost_Group, data=state_grouped,var.equal=TRUE) 

结果表明：拒绝原假设\(H0\)，两组之间差异有显著性。

另外一种写法：

low_frost <- as_tibble(state.x77) |> filter(Frost <= median(Frost)) 
high_frost <- as_tibble(state.x77) |> filter(Frost > median(Frost))
t.test(low_frost$Murder,high_frost$Murder)

    Welch Two Sample t-test

data:  low_frost$Murder and high_frost$Murder
t = 3.9964, df = 47.984, p-value = 0.0002206
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.814641 5.489359
sample estimates:
mean of x mean of y 
    9.204     5.552 

1. R语言中的t检验默认两组方差不等，并使用Welsh修正自由度

非独立样本t检验 (配对样本t检验)

非独立样本t检验又被称为配对样本t检验（Paired Samples t-Test）

如何判断样本独立性？

1. 样本来源：如果样本来自于不同的群体或条件下，且这些群体或条件之间没有直接联系，那么它们很可能是独立样本。

2. 测量方式：如果测量涉及到同一组对象在不同时间点或条件下的重复测量，那么这些样本很可能是非独立样本。

3. 数据关联性：如果样本数据之间存在配对关系，如同一个人的不同测量结果，那么这些样本是非独立的。

以下是独立样本和非独立样本的常见例子：

独立样本

1. 硬币投掷：每次投掷结果互不影响，视为独立。
2. 掷骰子：每次掷出点数不受前几次结果影响。
3. 问卷调查：对不同的人独立提问，假设每人回答不互相干扰。
4. 工厂机器生产零件：假设每件产品的质量不受前一件影响。
5. 温度测量：每天在不同的城市测量温度，城市之间温度数据无关联。

非独立样本

1. 不放回抽样：扑克牌中不放回抽取多张牌，后续抽牌会受到已抽出牌的影响。
2. 时间序列数据：例如，股票价格每天会受前一天价格影响。
3. 配对样本：同一人两次测试（如药物实验的前后效果对比），结果之间相关。
4. 家庭成员收入调查：同一家中成员的收入可能有关联。
5. 班级考试分数：班上同学成绩可能因共同的教学资源、学习氛围而相关。

MASS包中的sleep数据集包含了10个实验对象（ID）的睡眠时长(extra)，每个对象有两种处理方式（group）：服用安慰剂（1）和服用安眠药（2）。这种情况适合配对样本t检验。

配对样本t检验使用t.test(x,y,paired=TRUE)， 不能使用方程形式

with(sleep,
     t.test(extra[group == 1],
            extra[group == 2], paired = TRUE))

    Paired t-test

data:  extra[group == 1] and extra[group == 2]
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean difference 
          -1.58 

方差分析（ANOVA）

方差分析（Analysis of Variance，简称ANOVA）用于检验三个或更多组数据的均值是否存在显著差异。它是研究多个因素对一个连续变量影响的基本工具。

方差分析的原假设\(H0\)为组间差异无显著性，备择假设\(H1\)为至少2个组之间差异有显著性。方差分析使用F检验，即计算这两个变异的比率（F比率），以确定组间差异。

方差分析的数据也要满足以下基本假设：

1. 正态性：数据应该呈正态分布。
2. 方差齐性：各组的方差应该相等。
3. 独立性：观测值之间应该是相互独立的。

方差分析的种类

• 单因素组间方差分析(one-way ANOVA)

• 单因素组内方差分析（重复测量方差分析）（one-way repeated measures ANOVA）

• 多因素组间方差分析（MANOVA）

• 多因素组内方差分析（repeated measures MANOVA）

此外还有协方差分析（ANCOVA）: 含有协变量

使用aov()函数进行方差分析

R语言中使用aov(formula, data = NULL...) 函数进行方差分析，其中formula是方程式，data是数据框。

以下是一些典型的方程式¹

• y ~ A : 单因素ANOVA
• y ~ A + B : 单因素ANOVA+协变量（ANCOVA）
• y ~ A * B : 双因素ANOVA

注意以上方程中的 A,B各项自变量的顺序对分析结果有影响

1. 完整方程式解释见课本p205

方差分析举例

table(mtcars$cyl)

 4  6  8 
11  7 14 

mtcars数据集的cyl列包含3种气缸数量类型，如果想探究各种类型气缸数的汽车对油耗有无显著性差异，可以使用单因素方差分析。这里的“因素”就是气缸数。

# 注意：进行方差分析前确保自变量为因子类型！
mtcars$cyl = factor(mtcars$cyl)
fit <- aov(mpg ~ cyl, data = mtcars)
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)    
cyl          2  824.8   412.4    39.7 4.98e-09 ***
Residuals   29  301.3    10.4                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果表明4、6、8缸的各组汽车任意两组间油耗差异均有显著性，拒绝\(H0\)假设

组间成对比较

方差分析仅仅告诉了我们各组之间有差异显著性，但并没有说明各个组之间的差异性。使用TukeyHSD(fit)可以做到这一点

pairwise <- TukeyHSD(fit)
pairwise

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = mpg ~ cyl, data = mtcars)

$cyl
          diff        lwr        upr     p adj
6-4  -6.920779 -10.769350 -3.0722086 0.0003424
8-4 -11.563636 -14.770779 -8.3564942 0.0000000
8-6  -4.642857  -8.327583 -0.9581313 0.0112287

可以看到，如果取\(α=0.01\)的显著性水平，则6-4和8-4组之间有显著性差异

绘制图形

可以利用生成的组间对比结果绘图

library(ggplot2)
plotData <- pairwise$cyl
ggplot(data = plotData,mapping=aes(rownames(plotData),diff)) +
  geom_point(size=8,color="red") + 
  geom_errorbar(aes(ymin=lwr,ymax=upr,width=.2)) +
  geom_line(yintercept=0,color="red",linetype="dashed") +
  theme_bw() +
  coord_flip()

非参数假设检验

在实际的研究中，会经常碰见不满足参数检验的前提条件的情况，比如总体不服从正态分布，或者总体分布未知，这时就不能使用参数检验方法。

非参数检验是在总体分布未知或者知之甚少的情况下，通过样本数据对总体分布形态等特征进行推断的统计检验方法。由于此种检验方法在由样本数据推断总体的过程中不涉及总体分布的参数，不要求总体分布满足某些条件，使用条件较为宽松，因此被称为非参数假设检验

常用的非参数检验方法

• 非参数检验方法
  • 分类变量
    • 卡方检验
    • Fisher精确检验
    • CMH检验
  • 连续变量
    • 两组：Wilconxon符号秩检验（曼-惠特尼U检验，Mann-Whitney U test）
    • 多组：Kruskal-Wallis检验

卡方检验（Chi-square test）

卡方检验用于确定样本数据集中的两个分类变量之间是否存在显著关联。其基本原理是统计样本的实际观测值与理论推断值之间的偏离程度，卡方值越大，二者偏差程度越大；反之，二者偏差越小；若两个值完全相等时，卡方值就为0，表明理论值完全符合。

• 在没有其他的限定条件或说明时，卡方检验一般代指的是皮尔森卡方检验（Pearson’s Chi-squared test）。

• 进行卡方检验前需先构建 \(r × c\) 列联表。

• 在R语言中，卡方检验主要通过chisq.test()函数实现。这个函数可以对一个交叉表进行卡方检验，以确定两个分类变量之间是否独立。注意：chisq.test 函数目前只能处理二维表格。不能接受三维或更高维度的表格作为输入。

我们依然使用mtcars数据集进行演示，但卡方检验只适合于分类变量，我们需要对mpg进行分类变换：

mycars <- mtcars |>
  mutate(mpg_level=cut(mpg, breaks=c(-Inf, 20, 25, Inf), 
                          labels=c("Low", "Medium", "High")))

检查油耗和气缸数有无相关性：

# 构建列联表
table <- xtabs( ~ mpg_level + cyl ,data=mycars)
# 进行卡方检验
chisq.test(table)

    Pearson's Chi-squared test

data:  table
X-squared = 28.641, df = 4, p-value = 9.247e-06

结果表明：拒绝变量间相互独立的原假设\(H0\)，油耗和气缸数具有相关性。

检查油耗和操作方式有无相关性：

table <- xtabs(~ mpg_level + am ,data=mycars)
chisq.test(table)

    Pearson's Chi-squared test

data:  table
X-squared = 13.344, df = 2, p-value = 0.001266

结果表明：拒绝变量间相互独立的原假设\(H0\)，油耗和操作方式具有相关性。

检查引擎形状和操作方式有无相关性：

table <- xtabs(~vs + am ,data=mycars)
chisq.test(table)

    Pearson's Chi-squared test with Yates' continuity correction

data:  table
X-squared = 0.34754, df = 1, p-value = 0.5555

结果表明：无法拒绝变量间相互独立的原假设\(H0\)，引擎形状和操作方式无相关性。

Fisher精确检验（Fisher’s exact test）

Fisher精确检验是一种用于分析离散型数据的统计方法，特别适用于小样本量和低频事件的情况。特点是可以在小样本量下进行精确的推断，不依赖于正态分布假设。

检查油耗和气缸数有无相关性：

table <- xtabs(~mpg_level + cyl ,data=mycars)
fisher.test(table)

    Fisher's Exact Test for Count Data

data:  table
p-value = 3.361e-08
alternative hypothesis: two.sided

Cochran-Mantel-Haenszel检验(CMH test)

CMH检验突出的特点是，当存在第三个分类变量（称为混杂变量）时。这种检验可以用于分层样本分析。

CMH检验可以通过mantelhaen.test()函数来执行。这个函数可以处理2×2×k表格，其中k是分层变量的水平数（混杂变量）

table <- xtabs(~mpg_level + cyl + am ,data=mycars)
mantelhaen.test(table)

    Cochran-Mantel-Haenszel test

data:  table
Cochran-Mantel-Haenszel M^2 = 22.745, df = 4, p-value = 0.0001424

Wilconxon符号秩检验

Wilcoxon符号秩检验（Wilcoxon signed-rank test），又称为Mann-Whitney U检验。适用于连续型数值变量，用于确定两个相关样本或配对样本是否来自具有相同分布的总体。用于成对的数据集，可以检查两个样本的中位数是否存在显著差异，而不需要假设数据服从正态分布。通常用作单样本t检验或配对t检验的非参数替代方法。

wilcox.test(mpg ~  am , data=mtcars)

    Wilcoxon rank sum test with continuity correction

data:  mpg by am
W = 42, p-value = 0.001871
alternative hypothesis: true location shift is not equal to 0

Kruskal-Wallis检验

Kruskal-Wallis检验是一种非参数统计方法，用于确定三个或更多独立样本是否来自具有相同分布的总体。适用于比较两个以上独立样本的中位数是否存在显著差异。Kruskal-Wallis检验不依赖于数据的分布形态，因此它是单因素方差分析（ANOVA）的非参数替代方法。

kruskal.test(mpg ~ cyl, data = mtcars)

    Kruskal-Wallis rank sum test

data:  mpg by cyl
Kruskal-Wallis chi-squared = 25.746, df = 2, p-value = 2.566e-06

检测样本相关性的非参数检验

1. 相关系数

统计中的相关系数是衡量两个变量之间关系强度和方向的统计量。它量化了两个变量之间的线性或单调关系。相关系数的值介于-1（完全负相关）、0（没有相关性）和1（完全正相关）之间。

几种常用的相关系数：

• Pearson相关系数：用于衡量两个连续变量之间的线性关系。
• Spearman等级相关系数：用于衡量两个变量之间的单调关系，适用于非正态分布的数据或顺序变量。
• Kendall等级相关系数：也是一种衡量两个变量之间单调关系的非参数方法。
• Phi系数：衡量两个二元变量之间的相关性。可以看作是Pearson相关系数的特例，适用于2x2列联表。

计算相关系数

在R语言中，cov()和cor()函数是用于计算变量之间的协方差和相关系数的内置函数。

• cov()： 协方差函数。协方差是衡量两个变量联合变化趋势的度量，如果两个变量朝相同方向变化（一个增加，另一个也增加），协方差为正；如果两个变量朝相反方向变化（一个增加，另一个减少），协方差为负。
• cor()： 相关系数函数。相关系数是协方差的标准化形式，其值范围在-1到1之间

cov和cor用法类似：cor(x,use=...,method=...)，x为矩阵或数据框，use指定缺失值的处理方式，method即采用的相关系数的类型，可以是Pearson、Spearman或Kendall。

cov(iris[,1:4])

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    0.6856935  -0.0424340    1.2743154   0.5162707
Sepal.Width    -0.0424340   0.1899794   -0.3296564  -0.1216394
Petal.Length    1.2743154  -0.3296564    3.1162779   1.2956094
Petal.Width     0.5162707  -0.1216394    1.2956094   0.5810063

cor(iris[,1:4])

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

cor(iris[,1:4],method = "spearman")

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1667777    0.8818981   0.8342888
Sepal.Width    -0.1667777   1.0000000   -0.3096351  -0.2890317
Petal.Length    0.8818981  -0.3096351    1.0000000   0.9376668
Petal.Width     0.8342888  -0.2890317    0.9376668   1.0000000

2. 相关性的显著性检验

相关系数表可以清楚地看到变量两两之间的相关性，但需进一步判断相关性是否具有显著性。

R中使用cor.test(x,y,alternative=..., method=...)进行相关显著性判断。x和y为要检验相关性的变量，alternative则用来指定进行双侧检验或单侧检验（取值为”two.side”、“less”或”greater”)，而 method用以指定要计算的相关类型（“pearson”、“kendall”或”spearman”）

with(iris,cor.test(Sepal.Length,Petal.Width))

    Pearson's product-moment correlation

data:  Sepal.Length and Petal.Width
t = 17.296, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7568971 0.8648361
sample estimates:
      cor 
0.8179411 

结果表明，花萼长度和花瓣宽度显著相关（拒绝两者不相关的\(H0\)假设）

Part03 回归分析(Regression Analysis)

统计学中的回归分析是一种强大的预测和分析方法，用于研究一个或多个自变量（解释变量）对一个因变量（响应变量）的影响。回归分析可以帮助我们理解变量之间的关系，预测未来的趋势，以及做出数据驱动的决策。

本章介绍如何使用最小二乘法（OLS）进行回归分析

最小二乘法（OLS）

最小二乘法通过最小化误差的平方和来寻找数据的最佳函数匹配。在简单线性回归中，最小二乘法帮助我们找到一条直线，使得所有数据点到这条直线的垂直距离（即误差）的平方和最小。

用函数lm()拟合回归模型

lm()是R语言中拟合线性模型最基本的函数,其基本格式为：

lm(formula, data, ...)

其中，formula值要拟合的模型公式，data是数据框。模型公式的形式如下所示，~的左侧和右侧分别是应变量和自变量¹

\[ Y \sim X1 + X2 + ... + Xn \]

• 一元(简单)线性回归：\(Y \sim X1\)
• 多项式线性回归：\(Y \sim X + X^2 + X^3\)
• 多元线性回归：\(Y \sim X1 + X2 + X3\)

1. 符号含义见教材p169-p170

一元(简单)线性回归

# 探讨谋杀率和文盲率之间有无相关性
fit <- lm( Murder ~ Illiteracy , data = states)
summary(fit)

Call:
lm(formula = Murder ~ Illiteracy, data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5315 -2.0602 -0.2503  1.6916  6.9745 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.3968     0.8184   2.928   0.0052 ** 
Illiteracy    4.2575     0.6217   6.848 1.26e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.653 on 48 degrees of freedom
Multiple R-squared:  0.4942,    Adjusted R-squared:  0.4836 
F-statistic: 46.89 on 1 and 48 DF,  p-value: 1.258e-08

结果显示：

1. F检验的结果表明：文盲率对谋杀率有显著的正向影响（p<0.01）
2. Multiple R-squared和Adjusted R-squared这两个值，称为“拟合优度”和“修正的拟合优度”，是指回归方程对样本的拟合程度。这里我们可以看到，修正的拟合优度=0.4836 ，也就是大概拟合程度不到五成，表示拟合程度很一般。这个值越高越好，

从结果中可以构建回归方程：

\[ Murder = 2.3968 + 4.2575 × Illiteracy \] 画图表示为：

with(states,
     {
        plot(Illiteracy,Murder)
        abline(fit,col="red")
     })

ggplot有更为简便的内置统计函数：

ggplot(states, aes(Illiteracy,Murder)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red")

广义线性模型（glm）

广义线性模型（Generalized Linear Model，简称GLM）是线性回归模型的推广，它允许因变量（响应变量）的分布属于指数分布族，而不仅仅是正态分布。这种灵活性使得GLM可以应用于更广泛的数据类型和问题，包括分类问题、计数问题和比例问题等。

链接函数（Link function）是GLM中的一个关键概念，它将线性预测 \(η\) 与因变量 \(y\) 的均值 \(μ\) 联系起来。链接函数的选择取决于因变量的分布。例如，对于二项分布，常用的链接函数是logit函数¹。

GLM假设因变量 \(y\) 来自指数分布族，这包括正态分布、二项分布、泊松分布等。不同的分布对应不同的统计问题和链接函数。

1. 基因表达数据符合泊松分布，一般使用poisson(link = “log”)进行拟合

使用glm()拟合Logistic回归

当因变量是二值型或逻辑型变量时，可使用Logistic回归对自变量进行拟合。

我们使用AER包中著名的婚外情数据集Affairs举例

data("Affairs",package="AER")
summary(Affairs)

    affairs          gender         age         yearsmarried    children 
 Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  
 1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  
 Median : 0.000                Median :32.00   Median : 7.000            
 Mean   : 1.456                Mean   :32.49   Mean   : 8.178            
 3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            
 Max.   :12.000                Max.   :57.00   Max.   :15.000            
 religiousness     education       occupation        rating     
 Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  
 1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  
 Median :3.000   Median :16.00   Median :5.000   Median :4.000  
 Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  
 3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  
 Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000  

table(Affairs$affairs)

  0   1   2   3   7  12 
451  34  17  19  42  38 

当数据中不存在分类变量时，可根据需求手动进行创建

Affairs$has_affairs <-  ifelse(Affairs$affairs>0,1,0)
table(Affairs$has_affairs)

  0   1 
451 150 

# 创建逻辑回归模型
fit.full <- glm(has_affairs ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating, 
             data = Affairs, family = binomial)
summary(fit.full)

Call:
glm(formula = has_affairs ~ gender + age + yearsmarried + children + 
    religiousness + education + occupation + rating, family = binomial, 
    data = Affairs)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    1.37726    0.88776   1.551 0.120807    
gendermale     0.28029    0.23909   1.172 0.241083    
age           -0.04426    0.01825  -2.425 0.015301 *  
yearsmarried   0.09477    0.03221   2.942 0.003262 ** 
childrenyes    0.39767    0.29151   1.364 0.172508    
religiousness -0.32472    0.08975  -3.618 0.000297 ***
education      0.02105    0.05051   0.417 0.676851    
occupation     0.03092    0.07178   0.431 0.666630    
rating        -0.46845    0.09091  -5.153 2.56e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 675.38  on 600  degrees of freedom
Residual deviance: 609.51  on 592  degrees of freedom
AIC: 627.51

Number of Fisher Scoring iterations: 4

结果中可以看到，如果取alpha=0.05，则性别、子女、教育、职业对模型贡献不显著。可以去掉这些因素重新拟合

fit.part <- glm(has_affairs ~ age + yearsmarried + religiousness + rating, 
             data = Affairs, family = binomial)
summary(fit.part)

Call:
glm(formula = has_affairs ~ age + yearsmarried + religiousness + 
    rating, family = binomial, data = Affairs)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    1.93083    0.61032   3.164 0.001558 ** 
age           -0.03527    0.01736  -2.032 0.042127 *  
yearsmarried   0.10062    0.02921   3.445 0.000571 ***
religiousness -0.32902    0.08945  -3.678 0.000235 ***
rating        -0.46136    0.08884  -5.193 2.06e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 675.38  on 600  degrees of freedom
Residual deviance: 615.36  on 596  degrees of freedom
AIC: 625.36

Number of Fisher Scoring iterations: 4

可以使用卡方检验对两个模型进行评估

anova(fit.full,fit.part,test="Chisq")

Analysis of Deviance Table

Model 1: has_affairs ~ gender + age + yearsmarried + children + religiousness + 
    education + occupation + rating
Model 2: has_affairs ~ age + yearsmarried + religiousness + rating
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       592     609.51                     
2       596     615.36 -4  -5.8474   0.2108

\(p=0.21\)表明，两者差异不显著，因此选用fit.part即可：始终使用自变量最少的拟合模型。

Part04 主成分分析（PCA）

什么是主成分分析

主成分分析(PCA, principal component analysis)是一种数学降维方法,利用正交变换把一系列可能线性相关的变量转换为一组线性不相关的新变量。也称为主成分PC，用新变量在更小的维度下展示数据的特征。

PCA在生信中的应用举例为：比如转录组表达数据，其中有n个gene的表达量信息。这个n有时会非常大，会到上百万千万。若是直接去对比两个或多个数据，就显得非常困难。经过主成分分析后，就可以将这n个变量降维到r个量(r<n)。

这r个主成分会依据方差的大小进行排序，称作主成分（PC）1、主成分2、……主成分r。而每个主成分的方差在这一组变量中的总方差中所占的比例，即是主成分的贡献度。通常来说，我们仅考察贡献度前2或者前3的主成分，经过可视化后，即得到了二维或三维PCA散点图。

R中的PCA分析

主成份分析可通过求取协方差矩阵的本征向量实现，本征值最大的本征向量即为PC1, Rprincomp函数采样的就是这种方法。我们这里不去探讨内部细节，而是使用一个转录组数据说明PCA的分析方法：

load(url("http://roypub.biovr.cc/demo-rnaseq-data.rds"))
dataFilt_small[1:2,1:2]

                TCGA-3X-AAV9-01A-72R-A41I-07 TCGA-3X-AAVA-01A-11R-A41I-07
ENSG00000000003                         3350                         4252
ENSG00000000005                            1                            1

dim(dataFilt_small)

[1] 1000   44

这是一个标准的转录组数据，每个基因为一行，每个样本为一列。共有1000个gene，44个样本。

接下来我们使用R内置函数prcomp进行PCA分析。

使用PCA对基因降维

prcomp函数把每行当做一个因变量，每列代表一个自变量，每列是一个维度。此时应把基因作为列，所以在进行PCA之前先要将表达矩阵进行转置¹。

pca_gene <- prcomp(t(dataFilt_small))

结果返回一个列表，pca$rotation可查看每个主成分内容，pca$x就是原数据在各个PC上的投影

pca_gene$rotation[1:4,1:4]

                          PC1           PC2           PC3           PC4
ENSG00000000003 -3.327933e-03  1.277140e-02  1.636236e-02  2.214705e-02
ENSG00000000005  3.298104e-07 -7.257778e-06 -7.338082e-07  6.380562e-06
ENSG00000000419 -2.489401e-03  2.845017e-03 -7.098868e-04 -1.784246e-03
ENSG00000000457 -1.163847e-03  3.931079e-03 -2.506624e-03  2.500364e-03

pca_gene$x[1:4,1:4]

                                   PC1       PC2       PC3       PC4
TCGA-3X-AAV9-01A-72R-A41I-07 -52369.54  3380.958 -25368.73 -36744.85
TCGA-3X-AAVA-01A-11R-A41I-07 -10301.82 45656.946 -31615.26  15344.22
TCGA-3X-AAVB-01A-31R-A41I-07 -52227.44 21703.315 -25499.80 -21531.48
TCGA-3X-AAVC-01A-21R-A41I-07   7254.12 38119.208 -45418.20  22031.50

使用PC1和PC2绘制散点图

dt <- data.frame(pca_gene$x)
plot(dt$PC1,dt$PC2)

1. 1. 此处的例子仅为说明原理，实际使用时还要对数据进行对数变换、过滤、标准化、分组后再进行PCA分析。2. 除了R自带的prcomp函数外，也有很多PCA相关的第三方库和包。

使用PCA对样本降维

此时应把样本作为列，对于目前的数据来说，无需进行转置，直接分析即可。

pca_sample <- prcomp(dataFilt_small)
pca_sample$rotation[1:4,1:4]

                                    PC1        PC2          PC3       PC4
TCGA-3X-AAV9-01A-72R-A41I-07 -0.1616436 0.02995391  0.003011641 0.1190609
TCGA-3X-AAVA-01A-11R-A41I-07 -0.1266168 0.03108104 -0.091973080 0.1109210
TCGA-3X-AAVB-01A-31R-A41I-07 -0.1644853 0.03661047 -0.035808769 0.1175495
TCGA-3X-AAVC-01A-21R-A41I-07 -0.1045749 0.01945224 -0.072113740 0.1326800

pca_sample$x[1:4,1:4]

                       PC1       PC2       PC3       PC4
ENSG00000000003 -14571.600 -4182.899 -8809.916 -4433.327
ENSG00000000005  14981.850  2556.527  4028.525 -2063.358
ENSG00000000419   6798.633  1601.558  1401.469 -1046.312
ENSG00000000457   9613.625  2326.563  1444.987 -1201.904

dt <- data.frame(pca_sample$x)

plot(dt$PC1,dt$PC2)

使用PCA对样本降维

Part05 聚类分析

什么是聚类分析

聚类分析是一种无监督学习方法，用于将数据集中的对象分组，使得同一组内的对象之间的相似度高，而不同组之间的对象相似度低。在统计学中，聚类分析被广泛应用于市场研究、社交网络分析、图像分割、基因表达分析等领域。

聚类分析的主要方法

1. 层次聚类：通过构建一个多层次的嵌套聚类树来实现聚类，可以是凝聚的（从单个数据点开始，逐步合并）或分裂的（从一个包含所有数据点的聚类开始，逐步分割）。

2. K-Means聚类：也称为划分聚类。通过迭代地移动聚类中心（质心）来最小化每个点到其聚类中心的距离的平方和。

在R中进行聚类分析

1. 使用hclust完成层次聚类

# 加载数据集
data(iris)
# 选择数据集中的前两个特征进行聚类
iris_data <- iris[, 1:2]
# 计算距离矩阵
dist_matrix <- dist(iris)
# 执行层次聚类
hc_result <- hclust(dist_matrix)
# 绘制树状图
plot(hc_result)

# 使用ggdendro包绘制更美观的图形
library(ggplot2)
library(ggdendro)
ggdendrogram(hc_result) 

2. 使用kmeans完成划分聚类

# 执行K-Means聚类
set.seed(123) 
kmeans_result <- kmeans(iris_data, centers = 3, nstart = 25)
# 查看聚类结果
print(kmeans_result)

K-means clustering with 3 clusters of sizes 50, 53, 47

Cluster means:
  Sepal.Length Sepal.Width
1     5.006000    3.428000
2     5.773585    2.692453
3     6.812766    3.074468

Clustering vector:
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 2 3 2 3 2 3 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2
 [75] 3 3 3 3 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3 3 2 3 3 3 3
[112] 3 3 2 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 2 3 3 3 2 3 3 3 2 3 3 3 2 3
[149] 3 2

Within cluster sum of squares by cluster:
[1] 13.1290 11.3000 12.6217
 (between_SS / total_SS =  71.6 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

# 将聚类结果添加到iris数据框中
iris$cluster <- kmeans_result$cluster
# 可视化聚类结果
plot(iris_data, col = iris$cluster, pch = 19)
text(iris_data, labels = iris$Species[iris$cluster == 1], pos = 4, cex = 0.7)

附录（Appendix）

（book）R计算基因组学

https://compgenomr.github.io/book，开源书籍。内含大量R语言统计学应用知识