R语言医学数据分析实践-高级回归分析

多重线性回归分析是一种统计分析方法，用于建立自变量（解释变量）与因变量（响应变量）之间的线性关系模型。在多重线性回归中，可以通过多个自变量对因变量进行预测和解释，考虑各自变量对因变量的影响并对它们的系数进行估计。以下是多重线性回归的详细介绍。

多重线性回归模型的一般形式如下：

Y=β0+β1X1+β2X2+…+βpXp+ε

其中，Y是因变量，X1, X2, … , Xp是自变量，β0, β1, … , βp是对应系数，ε是误差项。系数表示当自变量增加一个单位时，因变量的变化量，即自变量与因变量变化的关系。通过最小化观测值与模型预测值之间的残差平方和，可以得到最佳的系数估计值，从而拟合回归模型。在建立多重线性回归模型后，需要对模型进行诊断，发现问题并进行改进。例如，发现离群值、共线性、异方差等问题时，需要采取相应措施。

在实际应用中，多重线性回归可以用于探索自变量对因变量的影响，预测因变量的取值，识别影响因变量的重要自变量等。通过适当的建模方法和对结果的解释，多重线性回归能有效地帮助我们理解变量之间的关系，进行预测和决策。在R语言中，可以使用lm()函数来拟合多重线性回归模型，并通过summary()函数查看模型的统计结果。执行多重线性回归的示例R代码如下：

#创建数据框
data <- data.frame(
  Age = c(22, 34, 35, 42, 40),        		#年龄
  Weight = c(62, 60, 75, 79, 99),    		#体重，单位：千克
  ExerciseFrequency = c(3, 1, 4, 6, 2), 	#每周运动频率
  BloodPressure = c(132, 112, 143, 159, 122)  	#血压，单位：毫米汞柱
)
#进行多重线性回归分析
model <- lm(BloodPressure ~ Age + Weight + ExerciseFrequency, data = data)
#查看模型
summary(model)

lm函数中的公式“BloodPressure ~ Age + Weight + ExerciseFrequency”表示我们正在拟合一个模型，其中血压（BloodPressure）是因变量，年龄（Age）、体重（Weight）和运动频率（ExerciseFrequency）是自变量；“data = data”指定了包含数据的数据框。summary函数将提供模型的详细摘要，包括每个自变量的回归系数、截距、R平方值、F统计量和p值等。这些统计量可以帮助我们了解模型的拟合情况和各个自变量对因变量的影响是否显著。

代码的运行结果如图5-5所示，该图展示了多重线性回归分析的结果。

图5-5

从图5-5中可以看到，模型的截距是104，年龄、体重和运动频率的系数分别是-0.09485、0.0276和9.57084。p值显示运动频率对血压的影响是显著的（p < 0.05），而年龄和体重的影响不显著。R平方值（Multiple R-squared）为0.997，表示模型解释了99.7%的血压变异。

Logistic（逻辑）回归分析是一种广义的线性回归分析模型，常用于数据挖掘、疾病自动诊断、经济预测等领域。Logistic回归根据给定的自变量数据集来估计事件的发生概率，由于结果是一个概率，因此因变量的范围在0和1之间。例如，它可以用于探讨引发疾病的危险因素，并根据危险因素预测疾病发生的概率等。以胃癌病情分析为例，选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群必定具有不同的体征与生活方式等。因此，因变量就为是否患有胃癌，其值为“是”或“否”，自变量可以有很多，如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。通过Logistic回归分析，可以得到自变量的权重，从而可以大致了解到底哪些因素是胃癌的危险因素，同时可以根据这些危险因素预测一个人患癌症的可能性。

Logistic回归用于预测二分类（如是/非、成功/失败等）或多分类问题，它将线性回归模型的输出映射为一个0~1的概率值，从而预测观测值属于某个类别的概率。训练好的Logistic回归模型可以用于预测新样本属于某一类别的概率。通常将概率大于0.5的归类为正类（1），概率小于或等于0.5的归类为负类（0）。

Logistic回归广泛应用于医学领域，可用于疾病风险预测、疾病预后预测等。在R语言中，Logistic回归分析可以通过调用广义线性回归模型函数glm()来实现，通常通过设置参数family为Binomial（二项式分布）来解决二元分类问题。在建立模型后，可以通过predict()函数进行预测，并使用性能评估指标对模型进行评估。以下是执行Logistic回归的示例R代码：

#创建数据框
data <- data.frame(
  Age = c(25, 30, 35, 40, 45, 50),            	#年龄
  Weight = c(63, 71, 89, 90, 110, 102),        	#体重，单位：千克
  Smokes = factor(c(0, 1, 1, 1, 0, 1)),        	#是否吸烟，0代表否，1代表是
  Hypertension = factor(c(0, 1, 0, 1, 1, 0))   	#是否患有高血压，0代表否，1代表是
)
#进行Logistic回归分析
model <- glm(Hypertension ~ Age + Weight + Smokes, data = data, family = binomial)
#查看模型
summary(model)

glm函数中的公式“Hypertension ~ Age + Weight + Smokes”表示我们正在拟合一个模型，其中是否患有高血压（Hypertension）是因变量，年龄（Age）、体重（Weight）和是否吸烟（Smokes）是自变量；“data = data”指定了包含数据的数据框；“family = binomial”指定了逻辑回归使用的分布族，对于二元分类问题，我们使用二项式分布。summary函数将提供模型的详细摘要，包括每个自变量的回归系数、截距、模型拟合优度指标（如伪R平方值）、似然比检验的p值等。这些统计量可以帮助我们了解模型的拟合情况以及各个自变量对因变量的影响是否显著。

