多元统计分析 期末总结

2. 多元数据的数学表达

矩阵及 R表示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 创建向量
c(1, 2, 3, 4)

# 创建矩阵
matrix()

# 矩阵转置
t()

# 矩阵相乘
A1 %*% A2

# 矩阵对角运算
diag()

# 矩阵求逆
solve()

# 矩阵的特征值 与 特征向量
eigen()

# 矩阵奇异值分解
svd()

# 矩阵维数
dim()

# 矩阵和
sum()

# 矩阵均值
mean()

数据框 及 R表示

1
2
3
4
5
6
7
8
9
10
11
12
13
data.frame()

# 按行组成
rbind()

# 按列组成
cbind()

# 按行显示(显示前5行)
head(d2.1)

# 数据框的应用
apply()

读数据

1
2
3
4
5
6
7
8
9
# 剪切板读取
read.table("clipboard", header = T)

# 文本文件读取
read.csv()

# Excel读取
library(openxlsx)
read.xlsx()

简单R分析

定量 变量的分析

1
2
3
4
5
6
7
# 直方图
# x:数值向量
# freq: bool: False: 频数直方图, True: 频率直方图
hist(x, freq = TRUE, ...)

# 散点图
plot(x, y, ...)

定性 变量的分析

1
2
3
4
5
6
7
8
9
10
11
# 二维列联表
table(年龄, 性别)

# 饼图
pie(table(投资结果))

# 条图
barplot(table(年龄), col = 1:7)

# 三维列联表
ftable()

3. 多元数据的直观表示

条图

箱尾图

星相图

脸谱图

调和曲线图

4. 多元相关与线性回归分析

1
2
3
4
5
# 矩阵散点图
pairs(x, ...)

# 多元数据散点图,等价于 pairs(d4.4, gap = 0)
plot(d4.4, gap = 0)
1
2
3
4
5
6
7
source("msaR.r")

# 矩阵相关系数效验
# x 为数值矩阵或数据框
msa.cor.test(d4.4)


5. 广义与一般线性模型

变量 的 取值类型

1. 因变量分类

2. 解释变量分类

模型 选择方式(5)

1. y 为 连续变量

2. y 为 0-1 变量

一般用 Logistic 回归模型

3. y 为 有序 变量

一般用 累积比数模型对数线性模型

解释变量 可以是 等级变量 或 因素

4. y 为 多分类 变量

宜用 对数线性模型多分类 Logistic 回归模型

解释变量x 既可是因素 又可是分类变量 或 等级变量

5. y 为 连续伴有删失 变量

一般用 Cox 比例风险模型

x 可是 因素 或 连续变量

广义线性模型

因变量为 非正态分布的线性模型 称为 广义线性模型,例如:Logistic 回归模型,对数线性模型,Cox 比例风险模型

1
2
3
4
5
6
7
# 广义线性模型函数 glm()

glm(formula, family = gaussian, data, ...)

formula 公式,拟合的模型
family 分布族:正态分布(gaussian)、二项分布(binomial)、泊松分布(poission)、伽马分布(gamma)
data 数据框

Logistic 模型

回归函数应限制在 [0, 1] 区间内的连续函数,而不能再沿用线性回归方程。

Logistic 函数(也称 Logit 变换) \[ y = f(x) = \frac{1}{1 + e^{-x}} = \frac{e^x}{1 + e^x} \]

Logistic 模型的估计

Logistic 模型 常用 极大似然估计法 估计参数,用 Newton - Raphson 迭代求解

1
2
3
4
5
6
7
d5.1 = read.xlsx("mvstats5.xlsx", "d5.1")

# 1. Logistic 回归模型
logit.glm <- glm(y ~ x1 + x2 + x3, family = binomial, data = d5.1)

# 2. Logistic 回归模型结果
summary(logit.glm)

有第一列系数(Coefficients) (Intercept, 截距)(Estimate, 估计)即可得到 Logistic 回归模型: \[ P = \frac { exp(0.597610 -1.496084 x_1 -0.001595 x_2 + 0.315865 x_3) } { 1 + exp(0.597610 -1.496084 x_1 -0.001595 x_2 + 0.315865 x_3) } \]\[ logit(P) = 0.597610 - 1.496084 x_1 -0.001595 x_2 + 0.315865 x_3 \] 在此模型中,由于参数 \(\beta_1, \beta_2\) 没有通过效验,可类似于线性模型,用 step() 做变量筛选。

1
2
3
4
5
# 逐步筛选法 变量选择
logit.step <- step(logit.glm, direction = "both")

# 逐步筛选法 变量选择 结果
summary(logit.step)

可得到新的回归方程: \[ P = \frac { exp(0.6190 -1.3728 x_1) } { 1 + (0.6190 -1.3728 x_1) } \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 对 视力正常 和 视力有问题 的驾驶员分别做预测,预测发送交通事故的 概率

# 预测 视力正常 的驾驶员的 Logistic 回归结果
pre1 <- predict(logit.step, data.frame(x1 = 1))

# 预测 视力正常 的驾驶员发生事故的概率
p1 <- exp(pre1) / (1 + exp(pre1))

# 预测 视力有问题 的驾驶员的 Logistic 回归结果
pre2 <- predict(logit.step, data.frame(x1 = 0))

# 预测 视力有问题 的驾驶员发生事故的 概率
p2 <- exp(pre2) / (1 + exp(pre2))

c(p1, p2)

可见:\(p_1 = 0.32, p_2 = 0.65\) ,说明视力有问题驾驶员发送交通事故的概率 是 视力正常的驾驶员的两倍以上

对数线性模型

1
2
# Poisson 分布族
fm <- glm(formula, family = poisson, data = data.frame)

\[ ln(E(y)) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \varepsilon \]

即: \[ E(y) = exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p) \]

