数据挖掘实验-Rstudio

数据挖掘实验-Rstudio

日期:22/4/29

Ps.孩子下载了Rstudio，打算有关R语言的实验都在Rstudio里做。

任务一：数据采集、抽取、预处理

例一：chengjidan.csv 实验

1、数据准备

chengjidan.csv文件要自己来搞，并且可以自己建立，那我就直接抄老师指导中的例子了。

注意1

我们在建立chengjidan.csv文件是不可直接将excel文件改后缀为csv，这样做的话文件导入后会是乱码，要进入Excel的另存为中选择csv类型。

2、打开Rstudio，导入文件并查看

首先，给出整个Rstudio的工作界面，左上角是R语言的文本编辑处，左下角是控制中心，当然，也可以直接在这里编写R语言，右上角有着环境中的变量，过程等，右下角有着文件，R语言的包等一些内容，现在还没有完全掌握。

直接在控制台输入导入文件的代码：

mydata <- read.csv("D:\\sjwjsy6\\chengjidan.csv",as.is = TRUE)

运行正确后，我们导入了chengjidan.csv文件，并进行查看。

3、对冗余属性和无关变量的处理

对于无关的属性列我们可以直接对列进行删除，这里我们以第2列和第9列为例进行删除

mydata2 <- mydata[-c(2,9)]

4、对缺失值的处理

首先，我们要知道整个表格中有几个缺失值：

> sum(is.na(mydata2)) # 判断mydata2中共有几个缺失值
[1] 2

喏，我们知道了mydata2中共有两个缺失值，接下来要找到它们在哪里

> mydata2[!complete.cases(mydata2),]  # 展示有缺失值的行的数据
     姓名 语文 数学 英语 历史 地理 生物
8  范彬彬   NA   91   83   81   88   80
12 堇思思   NA   75   80   70   84   71
> sub <- which(is.na(mydata2))   # 识别缺失值所在的行
> sub
[1] 23 27
> sub <- which(is.na(mydata2$语文))   # 识别缺失值所在的行
> sub
[1]  8 12

注意2

这里我 $语文 是在语文所在列找确实值，我原本认为不加这个它的输出也是对应列，但是通过上面的结果，我们知道绝不是这样。

它应该是按照列来判断，同时记录这是第几个数字，如果是空值输出这是第几个数据，这样他在碰到这两个空值时才会是 23 和 27（以上为个人猜测，应该是这样，但毕竟是猜的，还需要验证）

但现在我们知道了这两个空值的位置，要对整个mydata2表进行处理，分成有空值和没空值的两部分。

> mydata2_1 <- mydata2[-sub,];head(mydata2_1)
    姓名 语文 数学 英语 历史 地理 生物
1   冯楠   89 90.0   84   83   82   82
2   姜华   90 64.0   82   76   80   78
3   李平   86 69.0   82   72   80   74
4   郑磊   89  9.1   83   69   84   81
5 陈正真   87 87.0   82   80   84   84
6 王盼盼   89 69.0   84   71   80   62
> mydata2_2 <- mydata2[sub,];head(mydata2_2)
     姓名 语文 数学 英语 历史 地理 生物
8  范彬彬   NA   91   83   81   88   80
12 堇思思   NA   75   80   70   84   71

对于有空值的部分，我们可以按照自己的需求进行填充，这里我们按照均值来填充空值

> avg_chinese <- round(mean(mydata2_1$语文),0) #保留语文部分的整数
> mydata2_2$语文 <- rep(avg_chinese,2)    # 把缺失值填充为均值
> result <- rbind(mydata2_1,mydata2_2)   # 将两个部分合并
> result
     姓名 语文 数学 英语 历史 地理 生物
1    冯楠   89 90.0   84   83   82   82
2    姜华   90 64.0   82   76   80   78
3    李平   86 69.0   82   72   80   74
4    郑磊   89  9.1   83   69   84   81
5  陈正真   87 87.0   82   80   84   84
6  王盼盼   89 69.0   84   71   80   62
7    郑娜   94 96.0   85   90   90   85
9    杨娟   87 96.0   82   81   92   86
10   杨爽   88 87.0   82   71   88   72
11   朱清   91 60.0   80   60   84   75
13 叶文强   89 91.0   82   68   82   83
14   赵蕾   91 48.0   81   65   80   67
15 王自健   89 68.0   80   68   76   69
8  范彬彬   89 91.0   83   81   88   80
12 堇思思   89 75.0   80   70   84   71

到现在我们获得了对空值处理后的数据：result

5、对异常值的处理

我们可能期望我们要处理的数据都是完美正确的，有时可能因为输入时手滑了少摁几个0之类的，对于这部分异常值我们也要进行处理。

这里由于数据量不大，我们可以很明显的看到 郑磊 的数学成绩是9.1分，这明显是不合理的。当然在不能明显看出的时候我们也可以通过绘图来找到异常值。

> par(mfrow = c(1,2))   # 将绘图窗口划分为一行两列
> dotchart(result$数学) # 绘制单变量散点图
> boxplot(result$数学,horizontal = T) # 绘制箱线图

我们可以很明显的看出有一个数据点跑偏了

注意3

