上一次笔记刚刚谈过多元正态分布检验的STATA操作，方法相对单一，本次笔记继续这个这个话题，这次换用R。不是说R是靠程序命令来完成操作的吗？是的，不妨先来一段感受一下：
setwd("D:/Temp")
library(foreign)
multivnorm <- read.dta("Multivariate.dta")
U <-t(multivnorm[3:4])
library(mvnormtest)
mshapiro.test(U)
有没有头皮发麻焦躁不安？有没有恶心反胃生无可恋？我刚看的时候就这感觉。但无论是什么统计软件，在进行分析时无外乎这么几个步骤：按要求数据录入（输入或导入）、调用分析程序（菜单或命令）、统计结果的专业解读，R也是如此。而且不要忘了统计的关键在于选择合适的统计方法并对结果作出合理的解释，而不在于如何操作软件，统计软件仅仅是解决复杂统计计算的快速实现问题。我们无需记住或者会在R的操作台上敲出命令，我们只要知道用什么样的统计方法，统计软件能实现这个方法的过程或者命令在什么地方或者在哪里可以找到就足够了，然后照着葫芦画瓢就好。我们不生产命令，我们只是程序包的搬运工！R是一个开源项目，而且完全免费！完全免费！完全免费！重要的事情要讲三遍。R拥有数量庞大的被称为程序包（Package）的功能模块可以快速实现各种功能。包是R函数、数据、预编译代码以一种定义完善的格式组成的集合，放在一个被称为库（library）的目录中，除了一系列默认包，还可以通过下载安装各种功能的程序包。使用程序包前先要进行加载。R基本上囊括了在其他软件中尚不可用的、先进的统计方法，更新速度也非常快。单就多元正态性的检验，R中就有N多程序包可以实现，比如程序包mvnormtest【函数mshapiro.test】、程序包mvShapiroTest【mvShapiro.Test】、程序包MVN【函数mvn】、程序包mvnTest【函数AD.test、CM.test、DH.test、HZ.test、R.test、S2.test等】、程序包mvtnorm【函数pmvnorm】、程序包energy【函数mvnorm.e、mvnorm.test、mvnorm.etest】……安装程序包：以上有些包已经在默认存在库里了，但有些需要下载安装。以程序包mvnormtest的为例，下载安装如下：程序包>>装程序包，选择镜像下载地址，选择“mvnormtest”，确定。通过命令install.packages()也可实现上述过程，install.packages("程序包名称")可实现对某个包的安装。同样可安装其他包，包安装后即存在库中，以后使用时直接调用即可。
改变工作目录：文件>>改变工作目录…，选择目标文件夹即可。也可直接使用setwd("工作目录")，本例setwd("D:/Temp")，表示工作目录设置在D盘Temp文件夹。注意路径中采用的是反斜杠“/”而不是“\”。命令getwd()则会显示当前的工作目录。数据载入：可以直接键盘输入也可以导入其他文件。数据输入时要了解R中用于储存数据的对象类型，它包括标量、向量、矩阵、数组、数据框和列表，后面使用的数据包函数也会有固定的数据类型。直接输入：创建一个名为multivnorm的数据框，含有变量id（数值型）、group（字符型）、weight（数值型）和height（数值型），并进行数据录入。multivnorm <- data.frame(id=numeric(0),group=character(0),weight=numeric(0),height=numeric(0))
multivnorm <- edit(multivnorm)
注：id=numeric(0)表示创建一个名称为id但不含任何实际数据的变量。multivnorm <- edit(multivnorm)的等价简洁写法是fix(multivnorm)。
然后在打开的数据编辑器中进行数据录入即可。
文件导入：在进行数据导入时需要调用程序包foreign，你可以通过菜单：程序包>>加载程序包，选择foreign进行加载，也可以直接使用命令library(foreign)进行加载。我们在上一次的笔记中已经在STATA中对数据进行了录入并分析，我们可以把stata文件直接导入到R中。但要注意目前R3.6.1仅能导入stata12以下的数据，如果你的stata版本过高，需要把数据存为低版本的格式才可以被导入。将示例2的stata文件Multivariate放在R的工作目录D:\Temp下，然后导入到R中，命令如下：
library(foreign)
multivnorm <- read.dta("Multivariate.dta")
如果你想显示multivnorm的内容，可直接在命令行输入multivnorm即可，也可以用edit(multivnorm)对其进行编辑
多元正态分布分析（mshapiro.test {mvnormtest}）：不同的程序包中的函数对应的数据类型可能不同，在采用这些函数时，先要对数据进行变换，使之符合要求。不同的程序包及其函数的功能和用法可通过命令help.start()、help("函数")或？函数进行查找或查看，或者在帮助>>R函数帮助(文本)…、帮助>>搜寻帮助…等处进行查找。比如查询程序包mvnormtest可知其函数mshapiro.test可Performs the Shapiro-Wilk test for multivariate normality，应用格式为mshapiro.test(U)，U要求是数值型矩阵，且a matrix with number of columns (sample size) between 3 and 5000，并提供示例。因为在一般在录入时行表示观测/记录，列表示变量/字段，因此使用此函数需要先对行和列进行转置，转置函数为t("矩阵或数据框")。在进行分析前记得先要载入程序包mvnormtest。本例A、B两组需要分别进行检验，方括号[i,j]中表示第i行第j列，冒号表示范围，如[1:8,3:4]表示第1-8行第3-4列。A<-t(multivnorm[1:8,3:4])
B<-t(multivnorm[9:16,3:4])
library(mvnormtest)
mshapiro.test(A)
mshapiro.test(B)
结果：红色的为执行语句，蓝色为结果。本例显示A组多元正态分布SW检验W=0.912，P=0.366>0.05，呈二元正态分布；B组多元正态分布SW检验W=0.861，P=0.122>0.05，呈二元正态分布【解读参见留言】。
附完整语句命令【命令有重复，参见留言】：
setwd("D:/Temp")
library(foreign)
multivnorm <- read.dta("Multivariate.dta")
library(foreign)
multivnorm <- read.dta("Multivariate.dta")
A<-t(multivnorm[1:8,3:4])
B<-t(multivnorm[9:16,3:4])
library(mvnormtest)
mshapiro.test(A)
示例2的多元正态分布检验的R实现过程：此次采用R的mvn {MVN}。mvn {MVN}多变量正态检验功能强大，用法：mvn(data, subset = NULL, mvnTest = c("mardia", "hz", "royston", "dh", "energy"), covariance = TRUE, tol = 1e-25, alpha = 0.5, scale = FALSE, desc = TRUE, transform = "none", R = 1000, univariateTest = c("SW", "CVM", "Lillie", "SF", "AD"), univariatePlot = "none", multivariatePlot = "none", multivariateOutlierMethod = "none", bc = FALSE, bcType = "rounded", showOutliers = FALSE, showNewData = FALSE)。Data为矩阵或数据框，具体可参照帮助。
语句命令：
setwd("D:/Temp")
library(foreign)
multivnorm <- read.dta("RM.dta")
U<-multivnorm[1:10,2:6]
library(MVN)
mvn(U, mvnTest = c("dh"), multivariatePlot = "qq")
注：本例10行6列，列名为id、W0、W1、W2、W3、W4。分析时仅筛选列，行不进行筛选，U<-multivnorm[1:10,2:6]可以直接U<-multivnorm[2:6]或者U<-multivnorm[-1]都可以。
结果：本例采用Doornik-Hansen了多变量正态性检验，当然也可以在语句中修改命令，换成Marida, Royston, Henze-Zirkler's, E-Statistics等方法。结果显示E=6.79，P=0.75>0.05，数据满足多元正态分布。本例同时给出了多元正态性的QQ图。
使用不同的方法结果可能会不一致。
在不同的程序包中使用相同的方法，具体结果也可能会有差异，原因尚未探寻，但基本上结论是一致的。比如上面的RM数据我们通过mvn {MVN}可得到Doornik-Hansen和Henze-Zirkler's的检验结果分别是E=6.79，P=0.745>0.05和HZ=0.751，P=0.335>0.05。而在使用mvnTest程序包时，函数DH.test{mvnTest}和HZ.test{mvnTest}的结果分别为DH=6.361，P=0.784>0.05和HZ=0.691，P=0.616>0.05，对应语句如下：setwd("D:/Temp")
library(foreign)
multivnorm <- read.dta("RM.dta")
U<-multivnorm[1:10,2:6]
library(mvnTest)
library(mvtnorm)
DH.test(U,qqplot = TRUE)
HZ.test(U,qqplot = FALSE)
在使用STATA分析的结果则为DH=6.361，P=0.784>0.05和HZ=0.181，P=0.670>0.05，STATA的操作步骤参加多元正态分布的检验。
由于R检验多元正态分布的程序包很多，函数也很多，在此就不一一操作了，步骤其实很简单：导入数据，加载具有相应功能的程序包，按相应函数对分析数据的要求进行数据的处理，然后进行分析并对结果进行专业解读即可。END