内容提要

• 概率
• 参数估计
• 假设检验
• 方差分析
• 回归分析

概率

计算排列组合

• choose(n,k)
• factorial(n)

1

choose(5,2)

1

## [1] 10

1

factorial(5)

1

## [1] 120

生成随机数

1

rnorm(3, mean = c(-10,0,10), sd = 1)

1

## [1] -10.7509913   0.9316995  10.5312236

生成随机样本

1

sample(1:10, 5)

1

## [1] 9 2 4 7 3

1

sample(c("H","T"), 5, replace = T)

1

## [1] "H" "H" "H" "T" "T"

1

sample(c("H","T"), 10, replace = T, prob = c(0.2,0.8))

1

##  [1] "T" "T" "H" "T" "T" "T" "T" "H" "H" "T"

向量随机排列

1

sample(1:10)

1

##  [1]  6  7  2  1 10  9  3  5  4  8

使用summary汇总数据

summary对向量, 矩阵, 因子, 数据框都可以汇总

1

summary(1:100)

1
2

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   25.75   50.50   50.50   75.25  100.00

1

summary(mtcars[,1:3])

1
2
3
4
5
6
7

##       mpg             cyl             disp      
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8  
##  Median :19.20   Median :6.000   Median :196.3  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0

分类变量计算频数:table

1

table(rep(1:3, 30))

1
2
3

## 
##  1  2  3 
## 30 30 30

1

table(sample(c("m","f"), 19, replace = T))

1
2
3

## 
##  f  m 
## 10  9

两个分类变量的列联分析

使用table之后, 再进行summary

1
2
3

gender <- rep(c("m","f"), each = 10)
answer <- rep(c("Yes","No","Yes"), c(8,10,2))
summary(table(gender, answer))

1
2
3
4

## Number of cases in table: 20 
## Number of factors: 2 
## Test for independence of all factors:
## 	Chisq = 7.2, df = 1, p-value = 0.00729

参数估计

• 点估计: mean, var, sd
• 估计量的评选标准: 无偏性, 有效性, 相合性
• 区间估计
• 正态总体均值与方差的区间估计

单个正态总体 $N(\mu, \sigma^2)$

问题提法:

• 单个正态总体 $N(\mu, \sigma^2)$
• 置信水平 $1 - \alpha$
• 样本 $x_1, \cdots, x_n$
• 样本均值 $\bar{x}$, 样本方差 $S^2$
• 按$\sigma^2$是否已知的不同情况, 求总体均值 $\mu$ 的$1 - \alpha$置信区间

方差 $\sigma^2$ 已知

当方差 $\sigma^2$ 已知时, $\mu$ 的$1 - \alpha$置信区间为
$$\bar{x} \pm \frac{\sigma}{\sqrt n}z_{\alpha/2}$$

1
2

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 5)

练习: 如果方差已知为25, 给定上述100个随机数, 求该正态总体均值的90%置信区间

练习答案:

1
2
3

x1 <- mean(x) - 5/sqrt(length(x))*qnorm(0.95)
x2 <- mean(x) + 5/sqrt(length(x))*qnorm(0.95)
cat("(", x1,",", x2,")", sep = "")

1

## (0.3937645,2.038618)

方差 $\sigma^2$ 未知

当方差 $\sigma^2$ 未知时, $\mu$ 的$1 - \alpha$置信区间为
$$\bar{x} \pm \frac{S}{\sqrt n}t_{\alpha/2}(n-1)$$

1
2

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 5)

练习: 如果方差未知, 给定上述100个随机数, 求该正态总体均值的90%置信区间

练习答案

1
2
3
4

df <- length(x) - 1
x3 <- mean(x) - sd(x)/sqrt(length(x))*qt(0.95, df)
x4 <- mean(x) + sd(x)/sqrt(length(x))*qt(0.95, df)
cat("(", x3,",", x4,")", sep = "")

1

## (0.3823384,2.050044)

两个正态总体 $N(\mu_1, \sigma_1^2), N(\mu_2, \sigma_2^2)$

问题提法:

