一、判别分析
1.读取数据
setwd("C:\\Users\\91333\\Documents\\semiester5\\MultivariateStatisticalAnalysis\\DiscriminantAnalysis")
rain_data <- read.table(file = "HW1.csv", sep = ",", header = TRUE)
rain_data$Rain <- ifelse(rain_data$Rain == "yes", 1, 0)
2.线性判别
- 交叉验证：由于数据量为20，较小，故选择lda函数的"CV = T"，用留一法交叉验证。
- prior参数：为了使得模型的错判率尽量低，建立循环，使参数"prior"遍历不同的值。
error_rate1 <- rep(NA,9)
for (i in 1:9){
rain_modle1 <- lda(formula = Rain ~ .,
data = rain_data,
prior = c(i, 10-i) / 10,
CV = T)
error_rate1[i] <- sum(rain_modle1$class != rain_data$Rain)/nrow(rain_data)
}
error_rate1
[1] 0.55 0.65 0.40 0.40 0.30 0.35 0.50 0.45 0.45
which.min(error_rate1)
[1] 5
由运行结果可以看出，错判率随着i的增大呈二次变化，并且使得错判率最小的i为5，即prior = c(1, 1) / 2
3.二次判别
error_rate2 <- rep(NA,49)
for (i in 1:49){
rain_modle2 <- qda(formula = Rain ~ .,
data = rain_data,
prior = c(i, 50-i) / 50,
CV = T)
error_rate2[i] <- sum(rain_modle2$class != rain_data$Rain)/nrow(rain_data)
}
error_rate2
[1] 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.45 0.40 0.25 0.25 0.25 0.25 0.25
[15] 0.30 0.30 0.30 0.25 0.20 0.25 0.25 0.35 0.35 0.35 0.35 0.35 0.35 0.35
[29] 0.35 0.35 0.35 0.30 0.30 0.30 0.30 0.30 0.30 0.25 0.30 0.30 0.40 0.35
[43] 0.35 0.40 0.40 0.40 0.40 0.40 0.45
which(error_rate2 == min(error_rate2))
[1] 19
由于把i的循环值设置为1：9时，向量error_rate2的输出值为0.50 0.25 0.30 0.25 0.35 0.35 0.30 0.30 0.40，显然i取2和4时都出现了最错判率，二次判别的错判率变化要比线性判别的复杂，所以扩大i的循环范围进行探索与观察，发现，当把i的循环值设置为1：49时，出现了唯一的最小错判率0.20，并且当继续扩大i的循环范围时，情况不再变化，综上，二次判别的prior参数应该设置为c(19, 31) / 50。
4.选择模型
为了更好地在线性判别和二次判别中选择，我尝试了X1和X2的方差齐性检验，由于X1和X2的经验分布没有任何一种已知分布的显著特征，所以我选择了对分布没有任何要求的非参数的levene检验。
leveneTest(y = c(rain_data$X1, rain_data$X2), group = as.factor(rep(c(1,2),c(20,20))))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 2.8841 0.09763 .
38
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levene检验原假设为方差为齐性，检验结果时在原假设方差齐性成立的情况下，样本发生的概率为0.09763，虽然没有十分显著，但是结合线性判别和二次判别的最终错判率，我最终决定建立prior = c(19, 31) / 50的二次判别模型。
5.进行预测
rain_modle <- qda(formula = Rain ~ .,
data = rain_data,
prior = c(19, 31) / 50
)
tomorrow_predict <- predict(object = rain_modle, newdata = data.frame(X1 = 8.1, X2 = 2.0))
tomorrow_predict
$class
[1] 1 Levels: 0 1
$posterior
0 1
1 0.007994128 0.9920059
由上面的结果可知，预测的结果为明天下雨。
二、聚类分析
X <- data.frame(DXBZ = c(9.30,4.67,0.96,1.38,1.48,2.60,2.15,2.14,6.53,1.47,1.17,0.88,1.23,0.99,0.98,0.85,1.57,1.14,1.34,0.79,1.24,0.96,0.78,0.81,0.57,1.67,1.10,1.49,1.61,1.85), CZBZ = c(30.55,29.38,24.69,29.24,25.47,32.32,26.31,28.46,31.59,26.43,23.74,19.97,16.87,18.84,25.18,26.55,23.16,22.57,23.04,19.14,22.53,21.65,14.65,13.85,3.85,24.36,16.85,17.76,20.27,20.66), WMBZ = c(8.70,8.92,15.21,11.30,15.39,8.81,10.49,10.87,11.04,17.23,17.46,24.43,15.63,16.22,16.87,16.15,15.79,12.10,10.45,10.61,13.97,16.24,24.27,25.44,44.43,17.62,27.93,27.70,22.06,12.75), row.names = as.character(c("北京", "天津", "河北", "山西", "内蒙古", "辽宁", "吉林", "黑龙江", "上海", "江苏", "浙江", "安徽", "福建", "江西", "山东", "河南", "湖北", "湖南", "广东", "广西", "海南", "四川", "贵州", "云南", "西藏", "陕西", "甘肃", "青海", "宁夏", "新疆")))
1. 标准化并计算初始的距离矩阵
Province<-dist(scale(X))
2. 采用两种方法进行聚类：类平均法和mcquitty方法
hc1<-hclust(Province, "average")
hc2<-hclust(Province, "mcquitty")
3. 画出谱系图并且分别聚成五个类别
plot(hc1,hang=-1)
rect.hclust(hc1,k=5)
plot(hc2,hang=-1)
rect.hclust(hc2,k=5)
4. 输出聚类结果
re1<-cutree(hc1,k=5)
re2<-cutree(hc2,k=5)
sort(re1)
北京 天津 上海 河北 山西 内蒙古 辽宁 吉林 黑龙江 江苏
1 2 2 3 3 3 3 3 3 3
浙江 福建 江西 山东 河南 湖北 湖南 广东 广西 海南
3 3 3 3 3 3 3 3 3 3
四川 陕西 新疆 安徽 贵州 云南 甘肃 青海 宁夏 西藏
3 3 3 4 4 4 4 4 4 5
sort(re2)
sort(re2)
北京 天津 上海 河北 内蒙古 江苏 浙江 安徽 福建 江西
1 2 2 3 3 3 3 3 3 3
山东 河南 湖北 湖南 广东 广西 海南 四川 贵州 云南
3 3 3 3 3 3 3 3 3 3
陕西 甘肃 青海 宁夏 新疆 山西 辽宁 吉林 黑龙江 西藏
3 3 3 3 3 4 4 4 4 5
5. 对结果进行分析
在对变量：大学学历所占比例、初中学历所占比例、文盲所占比例的聚类分析中，我使用了类平均法和mcquitty方法聚类分析将我国30个省市自治区聚为5类，两种方法结果大致相同,结合实际，五类结果和我们的生活经验比较相符，以类平均法为例，在结国中，五类分别为：{北京}、{天津 上海}、{河北 山西 内蒙古 辽宁 吉林 黑龙江 江苏 浙江 福建 江西 山东 河南 湖北 湖南 广东 广西 海南 四川 陕西 新疆}、{安徽 贵州 云南 甘肃 青海 宁夏}、{西藏}。对应实际情况：北京是首府、西藏既是世界第三极又是少数民族地区情况特殊，都各为一类、天津和上海两大直辖市经济较为发达为一类、中西部经济较落后的六个省份为一类，其他大多数中东部地区省份为一类。
三、因子分析
1.构成相关矩阵
x <- c(1.000, 0.846, 0.805, 0.859, 0.473, 0.398, 0.301, 0.382,
0.846, 1.000, 0.881, 0.826, 0.376, 0.326, 0.277, 0.277,
0.805, 0.881, 1.000, 0.801, 0.380, 0.319, 0.237, 0.345,
0.859, 0.826, 0.801, 1.000, 0.436, 0.329, 0.327, 0.365,
0.473, 0.376, 0.380, 0.436, 1.000, 0.762, 0.730, 0.629,
0.398, 0.326, 0.319, 0.329, 0.762, 1.000, 0.583, 0.577,
0.301, 0.277, 0.237, 0.327, 0.730, 0.583, 1.000, 0.539,
0.382, 0.415, 0.345, 0.365, 0.629, 0.577, 0.539, 1.000)
names <- c("身高 x1", "手臂长 x2", "上肢长 x3", "下肢长 x4", "体重 x5", "颈围 x6", "胸围 x7", "胸宽 x8")
r <- matrix(x, nrow = 8, dimnames = list(names, names))
2.因子分析
题目中指出数据为305个观测，故n.obs = 305。
尝试factors = 1,2,3,4后发现factors = 3,4时，都有一列因子的因子载荷均为0，故不合适；factors = 1时Cumulative Var小于factors = 2，故也不合适，只剩下factors = 2。
factor1 <- factanal( covmat = r, factors = 2, n.obs = 305, rotation = "varimax")
3.结果分析
print(factor1, digits = 3, cut = .3)
Call:
factanal(factors = 2, covmat = r, n.obs = 305,
rotation = "varimax")
Uniquenesses:
身高 x1 手臂长 x2 上肢长 x3 下肢长 x4 体重 x5 颈围 x6 胸围 x7 胸宽 x8
0.166 0.110 0.166 0.196 0.099 0.360 0.414 0.538
Loadings:
Factor1 Factor2
身高 x1 0.869
手臂长 x2 0.929
上肢长 x3 0.896
下肢长 x4 0.862
体重 x5 0.917
颈围 x6 0.774
胸围 x7 0.752
胸宽 x8 0.643
Factor1 Factor2
SS loadings 3.333 2.618
Proportion Var 0.417 0.327
Cumulative Var 0.417 0.744
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 74.65 on 13 degrees of
freedom.
The p-value is 1.11e-10
1）假设检验
原假设为两个因子为充分的，p-value值为1.11e-10,故拒绝原假设，认为两个因子不充分。
2）方差贡献率
提取两个因子后，累计贡献率为0.744。
3）解释2个因子
根据上面的结果，身高、手臂长、上肢长、下肢长对factor1的因子载荷较大；体重、颈围、胸围和胸宽对factor2的因子载荷较大。
故可以把factor1称为纵向体长因子，把factor2称为横向体宽因子。
2347
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