在绘图时不要将左边Plots窗口缩小，因为我们绘制的图片是在这里显示的，太小了的话会表示图片过大的错误，如下：

在我们知道了异常值在哪里之后对于异常值的处理就要我们自己来判断了，对于本例，我们有理由怀疑是在信息录入的时候多点了个小数点（非常合理），所以我们只要把值改为 91 就可。

> result[4,3] <- 91
> result
     姓名 语文 数学 英语 历史 地理 生物
1    冯楠   89   90   84   83   82   82
2    姜华   90   64   82   76   80   78
3    李平   86   69   82   72   80   74
4    郑磊   89   91   83   69   84   81
5  陈正真   87   87   82   80   84   84
6  王盼盼   89   69   84   71   80   62
7    郑娜   94   96   85   90   90   85
9    杨娟   87   96   82   81   92   86
10   杨爽   88   87   82   71   88   72
11   朱清   91   60   80   60   84   75
13 叶文强   89   91   82   68   82   83
14   赵蕾   91   48   81   65   80   67
15 王自健   89   68   80   68   76   69
8  范彬彬   89   91   83   81   88   80
12 堇思思   89   75   80   70   84   71

例二：bank-additional-full.csv 实验

这个例子中的文件不需要我们自己建立，通过链接下载指定文件（bank-additional.zip），解压，其中的bank-additional-full.csv文件便是本例的数据。

http://archive.ics.uci.edu/ml/machine-learning-databases/00222/

1、导入csv文件

> mydata <- read.csv("E:\\bank-additional\\bank-additional\\bank-additional-full.csv",sep = ';')
> View(mydata)

2、了解数据

利用head()、str()、和summary()三个函数初步探索数据mydata的大体结构、每一列的数据类型及其数据类型

对数据处理如下，head() 是给我们展示前6个数据样本，str（）给我们展示了每个特征的数据类型并给出例子，summary（）对每个特征进行分析，以年龄为例，表明了年龄的范围在17到98之间，一等分点，中点，三等分点分别是32，38，47；年龄均值为40.02岁

> head(mydata)
  age       job marital   education default housing loan
1  56 housemaid married    basic.4y      no      no   no
2  57  services married high.school unknown      no   no
3  37  services married high.school      no     yes   no
4  40    admin. married    basic.6y      no      no   no
5  56  services married high.school      no      no  yes
6  45  services married    basic.9y unknown      no   no
    contact month day_of_week duration campaign pdays previous
1 telephone   may         mon      261        1   999        0
2 telephone   may         mon      149        1   999        0
3 telephone   may         mon      226        1   999        0
4 telephone   may         mon      151        1   999        0
5 telephone   may         mon      307        1   999        0
6 telephone   may         mon      198        1   999        0
     poutcome emp.var.rate cons.price.idx cons.conf.idx
1 nonexistent          1.1         93.994         -36.4
2 nonexistent          1.1         93.994         -36.4
3 nonexistent          1.1         93.994         -36.4
4 nonexistent          1.1         93.994         -36.4
5 nonexistent          1.1         93.994         -36.4
6 nonexistent          1.1         93.994         -36.4
  euribor3m nr.employed  y
1     4.857        5191 no
2     4.857        5191 no
3     4.857        5191 no
4     4.857        5191 no
5     4.857        5191 no
6     4.857        5191 no
> str(mydata)
'data.frame':	41188 obs. of  21 variables:
 $ age           : int  56 57 37 40 56 45 59 41 24 25 ...
 $ job           : chr  "housemaid" "services" "services" "admin." ...
 $ marital       : chr  "married" "married" "married" "married" ...
 $ education     : chr  "basic.4y" "high.school" "high.school" "basic.6y" ...
 $ default       : chr  "no" "unknown" "no" "no" ...
 $ housing       : chr  "no" "no" "yes" "no" ...
 $ loan          : chr  "no" "no" "no" "no" ...
 $ contact       : chr  "telephone" "telephone" "telephone" "telephone" ...
 $ month         : chr  "may" "may" "may" "may" ...
 $ day_of_week   : chr  "mon" "mon" "mon" "mon" ...
 $ duration      : int  261 149 226 151 307 198 139 217 380 50 ...
 $ campaign      : int  1 1 1 1 1 1 1 1 1 1 ...
 $ pdays         : int  999 999 999 999 999 999 999 999 999 999 ...
 $ previous      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ poutcome      : chr  "nonexistent" "nonexistent" "nonexistent" "nonexistent" ...
 $ emp.var.rate  : num  1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 ...
 $ cons.price.idx: num  94 94 94 94 94 ...
 $ cons.conf.idx : num  -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 -36.4 ...
 $ euribor3m     : num  4.86 4.86 4.86 4.86 4.86 ...
 $ nr.employed   : num  5191 5191 5191 5191 5191 ...
 $ y             : chr  "no" "no" "no" "no" ...
