一、多元正态的参数估计

1.1 样本均值

在R语言中，均值通常用函数mean()得到，但是mean()只能计算一维变量的样本均值，在面对多元随机变量的样本时，假设我们以数据框的形式保存样本，我们有以下方法可以得到样本均值：

对多元样本的每一个分量用mean()函数，可以用apply()或sapply()函数
以数据框类型保存的样本，可以用summary()函数返回各个变量的各项描述性数据，其中包括均值

例1.1：

计算有割草机家庭的样本均值向量

有割草机家庭		无割草机家庭
x1	x2	x1	x2
20.0	9.2	25.0	9.8
28.5	8.4	17.6	10.4
21.6	10.8	21.6	8.6
20.5	10.4	14.4	10.2
29.0	11.8	28.0	8.8
36.7	9.6	16.4	8.8
36.0	8.8	19.8	8.0
27.6	11.2	22.0	9.2
23.0	10.0	15.8	8.2
31.0	10.4	11.0	9.4
17.0	11.0	17.0	7.0
27.0	10.0	21.0	7.4

注：在输入数据时，通常会用一个新的变量(假设命名为Y)来表示每个观测所属的组，称为分组变量，这个变量在R中通常要转换成因子

> data4.1=read.csv('Table4-1.csv')
> head(data4.1)
###                y是分组变量，y=1表示有割草机家庭x1   x2 y
1  20.0  9.2 1
2  28.5  8.4 1
3  21.6 10.8 1
4  20.5 10.4 1
5  29.0 11.8 1
6  36.7  9.6 1
###
> apply(data4.1[data4.1$y==1,],2,mean)[1:2]    #用apply()函数运算
###x1       x2 
26.49167 10.13333 
###
> summary(data4.1[data4.1$y==1,])           #用summary()获取各分量的样本均值
###x1              x2              y    Min.   :17.00   Min.   : 8.40   Min.   :1  1st Qu.:21.32   1st Qu.: 9.50   1st Qu.:1  Median :27.30   Median :10.20   Median :1  Mean   :26.49   Mean   :10.13   Mean   :1   #该行为均值3rd Qu.:29.50   3rd Qu.:10.85   3rd Qu.:1  Max.   :36.70   Max.   :11.80   Max.   :1 
###

apply()用法：apply(A,margin,fun,...)

apply()函数用来对矩阵或数据框的每行或每列进行指定函数的运算。其中A为矩阵或数据框；margin指定对行或对列进行运算，当margin=1时对行进行运算，当margin=2时对列进行运算；fun是指定的函数

summary()用法：summary(object,...)

summary()多用于获取项目的摘要，包含部分信息。当object为数据框时，会返回各个变量的五数（最小值，下四分位数，中位数，上四分为数，最大值）和均值

1.2 样本协差阵

在R中，样本协差阵的获取非常简便，对数据框使用cov()函数即可

例1.2：

继上题，计算有割草机组的样本协差阵

> cov(data4.1[data4.1$y==1,][,1:2])
###x1        x2
x1 39.182652 -1.969697
x2 -1.969697  1.020606
###

cov()用法：cov(x,y=NULL,...)

当指定cov()的参数x和y，且两者都为一维向量时，会返回两个向量的样本协方差；而未指定参数y，且x为矩阵或数据框时，会返回以x每一列作为变量样本的协差阵

1.3 样本相关阵

获取样本相关阵的函数是cor()，其用法与cov()相同，两个一维向量返回相关系数；数据框返回相关阵

二、各类检验

2.1 正态性检验

正态性检验即检验样本是否来自正态总体的检验，原假设都为来自正态总体。正态性的检验方法有许多种，此处介绍小样本量(3~50)时所用的夏皮洛-威尔克检验。R中的夏皮洛-威尔克检验的函数为shapiro.test()

shapiro.test()一次只能对一维变量进行正态性检验，当面对多元随机变量的样本时，有以下方法

我们可以对其每一个分量都进行一次正态性检验，当所有分量都检验得出服从正态分布后，可以认为该多元随机变量服从多元正态分布
运用mvnormtest包内的mshapiro.test()函数进行多元正态性检验

