一、多元正態的參數估計

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()	獲取方差分析模型的檢驗結果