> summary(mydata)
      age            job              marital         
 Min.   :17.00   Length:41188       Length:41188      
 1st Qu.:32.00   Class :character   Class :character  
 Median :38.00   Mode  :character   Mode  :character  
 Mean   :40.02                                        
 3rd Qu.:47.00                                        
 Max.   :98.00                                        
  education           default            housing         
 Length:41188       Length:41188       Length:41188      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                           
                                                         
     loan             contact             month          
 Length:41188       Length:41188       Length:41188      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                      
                                                         
 day_of_week           duration         campaign     
 Length:41188       Min.   :   0.0   Min.   : 1.000  
 Class :character   1st Qu.: 102.0   1st Qu.: 1.000  
 Mode  :character   Median : 180.0   Median : 2.000  
                    Mean   : 258.3   Mean   : 2.568  
                    3rd Qu.: 319.0   3rd Qu.: 3.000  
                    Max.   :4918.0   Max.   :56.000  
     pdays          previous       poutcome        
 Min.   :  0.0   Min.   :0.000   Length:41188      
 1st Qu.:999.0   1st Qu.:0.000   Class :character  
 Median :999.0   Median :0.000   Mode  :character  
 Mean   :962.5   Mean   :0.173                     
 3rd Qu.:999.0   3rd Qu.:0.000                     
 Max.   :999.0   Max.   :7.000                     
  emp.var.rate      cons.price.idx  cons.conf.idx  
 Min.   :-3.40000   Min.   :92.20   Min.   :-50.8  
 1st Qu.:-1.80000   1st Qu.:93.08   1st Qu.:-42.7  
 Median : 1.10000   Median :93.75   Median :-41.8  
 Mean   : 0.08189   Mean   :93.58   Mean   :-40.5  
 3rd Qu.: 1.40000   3rd Qu.:93.99   3rd Qu.:-36.4  
 Max.   : 1.40000   Max.   :94.77   Max.   :-26.9  
   euribor3m      nr.employed        y            
 Min.   :0.634   Min.   :4964   Length:41188      
 1st Qu.:1.344   1st Qu.:5099   Class :character  
 Median :4.857   Median :5191   Mode  :character  
 Mean   :3.621   Mean   :5167                     
 3rd Qu.:4.961   3rd Qu.:5228                     
 Max.   :5.045   Max.   :5228   

3、缺失值处理

用prop.table()和table()函数探索job列、poutcome列和education列的缺失情况，我们可以看到在job列中，空值也就是unkown的占比为0.008，education列中unkown的占比为0.042，poutcome中的空值 nonexistent占比却高达0.863。缺失值之间的差异如此之大，迫使我们要采取不同的处理方法来进行缺失值的处理。

> prop.table(table(mydata$job))

       admin.   blue-collar  entrepreneur     housemaid    management       retired 
  0.253034865   0.224677090   0.035350102   0.025735651   0.070991551   0.041759736 
self-employed      services       student    technician    unemployed       unknown 
  0.034500340   0.096363018   0.021244052   0.163712732   0.024618821   0.008012042 
> prop.table(table(mydata$poutcome))

    failure nonexistent     success 
 0.10323395  0.86343110  0.03333495 
> prop.table(table(mydata$education))

           basic.4y            basic.6y            basic.9y         high.school 
       0.1013887540        0.0556472759        0.1467660484        0.2310138875 
         illiterate professional.course   university.degree             unknown 
       0.0004370205        0.1272943576        0.2954258522        0.0420268039 

job中的异常值数量少，如果为空我们就给它删掉

> mydata$job[which(mydata$job == 'unknown')] <- NA  #将unknown全不赋值NA
> mydata.clear1 <- na.omit(mydata)                  # 把NA的全删了

对于缺失值过多的poutcome列，我们直接生成一个没有poutcome列的新数据

mydata.clear2 <- subset(mydata.clear1,select = -poutcome)  # 直接把poutcome列删掉

对于education，我们直接对空值的部分赋值为最高学历：大学学历

> count <- names(which.min(table(mydata.clear2$education)))  # 找到最高学历
> mydata.clear2$education <- ifelse(mydata.clear2$education != 'unknown',
+                                   as.character(mydata.clear2$education),count) # 将学历的空值赋值为大学学历

注意4

对学历这样处理空值真的好吗？个人觉得不好，但说实话我也不知道该咋整。要不，先对其他数据进行数据挖掘形成预测模型，再对这些学历空值的数据进行预测？但又觉得这样有点本末倒置了，我对数据预处理就是要进行数据挖掘，这样先挖掘再处理空值，再挖掘总觉得有点不得劲，但又莫名的有点合理。但这样值不值得就两说了，为了一点数据而做一次挖掘，值吗？

4、异常值处理

对年龄我们随便修改几个数字（有异常值才能处理），比如我们把第一个人的年龄设为209，第二个人的年龄设为-3

> mydata.clear2[1,1]
[1] 56
> mydata.clear2[1,1]  <- 209
> mydata.clear2[2,1]  
[1] 57
> mydata.clear2[2,1]  <- -3

现在有了异常值，我们采用箱型图来查找异常值

> sp <- boxplot(mydata.clear2$age , boxwex = 0.7)
> xi <- 1.1
> sd.age <- sd(mydata.clear2$age)  # 计算年龄的标准差
> mn.age <- mean(mydata.clear2$age)  # 计算年龄均值
> points(xi,mn.age, col = 'red', pch = 18)  # 在（1.1 ，均值）处标点，红色菱形
> arrows(xi,mn.age - sd.age, xi , mn.age+sd.age,code = 3,
+        col = 'pink',angle = 75,length = .1) # 设置双向箭头图形
> text(rep(c(1.05,1.05),length=length(sp$out)),
+      labels = sp$out[order(sp$out)],sp$out[order(sp$out)] + rep(c(3,3,
+                                                                   130),length = length(sp$out)),col="red") # 给点进行标注