实现时可能会用到的函数有：

sapply()，对每个分量进行指定的检验
tapply()，对以分组变量指定的不同组别分别进行指定的检验

例2.1：

继上题，对不同类型家庭的随机向量数据进行正态性检验

> sapply(data4.1[,-3],shapiro.test)    #对各分量进行正态性检验，但是未分组
###x1                            x2                           
statistic 0.9654387                     0.9880936                    
p.value   0.5568611                     0.9897171                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]" 
###
> tapply(data4.1[,1],data4.1$y,shapiro.test)
###                            对分组后的x1进行正态性检验
$`0`Shapiro-Wilk normality testdata:  X[[i]]
W = 0.98551, p-value = 0.9971$`1`Shapiro-Wilk normality testdata:  X[[i]]
W = 0.95332, p-value = 0.6859###
> tapply(data4.1[,2],data4.1$y,shapiro.test)
###                            对分组后的x2进行正态性检验
$`0`Shapiro-Wilk normality testdata:  X[[i]]
W = 0.97557, p-value = 0.9596$`1`Shapiro-Wilk normality testdata:  X[[i]]
W = 0.98262, p-value = 0.992
###
##对有割草机家庭的随机向量数据进行正态性检验
> mshapiro.test(t(as.matrix(data4.1[data4.1$y==1,-3])))Shapiro-Wilk normality testdata:  Z
W = 0.96877, p-value = 0.8975##对无割草机家庭的随机向量数据进行正态性检验
> mshapiro.test(t(as.matrix(data4.1[data4.1$y==0,-3])))Shapiro-Wilk normality testdata:  Z
W = 0.98001, p-value = 0.9837

检验结果为均服从正态分布

sapply()：sapply(X,Fun,...)

sapply()用于对X的每个分量进行Fun函数运算，X应该是矩阵或数据框

tapply()：tapply(X,Index,Fun=NULL,...)

tapply()用于对以分组变量Index指示的每个组中对应的X的数据进行Fun函数运算

mshapiro.test()：mshapiro.test(U)

mshapiro.test()用于进行多元的夏皮洛-威尔克正态性检验，需要注意U只能是数据矩阵，当遇到用数据框存储的数据时要用as.matrix()转化为矩阵，且这个函数默认变量的数据按行排放，通常我们需要对矩阵再进行一次转置

另外，可以画出Q-Q图查看样本的正态性，常用的函数有qqnorm()和qqline()

qqnorm(x)，其中x为一维变量的样本，当画出的散点图越接近一条斜线，其正态性越强
qqline(x)，其中x为一维变量的样本，当画出散点的Q-Q图后，添加点所靠近的斜线，该斜线的斜率为标准差,截距为均值

2.2 均值向量的检验

一维的均值检验有很多，若样本服从正态分布我们可以用t.test()单个总体或双总体的t检验；若不服从正态分布，我们可以用wilcox.test()进行秩和检验，用法与t.test()类似；当遇到多个总体时，若各个变量的方差相差不大，我们可以用将各个变量的数据放到一列，然后用一个分组变量表示数据属于哪个变量，运用aov()进行方差分析，从而进行多总体的均值检验

当遇到多元随机变量的均值检验时，我们有以下方法：

对每个分量进行均值检验，通过正态性检验的用t检验，未通过正态性检验的用秩和检验
对通过多元正态检验的数据，运用ICSNP包中的HotellingsT2()函数进行均值向量的检验
多总体时，若协差阵齐性检验通过，可以用manova()进行多元方差分析

例2.2.1：

继上题，检验有割草机家庭和无割草机家庭的向量均值是否相同

