数据可视化是将数据以视觉表现展现的科学研究技术，用绘制图表的方式来表示数据及其之间关系的方法，让研究者可准确而快捷地了解数据的趋势、异常值和关联性等信息[1]。人类的眼睛对图案和色彩较为敏感，数据可视化技术是构建人与数据合作的重要桥梁。在大数据时代，数据可视化技术应用广泛，如在分析海量信息后提出决策分析即可表现出明显的优势[2]。多元线性回归分析（multiple linear regression analysis，MLRA），是研究多个自变量与一个因变量间是否存在线性变化关系的研究方法，通过建立回归模型来衡量不同自变量对因变量的影响能力，是医学研究中最常用统计学方法之一[3]。MLRA 首先依据研究的目的确定自变量（Xi）与因变量（Y）的函数关系，然后进行建模、参数估计和检验，最后进行模型的赋值运用[4]。MLRA 的本质是构建具有实际临床应用价值的模型。而确定模型能否构建和使用，在建模前需要满足变量间的线性关系和方差齐性等条件，在建模后还需进行显著性检验和残差分析等，通过数据可视化技术可直接观察模型是否构建合理，有无异常值和异常趋势等，较单纯公式计算的结果更直观。
R 语言软件是由新西兰奥克兰大学 Ross Ihaka 和 Robert Gentleman 开发，R 开发核心团队负责完善，目前众多研究领域广泛使用的一个免费开放源代码统计软件（官方网站：https://www.r-project.org/），常用于统计分析、绘图、数据挖掘等[5]。目前随着计算机技术的快速发展，基于 R 语言集成开发环境的 R Studio 软件，是一组可帮助我们直接便捷高效利用 R 来进行软件开发的专业集成工具，该软件可在各种 UNIX 平台，Windows 和 MacOS 上进行编译和运行，它在传统 R 语言统计基础上更提升了数据可视化功能，使曾经可能难以实现的制图过程如今可做到轻松求解[6]。本文结合具体实例介绍如何使用 R 语言代码，基于 R Studio 实现 MLRA 的数学建模、赋值应用及相关数据可视化，以期该软件有效提升 MLRA 在医学研究工作中的服务能力。
在安装 R Studio 软件（下载网址：http://www.R Studio.com/products/R Studio/）前我们须在电脑上安装 R 语言软件。本文所选用的数据为孙振球，徐勇勇主编的“十二五”规划教材《医学统计学》第十五章多元线性回归分析中的教学案例：27 例糖尿病患者的血糖及有关变量的测量结果[3]（表 1），数据需要提前录入 Excel 软件并保存为扩展名为 csv（逗号分隔）类型的文件，本例我们在 D 盘（亦根据个人电脑硬盘分区进行选择）建立名为 GLU.csv 的 Excel 表，对“列”名称进行修改，采用英文缩写便于代码书写和软件运行，具体为“胆固醇”定义为“CHOL”，“甘油三酯”定义为“TG”，“胰岛素”定义为“Insulin”，“糖化血红蛋白”定义为“HbA1c”，“空腹血糖”定义为“GLU”，并对应每列录入数据。以上步骤实现了数据清洗、整理，以便 R 软件的调用和读取。对以上数据进行分析数据格式化（表 2）。根据研究目的和内容，我们确定多个自变量为Xi（i=1，2，3…i）和它们对应的因变量为Y。对应法则f的解析式Y=f（Xi）表示[7]。
表1
27 例糖尿病患者血糖及可能相关变量
表2
多元回归分析数据格式
1 加载数据到 R Studio 中
在构建回归方程之前，我们首先需要为所以确定变量进行赋值。启动 R Studio 软件后，我们进行数据加载。
输入命令：
> fichier <- "D:/GLU.csv"
> mydata <-read.table（fichier,header=TRUE,sep=","）
选择 D 盘 GLU.csv 为目标文件，读取文件数据，header=TRUE 为逻辑值，以指定不同变量名称均由文件第一行给出，sep=“,”意为每一行的各个值之间都以“,”来进行分隔。
此外还需载入后续工作中用到的多个函数包（packages）。
输入命令：
> library（ggplot2）
> library（reshape）
> library（corrplot）
函数包是由广大第三方开发者根据一定的工作主题，编写一系列函数集成体系化的工具套装，并发布在 CRAN（The Comprehensive R Archive Network）工具仓库中。各研究者可根据当前工作需要，自由组合选用不同函数包中的函数。因此，函数包是 R 语言灵活性和活跃生命力的体现，为科研工作带来极大的便利。
正式工作开展前，需要将计划使用到的函数包用 library（）函数加载到工作环境中。本次演示中用到的第三方工具库包括：ggplot2、reshape 和 corrplot。
2 确认数据是否满足多元回归分析条件
2.1 全域关系
我们通过命令求解所有变量之间的关系，并绘制全域散点图观察，确认数据是否满足多元回归分析条件。
输入命令：
> attach（mydata）
> pairs（mydata［,2:6］）
> cor（mydata［,2:6］）
输出结果：见图 1。
图1
全域散点图
这个数据集在第一步中被命名为“mydata”，我们将该数据集通过连接函数 attach（）连接到搜索路径中，后续使用该数据时，不必引用数据集名称，可只引用每一列的列名，即字段名。pairs（）函数对每一对变量生成全域散点图，这类图形可将变量之间的关系形象化，再连接到统计函数 cor（）求解出上图对应的代表“斜率”的矩阵，为下文的高级运算提供重要基础数据。
2.2 线性拟合
观察上图似乎难以发现变量之间存在明显的线性关系，这就需要我们将散点进一步拟合寻找自变量与因变量是否可能存在线性关系，以确认观察值的独立性、因变量的正态性、变量间的线性关系是否满足进一步构建多元回归分析模型的条件。
输入命令：
> add.cor <- function（x,y）{
usr <- par（"usr"）; on.exit（par（usr））
par（usr = c（0, 1, 0, 1））
text（0.5,0.5,round（cor（x,y）,2））}
> add.abline <- function（x, y）{
points（x, y）
abline（lm（y～x）, col = 'red'）
}
> pairs（mydata［,2:6］,lower.panel = add.abline,upper.panel=add.cor）
输出结果：见图 2。
图2
全域散点线性拟合图
执行上述命令后生成全域散点线性拟合图，观察图中最下一行除 Insulin 外其余三项均与 GLU 呈现正相关，且散点具有一定的线性回归趋势，故本例数据可建立线性模型。
3 建立多元线性回归模型
在实例中我们设X1为胆固醇值（变量 CHOL），X2为甘油三酯值（变量 TG），X3为胰岛素值（变量 Insulin），X4为糖化血红蛋白百分比（变量 HbA1c），则解析式Y=f（Xi）可具体为多元线性回归方程[3]：
=b0+b1X1+b2X2+b3X3+b4X4，亦可改写为GLU（Y/X1，X2，X3，X4）=b0+b1CHOL+b2TG+b3Insulin+b4HbA1c，b0，b1，b2，b3和 b4为下文模型中偏回归系数β0，β1，β2，β3和β4参数的估计值；
表示在X1，X2，X3，X4取值时Y的平均值。
3.1 参数估计
公式中的回归参数b0，b1，b2，b3，b4都是我们不知道的，所以参数估计就是要通过自变量的不同数据来估计这些参数从而确定X和Y之间的关系。参数估计的核心目标是要尽可能计算出一条直线，并使这条直线上的每个点的Y值和实际数据Y值之差的平方和最小，即最小二乘法（least sum of squares，LS）：令（Y1实际-Y1预测）2+（Y2实际-Y2预测）2+……+（Yn实际-Yn预测）2值最小。本例中我们可通过 R 软件的 lm（）函数来求得模型中参数β0，β1，β2，β3，β4进行估计，从而得到多元线性回归方程b0，b1，b2，b3，b4参数值。
输入命令：
> model1 <- lm（GLU～CHOL+TG+Insulin+HbA1c）
> model1
输出结果：
R 软件的 lm（Y～X）函数为线性模型的估计命令，model1 为模型构建命令，相应数据为模型参数的估计值，具体为b0=5.9433，b1=0.1424，b2=0.3515，b3=-0.2706，b4=0.6382，带入数据后拟合的 MLRA 方程为：
=5.9433+0.1424X1+0.3515X2-0.2706X3+0.6382X4，如果我们想进一步得出图示上述模型中 27 位患者X1-4对Y的函数映射关系，则执行如下操作：
输入命令：
> datamel <- melt（mydata,id.vars = c（"ID"））
> ggplot（datamel, aes（x = ID））+
geom_line（aes（y = value, col = variable）,size = 1）
输出结果：见图 3。
图3
模型函数 X1-4 与 Y 映射关系
3.2 多重共线性的识别与处理
如一些自变量之间存在较强线性关系时，我们称它们之间存在多重共线性（multicollinearity），其共线性程度越强越容易引起建模的不良后果，如参数估计值的标准误过大、t检验结果过小、回归参数正负符号与客观实际不一致等。因此我们首先要进行共线性的识别，通过矩阵运算对比自变量之间的相关系数，如果其绝对值越接近 1 则认为二者共线性越强。
输入命令：
> corr <- cor（mydata［,2:6］）
> corr
输出结果：
输入命令：
> par（mfrow = c（1, 2））
> corrplot（corr,method="number"）+
corrplot（corr,method="circle"）
输出结果：见图 4。
图4
相关系数图
图 4 为矩阵运算的数值，其颜色参照右侧标尺，越接近深红色和深蓝色表明共线性越强，图 4（左）中除左上至右下的对角线表示自身比较外，其余数值均小于 0.65，说明自变量间共线性较弱，不需要对方程进行调整；图 4（右）为数值的可视化呈现使读者一目了然。
3.3 显著性检验
MLRA 的显著性检验是关于回归模型中因变量与自变量间函数关系的假设是否显著的检验，也就是探讨二者之间的线性关系是否密切的检验。同样是需要经过t检验，F检验，和决定系数R2相关检验[3]。在 R 软件中这三种检验的计算可同时采用函数 summary（）来完成。
输入命令：
> summary（model1）
输出结果：
从输出的结果可观察到标准误差系数表（coefficients）、残差的五数概括（residuals）、决定系数 R2（multiple R-squared，adjusted R-squared）和残余标准误差（residual standard error）等，回归方程参数的估计值在 Estimate 中也有列出，对应于H0：βi =0 和H1：βi≠0 采用t检验统计量的输出数值在t value 中列出。Std. Error 为回归系数估计值b（Estimate）的标准误，故t检验的公式与上述结果的对应为t value=Estimate/ Std. Error，依据t界值表[1]得 t0.05/2，22=2.074，t4>|t3|>2.074，P值小于 0.05，表明b3和b4显著，对同一资料不同自变量的t值是可进行相互比较的，t的绝对值越大就说明自变量对因变量Y的回归起的作用越大。本例说明t4对应的糖化血红蛋白值对空腹血糖的作用较其他变量更大，综合结果我们可进一步明确，对空腹血糖影响大小的顺序依次为糖化血红蛋白（HbA1c）、胰岛素（Insulin）、甘油三酯（CHOL）和总胆固醇（TG）；结果中我们还还可得到多元线性回归决定系数R2[3]，即 R2=SS回/SS总=133.7100/222.5510=0.6008，即输出数据中的 multiple R-squared，调整后的 R2为 adjusted R-squared=0.5282，这表明表示空腹血糖含量（Y）变化的 52.82% 可通过预测变量“CHOL”、“TG”、“Insulin”和“HbA1c”之间的线性关系来解释，故提示模型拟合度不良。
3.4 残差分析
残差即观察值与估计值（拟合值）之间的差，其蕴含着预设模型基本假设的重要信息。如果多元回归模型正确，则我们可将残差看作误差的观测值，且应符合模型的假设条件。我们利用残差提供的信息，进一步检测对应模型假设的合理性和其数据可靠性的方法，我们称之为残差分析，它是检查资料是否符合模型条件的有效工具，如果资料与模型的假设偏离较大，则采用最小二乘估计、假设检验及区间估计有可能出现问题。我们可对函数 lm（）结合函数 plot（）方法得到残差图，残差图包括两个子图形，即残差图—拟合值图和残差—正态 QQ 图。
输入命令：
> par（mfrow=c（1,2））
> plot（model1）
输出结果：见图 5。
图5
标准化残差图和残差正态 QQ 图
图 5（左）理想情况为数据点均匀分布在 y=0 两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征；图 5（右）理想情况为数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布；求解运算后输出图形具体情况为图 5（左）中散点分布不理想，红色线呈一个近似的“U”形，提示该模型中可能缺失了一个重要的未知变量；图 5（右）中，除残差最大和最小的 26 与 13，其余基本符合正态分布。
3.5 模型确立
多元线性回归模型的构建形式为：Y=β0+β1X1+β2X2+…+βmXm+e[1]，β0为截距，β1，β2，…βm即回归系数（partial regression coefficient），e即去除m个自变量对Y影响后的残差，因为在模型和函数关系构建之初，残差e是不可观测的，故在构建回归方程时不需要考虑残差。根据上文的显著性检验和残差分析，MLRA 方程的各自变量参数均无修正，则模型即等于上述 MLRA 方程
=5.9433+0.1424X1+0.3515X2−0.2706X3+0.6382X4。
4 根据实际情况运用模型进行预测
多元线性回归模型可应用于统计区间的预测，即采用多个自变量来计算因变量变化的统计区间。我们可把预设的自变量X代入回归方程，通过 R Studio 软件的函数 predict（）得到Y的个体预测值、Y个体值 95% 预测区间和Y总体均数的 95% 置信区间。如果我们仅针对某一患者已有自变量的观测值对Y进行预测，我们可直接将数据代入方程对应的 R 软件命令进行求解，如我们假设有一位患者，其化验单显示总胆固醇 10.50 mmol/L，甘油三酯 9.50 mmol/L，胰岛素 1.82 μU/ml，糖化血红蛋 9.85%，我们就可依据此进行估计空腹血糖的上述区间。
输入命令：
> newdata <- data.frame（CHOL=10.50,TG=9.50,Insulin=1.82,HbA1c=9.85）
> predict（model1,newdata）
> predict（model1,newdata,interval = "pred",level =0.95）
> predict（model1,newdata,interval = "conf",level =0.95）
输出结果：
我们输入 predict（model1，newdata）命令可计算出Y的个体预测值为 16.5717，为避免误差，我们可在“predict”方法中指定“interval”的参数，通过“level”参数指定置信水平为 0.95，对预测区间可使用“interval="pred"”，求解得（11.56282，21.58057）；Y总体均数的 95% 置信区间使用“interval="conf"”，求解得（13.79310，19.35029）。
如果我们欲对多位患者的自变量进行Y值预测，则应建立数据库以便 R 软件整体导入进行运算和分析。我们仍以上文 27 位患者观测指标的数据库为例，为避免数据重合和图像重合难以阅读观察，我们首先计算 Y 值的预测值和预测区间。数据导入及求解步骤如下。
输入命令：
> library（ggplot2）
> library（reshape2）
> myprediction<-predict（model1,interval="pred",level =0.95）
> head（myprediction,27）
输出结果：
上述结果是 myprediction 的赋值矩阵，我们需要把它变成数据框，然后和原始数据的Y进行合并，即 mydata 中的 GLU 值。
输入命令：
> myprediction <- data.frame（myprediction）
> myprediction <- cbind（y = mydata$GLU, myprediction）
> myprediction
输出结果：
上述表格是宽表格形式，如果进行可视化制图，我们需要将表格转换成长表格形式才能满足函数的读取要求。
输入命令：
> myprediction_long <- melt（myprediction,id.vars = 'x'）
> myprediction_long$x <- as.numeric（myprediction_long$x）
> ggplot（myprediction_long）+
geom_line（aes（x = x, y = value, color = variable）,size=1）
输出结果：见图 6。
图6
27 例患者空腹血糖（Y）实测值，预测值和预测区间上下限
如果我们将上述代码中的“myprediction”替换为“myconfidence”，“interval="pred"”改为“interval="conf"”，“预测区间上下限”改为“置信区间上下限”，则在进一步计算出置信区间上下限及绘图（图 7）。
图7
27 例患者空腹血糖（Y）实际值，预测值和置信区间上下限
5 讨论
多元线性回归分析（MLRA）的本质是为了估计或预测，如文中在建模后代入新患者的糖化血红蛋白观测值以推算空腹血糖预测值的统计区间，需要应用矩阵运算，且在建模前构建预测回归方程时应选择具有较高 R2值的方程以确保方程在引入自变量计算后的结果具有统计学意义（P<0.05）[3]。
MLRA 在医学领域有着广泛的应用，原因在于多因素混杂分析是医学研究中经常遇到的问题。我们已知的绝大多数疾病都是多种原因共同作用的结果，其预后也往往是由多因素所决定[8]，如本文例子中影响空腹血糖的因素不止上述几种，还可能有年龄、饮食习惯、工作环境和家族史等，在此基础上，我们在这些因素中还要明确不同因素的影响权重和可能存在的逻辑关系。这些因素可能影响着临床试验的各个方面，造成我们需要面临可能在难以保证基线相同的混杂情况下对不同的干预方法进行比较这样的难题，以上这些问题很多都可利用多元线性回归分析来处理。我们控制其混杂因素的一个简单而有效的方法就是将具体观测值通过函数的映射关系代入方程进行分析，但这就决定了 MLRA 计算出的预测区间不能说明实际的效果，仅可用于佐证之需。
目前，多元线性回归方程假设检验应用较为普遍的软件是 SPSS 软件、SAS 软件和 STATA 软件等。SPSS 软件虽然简便，但是受其模块化界面操作的限制，难以完成多种统计，SAS、STATA 软件虽然具有编程等强大功能，但其高昂的租赁年费也常令我们的工作面临困境。R 语言因其逻辑运算功能强大，统计制图多样，编程语言简明科学等优势，始终受到科研工作者的青睐，但 R 语言传统自带的环境操作起来较为不便，R Studio 作为 R 语言的高性能操作平台。它兼容了 R 语言数据处理、计算、制图以及通过全球不同地区 Comprehensive R Archive Network（CRAN）镜像站下载升级数据包以随时不断丰富数学计算函数模块等几乎全部功能外，简化了编程的难度，提升了操作的便捷性，且它还具有更加优良的数据可视化功能，更具有支持纯 R 脚本、脚本文档混排（rmarkdown）、脚本文档混排成书（bookdown）、交互式网络应用（shiny）等优势。其唯一可能存在的障碍即要求我们掌握 R 语言的编写，才能完成所需要的统计分析和图形绘制。故本文结合实例提供读者代码的具体编写方法，希望 R Studio 软件今后能为更多研究者所采用，为医学研究的发展做出贡献。
数据可视化是将数据以视觉表现展现的科学研究技术，用绘制图表的方式来表示数据及其之间关系的方法，让研究者可准确而快捷地了解数据的趋势、异常值和关联性等信息[1]。人类的眼睛对图案和色彩较为敏感，数据可视化技术是构建人与数据合作的重要桥梁。在大数据时代，数据可视化技术应用广泛，如在分析海量信息后提出决策分析即可表现出明显的优势[2]。多元线性回归分析（multiple linear regression analysis，MLRA），是研究多个自变量与一个因变量间是否存在线性变化关系的研究方法，通过建立回归模型来衡量不同自变量对因变量的影响能力，是医学研究中最常用统计学方法之一[3]。MLRA 首先依据研究的目的确定自变量（Xi）与因变量（Y）的函数关系，然后进行建模、参数估计和检验，最后进行模型的赋值运用[4]。MLRA 的本质是构建具有实际临床应用价值的模型。而确定模型能否构建和使用，在建模前需要满足变量间的线性关系和方差齐性等条件，在建模后还需进行显著性检验和残差分析等，通过数据可视化技术可直接观察模型是否构建合理，有无异常值和异常趋势等，较单纯公式计算的结果更直观。
R 语言软件是由新西兰奥克兰大学 Ross Ihaka 和 Robert Gentleman 开发，R 开发核心团队负责完善，目前众多研究领域广泛使用的一个免费开放源代码统计软件（官方网站：https://www.r-project.org/），常用于统计分析、绘图、数据挖掘等[5]。目前随着计算机技术的快速发展，基于 R 语言集成开发环境的 R Studio 软件，是一组可帮助我们直接便捷高效利用 R 来进行软件开发的专业集成工具，该软件可在各种 UNIX 平台，Windows 和 MacOS 上进行编译和运行，它在传统 R 语言统计基础上更提升了数据可视化功能，使曾经可能难以实现的制图过程如今可做到轻松求解[6]。本文结合具体实例介绍如何使用 R 语言代码，基于 R Studio 实现 MLRA 的数学建模、赋值应用及相关数据可视化，以期该软件有效提升 MLRA 在医学研究工作中的服务能力。
在安装 R Studio 软件（下载网址：http://www.R Studio.com/products/R Studio/）前我们须在电脑上安装 R 语言软件。本文所选用的数据为孙振球，徐勇勇主编的“十二五”规划教材《医学统计学》第十五章多元线性回归分析中的教学案例：27 例糖尿病患者的血糖及有关变量的测量结果[3]（表 1），数据需要提前录入 Excel 软件并保存为扩展名为 csv（逗号分隔）类型的文件，本例我们在 D 盘（亦根据个人电脑硬盘分区进行选择）建立名为 GLU.csv 的 Excel 表，对“列”名称进行修改，采用英文缩写便于代码书写和软件运行，具体为“胆固醇”定义为“CHOL”，“甘油三酯”定义为“TG”，“胰岛素”定义为“Insulin”，“糖化血红蛋白”定义为“HbA1c”，“空腹血糖”定义为“GLU”，并对应每列录入数据。以上步骤实现了数据清洗、整理，以便 R 软件的调用和读取。对以上数据进行分析数据格式化（表 2）。根据研究目的和内容，我们确定多个自变量为Xi（i=1，2，3…i）和它们对应的因变量为Y。对应法则f的解析式Y=f（Xi）表示[7]。
表1
27 例糖尿病患者血糖及可能相关变量
表2
多元回归分析数据格式
1 加载数据到 R Studio 中
在构建回归方程之前，我们首先需要为所以确定变量进行赋值。启动 R Studio 软件后，我们进行数据加载。
输入命令：
> fichier <- "D:/GLU.csv"
> mydata <-read.table（fichier,header=TRUE,sep=","）
选择 D 盘 GLU.csv 为目标文件，读取文件数据，header=TRUE 为逻辑值，以指定不同变量名称均由文件第一行给出，sep=“,”意为每一行的各个值之间都以“,”来进行分隔。
此外还需载入后续工作中用到的多个函数包（packages）。
输入命令：
> library（ggplot2）
> library（reshape）
> library（corrplot）
函数包是由广大第三方开发者根据一定的工作主题，编写一系列函数集成体系化的工具套装，并发布在 CRAN（The Comprehensive R Archive Network）工具仓库中。各研究者可根据当前工作需要，自由组合选用不同函数包中的函数。因此，函数包是 R 语言灵活性和活跃生命力的体现，为科研工作带来极大的便利。
正式工作开展前，需要将计划使用到的函数包用 library（）函数加载到工作环境中。本次演示中用到的第三方工具库包括：ggplot2、reshape 和 corrplot。
2 确认数据是否满足多元回归分析条件
2.1 全域关系
我们通过命令求解所有变量之间的关系，并绘制全域散点图观察，确认数据是否满足多元回归分析条件。
输入命令：
> attach（mydata）
> pairs（mydata［,2:6］）
> cor（mydata［,2:6］）
输出结果：见图 1。
图1
全域散点图
这个数据集在第一步中被命名为“mydata”，我们将该数据集通过连接函数 attach（）连接到搜索路径中，后续使用该数据时，不必引用数据集名称，可只引用每一列的列名，即字段名。pairs（）函数对每一对变量生成全域散点图，这类图形可将变量之间的关系形象化，再连接到统计函数 cor（）求解出上图对应的代表“斜率”的矩阵，为下文的高级运算提供重要基础数据。
2.2 线性拟合
观察上图似乎难以发现变量之间存在明显的线性关系，这就需要我们将散点进一步拟合寻找自变量与因变量是否可能存在线性关系，以确认观察值的独立性、因变量的正态性、变量间的线性关系是否满足进一步构建多元回归分析模型的条件。
输入命令：
> add.cor <- function（x,y）{
usr <- par（"usr"）; on.exit（par（usr））
par（usr = c（0, 1, 0, 1））
text（0.5,0.5,round（cor（x,y）,2））}
> add.abline <- function（x, y）{
points（x, y）
abline（lm（y～x）, col = 'red'）
}
> pairs（mydata［,2:6］,lower.panel = add.abline,upper.panel=add.cor）
输出结果：见图 2。
图2
全域散点线性拟合图
执行上述命令后生成全域散点线性拟合图，观察图中最下一行除 Insulin 外其余三项均与 GLU 呈现正相关，且散点具有一定的线性回归趋势，故本例数据可建立线性模型。
3 建立多元线性回归模型
在实例中我们设X1为胆固醇值（变量 CHOL），X2为甘油三酯值（变量 TG），X3为胰岛素值（变量 Insulin），X4为糖化血红蛋白百分比（变量 HbA1c），则解析式Y=f（Xi）可具体为多元线性回归方程[3]：
=b0+b1X1+b2X2+b3X3+b4X4，亦可改写为GLU（Y/X1，X2，X3，X4）=b0+b1CHOL+b2TG+b3Insulin+b4HbA1c，b0，b1，b2，b3和 b4为下文模型中偏回归系数β0，β1，β2，β3和β4参数的估计值；
表示在X1，X2，X3，X4取值时Y的平均值。
3.1 参数估计
公式中的回归参数b0，b1，b2，b3，b4都是我们不知道的，所以参数估计就是要通过自变量的不同数据来估计这些参数从而确定X和Y之间的关系。参数估计的核心目标是要尽可能计算出一条直线，并使这条直线上的每个点的Y值和实际数据Y值之差的平方和最小，即最小二乘法（least sum of squares，LS）：令（Y1实际-Y1预测）2+（Y2实际-Y2预测）2+……+（Yn实际-Yn预测）2值最小。本例中我们可通过 R 软件的 lm（）函数来求得模型中参数β0，β1，β2，β3，β4进行估计，从而得到多元线性回归方程b0，b1，b2，b3，b4参数值。
输入命令：
> model1 <- lm（GLU～CHOL+TG+Insulin+HbA1c）
> model1
输出结果：
R 软件的 lm（Y～X）函数为线性模型的估计命令，model1 为模型构建命令，相应数据为模型参数的估计值，具体为b0=5.9433，b1=0.1424，b2=0.3515，b3=-0.2706，b4=0.6382，带入数据后拟合的 MLRA 方程为：
=5.9433+0.1424X1+0.3515X2-0.2706X3+0.6382X4，如果我们想进一步得出图示上述模型中 27 位患者X1-4对Y的函数映射关系，则执行如下操作：
输入命令：
> datamel <- melt（mydata,id.vars = c（"ID"））
> ggplot（datamel, aes（x = ID））+
geom_line（aes（y = value, col = variable）,size = 1）
输出结果：见图 3。
图3
模型函数 X1-4 与 Y 映射关系
3.2 多重共线性的识别与处理
如一些自变量之间存在较强线性关系时，我们称它们之间存在多重共线性（multicollinearity），其共线性程度越强越容易引起建模的不良后果，如参数估计值的标准误过大、t检验结果过小、回归参数正负符号与客观实际不一致等。因此我们首先要进行共线性的识别，通过矩阵运算对比自变量之间的相关系数，如果其绝对值越接近 1 则认为二者共线性越强。
输入命令：
> corr <- cor（mydata［,2:6］）
> corr
输出结果：
输入命令：
> par（mfrow = c（1, 2））
> corrplot（corr,method="number"）+
corrplot（corr,method="circle"）
输出结果：见图 4。
图4
相关系数图
图 4 为矩阵运算的数值，其颜色参照右侧标尺，越接近深红色和深蓝色表明共线性越强，图 4（左）中除左上至右下的对角线表示自身比较外，其余数值均小于 0.65，说明自变量间共线性较弱，不需要对方程进行调整；图 4（右）为数值的可视化呈现使读者一目了然。
3.3 显著性检验
MLRA 的显著性检验是关于回归模型中因变量与自变量间函数关系的假设是否显著的检验，也就是探讨二者之间的线性关系是否密切的检验。同样是需要经过t检验，F检验，和决定系数R2相关检验[3]。在 R 软件中这三种检验的计算可同时采用函数 summary（）来完成。
输入命令：
> summary（model1）
输出结果：
从输出的结果可观察到标准误差系数表（coefficients）、残差的五数概括（residuals）、决定系数 R2（multiple R-squared，adjusted R-squared）和残余标准误差（residual standard error）等，回归方程参数的估计值在 Estimate 中也有列出，对应于H0：βi =0 和H1：βi≠0 采用t检验统计量的输出数值在t value 中列出。Std. Error 为回归系数估计值b（Estimate）的标准误，故t检验的公式与上述结果的对应为t value=Estimate/ Std. Error，依据t界值表[1]得 t0.05/2，22=2.074，t4>|t3|>2.074，P值小于 0.05，表明b3和b4显著，对同一资料不同自变量的t值是可进行相互比较的，t的绝对值越大就说明自变量对因变量Y的回归起的作用越大。本例说明t4对应的糖化血红蛋白值对空腹血糖的作用较其他变量更大，综合结果我们可进一步明确，对空腹血糖影响大小的顺序依次为糖化血红蛋白（HbA1c）、胰岛素（Insulin）、甘油三酯（CHOL）和总胆固醇（TG）；结果中我们还还可得到多元线性回归决定系数R2[3]，即 R2=SS回/SS总=133.7100/222.5510=0.6008，即输出数据中的 multiple R-squared，调整后的 R2为 adjusted R-squared=0.5282，这表明表示空腹血糖含量（Y）变化的 52.82% 可通过预测变量“CHOL”、“TG”、“Insulin”和“HbA1c”之间的线性关系来解释，故提示模型拟合度不良。
3.4 残差分析
残差即观察值与估计值（拟合值）之间的差，其蕴含着预设模型基本假设的重要信息。如果多元回归模型正确，则我们可将残差看作误差的观测值，且应符合模型的假设条件。我们利用残差提供的信息，进一步检测对应模型假设的合理性和其数据可靠性的方法，我们称之为残差分析，它是检查资料是否符合模型条件的有效工具，如果资料与模型的假设偏离较大，则采用最小二乘估计、假设检验及区间估计有可能出现问题。我们可对函数 lm（）结合函数 plot（）方法得到残差图，残差图包括两个子图形，即残差图—拟合值图和残差—正态 QQ 图。
输入命令：
> par（mfrow=c（1,2））
> plot（model1）
输出结果：见图 5。
图5
标准化残差图和残差正态 QQ 图
图 5（左）理想情况为数据点均匀分布在 y=0 两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征；图 5（右）理想情况为数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布；求解运算后输出图形具体情况为图 5（左）中散点分布不理想，红色线呈一个近似的“U”形，提示该模型中可能缺失了一个重要的未知变量；图 5（右）中，除残差最大和最小的 26 与 13，其余基本符合正态分布。
3.5 模型确立
多元线性回归模型的构建形式为：Y=β0+β1X1+β2X2+…+βmXm+e[1]，β0为截距，β1，β2，…βm即回归系数（partial regression coefficient），e即去除m个自变量对Y影响后的残差，因为在模型和函数关系构建之初，残差e是不可观测的，故在构建回归方程时不需要考虑残差。根据上文的显著性检验和残差分析，MLRA 方程的各自变量参数均无修正，则模型即等于上述 MLRA 方程
=5.9433+0.1424X1+0.3515X2−0.2706X3+0.6382X4。
4 根据实际情况运用模型进行预测
多元线性回归模型可应用于统计区间的预测，即采用多个自变量来计算因变量变化的统计区间。我们可把预设的自变量X代入回归方程，通过 R Studio 软件的函数 predict（）得到Y的个体预测值、Y个体值 95% 预测区间和Y总体均数的 95% 置信区间。如果我们仅针对某一患者已有自变量的观测值对Y进行预测，我们可直接将数据代入方程对应的 R 软件命令进行求解，如我们假设有一位患者，其化验单显示总胆固醇 10.50 mmol/L，甘油三酯 9.50 mmol/L，胰岛素 1.82 μU/ml，糖化血红蛋 9.85%，我们就可依据此进行估计空腹血糖的上述区间。
输入命令：
> newdata <- data.frame（CHOL=10.50,TG=9.50,Insulin=1.82,HbA1c=9.85）
> predict（model1,newdata）
> predict（model1,newdata,interval = "pred",level =0.95）
> predict（model1,newdata,interval = "conf",level =0.95）
输出结果：
我们输入 predict（model1，newdata）命令可计算出Y的个体预测值为 16.5717，为避免误差，我们可在“predict”方法中指定“interval”的参数，通过“level”参数指定置信水平为 0.95，对预测区间可使用“interval="pred"”，求解得（11.56282，21.58057）；Y总体均数的 95% 置信区间使用“interval="conf"”，求解得（13.79310，19.35029）。
如果我们欲对多位患者的自变量进行Y值预测，则应建立数据库以便 R 软件整体导入进行运算和分析。我们仍以上文 27 位患者观测指标的数据库为例，为避免数据重合和图像重合难以阅读观察，我们首先计算 Y 值的预测值和预测区间。数据导入及求解步骤如下。
输入命令：
> library（ggplot2）
> library（reshape2）
> myprediction<-predict（model1,interval="pred",level =0.95）
> head（myprediction,27）
输出结果：
上述结果是 myprediction 的赋值矩阵，我们需要把它变成数据框，然后和原始数据的Y进行合并，即 mydata 中的 GLU 值。
输入命令：
> myprediction <- data.frame（myprediction）
> myprediction <- cbind（y = mydata$GLU, myprediction）
> myprediction
输出结果：
上述表格是宽表格形式，如果进行可视化制图，我们需要将表格转换成长表格形式才能满足函数的读取要求。
输入命令：
> myprediction_long <- melt（myprediction,id.vars = 'x'）
> myprediction_long$x <- as.numeric（myprediction_long$x）
> ggplot（myprediction_long）+
geom_line（aes（x = x, y = value, color = variable）,size=1）
输出结果：见图 6。
图6
27 例患者空腹血糖（Y）实测值，预测值和预测区间上下限
如果我们将上述代码中的“myprediction”替换为“myconfidence”，“interval="pred"”改为“interval="conf"”，“预测区间上下限”改为“置信区间上下限”，则在进一步计算出置信区间上下限及绘图（图 7）。
图7
27 例患者空腹血糖（Y）实际值，预测值和置信区间上下限
5 讨论
多元线性回归分析（MLRA）的本质是为了估计或预测，如文中在建模后代入新患者的糖化血红蛋白观测值以推算空腹血糖预测值的统计区间，需要应用矩阵运算，且在建模前构建预测回归方程时应选择具有较高 R2值的方程以确保方程在引入自变量计算后的结果具有统计学意义（P<0.05）[3]。
MLRA 在医学领域有着广泛的应用，原因在于多因素混杂分析是医学研究中经常遇到的问题。我们已知的绝大多数疾病都是多种原因共同作用的结果，其预后也往往是由多因素所决定[8]，如本文例子中影响空腹血糖的因素不止上述几种，还可能有年龄、饮食习惯、工作环境和家族史等，在此基础上，我们在这些因素中还要明确不同因素的影响权重和可能存在的逻辑关系。这些因素可能影响着临床试验的各个方面，造成我们需要面临可能在难以保证基线相同的混杂情况下对不同的干预方法进行比较这样的难题，以上这些问题很多都可利用多元线性回归分析来处理。我们控制其混杂因素的一个简单而有效的方法就是将具体观测值通过函数的映射关系代入方程进行分析，但这就决定了 MLRA 计算出的预测区间不能说明实际的效果，仅可用于佐证之需。
目前，多元线性回归方程假设检验应用较为普遍的软件是 SPSS 软件、SAS 软件和 STATA 软件等。SPSS 软件虽然简便，但是受其模块化界面操作的限制，难以完成多种统计，SAS、STATA 软件虽然具有编程等强大功能，但其高昂的租赁年费也常令我们的工作面临困境。R 语言因其逻辑运算功能强大，统计制图多样，编程语言简明科学等优势，始终受到科研工作者的青睐，但 R 语言传统自带的环境操作起来较为不便，R Studio 作为 R 语言的高性能操作平台。它兼容了 R 语言数据处理、计算、制图以及通过全球不同地区 Comprehensive R Archive Network（CRAN）镜像站下载升级数据包以随时不断丰富数学计算函数模块等几乎全部功能外，简化了编程的难度，提升了操作的便捷性，且它还具有更加优良的数据可视化功能，更具有支持纯 R 脚本、脚本文档混排（rmarkdown）、脚本文档混排成书（bookdown）、交互式网络应用（shiny）等优势。其唯一可能存在的障碍即要求我们掌握 R 语言的编写，才能完成所需要的统计分析和图形绘制。故本文结合实例提供读者代码的具体编写方法，希望 R Studio 软件今后能为更多研究者所采用，为医学研究的发展做出贡献。
