书名：R语言医学多元统计分析
ISBN：978-7-115-62011-8
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
编 者 赵 军 戴静毅
责任编辑 吴晋瑜
人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
网址 http://www.ptpress.com.cn
读者服务热线：(010)81055410
反盗版热线：(010)81055315
随着医学研究和信息技术的快速发展，多元数据分析方法广泛应用于医学各个领域。R是一款优秀的开源软件，有着强大的统计计算与可视化功能。本书使用R语言，结合精选的医学实例介绍常用多元统计分析方法。
统计分析方法只有在实际应用中才能得到最直接、最生动的验证。本书强调实战和应用，尽量淡化统计公式的推导和计算过程。通过本书的学习，读者能熟练使用R语言及相关包实现多元统计计算，还能更深入地理解多元数据分析方法。
本书可作为医学院校高年级本科生或研究生的多元统计分析课程教材，亦可作为其他专业读者和科研工作者从事科研活动的参考资料。全书附有代码和数据集，每章后都有习题，书后附有习题参考答案，可供读者自学使用。
生物医学现象变化万端，因果关系错综复杂。多变量、大样本已经成为医学数据的常态，多元数据分析方法在医学研究和实践中的应用越来越广泛。复杂的多元统计分析计算离不开软件。R是一款免费的开源软件，具有强大的统计计算与可视化功能。本书使用R语言，结合精选的医学实例介绍常用的多元统计分析方法。
统计分析方法只有在实际的应用中才能得到最直接、最生动的验证。本着“让非统计专业读者易理解”的原则，本书强调实战和应用，着重介绍多元统计分析的思路和方法、R语言实现和结果解释，尽量淡化统计公式的推导和计算过程。全书共分11章，分别为绪论、多元数据可视化、多元数据的组间比较、聚类分析、判别分析、主成分分析、因子分析、结构方程模型、典型相关分析、偏最小二乘回归分析、对应分析，基本涵盖了医学研究中常用的多元统计方法。为方便读者学习，书末附有多元统计分析中用到的矩阵运算的R语言实现。本书在介绍经典的统计方法的同时，注意吸收与医学科研实践密切相关的前沿方法以及相关R包的使用。
本书假定读者有一定的统计学基础，了解R语言的基本用法。书中配有大量的案例解析和程序示例，以及使用R绘制的图形，所有代码均在R 4.1.2环境下运行通过。书中每一章都配有习题，书末附有习题参考答案。书中示例和习题的数据集和源程序文件可以从异步社区（https://www.epubit.com）下载。书中所有R语言的函数均会带上小括号，以便同普通文本区分开来。书中代码和输出部分以浅灰色背景呈现，采用Courier New字体。除了安装R后自带的核心包，本书还用到了其他一些R包，这些包都可以从R的综合网站CRAN自由获取。在R控制台输入下面的命令可以一次性安装这些包（按照在书中出现的顺序）：
> install.packages(c("StatMatch", "philentropy", "vcd", "MVN", "car",
+ "ggplot2", "corrplot", "GGally", "aplpack", "fmsb", "mclust",
+ "ggpubr", "MSG" ,"ICSNP", "tidyr", "profileR", "biotools",
+ "mclust", "ggpubr", "NbClust", "cluster", "class", "klaR", "rpart",
+ "randomForest", "caret", "FactoMineR", "factoextra", "pls", "lavaan",
+ "psych", "semPlot", "CCA", "CCP", "gplots"))本书适合临床医学、公共卫生及其他医学相关专业的高年级本科生或研究生使用，亦可作为其他专业读者和科研工作者进行数据分析的参考书。读者可以从头至尾逐章学习，也可以根据自己遇到的实际问题有选择地在相应章节找到解决方案。
本书参阅了许多国内外教材和资料，并引用了部分示例数据，在此向相关作者表示衷心的感谢。此外，特别感谢人民邮电出版社的王峰松编辑和吴晋瑜编辑在本书出版过程中给予的支持和协助。
由于作者水平有限，书中难免有不妥和疏漏之处，欢迎读者提出批评、意见和建议，我的电子邮箱地址是zhaojun@hbmu.edu.cn。在医学大数据时代，让我们抓住机遇，共同努力与进步！
赵 军
于湖北十堰
本书由异步社区出品，社区（http://www.epubit.com）为您提供相关资源和后续服务。
本书提供以下资源：
● 书中示例和习题解答的源程序文件；
● 书中彩图文件。
要获得以上资源，您可以扫描下方二维码，根据指引领取。
注意：为保证购书读者的权益，该操作会给出相关提示，要求输入提取码进行验证。
作者和编辑尽最大努力来确保书中内容的准确性，但难免会存在疏漏。欢迎您将发现的问题反馈给我们，帮助我们提升图书的质量。
当您发现错误时，请登录异步社区（https://www.epubit.com/），按书名搜索，进入本书页面，单击“发表勘误”，输入勘误信息，单击“提交勘误”按钮即可（见下图）。本书的作者和编辑会对您提交的勘误进行审核，确认并接受后，将赠予您异步社区的100积分（积分可用于在异步社区兑换优惠券、样书或奖品）。
我们的联系邮箱是contact@epubit.com.cn。
如果您对本书有任何疑问或建议，请您发邮件给我们，并请在邮件标题中注明本书书名，以便我们更高效地做出反馈。
如果您有兴趣出版图书、录制教学视频，或者参与图书翻译、技术审校等工作，可以发邮件给本书的责任编辑（sunzhesi@ptpress.com.cn）。
如果您所在的学校、培训机构或企业，想批量购买本书或异步社区出版的其他图书，也可以发邮件给我们。
如果您在网上发现有针对异步社区出品图书的各种形式的盗版行为，包括对图书全部或部分内容的非授权传播，请您将怀疑有侵权行为的链接发邮件给我们。您的这一举动是对作者权益的保护，也是我们持续为您提供有价值的内容的动力之源。
“异步社区”（www.epubit.com）是由人民邮电出版社创办的IT专业图书社区，于2015年 8 月上线运营，致力于优质内容的出版和分享，为读者提供高品质的学习内容，为作译者提供专业的出版服务，实现作者与读者在线交流互动，以及传统出版与数字出版的融合发展。
“异步图书”是异步社区策划出版的精品IT图书的品牌，依托于人民邮电出版社在计算机图书领域30余年的发展与积淀。异步图书面向IT行业以及各行业使用IT技术的用户。
多元统计分析（multivariate statistical analysis），简称多元分析，是研究客观事物中多个指标（或多个因素）之间相互依赖的统计规律性的一门学科，是数理统计学的一个重要分支。多元分析能够在多个对象和多个指标互相关联的情况下分析它们的统计规律，广泛运用于医学研究。
医学多元统计分析中的“多元”一般是指研究的结局指标（因变量）有多个。在大多数医学研究中，对每个研究对象的观测结果往往不止一个，需要用很多个反应变量来表示。例如，血脂记录有胆固醇、甘油三酯、磷脂、非酯化脂肪酸等指标；研究儿童的生长发育通常需要测量身高、体重、胸围、肺活量等指标。这种有多个反应变量的数据称为多变量数据，也称多元数据（multivariate data）。多元数据通常可以表示成表1-1中的表格形式。
表1-1 多元数据的表格形式
| 样品编号（i） |
观测指标（j） |
|||
|---|---|---|---|---|
| X1 |
X2 |
… |
Xm |
|
| 1 |
|
|
… |
|
| 2 |
|
|
… |
|
|
|
|
|
|
|
| n |
|
|
… |
|
将表1-1中的数据用一个n行m列的矩阵来表达，就是一个多元数据矩阵。
这个矩阵可简写成
。
例如，表1-2是36例肝硬化患者的部分资料，包括性别、年龄组和4项临床检测指标。其中，临床检测指标包括纤维蛋白原（FIB，g/dL）、凝血酶原时间（PT，s）、凝血酶原活动度（PTA，%）、血清胆碱酯酶（CHE，U/L）。为使变量服从正态分布，对PT和CHE取了自然对数。
表1-2 36例肝硬化患者的部分资料
| 编号 |
性别 |
年龄组 |
FIB |
lnPT |
PTA |
lnCHE |
|---|---|---|---|---|---|---|
| 1 |
male |
<40 |
2.80 |
2.53 |
110 |
8.76 |
| 2 |
male |
40—59 |
3.02 |
2.57 |
103 |
8.37 |
| 3 |
female |
40—59 |
2.45 |
2.57 |
101 |
8.02 |
| 4 |
male |
40—59 |
2.59 |
2.57 |
81 |
7.43 |
| 5 |
male |
40—59 |
3.52 |
2.66 |
85 |
8.29 |
| 6 |
male |
40—59 |
2.50 |
2.51 |
80 |
8.52 |
| 7 |
female |
40—59 |
2.49 |
2.62 |
89 |
8.23 |
| 8 |
male |
<40 |
3.39 |
2.56 |
104 |
8.75 |
| 9 |
male |
40—59 |
2.35 |
3.03 |
44 |
7.53 |
| 10 |
male |
40—59 |
3.00 |
2.61 |
92 |
8.13 |
| 11 |
female |
<40 |
2.37 |
2.42 |
107 |
8.71 |
| 12 |
female |
60+ |
2.25 |
2.74 |
73 |
8.41 |
| 13 |
male |
40—59 |
2.22 |
2.53 |
78 |
8.15 |
| 14 |
female |
40—59 |
1.83 |
2.44 |
89 |
8.43 |
| 15 |
male |
<40 |
2.17 |
2.59 |
70 |
7.50 |
| 16 |
male |
40—59 |
0.78 |
3.25 |
32 |
7.55 |
| 17 |
male |
40—59 |
1.57 |
3.06 |
42 |
7.17 |
| 18 |
male |
40—59 |
2.22 |
2.66 |
85 |
8.52 |
| 19 |
female |
40—59 |
1.85 |
2.81 |
63 |
8.09 |
| 20 |
male |
40—59 |
1.24 |
2.97 |
49 |
7.05 |
| 21 |
male |
40—59 |
2.08 |
2.70 |
63 |
7.81 |
| 22 |
male |
40—59 |
2.38 |
2.71 |
75 |
8.68 |
| 23 |
female |
60+ |
3.33 |
2.62 |
94 |
8.72 |
| 24 |
male |
40—59 |
1.44 |
2.75 |
68 |
7.48 |
| 25 |
female |
40—59 |
2.84 |
2.59 |
97 |
8.78 |
| 26 |
male |
60+ |
2.48 |
2.73 |
74 |
7.72 |
| 27 |
male |
40—59 |
2.83 |
2.60 |
91 |
8.59 |
| 28 |
male |
40—59 |
2.34 |
2.79 |
67 |
8.68 |
| 29 |
female |
60+ |
2.40 |
2.39 |
113 |
8.79 |
| 30 |
male |
60+ |
3.91 |
2.74 |
72 |
7.83 |
| 31 |
male |
60+ |
2.32 |
2.76 |
71 |
8.62 |
| 32 |
female |
40—59 |
3.60 |
2.53 |
114 |
7.82 |
| 33 |
female |
<40 |
0.96 |
3.08 |
36 |
7.47 |
| 34 |
male |
60+ |
2.59 |
2.70 |
78 |
7.53 |
| 35 |
female |
60+ |
1.71 |
2.92 |
54 |
7.99 |
| 36 |
female |
40—59 |
5.10 |
2.30 |
128 |
8.97 |
读入表1-2中的数据并查看变量的类型：
> cirr <- read.csv('cirrhosis.csv')
> str(cirr)
'data.frame': 36 obs. of 6 variables:
$ sex : chr "male" "male" "female" "male" ...
$ agegrp: chr "<40" "40-59" "40-59" "40-59" ...
$ FIB : num 2.8 3.02 2.45 2.59 3.52 2.5 2.49 3.39 2.35 3 ...
$ lnPT : num 2.53 2.57 2.57 2.57 2.66 2.51 2.62 2.56 3.03 2.61 ...
$ PTA : int 110 103 101 81 85 80 89 104 44 92 ...
$ lnCHE : num 8.76 8.37 8.02 7.43 8.29 8.52 8.23 8.75 7.53 8.13 ...R4.0.0版本以后，函数read.csv( )的参数stringsAsFactors默认为FALSE，因此上面读入的两个变量sex和agegrp为字符型。下面使用函数factor( )将它们转化为因子（factor）型：
> cirr$sex <- factor(cirr$sex)
> cirr$agegrp <- factor(cirr$agegrp)
> str(cirr)
'data.frame': 36 obs. of 6 variables:
$ sex : Factor w/ 2 levels "female","male": 2 2 1 2 2 2 1 2 2 2 ...
$ agegrp: Factor w/ 3 levels "<40","40-59",..: 1 2 2 2 2 2 2 1 2 2 ...
$ FIB : num 2.8 3.02 2.45 2.59 3.52 2.5 2.49 3.39 2.35 3 ...
$ lnPT : num 2.53 2.57 2.57 2.57 2.66 2.51 2.62 2.56 3.03 2.61 ...
$ PTA : int 110 103 101 81 85 80 89 104 44 92 ...
$ lnCHE : num 8.76 8.37 8.02 7.43 8.29 8.52 8.23 8.75 7.53 8.13 ...对于单个变量，常用的描述性统计量有均值、方差、标准差等。对于多元数据，各个变量之间往往存在相互联系，它们之间的作用也会相互影响。因此，在分析多元数据时，我们还需要考虑各个变量之间的相互关联。多元分析中的描述性统计量主要有均值向量、协方差矩阵、相关系数矩阵等。与一元分析类似，多元分析中的统计量也是从样本计算得到的。
样本的均值向量（means vector）处于样本数据的“中心”，由各个指标的均值组成。例如，使用函数colMeans()计算表1-2中4项检测指标的样本均值向量：
> bio <- cirr[, 3:6]
> colMeans(bio)
FIB lnPT PTA lnCHE
2.470000 2.678056 79.777778 8.133056在一元分析中，用方差描述变量的离散程度；而在多元分析中，除了计算变量自身的方差，还需计算变量之间的协方差。两个变量
和
的样本协方差计算公式为
（1.1）
其中，
为样本量。当
时，
就是的方差。
将各指标的方差、协方差用矩阵的形式表示就得到方差-协方差矩阵，简称协方差矩阵（covariance matrix）。对于包含
个变量的随机向量，其样本协方差矩阵可以表示为
显然，协方差矩阵是一个对称矩阵。
对于表1-2中的数据，4项检测指标的样本协方差矩阵可以用函数var( )计算得到：
> var(bio)
FIB lnPT PTA lnCHE
FIB 0.6951200 -0.10213714 14.321714 0.22538286
lnPT -0.1021371 0.03695325 -4.075016 -0.06602817
PTA 14.3217143 -4.07501587 530.177778 8.43441270
lnCHE 0.2253829 -0.06602817 8.434413 0.27272468相关系数常用于描述两个连续型变量之间的关系，其符号（±）表明相关关系的方向（正相关或负相关），其绝对值的大小反映关系的强弱。两个变量和的样本相关系数计算公式为
（1.2）
其中，
为样本量。相关系数的取值在−1与1之间。
将各个指标之间的相关系数用矩阵的形式表示就得到相关系数矩阵（correlation coefficient matrix）。样本相关系数矩阵通常用R表示：
与协方差矩阵类似，相关系数矩阵也是一个对称矩阵。因为变量自身的相关系数为1，所以R的对角线上的元素均为1。
样本相关系数矩阵可以用函数cor( )计算得到，例如：
> cor(bio)
FIB lnPT PTA lnCHE
FIB 1.0000000 -0.6372759 0.7460267 0.5176411
lnPT -0.6372759 1.0000000 -0.9206450 -0.6577195
PTA 0.7460267 -0.9206450 1.0000000 0.7014260
lnCHE 0.5176411 -0.6577195 0.7014260 1.0000000实际上，如果对每个变量作标准化变换（减去其均值，除以其标准差），那么标准化后的变量的协方差矩阵就等于原变量的相关系数矩阵。标准化可以借助函数scale( )实现：
> var(scale(bio))
FIB lnPT PTA lnCHE
FIB 1.0000000 -0.6372759 0.7460267 0.5176411
lnPT -0.6372759 1.0000000 -0.9206450 -0.6577195
PTA 0.7460267 -0.9206450 1.0000000 0.7014260
lnCHE 0.5176411 -0.6577195 0.7014260 1.0000000在分类问题中，一般会把特征接近的事物归为一类，把特征不同的事物归为不同的类。因此，首先需要建立定量指标来刻画事物之间的相近程度或相似程度。这类指标就是距离、相异系数、相似程度和列联系数。它们是学习多元统计分析方法的基础。
一个研究对象常常需要用多个变量来刻画其特征。如果n个对象需要用p个数值型变量描述，那么可以把这n个对象看成p维空间中的n个点。很自然地，两个对象之间的相似程度可用p维空间中两点之间的距离来度量。令dij表示对象Xi与Xj的距离，常用的距离定义有以下几种。
（1）欧氏距离
（1.3）
（2）绝对值距离（曼哈顿距离）
（1.4）
（3）切比雪夫距离（棋盘距离）
（1.5）
（4）闵氏距离
（1.6）
（5）兰氏距离（堪培拉距离）
（1.7）
其中，欧氏距离是人们较为熟悉的，也是使用最多的距离。stats包里的函数dist( )可以用于计算上述各种距离，函数中的参数method用于设定计算距离的方法。欧氏距离、绝对值距离、切比雪夫距离、闵氏距离、兰氏距离的参数method的取值分别为euclidean、manhattan、maximum、minkowski、canberra，默认为“euclidean”，即欧氏距离。例如，使用表1-2中的4项临床检测指标计算前5个研究对象的距离：
> dist(bio[1:5,])
1 2 3 4
2 7.014421
3 9.037240 2.108886
4 29.031269 22.024271 20.009190
5 25.015119 18.007346 16.038264 4.196737上面的命令选取了数据框bio的前5行，用函数dist( )计算了5个对象两两之间的欧氏距离。
距离矩阵是一个对称矩阵，函数dist( )的默认输出只显示距离矩阵的下三角，我们可以通过设置参数diag和upper为TRUE显示完整的距离矩阵：
> dist(bio[1:5,], diag = TRUE, upper = TRUE)
1 2 3 4 5
1 0.000000 7.014421 9.037240 29.031269 25.015119
2 7.014421 0.000000 2.108886 22.024271 18.007346
3 9.037240 2.108886 0.000000 20.009190 16.038264
4 29.031269 22.024271 20.009190 0.000000 4.196737
5 25.015119 18.007346 16.038264 4.196737 0.000000需要注意的是，在使用欧氏距离时，变量的量纲不能相差太大。当变量的量纲不同，测量值的变异相差悬殊时，需要先对数据进行标准化处理，然后用标准化后的数据计算距离。
由1.2.1节中算得的均值向量可以看出，各指标之间由于测量单位的不同导致数值的差异较大。如果直接计算行与行之间的距离，会出现“大数吃小数”现象，即PTA所在列的权重远大于其他各列。因此，这里需要用函数scale( )将数据标准化后再计算距离矩阵：
> bio.scale <- scale(bio)
> dist(bio.scale[1:5,])
1 2 3 4
2 0.8735292
3 1.5427868 0.9613100
4 2.8599002 2.1020932 1.4349353
5 1.7865971 1.1015574 1.6175318 2.0507346由于两个变量的量纲不同，数据的变异也相差较大，上面得到的两个距离矩阵有很大的不同。此外，在计算距离矩阵时，还应尽可能地避免变量间的多重相关。多重相关所造成的信息重叠，会片面强调某些变量的重要性。
鉴于上述距离的不足，一种改进的距离就是马氏距离：
（1.8）
其中，
为总体的协方差矩阵，实际中常用样本协方差矩阵估计。马氏距离不受量纲的影响。当变量之间彼此完全不相关时，为单位矩阵，此时马氏距离就是欧氏距离。
在R中，stats包的函数mahalanobis( )可用于计算样本点到某个中心点的马氏距离。需要注意的是，该函数的返回值是上面定义的马氏距离的平方。
StatMatch包的函数mahalanobis.dist( )可以用来计算数据框中各样品之间的马氏距离。使用该包之前需要先安装。
> library(StatMatch)
> mahalanobis.dist(bio)[1:5, 1:5]
1 2 3 4 5
1 0.000000 1.146037 1.717457 3.972519 3.052284
2 1.146037 0.000000 1.274836 3.059483 2.248249
3 1.717457 1.274836 0.000000 2.884375 3.273679
4 3.972519 3.059483 2.884375 0.000000 2.992607
5 3.052284 2.248249 3.273679 2.992607 0.000000上面计算对象之间的马氏距离，其中的总体协方差矩阵用样本协方差矩阵估计。用户也可以使用函数中的参数vc自定义协方差矩阵。
如果研究对象的特征是用分类变量来刻画的，我们可以用相异系数（dissimilarity coefficient）度量对象之间的相异性。两个对象i和j之间的相异系数根据对象之间的不匹配率来衡量，计算公式为
（1.9）
其中，k是匹配的数目（即对象i和j取值相同的属性数），p是刻画对象的属性总数。
将研究对象彼此之间的相异系数用矩阵的形式表示就得到相异系数矩阵（dissimilarity coefficient matrix）。StatMatch包的函数gower.dist( )可以用来计算对象之间的相异系数矩阵。例如，在数据集cirr中，如果只用性别描述对象的特征（此时p = 1，k为0或者1），那么对象之间的相异系数矩阵为：
> g1 <- gower.dist(cirr[, 1])与距离矩阵类似，相异系数矩阵也是一个对角线上全为0的方阵。
> dim(g1)
[1] 36 36
> diag(g1)
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0由于该矩阵较大，下面仅显示前5个对象的相异系数矩阵：
> g1[1:5, 1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 1 0 0
[2,] 0 0 1 0 0
[3,] 1 1 0 1 1
[4,] 0 0 1 0 0
[5,] 0 0 1 0 0可以看出，性别相同的对象之间的相异系数为0，而性别不同的对象之间的相异系数为1。
进一步地，用两个分类变量（性别、年龄组）描述对象的特征（此时p = 2），相异系数矩阵为：
> g2 <- gower.dist(cirr[, 1:2])
> g2[1:5, 1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.0 0.5 1.0 0.5 0.5
[2,] 0.5 0.0 0.5 0.0 0.0
[3,] 1.0 0.5 0.0 0.5 0.5
[4,] 0.5 0.0 0.5 0.0 0.0
[5,] 0.5 0.0 0.5 0.0 0.0从上面的相异系数矩阵可以看出，性别和年龄组都相同的对象之间的相异系数为0（如对象2与对象4），性别和年龄组有其中一个相同的对象之间的相异系数为0.5（如对象1与对象2），而性别和年龄组均不同的对象之间的相异系数为1（如对象1与对象3）。
实际上，函数 gower.dist( )不仅可以用于计算基于分类变量的相异系数，还可以依据Gower相异系数的定义（Gower，1971）扩展到混合类型（包括逻辑型、因子型、字符型、有序因子型、数值型）的数据。Gower相异系数也称Gower距离，它依据不同类型的变量给出了不同的计算方法，并将计算结果映射到共同的值域区间[0, 1]上，值越大表示相异程度越大。例如：
> g3 <- gower.dist(cirr)
> g3[1:5, 1:5]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.0000000 0.2315684 0.4387740 0.3554576 0.3086653
[2,] 0.2315684 0.0000000 0.2241809 0.1408645 0.0770969
[3,] 0.4387740 0.2241809 0.0000000 0.2608194 0.2796814
[4,] 0.3554576 0.1408645 0.2608194 0.0000000 0.1408094
[5,] 0.3086653 0.0770969 0.2796814 0.1408094 0.0000000从上面的输出结果可以看到，在前5个研究对象中，第1个和第3个对象之间的差异最大，而第2个和第5个对象之间的差异最小。
在相异系数矩阵g3中，最大的值为：
> max(g3)
[1] 0.7910371除去对角线上的0，最小的值为：
> min(g3[g3 != 0])
[1] 0.0349978如果想找到整个数据集cirr中差异最大和最小的对象，可以使用下面的命令：
> which(g3 == max(g3), arr.ind = TRUE)
row col
[1,] 29 16
[2,] 16 29
> cirr[c(16, 29), ]
sex agegrp FIB lnPT PTA lnCHE
16 male 40-59 0.78 3.05 32 7.55
29 female 60+ 2.40 2.39 113 8.79
> which(g3 == min(g3[g3 != 0]), arr.ind = TRUE)
row col
[1,] 34 26
[2,] 26 34
> cirr[c(26, 34), ]
sex agegrp FIB lnPT PTA lnCHE
26 male 60+ 2.48 2.73 74 7.72
34 male 60+ 2.59 2.70 78 7.53cluster包的函数daisy( )也可以用来计算Gower相异系数，在聚类分析中经常用到。
样品间的亲疏关系用距离度量，而变量间的相似程度用相似系数度量。常用的相似系数有相关系数和夹角余弦。
Pearson相关系数是最常用的一种相关系数，其计算公式见式（1.2）。函数cor( )可用于计算Pearson样本相关系数，例如：
> cor(bio$FIB, bio$PTA)
[1] 0.7460267需要注意的是，计算Pearson相关系数时，要求变量服从二元正态分布。函数cor( )中的参数method默认为“pearson”，即默认计算Pearson相关系数。如果变量不服从正态分布，可以将method设为“spearman”计算Spearman秩相关系数。
夹角余弦使用向量空间中两个向量的夹角的余弦值来衡量它们之间的相似程度。两变量与的夹角余弦定义为
（1.10）
R的基本包中没有计算夹角余弦的函数，不过我们很容易根据上面的公式计算。例如，表1-2中FIB与PTA两个变量的夹角余弦为：
> x <- bio$FIB
> y <- bio$PTA
> sum(x * y) / sqrt(sum(x ^ 2) * sum(y ^ 2))
[1] 0.9770776实际上，
就是向量
与
的夹角余弦。若将原始数据标准化，则夹角余弦与Pearson样本相关系数等价。
> x.scale <- scale(x)
> y.scale <- scale(y)
> sum(x.scale * y.scale) / sqrt(sum(x.scale ^ 2) * sum(y.scale ^ 2))
[1] 0.7460267Philentropy 包实现了 46 种不同的距离算法、相似性度量以及信息熵的度量。它为分类问题、信息理论和机器学习算法等提供了核心的计算框架。使用该包之前需要先安装。
> library(philentropy)使用函数getDistMethods( )查看可以使用的距离算法：
> getDistMethods()
[1] "euclidean" "manhattan" "minkowski"
[4] "chebyshev" "sorensen" "gower"
[7] "soergel" "kulczynski_d" "canberra"
[10] "lorentzian" "intersection" "non-intersection"
[13] "wavehedges" "czekanowski" "motyka"
[16] "kulczynski_s" "tanimoto" "ruzicka"
[19] "inner_product" "harmonic_mean" "cosine"
[22] "hassebrook" "jaccard" "dice"
[25] "fidelity" "bhattacharyya" "hellinger"
[28] "matusita" "squared_chord" "squared_euclidean"
[31] "pearson" "neyman" "squared_chi"
[34] "prob_symm" "divergence" "clark"
[37] "additive_symm" "kullback-leibler" "jeffreys"
[40] "k_divergence" "topsoe" "jensen-shannon"
[43] "jensen_difference" "taneja" "kumar-johnson"
[46] "avg" 函数distance( )是philentropy包中的核心函数，用于计算各种距离。其中的参数method可以设置为上面46种算法之一。参数method的默认值为“euclidean”，即计算欧氏距离。例如，1.3.1节中计算的欧氏距离也可以用下面的命令实现：
> distance(bio[1:5,])
Metric: 'euclidean'; comparing: 5 vectors.
v1 v2 v3 v4 v5
v1 0.000000 7.014421 9.037240 29.031269 25.015119
v2 7.014421 0.000000 2.108886 22.024271 18.007346
v3 9.037240 2.108886 0.000000 20.009190 16.038264
v4 29.031269 22.024271 20.009190 0.000000 4.196737
v5 25.015119 18.007346 16.038264 4.196737 0.000000再如，本小节前面两个变量的夹角余弦也可以用下面的命令实现：
> distance(rbind(x, y), method = 'cosine')
Metric: 'cosine'; comparing: 2 vectors.
cosine
0.9770776上面的命令用函数rbind( )将两个变量x和y按行排列。这是因为函数distance( )计算的是行与行之间的距离或相关系数。在通常情况下，矩阵或者数据框的行表示记录、列表示变量。因此，在计算变量之间的相关系数时，需要将矩阵或者数据框用函数t( )作行与列的转置。
> distance(t(scale(bio)), method = 'cosine')
Metric: 'cosine'; comparing: 4 vectors.
v1 v2 v3 v4
v1 1.0000000 -0.6372759 0.7460267 0.5176411
v2 -0.6372759 1.0000000 -0.9206450 -0.6577195
v3 0.7460267 -0.9206450 1.0000000 0.7014260
v4 0.5176411 -0.6577195 0.7014260 1.0000000上面的命令将数据框bio标准化后计算了各个变量之间的夹角余弦，此时的夹角余弦等价于Pearson相关系数。
> cor(bio)
FIB lnPT PTA lnCHE
FIB 1.0000000 -0.6372759 0.7460267 0.5176411
lnPT -0.6372759 1.0000000 -0.9206450 -0.6577195
PTA 0.7460267 -0.9206450 1.0000000 0.7014260
lnCHE 0.5176411 -0.6577195 0.7014260 1.0000000对于分类变量，常用列联系数（contingency coefficient）表示列联表资料的关联强度。常用的列联系数有Phi系数、Pearson列联系数和克莱姆V系数。
Phi系数只适用于四格表资料，其计算公式为
（1.11）
Pearson列联系数是Phi系数的校正和推广，可以用于多维列联表资料，其计算公式为
（1.12）
克莱姆V系数是在列联表行列数不同时对Pearson列联系数的修正方法，其计算公式为
（1.13）
式（1.11）、式（1.12）和式（1.13）中，
为列联表的值，n为样本量，r和c分别是列联表的行数和列数。
上面3种列联系数的取值范围都在0和1之间。值越接近0，行变量和列变量关系越不密切；值越接近1，行变量和列变量关系越密切。vcd包里的函数assocstats( )可以用来计算上述3种列联系数，例如：
> mytable <- table(cirr$sex, cirr$agegrp)
> mytable
<40 40-59 60+
female 2 7 4
male 3 16 4
> library(vcd)
> assocstats(mytable)
X^2 df P(> X^2)
Likelihood Ratio 1.0043 2 0.60522
Pearson 1.0229 2 0.59963
Phi-Coefficient : NA
Contingency Coeff.: 0.166
Cramer's V : 0.169结果显示，值为1.0229，Pearson列联系数和克莱姆V系数分别为0.166和0.169。这表明在该数据集里，性别和年龄组之间的关联比较弱。由于上面的列联表是2行3列的，Phi系数在这里不适用，所以给出的结果是NA。
正态分布在统计学中有举足轻重的地位。一方面是因为它可以用来描述很多自然现象，另一方面它又是许多其他理论分布的渐近或者极限。多元正态分布可以看作一元正态分布的推广。很多多元统计分析的理论都是建立在多元正态总体的基础上。
在概率论中，一元正态分布的密度函数为
（1.14）
为了便于推广，我们将上式改写为
（1.15）
其中，
表示
的转置。这里，因为
和
都是一维的数，所以与相同。
可以看作以标准差为测量单位的到的距离的平方。
将式（1.15）中的一元距离推广到多元距离
，并对前面的系数进行相应的变换，就可以得到多元正态密度函数。
设p元随机向量X = 
的概率密度函数为
（1.16）
则称X服从p元正态分布，记作
，也称X为p元正态变量。其中，
为协方差矩阵
的行列式。
MASS包的mvrnorm( )函数可以生成服从多元正态分布的随机数向量。例如，下面的命令生成了一个均值向量为（2, 3, 5）、协方差矩阵为3阶单位矩阵、长度为50的服从多元正态分布的随机数向量。为了让结果具有重复性，我们设定了生成随机数的种子。
> mu <- c(2, 3, 5) # 设定均值向量
> sigma <- diag(3) # 设定协方差矩阵
> n <- 50 # 设定样本量
> set.seed(1234) # 设定随机数种子
> library(MASS)
> m <- mvrnorm(n, mu, sigma)
> head(m)
[,1] [,2] [,3]
[1,] 1.289593 3.253319 4.439524
[2,] 2.256884 2.971453 4.769823
[3,] 1.753308 2.957130 6.558708
[4,] 1.652457 4.368602 5.070508
[5,] 1.048381 2.774229 5.129288
[6,] 1.954972 4.516471 6.715065要检验一个随机向量是否服从多元正态分布，可以借助MVN包里的函数mvn( )。例如，检验1.4.1节生成的随机向量是否服从多元正态分布：
> library(MVN)
> mvn(m)
$multivariateNormality
Test HZ p value MVN
1 Henze-Zirkler 0.5421244 0.7899205 YES
$univariateNormality
Test Variable Statistic p value Normality
1 Anderson-Darling Column1 0.3126 0.5379 YES
2 Anderson-Darling Column2 0.1898 0.8948 YES
3 Anderson-Darling Column3 0.2100 0.8529 YES
$Descriptives
n Mean Std.Dev Median Min Max 25th
1 50 1.746100 0.9893339 1.695130 -0.05324722 4.100109 1.049417
2 50 3.146408 0.9054472 3.152579 0.69083112 5.187333 2.638703
3 50 5.034404 0.9258700 4.927360 3.03338284 7.168956 4.440683
75th Skew Kurtosis
1 2.290086 0.41177454 -0.38082462
2 3.629436 -0.04353318 0.04864663
3 5.698177 0.16269192 -0.52369419函数mvn( )中，参数mvnTest的默认值为“hz”，即对随机向量做Henze-Zirkler多元正态性检验；参数univariateTest默认为“AD”，即进行单变量的Anderson-Darling正态分布拟合优度检验。根据上面的输出结果可以认为，各列变量均服从单变量的正态分布，且 3个变量服从多元正态分布（P = 0.79）。
需要强调的是，单变量的正态性是多变量多元正态性的必要非充分条件。即使每个变量的分布呈正态分布，所有变量的组合也未必呈多元正态分布。例如，对于表1-2中的4项临床检测指标，检测它们的正态性：
> mvn(bio, univariateTest = "SW")
$multivariateNormality
Test HZ p value MVN
1 Henze-Zirkler 1.047387 0.00628909 NO
$univariateNormality
Test Variable Statistic p value Normality
1 Shapiro-Wilk FIB 0.9535 0.1348 YES
2 Shapiro-Wilk lnPT 0.9480 0.0903 YES
3 Shapiro-Wilk PTA 0.9853 0.9041 YES
4 Shapiro-Wilk lnCHE 0.9418 0.0579 YES
$Descriptives
n Mean Std.Dev Median Min Max 25th 75th
FIB 36 2.470000 0.8337386 2.39 0.78 5.10 2.1475 2.8325
lnPT 36 2.678056 0.1922323 2.64 2.30 3.08 2.5675 2.7525
PTA 36 79.777778 23.0255896 79.00 32.00 128.00 67.7500 94.7500
lnCHE 36 8.133056 0.5222305 8.19 7.15 8.97 7.6775 8.5975
Skew Kurtosis
FIB 0.6046513 1.2879970
lnPT 0.5042816 -0.3325948
PTA -0.1565475 -0.5762517
lnCHE -0.2425023 -1.2648804在上面的命令中，我们将参数univariateTest设为“SW”，即进行单变量的Shapiro-Wilk正态性检验。检验结果表明，所有P值均大于0.05，即每个变量都服从正态分布。但是，Henze-Zirkler检验的结果表明这4个指标不服从多元正态分布（P = 0.006 < 0.05）。输出结果还给出了4项指标的常用描述性统计量。
在式（1.16）中，如果p = 2，则称X服从二元正态分布。设
和
的均值分别为
和
，方差分别为
和
，和的相关系数为
，则X的概率密度函数表达式为
（1.17）
为了直观地展示二元正态分布，我们用下面的命令绘制了均值向量为（0, 0）、方差都为1、相关系数为0的二元正态分布的三维分布曲面。
> mu1 <- 0; mu2 <- 0 # 设定均值
> s1 <- 1; s2 <- 1 # 设定方差
> rho <- 0 # 设定相关系数
> # 定义密度函数
> dens <- function(x, y){
+ (2 * pi * s1 * s2 * sqrt(1 - rho ^ 2)) ^ -1 *
+ exp(-0.5 * (1 - rho ^ 2) ^ -1 *
+ ((x - mu1) ^ 2 / s1 ^ 2 -
+ 2 * rho * (x - mu1) * (y - mu2) / (s1 * s2) +
+ (y - mu2) ^ 2 / s2 ^ 2))
+ }
> x <- seq(-3, 3, length = 50)
> y <- seq(-3, 3, length = 50)
> z <- outer(x, y, dens) # 计算z = f(x, y), 并且输出给z
> persp(x, y, z, theta = 55, phi = 25) # 绘制二元正态分布曲面绘图结果如图1-1所示。读者可以重新设定不同的参数绘制相应的二元正态分布曲面。
图1-1 二元正态分布曲面（
,
,
）
对于二元正态分布，在某一个变量的取值固定时，另外一个变量服从（单变量的）正态分布。因此，从图1-1的剖面可以看到很多条正态分布曲线。
从一元正态分布很容易得到变量的95%参考值范围。一般来说，由于多元正态分布确定的多元参考值范围考虑了多个指标之间的相关性，因此要比单个指标的参考值范围的简单联合更合理。对于服从二元正态分布的变量，它们的参考值范围在几何上是一个椭圆，称为置信椭圆（confidence ellipse）。car包中的函数dataEllipse( )可以用于绘制置信椭圆，例如：
# install.packages("car")
> library(car)
> ellipse <- dataEllipse(cirr$FIB, cirr$PTA, levels = 0.95,
+ lwd = 1.5, center.cex = 1,
+ xlim = c(0, 5), ylim = c(10, 140))绘图结果如图1-2（a）所示。该置信椭圆表示，根据样本数据，肝硬化患者的FIB和PTA的检测值落在该区域的约占95%。
ggplot2包也可以用于绘制置信椭圆，绘图结果如图1-2（b）所示。
# install.packages("ggplot2")
> library(ggplot2)
> ggplot(cirr, aes(FIB, PTA)) + geom_point(color = 'red') +
+ stat_ellipse(type = 'norm')（a）
（b）
图1-2 二维相关数据FIB和PTA的95%参考值范围
使用car包绘制置信椭圆的好处是可以提取椭圆上各点的坐标。查看对象ellipse的前6行，结果如下。
> head(ellipse)
x y
[1,] 4.060018 138.6389
[2,] 4.222379 138.1928
[3,] 4.358175 136.8611
[4,] 4.465349 134.6641
[5,] 4.542276 131.6351
[6,] 4.587789 127.8200多元统计分析的内容既包括一元统计分析方法的推广，也包括多个随机变量本身特有的性质和应用。基础医学、临床医学和公共卫生实践为多元统计分析的应用提供了广阔的空间。多元统计分析方法在医学上的应用大致可以归纳为下面5个方面。
（1）多元数据的统计描述。
（2）多变量的参数估计与假设检验。主要是以多元正态分布的均值向量和协方差矩阵为代表的参数估计和假设检验。
（3）降维。将相互依赖的多个变量通过变量变换等方式转换成互不相关的变量，或把高维空间的数据投影到低维空间。主成分分析、因子分析、对应分析等就利用了这样的方法。
（4）分类与判别。根据测量指标将一些“相似”的对象和变量按照某些规则归为一类。聚类分析、判别分析等就是用于构造分类的方法。
（5）研究多个变量之间的关系。
1-1 为了解某地9岁男童的生长发育情况，某医师在该地随机选取并测量了20名9岁男童的身高、体重和胸围，数据见表1-3。
（1）试计算样本的均值向量、协方差矩阵和相关系数矩阵。
（2）使用马氏距离计算20名男童生长发育之间的相似度，并找出生长发育差别最大的两个。
表1-3 某地20名9岁男童生长发育数据
| 编号 |
身高/cm |
体重/kg |
胸围/cm |
编号 |
身高/cm |
体重/kg |
胸围/cm |
|
|---|---|---|---|---|---|---|---|---|
| 1 |
136 |
30.5 |
67.5 |
11 |
128 |
28.7 |
58.5 |
|
| 2 |
140 |
31.5 |
69.5 |
12 |
130 |
29.6 |
59.9 |
|
| 3 |
138 |
25.6 |
67.8 |
13 |
128 |
25.4 |
55.7 |
|
| 4 |
136 |
32.5 |
62.5 |
14 |
136 |
26.9 |
64.8 |
|
| 5 |
134 |
33.2 |
63.6 |
15 |
131 |
28.8 |
63.5 |
|
| 6 |
129 |
29.8 |
59.6 |
16 |
146 |
42.5 |
70.5 |
|
| 7 |
142 |
35.6 |
69.8 |
17 |
145 |
41.5 |
71.2 |
|
| 8 |
136 |
34.2 |
64.8 |
18 |
142 |
39.8 |
71.9 |
|
| 9 |
126 |
27.8 |
59.7 |
19 |
139 |
38.9 |
69.8 |
|
| 10 |
140 |
34.6 |
65.4 |
20 |
135 |
33.8 |
60.4 |
1-2 依据表1-3中的数据，绘制二元变量身高和体重的95%置信椭圆。