#上题已得出各类家庭的数据均通过正态检验> t.test(x1~y,data=data4.1)            #对x1分量进行t检验Welch Two Sample t-testdata:  x1 by y
t = -3.2508, df = 20.458, p-value = 0.003919
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:-12.073229  -2.643437
sample estimates:
mean in group 0 mean in group 1 19.13333        26.49167 > t.test(x2~y,data=data4.1)            #对x2分量进行t检验Welch Two Sample t-testdata:  x2 by y
t = -3.1203, df = 21.956, p-value = 0.004991
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:-2.1918725 -0.4414609
sample estimates:
mean in group 0 mean in group 1 8.816667       10.133333 > library(ICSNP)
> attach(data4.1)
> HotellingsT2(cbind(x1,x2)~y)            #霍特林T方检验，用于多元正态向量Hotelling's two sample T2-testdata:  cbind(x1, x2) by y
T.2 = 12.257, df1 = 2, df2 = 21, p-value = 0.000297
alternative hypothesis: true location difference is not equal to c(0,0)

对各分量进行t检验的结果是：有割草机家庭与无割草机家庭的两个分量都不相等

对向量整体进行霍特林T方检验的结果是：有割草机家庭与无割草机家庭的向量均值不相等

该题也可以用多元方差分析，只是水平数为2，其实也是通过对每个分量进行检验

例2.2.2：

现有如下表所示各省统计数据，试检验它们的均值向量是否等于 (1081,1822,115,179)

序号	省份	工资性收入	家庭性收入	财产性收入	转移性收入
1	北京	4524.25	1778.33	588.04	455.64
2	天津	2720.85	2626.46	152.88	79.64
3	河北	1293.50	1988.58	93.74	105.81
4	山西	1177.94	1563.52	62.70	86.49
5	内蒙古	504.46	2223.26	73.05	188.10
6	辽宁	1212.20	2163.49	113.24	201.28

注：以上为部分表格，表格全部内容在此不展示

> data3=read.csv('Table_0.csv',encoding='UTF-8')    #读取数据
> mu_bar=c(1081,1882,115,179)
> rownames(data3)=data3[,1]            #将省名赋给数据框行名
> data3=data3[,-1]                     #去除省名一列
###
假设通过了正态性检验
###
> HotellingsT2(data3,mu=mu_bar)Hotelling's one sample T2-testdata:  data3
T.2 = 1.8443, df1 = 4, df2 = 27, p-value = 0.1494
alternative hypothesis: true location is not equal to c(1081,1882,115,179)

检验结果为各省统计数据的均值向量等于(1081,1822,115,179)

例2.2.3：

在数据New drug.xls中，各变量的意义为drug（药），取值1表示对病人给以新药，取值2表示对病人给以安慰剂，resp1-resp3是治疗后病人三个时点的呼吸状况，pulse1-pulse3是病人三个时点的脉搏。试分析这两方法的各次重复测定均值向量是否有显著差异？

drug	resp1	resp2	resp3	pulse1	pulse2	pulse3
1	3.4	3.3	3.3	2.2	2.1	2.1
1	3.4	3.4	3.3	2.2	2.1	2.2
1	3.3	3.4	3.4	2.3	2.4	2.3
2	3.3	3.3	3.3	2.8	2.9	2.7
2	3.2	3.3	3.4	2.6	2.7	2.7
2	3.2	3.2	3.2	2.7	2.9	2.7

注：以上为部分表格，表格全部内容在此不展示

题目要求检验不用用药的组之间，向量(resp1,resp2,resp3,pulse1,pulse2,pulse3)的均值是否相等。因为drug只有2个水平，可以对每个分量进行t检验，但是分量比较多会比较麻烦；也可以用多元方差分析，查看结果也是对每个分量的检验，不过需要先进行协差阵检验；用霍特林T2检验会比较简单。