代码的运行结果如图5-6所示，该图展示了Logistic回归的结果。

图5-6

从图5-6中可以看到，年龄、体重和是否吸烟的系数在统计上是不显著的（p＞0.05），这表明年龄、体重和吸烟与患高血压之间不存在显著的关联；伪R平方值则可以帮助我们了解模型对因变量变异的解释程度。

为了说明如何在临床研究中运用回归分析，笔者从REGICOR研究中选取了一部分数据用于回归分析实战。REGICOR是一个针对来自西班牙东北部的参与者进行的横断面研究，包括人口统计学信息（年龄、性别、身高、体重、腰围等）、血脂特征（总胆固醇和胆固醇、甘油三酯等）、问卷调查信息（体格、活动、生活质量等）等。此外，心血管事件和死亡信息来自医院和官方登记处。

本研究旨在探究心血管事件发生的危险因素，以有无心血管事件发生作为因变量进行回归分析。R代码如下：

#载入regicor数据集
library(compareGroups)
data("regicor")
regicor<-na.omit(regicor)	#删除缺失值
#为避免繁杂的计算，编写循环函数进行批量单因素回归分析
Uni_glm_model=function(x){ 
  FML=as.formula(paste0("cv~",x)) 	#构建分析
  glm1<- glm(FML,family = binomial,data = regicor) 	#单因素分析
  glm2=summary(glm1) 	#处理分析结果
  OR=round(exp(coef(glm1)),2) 	#提取风险率
  SE=glm2$coefficients[,2]
  CI5=round(exp(coef(glm1)-1.96*SE),2) 	#计算风险率范围
  CI95=round(exp(coef(glm1)+1.96*SE),2)
  CI=paste0(CI5,"-",CI95)
  P=signif(glm2$coefficients[,4],3) 	#提取p值，保留2位小数
  Uni_glm_model <- data.frame("characteristics"=x,
                              "OR"=OR,
                              "CI"=CI,
                              "p"=P)[-1,]
  return(Uni_glm_model)
}
#提取需要纳入分析的变量
variable.names=colnames(regicor)[c(2:5,10:13)]
#应用函数
Uni_glm=lapply(variable.names,Uni_glm_model)
#把列表转换为数据框
Uni_glm<-do.call(rbind, Uni_glm)
Uni_glm$characteristics<-rownames(Uni_glm)
rownames(Uni_glm) <- NULL
#输出结果
print(Uni_glm)

代码的运行结果如图5-7所示。图中以有无心血管事件发生为分组变量，统计了各个变量的组间差异。

图5-7

为了进一步找出独立危险因素，将在单因素回归中具有统计学意义的变量纳入多因素回归模型。R示例代码如下：

#将单因素回归有统计学意义的变量(p<0.05)纳入多因素回归
milti_glm_model<- glm(cv~age+smoker+triglyc+ldl,
                      family = binomial,
                      data = regicor)
summary(milti_glm_model)
library(broom)
Muti_uni_result <- data.frame(exp(coef(milti_glm_model)),
                              exp(confint(milti_glm_model)),
                              tidy(milti_glm_model)$p.value)
names(Muti_uni_result) <- c("OR","CI5","CI95","P")
Muti_uni_result

代码运行结果如图5-8所示，该图展示了多因素回归分析的结果：OR（odds ratio，比值比）值、P值、CI95。从结果中可以看出，年龄、吸烟和高甘油三酯是心血管事件发生的独立危险因素。

图5-8

此外，当我们的研究目的是评估两个变量之间的相关性并消除其他混杂因素的影响时，也可以通过同样的代码输出不同调整混杂因素模型的回归分析结果。如果我们想进一步分析吸烟和死亡事件发生的相关性，并计算未调整混杂因素、调整年龄和性别的影响以及调整年龄、性别、胆固醇和甘油三酯的影响的结果，可以使用如下R代码：

#评估吸烟和胆固醇水平的关系
#未调整混杂因素model1
smoke_chol_model1<- glm(death~smoker,family = binomial,data = regicor)
smoke_chol_result1 <- data.frame(exp(coef(smoke_chol_model1)),
                                 exp(confint(smoke_chol_model1)),
                                 tidy(smoke_chol_model1)$p.value)                 
names(smoke_chol_result1) <- c("OR","CI5","CI95","P")
smoke_chol_result1

#调整年龄和性别的影响model2
smoke_chol_model2<- glm(death~smoker+age+sex,family = binomial,
                        data = regicor)
smoke_chol_result2 <- data.frame(exp(coef(smoke_chol_model2)),
                                 exp(confint(smoke_chol_model2)),
                                 tidy(smoke_chol_model2)$p.value)                 
names(smoke_chol_result2) <- c("OR","CI5","CI95","P")
smoke_chol_result2

#调整年龄、性别、甘油三酯和胆固醇的影响model3
smoke_chol_model3<- glm(death~smoker+age+sex+chol+triglyc,family = binomial,data = regicor)
smoke_chol_result3 <- data.frame(exp(coef(smoke_chol_model3)),
                                 exp(confint(smoke_chol_model3)),
                                 tidy(smoke_chol_model3)$p.value)                 
names(smoke_chol_result3) <- c("OR","CI5","CI95","P")
smoke_chol_result3

运行结果如图5-9所示，该图展示了多因素回归分析中三种不同调整混杂因素模型计算的吸烟和死亡事件相关性的结果（OR值、P值、CI95%），整理成表格如图5-10所示。

图5-9

图5-10