透过箱型图我们便可以得知异常值的数值，就是这密集的我有点看不清了已经。

既然知道了年龄的异常值，那么在经过人工或其他方式确定正确值之后便可以对异常值进行修改

> mydata.clear2$age[which(mydata.clear2$age == -3)] <- 30
> mydata.clear2$age[which(mydata.clear2$age == 209)] <- 47

5、数据分段

接下来需要将age列进行分段处理，生成的新列newage要求小于18的修改为1，大于等于18且小于30的修改为2，大于等于30且小于60的修改为3，大于等于60的修改为4。

> mydata.clear3 <- within(mydata.clear2,{
+     newage = NA              # 给新列赋值为空
+     newage[age < 18] = 1     # 按条件给新列进行赋值
+     newage[age >= 18 & age < 30] = 2
+     newage[age >= 30 & age < 60] = 3
+     newage[age > 60] = 4
+ })
> table(mydata.clear3$newage) #统计

    1     2     3     4 
    5  5639 34050   889 

6、分层抽样

按照mydata.clear3中的y列进行分层抽样，抽取的样本赋值给newdata，抽取mydata.clear3的80%

首先，安装包：creat

# 导入包
library(caret)
# 分层抽样
index <- createDataPartition(mydata.clear3,p=0.8,list=F)
newdata <- mydata.clear3[index,]
prop.table(table(mydata.clear3$y))

       no       yes 
0.8873458 0.1126542 
prop.table(table(newdata$y))

       no       yes 
0.8873252 0.1126748 
dim(newdata)[1]/dim(mydata.clear3)[1]
[1]0.8000147

任务二：R语言基本操作

例一 ：R语言基本运算符号和常用函数的使用

1、基本运算符号

> 2+5             # 加法
[1] 7
> 5 %/% 2         # 整除
[1] 2
> 5 / 2           # 除法
[1] 2.5
> 4.5 %% 2          # 求模
[1] 0.5
> 2>5             # 判断
[1] FALSE
> 3^2             # 幂
[1] 9
> # 逻辑运算
> TRUE | 2        # 逻辑运算或
[1] TRUE
> FALSE & 2       # 逻辑运算与
[1] FALSE
> xor(TRUE,FALSE) # 逻辑运算异或
[1] TRUE
> !TRUE           # 逻辑运算非
[1] FALSE
> # 集合
> x <- c(2,5,3,6)
> y <- c(3,5,7,8)
> union(x,y) # 求并集
[1] 2 5 3 6 7 8
> intersect(x,y) # 求交集
[1] 5 3
> setequal(x,y)  # 判断x与y是否相等
[1] FALSE
> combn(x,3)     # x 中的元素每次取n个的所有组合
     [,1] [,2] [,3] [,4]
[1,]    2    2    2    5
[2,]    5    5    3    3
[3,]    3    6    6    6

2、常用函数

> x <- 5:2
> sum(x)           # 对x中的元素求和
[1] 14
> prod(2:5)        # 5的阶乘
[1] 120
> which.max(x)     # 返回x中的最大值的下标
[1] 1
> mean(x)          # x中所有元素的均值
[1] 3.5
> var(x)           # x中元素的方差（用n-1做分母）
[1] 1.666667
> rev(x)           # 对x中的元素取逆序
[1] 2 3 4 5
> sort(x)          # 将x中的元素取升序
[1] 2 3 4 5
> cumsum(x)        # 求累计和，返回一个向量 它的第i个元素是从x[1] 到 x[i] 的和
[1]  5  9 12 14
> cummin(x)        # 求累计最小值（从左到右）
[1] 5 4 3 2
> cummax(x)        # 求累计最大值（从左到右）
[1] 5 5 5 5

3、正态分布

> # 统计函数--标准正态分布
> x <- pretty(c(-4,4),30)
> y <- dnorm(x)
> plot(x,y,type="l",xlab="NormalDeviate",ylab="Density",yaxs="i")

4、矩阵

> # 数组与矩阵相关函数
> a <- array(1:20,dim=c(2,2,5))
> a
, , 1

     [,1] [,2]
[1,]    1    3
[2,]    2    4

, , 2

     [,1] [,2]
[1,]    5    7
[2,]    6    8

, , 3

     [,1] [,2]
[1,]    9   11
[2,]   10   12

, , 4

     [,1] [,2]
[1,]   13   15
[2,]   14   16

, , 5

     [,1] [,2]
[1,]   17   19
[2,]   18   20

> 
> b <- rbind(a[1,,],a[2,,])  # 合成矩阵
> b
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    3    7   11   15   19
[3,]    2    6   10   14   18
[4,]    4    8   12   16   20
> 
> aperm(b) # 数组转置
     [,1] [,2] [,3] [,4]
[1,]    1    3    2    4
[2,]    5    7    6    8
[3,]    9   11   10   12
[4,]   13   15   14   16
[5,]   17   19   18   20
> dim(b)   # 对象的维向量
[1] 4 5