• 两个正态总体 $N(\mu_1, \sigma_1^2), N(\mu_2, \sigma_2^2)$
• 置信水平 $1 - \alpha$
• 样本$x_1, \cdots, x_{n_1}$和$y_1, \cdots, y_{n_2}$分别来自两个总体
• 样本均值 $\bar{x}, \bar{y}$, 样本方差 $S_1^2, S_2^2$
• 按$\sigma^2$的不同情况, 求总体均值 $\mu_1 - \mu_2$ 的$1 - \alpha$置信区间

方差 $\sigma_1^2, \sigma_2^2$ 已知

当方差 $\sigma_1^2, \sigma_2^2$ 已知时, $\mu_1 - \mu_2$ 的$1 - \alpha$置信区间为
$$\bar{x}-\bar{y} \pm z_{\alpha/2}{\sqrt {\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}}$$

1
2
3

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 3)
y <- rnorm(80, mean = 3, sd = 5)

练习: 如果上述两个总体方差已知分别为9和25, 给定上述两组随机数, 求这两个正态总体均值差的90%置信区间

练习答案:

1
2
3
4
5

mean.diff <- mean(x) - mean(y)
new.sd <- sqrt(9/length(x) + 25/length(y))
x5 <- mean.diff - qnorm(0.95)*new.sd
x6 <- mean.diff + qnorm(0.95)*new.sd
cat("(", x5,",", x6,")", sep = "")

1

## (-2.813356,-0.7262706)

方差 $\sigma_1^2, \sigma_2^2$ 未知却相等

当方差 $\sigma_1^2, \sigma_2^2$ 未知但 $\sigma_1^2 = \sigma_2^2 = \sigma^2$ 时, $\mu_1 - \mu_2$ 的$1 - \alpha$置信区间为
$$\bar{x}-\bar{y} \pm t_{\alpha/2}(n_1 + n_2 - 2)S_w{\sqrt {\frac{1}{n_1}+\frac{1}{n_2}}}$$
其中
$$S_w = \frac{(n_1 -1)S_1^2 + (n_2 -1)S_2^2}{n_1 + n_2 - 2}$$

用模拟验证置信区间的含义

1
2
3
4
5
6
7
8
9
10

set.seed(1234)
n <- 100
x1 <- numeric(0)
x2 <- numeric(0)
for(i in 1:n){
  x <- rnorm(n, mean = 2, sd = 5)
  x1[i] <- mean(x) - 5/sqrt(length(x))*qnorm(0.95)
  x2[i] <- mean(x) + 5/sqrt(length(x))*qnorm(0.95)
}
sum(x1 > 2) + sum(x2 < 2)

1

## [1] 10

用模拟图形验证置信区间

1
2
3
4
5
6
7
8
9

plot(0:(n+1),0:(n+1),type = "n", ylim = c(-1,5))
for(i in 1:n){
  segments(i,x1[i],i,x2[i])
  if(x1[i] > 2){points(i,x1[i],pch = 20, col = "red", cex = 2)
    }else points(i,x1[i],pch = 20)
  if(x2[i] < 2){points(i,x2[i],pch = 20, col = "red", cex = 2)
    }else points(i,x2[i],pch = 20)  
}
segments(0,2,n+1,2)

假设检验

• 原假设$H_0$, 备择假设$H_1$
• 单侧假设, 双侧假设
• 检验统计量
• 显著性水平(置信水平)
• 临界值
• 拒绝域

单个正态总体 $N(\mu, \sigma^2)$均值$\mu$的检验

问题提法:

• 单个正态总体, 抽取$n$个样本, 检验总体均值是否为某个常数$\mu_0$
• 总体方差$\sigma^2$已知时, 检验统计量$\frac{\bar{x}-\mu_0}{\sigma/\sqrt{n}}$服从正态分布
• 总体方差$\sigma^2$未知时, 检验统计量$\frac{\bar{x}-\mu_0}{s/\sqrt{n}}$服从$t$分布

实际中, 正态总体方差经常为未知, 所以t检验用的更多

t.test 进行单正态总体 t检验

检验以下样本所属总体的均值是否为2, 显著性水平为$\alpha = 0.1$

1
2
3

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 5)
t.test(x, mu = 2, conf.level = 0.9)

1
2
3
4
5
6
7
8
9
10
11