> data_drug=read.csv('new drug.csv',encoding='UTF-8')
> names(data_drug)[1]='drug'        #UTF-8格式的csv文件读取后，第一列的名字会有变动，此处改回
> attach(data_drug)
###
对每个变量进行正态性检验后
得知随机向量不服从多元正态分布
因此不能用t检验和霍特林T方检验，不过可以对每个分量进行秩和检验
假设数据通过了协差阵检验
接下来进行多元方差分析
###
> modle_drug=manova(cbind(resp1,resp2,resp3,pulse1,pulse2,pulse3)~drug,data=data_drug)
> summary.aov(modle_drug)Response resp1 :Df   Sum Sq  Mean Sq F value   Pr(>F)   
drug         1 0.040833 0.040833  14.412 0.003507 **
Residuals   10 0.028333 0.002833                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response resp2 :Df   Sum Sq  Mean Sq F value   Pr(>F)   
drug         1 0.040833 0.040833  14.412 0.003507 **
Residuals   10 0.028333 0.002833                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response resp3 :Df   Sum Sq   Mean Sq F value Pr(>F)
drug         1 0.020833 0.0208333  3.0488 0.1114
Residuals   10 0.068333 0.0068333               Response pulse1 :Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 0.65333 0.65333      70 7.936e-06 ***
Residuals   10 0.09333 0.00933                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response pulse2 :Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 1.08000 1.08000  79.024 4.623e-06 ***
Residuals   10 0.13667 0.01367                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Response pulse3 :Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 0.75000 0.75000  64.286 1.155e-05 ***
Residuals   10 0.11667 0.01167                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

只有resp3的检验结果为相同，其他都为不同，所以认为两方法的各次重复测定均值向量有显著差异

t.test()用法:

t.test(x,y=NULL,alternative=c("two.sided","less","greater"),mu=0,...)或

t.test(formula,data,...)

指定x未指定y时，进行单总体的t检验，可以指定mu检验其是否与mu相同

alternative指定双边检验或左尾检验或右尾检验

也可以用formula类的参数，即x~y的类型，y是分组变量，会对不同组的x进行双总体t检验

wilcox.test()用法：

wilcox.test()用法与t.test()相同，此处不赘述

HotellingsT2()用法:

HotellingsT2(X,Y=NULL,mu=NULL,test="f",...)或

HotellingsT2(formula,...)

X与Y为矩阵或数据框，未指定Y时进行单总体的检验，可以指定mu检验其是否与mu相同

test参数指定近似统计量，默认为f，即F近似，可以指定"chi"，即卡方近似

可以用formula类参数，与先前用法相同，但是HotellingsT2()没有data参数

manova()用法:

manova(formula,data,...)

manova()的formula参数用法aov()类似，manova()返回的是多元方差分析的模型，将其赋给某个变量，然后用aov.summary()函数可以看每个变量的检验

2.3 协差阵检验

在进行多元方差分析前需要进行协差阵齐性检验，协差阵检验可以用heplots包内的boxM()函数。

例2.3：

继有无割草机家庭数据，检验两组的协差阵是否有差异

> boxM(data4.1[,-3],group=data4.1[,3])Box's M-test for Homogeneity of Covariance Matricesdata:  data4.1[, -3]
Chi-Sq (approx.) = 0.99346, df = 3, p-value = 0.8028

检验结果为两组协差阵相同

boxM()用法：

boxM(formula,data,...)或

boxM(Y,group,...)

formula类参数的用法与之前的函数相同

Y是数据矩阵或数据框，group是指定的分组变量

boxM()函数进行的是协差阵齐性检验，在分组变量的水平数大于2时也可以使用

三、小结

总结本文提到的函数和应用场景

参数估计		正态性检验
函数	应用场景	函数	应用场景
mean()	计算一维变量的样本均值	shapiro.test()	小样本正态性检验
apply()	对矩阵或数据框的行或列进行运算	mshapiro.test()	多元小样本正态性检验
sapply()	对矩阵或数据框的每个变量进行运算	sapply()	对每个变量进行指定运算或检验
summary()	对数据框使用时返回每个变量的统计描述	tapply()	对以分组变量指定的不同组别分别进行指定的运算或检验
cov()	获取协方差或协差阵	qqnorm()	画Q-Q图
cor()	获取相关系数或相关阵	qqline()	在Q-Q图中添加正态标准线

均值向量检验		协差阵检验
函数	应用场景	函数	应用场景
t.test()	正态样本的单双总体均值检验	boxM()	协差阵齐性检验
wilcox.test()	非正态样本的单双总体均值检验
HotellingsT2()	多元正态样本的单双总体均值检验
aov()	方差齐性情况下的方差分析
manova()	协差阵齐性下的多元方差分析
aov.summary()	获取方差分析模型的检验结果