例二 ：R语言可视化绘图

1、数据准备

> # 生成20个1到5以内的随机数
> m <- round(runif(20,1,5))
> m
 [1] 5 2 1 4 2 5 3 2 2 1 3 4 3 2 1 3 2 3 2 4
> 
> # 统计
> x = table(m)
> x
m
1 2 3 4 5 
3 7 5 3 2 

2、条形图

> # 绘制条形图
> barplot(x,main="统计分布图",xlab="数值",ylab="个数")

3、饼图

> # 绘制饼图
> pie(x,col=c("purple","violetred1","green3","cornsilk"),main="统计分布图")

注意5：

这里在给颜色的时候只给出了四个颜色，但应该给出5种颜色更合理。颜色不仅可以按照本例中直接用英文的方式给出，也可按照下面的方式给出。

cols = c("#ED1C24","#22B14C","#FFC90E","#3f48CC","#EB36E0")
pie(x,col=cols,main="统计分布图")

4、直方图

> # 绘制直方图
> hist(m,main="统计分布图",xlab="数值",ylab="个数")

5、核密度图

> # 绘制核密度图
> plot(density(m))

6、点图

> # 绘制点图
> dotchart(m,col="purple",main="统计分布图",xlab="数值",ylab="个数")

例三 :R语言数据分析

1、描述性统计分析

已知15位学生的体重（单位：千克），利用数据进行描述性统计分析。

75.0 64.0 47.4 66.9 62.2 62.2 58.7 63.5 66.6 64.0 57.0 69.0 56.9 50.0 72.0

> # 数据准备
> w <- c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)

计算描述性统计量

> # 通过summary(w)计算描述性统计量
> summary(w)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  47.40   57.85   63.50   62.36   66.75   75.00 

> # 通过sapply()计算描述性统计量
> mystats <- function(x, na.omit=FALSE){
+   if(na.omit)
+       x<-x[!is.na(x)]
+     m<-mean(x)
+     n<-length(x)
+     s<-sd(x)
+     skew<-sum((x-m)^3/s^3)/n
+     kurt<-sum((x-m)^4/s^4)/n-3
+     return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
+ }
> sapply(w,mystats)
         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
n           1    1  1.0  1.0  1.0  1.0  1.0  1.0  1.0     1     1     1
mean       75   64 47.4 66.9 62.2 62.2 58.7 63.5 66.6    64    57    69
stdev      NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
skew       NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
kurtosis   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
         [,13] [,14] [,15]
n          1.0     1     1
mean      56.9    50    72
stdev       NA    NA    NA
skew        NA    NA    NA
kurtosis    NA    NA    NA

> # 通过psych包中的describe()计算描述性统计量
> library(psych)
> describe(w)
   vars  n  mean   sd median trimmed  mad  min max range  skew kurtosis
X1    1 15 62.36 7.51   63.5   62.54 7.12 47.4  75  27.6 -0.35    -0.65
     se
X1 1.94

2、多元统计分析

数据准备：已知有数据集

25.6 22.2 28.0 29.8 24.4 30.0 29.0 27.5 25.0 27.7 23.0 32.2 28.8 28.0 31.5 25.9 20.6 21.2 22.0 21.2

数据集用5个因子水平测量，问是否存在差异。

> # 输入数据
> x <- c(25.6,22.2,28.0,29.8,24.4,30.0,29.0,27.5,25.0,27.7,23.0,32.2,28.8,28.0,31.5,25.9,20.6,21.2,22.0,21.2)
> 
> # 对数据进行格式转换
> # gl指定因子，5是水平，4是重复次数
> b <- data.frame(x,a=gl(5,4,20))
> b
      x a
1  25.6 1
2  22.2 1
3  28.0 1
4  29.8 1
5  24.4 2
6  30.0 2
7  29.0 2
8  27.5 2
9  25.0 3
10 27.7 3
11 23.0 3
12 32.2 3
13 28.8 4
14 28.0 4
15 31.5 4
16 25.9 4
17 20.6 5
18 21.2 5
19 22.0 5
20 21.2 5

考察方差齐次性（用bartlett检验）

> # 考察方差齐次性
> bartlett.test(x~a,data=b)

	Bartlett test of homogeneity of variances

data:  x by a
Bartlett's K-squared = 7.0966, df = 4, p-value =
0.1309

进行方差分析

> # 方差分析
> m1 <- aov(x~a,data=b)
> summary(m1)
            Df Sum Sq Mean Sq F value Pr(>F)  
a            4  132.0   32.99   4.306 0.0162 *
Residuals   15  114.9    7.66                 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> # 考察具体的差异（多重比较）
> TukeyHSD(m1)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = x ~ a, data = b)

$a
      diff        lwr        upr     p adj
2-1  1.325  -4.718582  7.3685818 0.9584566
3-1  0.575  -5.468582  6.6185818 0.9981815
4-1  2.150  -3.893582  8.1935818 0.8046644
5-1 -5.150 -11.193582  0.8935818 0.1140537
3-2 -0.750  -6.793582  5.2935818 0.9949181
4-2  0.825  -5.218582  6.8685818 0.9926905
5-2 -6.475 -12.518582 -0.4314182 0.0330240
4-3  1.575  -4.468582  7.6185818 0.9251337
5-3 -5.725 -11.768582  0.3185818 0.0675152
5-4 -7.300 -13.343582 -1.2564182 0.0146983