1
2
3
4
5
6
7
# 模型的 效验过程

# 多元 对数线性模型
log.glm <- glm(y ~ factor(x1) + factor(x2), family = poisson(link = log), data = d5.2)

# 多元 对数线性模型 结果
summary(log.glm)

从检验结果看出,\(p_1 = 0.0031 < 0.01, p_2 < 0.01\) ,说明收入和满意程度对产品有重要影响

一般线性模型

完全随机设计模型

1
2
3
# 完全随机设计模型 方差分析
# 方差分析(Analysis of Variance, anova)
anova(lm(Y ~ factor(A), data = d5.3))

P < 0.05,说明各机器 生产的薄板厚度有显著差异

随机区组设计模型

1
2
# 随机单位组设计模型 方差分析
anova(lm(Y ~ factor(A) + factor(B), data = d5.4))

\(P_A > 0.05\) ,说明各种燃料A 对火箭射程 无显著影响

\(P_B > 0.05\) ,说明各种推进器B 对火箭射程也 无显著影响

析因设计模型

1
anova(lm(Y ~ factor(A) * factor(B), data = d5.5))

\(P_A < 0.05\) ,说明 不同方法对回收率 有显著影响

\(P_B < 0.05\),说明不同化合物对回收率 有显著影响

\(P_{AB} < 0.05\),说明 方法和化合物之间 的交互作用 对回收率 有显著影响

正交设计模型

1
anova(lm(Y ~ A + B + A * B + C + D), data = d5.6)

\(P_A > 0.05\),说明反应温度A 对农药收率 无显著影响

...

\(P_{AB} < 0.05\),说明反应温度A和 反应时间B 之间的交互作用 对农药的收率 有显著影响

6. 判别分析

Fisher 判别法

线性判别函数(Linear Discriminatory Function, LDF)为: \[ Y = a_1X_1 + a_2X_2 + ... + a_pX_p = a'X \]

判别系数:\(a_i(i=1,2,...,p)\)

指标:\(X_1,X_2,...,X_p\)

判别函数根据指标 之值区分各样品 应归属与哪一类

1
2
3
4
# Fisher 线性判别

library(MASS)
ld = lda(G ~ x1 + x2);ld
1
2
3
lp = predict(ld)  # 根据线性判别模型预测所属类别
G1 = lp$class # 预测的所属类别结果
data.frame(G, G1) # 显示结果
1
tab1 = table(G, G1);tab1 # 判别矩阵
1
sum(diag(prop.table(tab1)))  # 判对率

非线性判别模型

异方差阵:曲线判别

1
2
3
library(MASS)

qd = qda(G ~ x1 + x2);qd
1
2
3
4
5
6
7
8
# 根据 qd 模型预测所属类别
qp = predict(qd)

# 预测的所属类别结果
G2 = qp$class

# 显示结果
data.frame(G, G1, G2)
1
tab2 = table(G, G2);tab2 # 判别矩阵
1
sum(diag(tab2)) / sum(tab2) # 判对率
1
2
# 判定
predict(qd, data.frame(x1 = 8.1, x2 = 2.0))

若假定协方差矩阵相等,进行 线性判别分析 的回带效果 要好于非线性判别

Bayes 判别法

1
2
# 先验概率相等的 Bayes 判别模型
ld41 = lda(G3 ~ Q + C + P, prior = c(1, 1, 1)/3);lda41
1
2
# 先验概率不相等的 Bayes 判别模型
ld42 = lda(G3 ~ Q + C + P, prior = c(5, 8, 7)/20);lda42

7. 聚类分析

简介

聚类分析类型:

  • Q 型:对样品 聚类
  • R型:对变量 聚类

聚类方法可分为:

  1. 系统聚类法
  2. 快速聚类法

聚类统计量(距离公式)

1
2
# 默认 Euclidean 聚类
D = dist(df_data)
1
2
# 添加上三角距离
dist(df_data, upper = True)
1
dist(df_data, upper = True, diag = True)

系统聚类法

1
2
3
4
5
6
# 默认最长法 complete
hc <- hclust(D)

data.frame(hc$merge, hc$height)

plot(hc)
1
2
3
4
5
6
7
8
hc = hclust(D, "single")

names(hc)

data.frame(hc$merge, hc$height)

# 系统聚类图 Dendrogram
plot(hc)
1
2
3
4
5
6
7
plot(hc)

# 分2类框
rect.hclust(hc, 2)

# 分2类
cutree(hc, 2)

快速聚类法

1
kmeans(x, centers, ...)
1
2
3
4
5
6
7
8
9
# K-means 聚类
km = kmeans(X, 2)

# 分类结果
km$cluster

plot(X, pch = km$cluster)

#plot(X, pch = kc, col = kc)

注意

KMeans 原理

K-Means是典型的聚类算法,K-Means算法中的k表示的是聚类为k个簇,means代表取每一个聚类中数据值的均值作为该簇的中心,或者称为质心,即用每一个的类的质心对该簇进行描述。

K-Means步骤

  1. 创建k个点作为起始质心。
  2. 计算每一个数据点到k个质心的距离。把这个点归到距离最近的哪个质心。
  3. 根据每个质心所聚集的点,重新更新质心的位置。
  4. 重复2,3,直到前后两次质心的位置的变化小于一个阈值。

降维方法

  1. PCA
  2. Fisher 线性判别

降维原理

标签两类->Logistic回归

Q&A

补充

参考