## 
## 	One Sample t-test
## 
## data:  x
## t = -1.5607, df = 99, p-value = 0.1218
## alternative hypothesis: true mean is not equal to 2
## 90 percent confidence interval:
##  0.3823384 2.0500441
## sample estimates:
## mean of x 
##  1.216191

两个正态总体 $N(\mu_1, \sigma_1^2), N(\mu_2, \sigma_2^2)$均值差的检验

问题提法:

• 两正态总体, 分别抽取$n_1$和$n_2$个样本, 检验均值差是否为某个常数
• 总体方差$\sigma_1^2, \sigma_2^2$已知时, 检验统计量服从正态分布
• 总体方差未知但相等($=\sigma^2$)时, 检验统计量服从$t$分布

实际中, 正态总体方差经常为未知, 所以t检验用的更多

t.test 进行两正态总体均值差 t检验

1
2
3
4

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 3)
y <- rnorm(80, mean = 3, sd = 3)
t.test(x,y, var.equal = T)

1
2
3
4
5
6
7
8
9
10
11

## 
## 	Two Sample t-test
## 
## data:  x and y
## t = -3.7543, df = 178, p-value = 0.0002352
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.5172894 -0.7827147
## sample estimates:
## mean of x mean of y 
##  1.529715  3.179717

成对数据配对 t检验

问题提法:

• 一般出现于药物测试等数据中(用药前, 用药后)
• 假设配对数据做差之后满足正态分布$N(\mu, \sigma^2)$(但方差未知)
• 检验药效如何($\mu$和0的大小关系)
• 总体方差未知, 检验统计量服从$t$分布

t.test 进行 成对数据配对 t检验

1
2
3
4

set.seed(1234)
x <- rnorm(50, mean = 130, sd = 10)
y <- rnorm(50, mean = 120, sd = 8)
t.test(x,y,paired = T)

1
2
3
4
5
6
7
8
9
10
11

## 
## 	Paired t-test
## 
## data:  x and y
## t = 2.4616, df = 49, p-value = 0.0174
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.7994091 7.9070582
## sample estimates:
## mean of the differences 
##                4.353234

t.test 的其他参数

• ?t.test
• 单侧双侧: alternative = c(“two.sided”, “less”, “greater”)
• var.equal = FALSE, 两总体时方差是否按相等处理
• conf.level = 0.95

置信区间和假设检验中忽略了的内容

• 单个总体方差的置信区间和假设检验
• 两个总体方差比的置信区间和假设检验

以上内容可自行编程解决

置信区间和假设检验的关系

置信区间上下限可以从假设检验的结果中获得, 如上例的90%置信区间

1
2
3
4
5
6

set.seed(1234)
x <- rnorm(100, mean = 2, sd = 5)
df <- length(x) - 1
x3 <- mean(x) - sd(x)/sqrt(length(x))*qt(0.95, df)
x4 <- mean(x) + sd(x)/sqrt(length(x))*qt(0.95, df)
cat("(", x3,",", x4,")", sep = "")

1

## (0.3823384,2.050044)

置信区间上下限和假设检验的结果

1

t.test(x, conf.level = 0.9)

1
2
3
4
5
6
7
8
9
10
11

## 
## 	One Sample t-test
## 
## data:  x
## t = 2.4217, df = 99, p-value = 0.01727
## alternative hypothesis: true mean is not equal to 0
## 90 percent confidence interval:
##  0.3823384 2.0500441
## sample estimates:
## mean of x 
##  1.216191

1

t.test(x, conf.level = 0.9)$conf.int[1]

1

## [1] 0.3823384

1

t.test(x, conf.level = 0.9)$conf.int[2]

1

## [1] 2.050044

方差分析

• 单因素方差分析
• 双因素方差分析

单因素方差分析

问题和解法:

• 多个正态总体的均值是否相等的比较
• $H_0:$ 相等, $H_1:$ 不相等
• 平方和分解: 组内平方和+组间平方和
• 检验统计量服从$F$分布
• 当$F$值过大时, 拒绝$H_0$, 认为各组均值不全相等

单因素方差分析的数据形式

1
2
3
4
5
6
7

set.seed(1234)
n1 <- 30; n2 <- 40; n3 <- 35
x1 <- rnorm(n1, mean = 1)
x2 <- rnorm(n2, mean = 2)
x3 <- rnorm(n3, mean = 0)
group <- c(rep(1, n1),rep(2, n2),rep(3, n3))
anova.data <- data.frame(x = c(x1,x2,x3), group)