可以看出，5-1，5-2，5-3，5-4之间的差异都很明显。

3、判别分析

选用kknn软件包中的miete数据集进行算法演示。miete数据集记录了1994年慕尼黑的住房佣金标准中一些有趣变量，比如房子的面积、是否有浴室、是否有中央供暖、是否供应热水等等，这些因素都影响着佣金的高低。

> # 导入
> library(kknn)
> data(miete) 
> head(miete) # 展示6个样本数据
       nm wfl     bj bad0 zh ww0 badkach fenster kueche
1  693.29  50 1971.5    0  1   0       0       0      0
2  736.60  70 1971.5    0  1   0       0       0      0
3  732.23  50 1971.5    0  1   0       0       0      0
4 1295.14  55 1893.0    0  1   0       0       0      0
5  394.97  46 1957.0    0  0   1       0       0      0
6 1285.64  94 1971.5    0  1   0       1       0      0
  mvdauer bjkat wflkat      nmqm rooms nmkat adr wohn
1       2     4      1 13.865800     1     3   2    2
2      26     4      2 10.522857     3     3   2    2
3       1     4      1 14.644600     1     3   2    2
4       0     1      2 23.548000     3     5   2    2
5      27     3      1  8.586304     3     1   2    2
6       2     4      3 13.677021     4     5   2    2
> dim(miete)  #查看样本和变量个数
[1] 1082   17

一共有1082条数据样本，样本维度为17

> # 查看变量信息
> summary(miete)
       nm              wfl               bj       bad0    
 Min.   : 127.1   Min.   : 20.00   Min.   :1800   0:1051  
 1st Qu.: 543.6   1st Qu.: 50.25   1st Qu.:1934   1:  31  
 Median : 746.0   Median : 67.00   Median :1957           
 Mean   : 830.3   Mean   : 69.13   Mean   :1947           
 3rd Qu.:1030.0   3rd Qu.: 84.00   3rd Qu.:1972           
 Max.   :3130.0   Max.   :250.00   Max.   :1992           
 zh      ww0      badkach fenster  kueche     mvdauer     
 0:202   0:1022   0:446   0:1024   0:980   Min.   : 0.00  
 1:880   1:  60   1:636   1:  58   1:102   1st Qu.: 2.00  
                                           Median : 6.00  
                                           Mean   :10.63  
                                           3rd Qu.:17.00  
                                           Max.   :82.00  
 bjkat   wflkat       nmqm            rooms       nmkat  
 1:218   1:271   Min.   : 1.573   Min.   :1.000   1:219  
 2:154   2:513   1st Qu.: 8.864   1st Qu.:2.000   2:230  
 3:341   3:298   Median :12.041   Median :3.000   3:210  
 4:226           Mean   :12.647   Mean   :2.635   4:208  
 5: 79           3rd Qu.:16.135   3rd Qu.:3.000   5:215  
 6: 64           Max.   :35.245   Max.   :9.000          
 adr      wohn   
 1:  25   1: 90  
 2:1035   2:673  
 3:  22   3:319  

可以选择nmkat(净租金)作为待判别变量：一是，由于该变量在含义上容易受其他变量影响，为被解释变量；二是，nmkat自身就有5个等级类别，其相应的样本量依次为 219、230、210、208、215，即每一类的样本量都为200个左右，分布比较均匀。

为了提高判别效果，我们考虑采用分层抽样的方式，由于miete数据集中，待判别变量nmkat的5个等级分布比较均匀，因此采用5个等级按等量抽取样本。（如果分布不均匀，则采用按比例抽取样本）

> #以nmkat变量的5个等级划分层次，进行分层抽样  
> sub_train = strata(miete,stratanames = "nmkat",size=rep(n,5),method="srswor")  
> head(sub_train)
   nmkat ID_unit      Prob Stratum
1      3       1 0.6857143       1
3      3       3 0.6857143       1
8      3       8 0.6857143       1
16     3      16 0.6857143       1
20     3      20 0.6857143       1
22     3      22 0.6857143       1

这样我们就完成了分层抽样，对每一个按照 nmkat 进行划分的分层都从中无放回的抽出 144 个样本数据，所以新生成的sub_train中一共有720条数据。

接下来是形成训练集

> # 获取如上ID_unit所对应的样本构成训练集，并删除变量1、3、12  
> data_train = getdata(miete[ c(sub_train$ID_unit),],1) 
> data_train = data_train[c(-1,-3,-12)]
> dim(data_train)     #显示训练集的维度
[1] 720  15

注意6：

老师，我觉得知道里的这里是有错误的：

第一是我认为我们从1082条数据里抽出720条数据是用来做测试集的，并且在原本的代码中测试集里样本数为0根本就没有意义。

第二是我认为原本的代码有错误，从源数据集中形成训练集，那么训练集中的样本数绝不可能超过源数据集的数据量。

以上是学生的想法，如有错误之处，望指正。

> # 查看
> head(data_train)
       nm     bj bad0 zh ww0 badkach fenster kueche mvdauer wflkat
1  693.29 1971.5    0  1   0       0       0      0       2      1
2  736.60 1971.5    0  1   0       0       0      0      26      2
3  732.23 1971.5    0  1   0       0       0      0       1      1
16 812.80 1971.5    0  1   0       1       0      0      12      3
22 765.39 1971.5    0  1   0       1       0      0       3      3
27 759.20 1957.0    0  1   0       1       0      0       0      2
        nmqm rooms nmkat adr wohn
1  13.865800     1     3   2    2
2  10.522857     3     3   2    2
3  14.644600     1     3   2    2
16  8.739785     4     3   2    2
22  8.230000     4     3   2    1
27  9.989474     3     3   2    2

接下来进行线性判别

> library(MASS) 
> fit_lda1 = lda(nmkat ~.  , data_train) #以公式格式进行线性判别 
> names(fit_lda1) 
 [1] "prior"   "counts"  "means"   "scaling" "lev"     "svd"    
 [7] "N"       "call"    "terms"   "xlevels"

> fit_lda1$prior 
  1   2   3   4   5 
0.2 0.2 0.2 0.2 0.2 

可以看到，各类别的先验概率在5个等级中都为0.2，之和为1，即它们都相等，这与它们分别均匀对应。

> fit_lda1$means
         nm       bj       bad01       zh1        ww01
1  384.1506 1938.753 0.069444444 0.5833333 0.138888889
2  596.8247 1949.931 0.027777778 0.8541667 0.034722222
3  757.2785 1948.049 0.020833333 0.8472222 0.034722222
4  977.7442 1950.267 0.013888889 0.8958333 0.034722222
5 1475.5799 1951.434 0.006944444 0.9583333 0.006944444
   badkach1   fenster1    kueche1   mvdauer    wflkat.L
1 0.3750000 0.06250000 0.04166667 13.486111 -0.27498597
2 0.5694444 0.07638889 0.05555556 10.590278 -0.13258252
3 0.5555556 0.05555556 0.05555556 12.222222 -0.07365696
4 0.6388889 0.02777778 0.11111111 10.312500  0.14240345
5 0.8263889 0.02777778 0.21527778  4.861111  0.39283710
     wflkat.Q      nmqm    rooms        adr.L
1 -0.05103104  8.481306 2.104167 -0.014731391
2 -0.19561897 11.375869 2.333333 -0.034373246
3 -0.29768105 12.946274 2.506944 -0.004910464
4 -0.31469139 14.479111 2.819444  0.014731391
5 -0.05103104 16.887573 3.361111  0.044194174
       adr.Q     wohn.L     wohn.Q
1 -0.7739707 0.08347788 -0.3997431
2 -0.7569604 0.10803020 -0.3572173
3 -0.7739707 0.12767206 -0.3912379
4 -0.7739707 0.18659762 -0.3912379
5 -0.7059293 0.30935922 -0.1786086

从上面的结果中，可以看到一些很能反映现实情况的数据特征。比如，住房面积wfl变量，它明显随着租金nmkat的升高而逐步提高。这与我们的常识“房子的面积越大，租金就越贵”是十分吻合的。

4、聚类分析

这里我们选用R语言自带包USArrests，美国50个洲的4个犯罪指标，利用K-means进行聚类分析。

首先，导入

> head(USArrests)
           Murder Assault UrbanPop Rape
Alabama      13.2     236       58 21.2
Alaska       10.0     263       48 44.5
Arizona       8.1     294       80 31.0
Arkansas      8.8     190       50 19.5
California    9.0     276       91 40.6
Colorado      7.9     204       78 38.7

进行聚类分析

> # 利用K-means进行聚类分析
> kmeans.result=kmeans(USArrests,3)
> kmeans.result
K-means clustering with 3 clusters of sizes 14, 16, 20

Cluster means:
     Murder  Assault UrbanPop     Rape
1  8.214286 173.2857 70.64286 22.84286
2 11.812500 272.5625 68.31250 28.37500
3  4.270000  87.5500 59.75000 14.39000

Clustering vector:
       Alabama         Alaska        Arizona 
             2              2              2 
      Arkansas     California       Colorado 
             1              2              1 
   Connecticut       Delaware        Florida 
             3              2              2 
       Georgia         Hawaii          Idaho 
             1              3              3 
      Illinois        Indiana           Iowa 
             2              3              3 
        Kansas       Kentucky      Louisiana 
             3              3              2 
         Maine       Maryland  Massachusetts 
             3              2              1 
      Michigan      Minnesota    Mississippi 
             2              3              2 
      Missouri        Montana       Nebraska 
             1              3              3 
        Nevada  New Hampshire     New Jersey 
             2              3              1 
    New Mexico       New York North Carolina 
             2              2              2 
  North Dakota           Ohio       Oklahoma 
             3              3              1 
        Oregon   Pennsylvania   Rhode Island 
             1              3              1 
South Carolina   South Dakota      Tennessee 
             2              3              1 
         Texas           Utah        Vermont 
             1              3              3 
      Virginia     Washington  West Virginia 
             1              1              3 
     Wisconsin        Wyoming 
             3              1 