箱线图查看数据均值是否有差异

1

boxplot(x ~ group, anova.data)

使用oneway.test函数进行方差分析

1

oneway.test(x ~ group, anova.data)

1
2
3
4
5

## 
## 	One-way analysis of means (not assuming equal variances)
## 
## data:  x and group
## F = 26.281, num df = 2.000, denom df = 66.962, p-value = 3.762e-09

双因素和单因素的区别只在于formula(见线性回归部分)

线性回归

• 一元线性回归
• 多元线性回归
• 模型检验(模型, 单个系数, 残差)

常见的线性回归模型

• $y = \alpha + \beta x + \epsilon$
• $y = \beta x + \epsilon$
• $y = \alpha + \beta_1 x_1 + \beta_2 x_2 + \epsilon$
• $y = \alpha + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1x_2 + \epsilon$

线性模型——最小二乘估计

线性模型 $y = 3 + 2x + e$ 的数据准备：

1
2
3
4
5

set.seed(1234)
x <- runif(50, 1, 3)
e <-  rnorm(50)
y <-  3 + 2 * x + e
lm.data <- data.frame(x, y, e)

注意：一定要将所有变量的数据放到数据框中, 并命名

一元线性回归

对上述模型进行线性回归分析：

1

lm(y ~ x, data = lm.data)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = y ~ x, data = lm.data)
## 
## Coefficients:
## (Intercept)            x  
##       2.572        2.104

公式的基本用法

• ?lm, ?formula
• y ~ x 表示 y 是因变量, x 是自变量
• y 对 x 做线性回归分析, 默认是带有截距项的
• y ~ x + 1, 结论和 y ~ x 是一样的, 只是将截距项显式表示了出来
• y ~ x - 1, 是不带有截距项的回归
• 如果有两个自变量 x1 和 x2, 则写为 y ~ x1 + x2
• 交互项的公式见后面, 其他formula用法见帮助文档

用图形展示一元线性回归

1
2
3
4

linear.model <- lm(y ~ x, data = lm.data)
plot(x,y)
abline(linear.model$coefficients)
abline(3,2, col = "red")

估计值和真实值的差别

不带截距项的一元线性回归

1
2
3

y <- 2 * x + e
no.intercept <- data.frame(x, y, e)
lm(y ~ x - 1, no.intercept)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = y ~ x - 1, data = no.intercept)
## 
## Coefficients:
##     x  
## 1.898

多元线性回归

数据准备:

1
2
3
4
5

set.seed(1234)
x <- runif(50, 1, 3)
y <-  rnorm(50, 3, 2)
z <-  10 + 2*x + 5*y + e
multilinear <- data.frame(x, y, z, e)

多元线性回归

1

lm(z ~ x + y, multilinear)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = z ~ x + y, data = multilinear)
## 
## Coefficients:
## (Intercept)            x            y  
##         8.5          2.0          5.5

获取回归统计量

将线性模型赋值给一个变量, 然后对其运行以下函数

• anova 方差分析表
• coefficients(coef) 回归系数
• confint 给出回归系数的置信区间
• deviance 残差平方和
• fitted 给出拟合y值的向量
• residuals(resid) 给出模型残差
• summary 重要统计量

回归统计量提取

1
2

m <- lm(y ~ x, data = lm.data)
coef(m);confint(m);

1
2

## (Intercept)           x 
##    2.571960    2.103918

1
2
3

##                2.5 %   97.5 %
## (Intercept) 1.369816 3.774104
## x           1.505274 2.702562

请自行尝试其他函数, 尤其是 summary

运行有交互项的线性回归

• 带有交互项的回归模型: $y = \beta_0 + \beta_1 u + \beta_2 v + \beta_3 uv + \epsilon$
• lm(formula) 其中formula可以是以下三种等价情况
  • y ~ u*v
  • y ~ u + v + u*v
  • y ~ (u + v)^2

交互项举例

1
2
3
4
5

set.seed(1234)
x <- runif(100, 2, 5)
y <- runif(100, 3, 6)
z <- 1 + 2*x + 3*y + 5*x*y + rnorm(100)
lm(z ~ x*y)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = z ~ x * y)
## 
## Coefficients:
## (Intercept)            x            y          x:y  
##      -1.899        2.950        3.813        4.740

选择最合适的回归变量

如果数据可选择的变量太多, 则需要筛选一部分作为回归变量

1
2
3
4
5
6
7

full.model <- lm(y ~ x1 + x2 + x3 + x4)
reduced.model <- step(full.model, direction = "backward")

min.model <- lm(y ~ 1)
fwd.model <- step(min.model, 
                  direction = "forward", 
                  scope = ( ~ x1 + x2 + x3 + x4))

在回归公式中使用表达式

• 问题1: $y = \beta_0 + \beta_1(u + v) + \epsilon$ 如何建立回归模型?
• 问题2: $y = \beta_0 + \beta_1 u + \beta_2 u^2 + \epsilon$ 如何建立回归模型?
• 对问题1: lm(y ~ u + v) 可以吗
• 对问题2: lm(y ~ u + u^2) 可以吗

正确方法:

• 问题1: lm(y ~ I(u + v))
• 问题2: lm(y ~ u + I(u^2))

举例

1
2
3
4

set.seed(1234)
x <- runif(100, 1, 5)
y <- 3 + x + x^2 + rnorm(100, sd = 0.5)
lm(y ~ x + x^2)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = y ~ x + x^2)
## 
## Coefficients:
## (Intercept)            x  
##      -4.055        6.777

1

lm(y ~ x + I(x^2))

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = y ~ x + I(x^2))
## 
## Coefficients:
## (Intercept)            x       I(x^2)  
##      3.4047       0.7102       1.0492

多项式回归

$y = \beta_0 + \beta_1 x_1 + \beta_2 x_1^2 + \beta_3 x_1^3 + \epsilon$

• 错误写法: lm(y ~ x + x^2 + x^3)
• 正确写法: lm(y ~ poly(x,3,raw = TRUE))
• 也可以: lm(y ~ x + I(x^2) + I(x^3))

数据转换

对下图, 应如何拟合? 是线性吗?

以上直接用线性拟合显然不好

1
2

plot(x,y)
abline(coef(lm(y ~ x)))

对以上数据尝试取对数

1
2

z <- log(y)
plot(x,z)

也可以在公式中直接对数据变换

1
2
3
4

set.seed(1234)
x <- runif(100, -1, 2)
y <- exp(1 + x + rnorm(100, sd = 0.1))
lm(log(y) ~ x)

1
2
3
4
5
6
7

## 
## Call:
## lm(formula = log(y) ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      1.0084       0.9985

绘制残差图

1
2

m <- lm(log(y) ~ x)
plot(m, which = 1)

线性回归模型诊断

1
2

par(mfrow = c(2,2))
plot(m)

用线性模型进行预测

• 建立模型
• 将需要预测的自变量放在一个数据框里
• 使用predict函数进行预测

模型预测举例

数据准备

1
2
3
4
5
6
7

set.seed(1234)
x <- seq(1,10,by = 0.5)
y <- 1 + 2 * x + rnorm(length(x))
data <- data.frame(x,y)
original <- lm(y ~ x, data)
preds <- data.frame(x = c(1.7, 2.8, 5.4))
predict(original, newdata = preds)

1
2

##         1         2         3 
##  4.068607  6.251350 11.410563

预测效果

1
2
3
4

new.y <- predict(original, newdata = preds)
plot(data)
abline(coef(original), lwd = 2)
points(c(1.7, 2.8, 5.4), new.y, pch = 20, col = "red", cex = 3)

正态性检验

检验单个变量是否服从正态分布, 使用 shapiro.test

1
2
3

set.seed(123)
x <- rnorm(100, mean = 5, sd = 3)
shapiro.test(x)

1
2
3
4
5

## 
## 	Shapiro-Wilk normality test
## 
## data:  x
## W = 0.99388, p-value = 0.9349

如果 p-value 太小，不是正态分布

使用qq图查看正态性

1
2

qqnorm(x)
qqline(x)

其他内容

• Logistic回归(glm)
• 聚类分析(dist, hclust, cutree)
• 主成分分析(prcomp)
• 因子分析(factanal)