Within cluster sum of squares by cluster:
[1]  9136.643 19563.863 19263.760
 (between_SS / total_SS =  86.5 %)

Available components:

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

> # 画出"Assault"和 "Rape"的K-means聚类分析图
> plot(USArrests[c("Assault","Rape")],col=kmeans.result$cluster)
> points(kmeans.result$centers[ c("Assault" ,"Rape")] ,col = 1:3,pch = 8,cex=2)

由上图可知，"Assault"和 "Rape"关系被分为三类。

5、主成分分析

我们对128个成年男子身材测量16项指标：身高、坐高、胸围、头高、裤长、下档、手长、领围、前胸、后背、肩厚、肩宽、袖长、肋围、腰围、腿肚，分别为X1，X2，X3，…相关矩阵如下表，从相关矩阵R出发进行主成分分析，对16项指标进行分类。

> x<-c(1.00,
+      0.79, 1.00,
+      0.36, 0.31, 1.00,
+      0.96, 0.74, 0.38, 1.00,
+      0.89, 0.58, 0.31, 0.90, 1.00,
+      0.79, 0.58, 0.30, 0.78, 0.79, 1.00,
+      0.76, 0.55, 0.35, 0.75, 0.74, 0.73, 1.00,
+      0.26, 0.19, 0.58, 0.25, 0.25, 0.18, 0.24, 1.00,
+      0.21, 0.07, 0.28, 0.20, 0.18, 0.18, 0.29,-0.04, 1.00,
+      0.26, 0.16, 0.33, 0.22, 0.23, 0.23, 0.25, 0.49,-0.34, 1.00,
+      0.07, 0.21, 0.38, 0.08,-0.02, 0.00, 0.10, 0.44,-0.16, 0.23,1.00,
+      0.52, 0.41, 0.35, 0.53, 0.48, 0.38, 0.44, 0.30,-0.05, 0.50,0.24, 1.00,
+      0.77, 0.47, 0.41, 0.79, 0.79, 0.69, 0.67, 0.32, 0.23, 0.31,0.10, 0.62, 1.00,
+      0.25, 0.17, 0.64, 0.27, 0.27, 0.14, 0.16, 0.51, 0.21, 0.15,0.31, 0.17, 0.26, 1.00,
+      0.51, 0.35, 0.58, 0.57, 0.51, 0.26, 0.38, 0.51, 0.15, 0.29,0.28, 0.41, 0.50, 0.63, 1.00,
+      0.21, 0.16, 0.51, 0.26, 0.23, 0.00, 0.12, 0.38, 0.18, 0.14,0.31, 0.18, 0.24, 0.50, 0.65, 1.00)

> # 输入变量名称并将矩阵生成相关矩阵
> names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9","X10", "X11", "X12", "X13", "X14", "X15", "X16")
> R<-matrix(0, nrow=16, ncol=16, dimnames=list(names, names))
> for (i in 1:16){for (j in 1:i){R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]}}

> # 作主成分分析
> pr<-princomp(covmat=R); load<-loadings(pr)
> 
> # 画散点图
> plot(load[,1:2]); text(load[,1], load[,2], adj=c(-0.4, 0.3))

, 0.28, 0.20, 0.18, 0.18, 0.29,-0.04, 1.00,

•  0.26, 0.16, 0.33, 0.22, 0.23, 0.23, 0.25, 0.49,-0.34, 1.00,
  

•  0.07, 0.21, 0.38, 0.08,-0.02, 0.00, 0.10, 0.44,-0.16, 0.23,1.00,
  

•  0.52, 0.41, 0.35, 0.53, 0.48, 0.38, 0.44, 0.30,-0.05, 0.50,0.24, 1.00,
  

•  0.77, 0.47, 0.41, 0.79, 0.79, 0.69, 0.67, 0.32, 0.23, 0.31,0.10, 0.62, 1.00,
  

•  0.25, 0.17, 0.64, 0.27, 0.27, 0.14, 0.16, 0.51, 0.21, 0.15,0.31, 0.17, 0.26, 1.00,
  

•  0.51, 0.35, 0.58, 0.57, 0.51, 0.26, 0.38, 0.51, 0.15, 0.29,0.28, 0.41, 0.50, 0.63, 1.00,
  

•  0.21, 0.16, 0.51, 0.26, 0.23, 0.00, 0.12, 0.38, 0.18, 0.14,0.31, 0.18, 0.24, 0.50, 0.65, 1.00)
  

```R
> # 输入变量名称并将矩阵生成相关矩阵
> names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9","X10", "X11", "X12", "X13", "X14", "X15", "X16")
> R<-matrix(0, nrow=16, ncol=16, dimnames=list(names, names))
> for (i in 1:16){for (j in 1:i){R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]}}

> # 作主成分分析
> pr<-princomp(covmat=R); load<-loadings(pr)
> 
> # 画散点图
> plot(load[,1:2]); text(load[,1], load[,2], adj=c(-0.4, 0.3))

[外链图片转存中…(img-Ro96xItl-1651563236695)]

上图右上角的点X1,X2,X4,X5,X6,X7,X13可以看作一类为”长”类,左下角的点X3,X8,X11,X14,X15,X16可以看作一类为”围”类,中间的点X9,X10,X12可以看作一类为体型特征指标。