5. 基础统计方法——关联性分析/简单线性回归

最新推荐文章于 2023-12-30 14:36:57 发布

Dr_long1996

最新推荐文章于 2023-12-30 14:36:57 发布

阅读量2.6k

点赞数 2

分类专栏：临床统计分析 R语言文章标签： r语言

本文链接：https://blog.csdn.net/Dr_long1996/article/details/131054609

版权

R语言同时被 2 个专栏收录

16 篇文章 16 订阅

订阅专栏

临床统计分析

13 篇文章 3 订阅

订阅专栏

基本概念：

两个连续随机变量间的线性练习称为线性相关，两个分类变量间的联系称为关联。

相关种类：

正相关：一般若Y有随着X增大而线性上升的趋势，称为正相关；
负相关：若Y有随着X增大而线性下降的趋势，称为负相关；
零相关：若Y或X不随另一个变量的改变而改变，则称为零相关；
非线性相关：若散点呈曲线形状，则变量可能呈曲线关系，不宜做线性相关

关联强度的评价指标：

两连续型随机变量间联系的强度用相关系数（correlation coefficient）来描述。

两分类变量间联系的强度用关联系数（association coefficient）来描述。

临床含义：相关系数最大为1，最小为-1，联系越密切则相关系数接近于1或-1，没有线性联系的相关系数为0。相关系数为正，称为正相关，相关系数为负，称为负相关，相关系数为0称为零相关。

连续型随机变量间的相关分析

pearson相关是要求两个变量都服从正态分布，spearman相关系数和kendall's Tau相关系数都是非参数的相关度量，对数据分布没有特殊要求，或者等级资料。

1. 相关系数的计算

案例：使用数据集birthwt，将其中的连续数值变量提取出来形成数据框cont.vars，计算变量age，lwt，bwt的相关性。

data(birthwt,package = "MASS")
cont.vars <- dplyr::select(birthwt,age,lwt,bwt)#将连续性变量单独提取出来

根据中心极限定理：当每组样本量>=50时，可以认为样本的均数近似服从正态分布。所以可以选择person相关，使用函数cor()计算样本的相关系数，method=“pearson”或“spearman”参数。

cor(cont.vars,method = "pearson")
cor(cont.vars,method = "spearman")

函数cor()计算相关系数，cov()和var()计算协变量。容易混淆！！！！

2. 相关系数的假设检验

方法一：base包中的cor.test()函数，但是只能每次检验两个连续变量的相关性。

cor.test(cont.vars$lwt,cont.vars$bwt,method = "spearman")
cor.test(cont.vars$lwt,cont.vars$bwt,method = "pearson")

从spearman秩相关的假设检验可以看到相关系数 $\rho$ =0.2488882，p=0.00055<0.001，因此拒绝H0，接受H1，可以认为，母亲最后一次月经时的体重（lwt）与新生儿出生时体重（bwt）之间存在正相关关系。

从pearson积矩相关系数的假设检验可以看到相关系数 $\rho$ =0.1857，相关系数95%置信区间为（0.044，0.31998），p=0.0105<0.05，按照α=0.05水准，拒绝H0，接受H1，即认为两变量间线性相关有统计学差异。可以认为，母亲最后一次月经时的体重（lwt）与新生儿出生时体重（bwt）之间存在正相关关系。

方法二：cor.test只能检验一个相关系数，psych包中corr.test()可以计算多个相关系数矩阵并进行检验

library(psych)
corr.test(cont.vars)
corr.test(cont.vars,method = "spearman")

首先输出一个相关系数矩阵，矩阵的相关性检验p值。默认方法method=pearson，还可以选择spearman和kendall

置信区间：

cor.test会展示自信区间，但是corr.test没有，那么如何操作呢？

print(corr.test(cont.vars),short = F)#展示自信区间
print(corr.test(cont.vars,method="spearman"),short = F)#展示自信区间

3. 偏相关

偏相关是指控制一个或者多个数值型变量时，另外两个变量之间的相关性。使用ggm包中pcor计算偏相关系数。

library(ggm)
names(cont.vars)
#[1] "age" "lwt" "bwt"
r <- pcor(c(2,3,1),cov(cont.vars))
#pcor()第一个参数中的向量的前两个数字代表要计算相关规系数的变量的下标，其余的数字待变条件变量（控制变量）的下标
nrow(cont.vars)
#[1] 189
pcor.test(r,q=1,n=189)
#q表示条件变量的个数,这里只有一个变量(age)，n表示样本数

上面的结果表明，在控制了年龄的影响时，母亲体重和新生儿体重的相关系数有统计学意义。

4. 相关的解释——应注意的问题

（1）如果从散点图可见两随机变量呈线性相关的趋势，假设检验也拒绝总体相关系数为0的假设，则我们可以推断两变量是“线性相关”的，但我们不能因此推断两变量在生物学上有任何的联系，甚至认为呈因果关系。

（2）如果经检验后，不拒绝总体相关系数 $\rho$ =0，我们不能轻易下结论认为两个变量“无关”，应注意样本含量是否够大，其次散点图是不是呈现曲线关系，是否对应资料进行分层分析等。

（3）异常值的处理：校对原始数据，其次分别保留和剔除后进行一次分析。

分类变量间的关联分析

关联分析思想：对分类变量间联系，可作关联分析，即对两个分类变量交叉分类，对计数所得的频数资料（列联表）作关于两种属性独立的 $\chi ^{2}$ 检验。

$\chi ^{2}$ 检验的基本方法可以参考：5.基本统计方法-分类变量的组间比较_Dr_long1996的博客-CSDN博客

1. 交叉分类2 X 2表的关联分析

案例：收集了3154例冠心病患者的观察资料，研究者将观察对象的行为分为A型和B型。分析行为类型和冠心病的发病是否有关联性。

my.matrix <- matrix(c(178,79,1411,1486),nrow = 2);my.matrix
row.names(my.matrix)=c("A","B")#行为类型
colnames(my.matrix)=c("+","-")#是否患有冠心病
chisq.test(my.matrix)
chisq.test(my.matrix)$expected
chisq.test(my.matrix,correct = F)

如果一种属性的概率分布与另一种属性的概率分布无关，则称为两种属性相互独立（independence），否则称这两种属性之间存在关联性。

$\chi ^{^{2}}$ =39.90，自由度df=1，p值<0.01，说明行为类型与冠心病之间存在着关联性。

关联系数（association coefficient，r）

如果两个数据不独立，那么就要进一步了解相关性的强弱。使用vcd包中assocstats()来计算列联表的phi系数，列联系数和Cramer’s V系数，其中Phi系数只适用四格表

行为类型与冠心病之间存在关联，但是关联程度如何？

assocstats(my.matrix)

从上述结果可以看到关联系数r=0.112。

对于2 X 2的列联表而言，关联系数r介于0和 $\sqrt{0.5}$ 之间，其数值越大，关联程度越高；总体的关联系数是否等于0的假设检验与上述关于两种属性的独立性检验等价。

2. 2 $\times$ 2配对资料的关联性分析

案例：对103名患者进行影像学检查（A）和生化检查（B），结果均分为疾病（+）和正常（-）两类，现欲分析A、B两法检查结果的关联性。

my.matrix <- matrix(c(50,8,15,30),nrow = 2);my.matrix
row.names(my.matrix)=c("+","-")
colnames(my.matrix)=c("+","-")
mcnemar.test(my.matrix)

使用McNemar检验解决了两种方法的阳性率是否相等的问题，但若要了解两种方法之间的关联，需要作两种属性的关联分析。

chisq.test(my.matrix,correct = F)

assocstats(my.matrix)

本案例中 $\chi ^{2}$ =30.43，p<0.05，认为两种之间存在关联，关联系数r=0.478。与Mcn检验的p=0.2109，尚不能认为两种方法之间有差异。本质含义是相同的，但是理解的思路却不同。

配对设计的此类资料类型既可以进行频率的比较（McNemar检验），又可以进行关联分析（chisq.test()），资料的整理完全一致，但由于分析目的不同，计算方法也不同：进行频率比较时，McNemar检验仅对配对中不一致的数据进行检验。

3. 多分类资料的关联分析

案例：有研究表明，不同国籍人的血型分布不同的，现有2500不同国籍人的血型分布资料，请问国籍与血型的分布是否有关？

my.matrix <- matrix(c(450,300,190,
                      410,250,250,
                      100,350,40,
                      40,100,20),nrow = 3);my.matrix
row.names(my.matrix)=c("美国人","中国人","挪威人")
colnames(my.matrix)=c("O","A","B","AB")
addmargins(my.matrix)

现在的“关联”和以前的“比较”是两类不同的问题和设计。多组频率分布资料的比较，其设计是，首先确定三种国籍人群，调查其血型分布资料，分析三种国籍的血型分布是否不同。本例2500例患者是一次调查的结果，可以看作是总体种一份随机抽样，按照两个属性交叉分类，统计频数，分析的是两种属性之间的是否有关联。

chisq.test(my.matrix)

$\chi ^{^{2}}$ =332.97，P<0.001，拒绝H0，说明血型与国籍之间有关联性。

assocstats(my.matrix)

关联系数r=0.343

关联性分析与多个频率比较的 $\chi ^{2}$ 检验一样，要求不能有超过1/5以上的理论频数小于5，或者不能有1个理论频数小于1.

4. 两个有序多分类变量

spearman秩相关法进行分析或者Kendall’s Tau-b相关性分析

案例：欲研究年龄(Age)与冠状动脉粥样硬化等级(Grade)之间的关系，抽样调查了283例年龄≥30岁的居民，收集其年龄和冠状动脉粥样硬化等级数据。年龄按30-39岁、40-49岁、50-59岁、≥60岁统计频数，冠状动脉粥样硬化按照0级、1级、2级、3级统计频数。

age=c(rep(1,70+23+4+3),rep(2,27+25+9+4),rep(3,16+23+13+8),rep(4,9+20+15+14))
grade=c(rep(1,70),rep(2,23),rep(3,4),rep(4,3),
    rep(1,27),rep(2,25),rep(3,9),rep(4,4),
    rep(1,16),rep(2,23),rep(3,13),rep(4,8),
    rep(1,9),rep(2,20),rep(3,15),rep(4,14))
mydata=data.frame(age,grade)
mydata$age=factor(mydata$age,
                  levels = c(1,2,3,4),
                  labels = c("30-39","40-49","50-59","60-69"))
mydata$grade=factor(mydata$grade,
                    levels = c(1,2,3,4),
                    labels = c("0级","1级","2级","3级"),ordered = T)

将原始数据制作成为列联表，然后进行关联分析。

#创建列联表
mytable <- table(mydata$age,mydata$grade);mytable
#统计描述——计算构成比
prop.table(mytable,margin = 1) #按照行计算百分比
prop.table(mytable,margin = 2)#按照列计算百分比
#关联性分析
#install.packages("DescTools")
library(DescTools) #调用包“DescTools”
KendallTauB(mytable,conf.level=0.95) #计算Tau-b系数及其95%CI

年龄与动脉粥样硬化等级间的Kendall’s Tau-b相关系数为0.410(95%CI：0.33-0.49)，即年龄与动脉粥样硬化的等级有相关性，且为正相关

注意：kendallTauB()主要用于列联表频数格式的数据，如果原始的长数据格式还可以直接用cor()函数，设置参数method=“kendall”实现上述关联分析。

cor.test(age,grade,use="everything",method = "kendall")

关联分析总结：

两个连续变量间呈线性相关时，可以使用Pearson相关分析，不满足Pearson相关分析的适用条件时，可以使用Spearman相关系数来描述。
两个连续变量既可以使用Pearson相关分析，也可以使用Kendall's tau-b等级相关系数描述，但后者更多适用于两个分类变量均为有序分类的情况(也可以用于有序分类变量+连续变量)。

关联分析的可视化

1. 连续型变量散点图

利用连续型变量部分的数据进行分析：

data(birthwt,package = "MASS")
cont.vars <- dplyr::select(birthwt,age,lwt,bwt)#将连续性变量单独提取出来

方法一：pairs()绘制散点图

pairs(cont.vars)

方法二：car包中的函数scatterplotMatrix用于创建点图矩阵

library(car)
scatterplotMatrix(cont.vars)

方法三：其他方式展示两变量之间相关性

#install.packages("corrplot")
library(corrplot)
#install.packages("corrgram")
library(corrgram)
corrplot(cor(cont.vars))#颜色深浅表示相关性强弱，其中蓝色表示正相关，红色负相关
corrgram(cont.vars,upper.panel = panel.pie)#upper.panel=展示的图形

利用相关系数绘制的圈图展示两变量之间的相关性。

2. 分类变量马赛克图

案例：使用R自带数据集Arthritis，计算治疗方式（安慰剂和标准治疗）于疗效之间的关联程度，原始数据表格如下：

library(vcd)
mytable <- table(Arthritis$Treatment,Arthritis$Improved);mytable

将数据统计形成一个列联表，分类变量:用关联图或者马赛克图展示

方法一：关联图

assocplot(mytable)

从上图可以看到安慰剂治疗于疗效none正相关，与显著疗效负关联，治疗组同理。

方法二：马赛克图

mosaic(~Treatment+Improved,data=Arthritis)#两变量的马赛克图
mosaic(~Sex+Treatment+Improved,data=Arthritis)#多变量的马赛克图

如果要在一张图中展示数值型变量和分类变量的关联，可以使用GGally包中的ggpair

#install.packages("GGally")
library(GGally)
library(ggplot2)
dat <- dplyr::select(birthwt,age,lwt,bwt,smoke,race)
ggpairs(dat)

简单线性回归

简单线性回归：一个变量Y（反应变量）随着另一个变量X（解释变量）的变化而变化，且呈直线变化趋势

线性回归分析的要求数据满足线性、独立性、正态性、等方差的前提

线性：反应变量Y与自变量X呈直线变化趋势，通过散点图来考察两变量是否呈线性趋势

独立性：两个个体的取值是相互独立的，一个个体的取值不受另一个个体的影响，多于非独立样本，可以选择重复测量资料的方差分析、多水平模型。

正态性：在给定x取值时，y的取值服从正太分布，也即y的残差服从正态分布。通过残差图或者正态概率图考察

等方差：对于不同的x值，y的总体变异相同

案例：某种疾病患病后会影响尿肌酐的表达，收集了一部分患病和正常患者的尿肌酐，欲研究尿肌酐是否会随着年龄的改变有变化，数据如下：

UCR <- data.frame(age=c(13,11,9,6,8,10,12,7,10,9,11,12,15,16,8,7,10,15),
                  ucr=c(3.54,3.01,3.09,2.48,2.56,3.36,3.18,2.65,3.01,2.83,2.92,3.09,3.98,3.89,2.21,2.39,2.74,3.36),
                  group=c(rep(0,8),rep(1,10)))

在数据栏目中可以查看全部数据情况，数据集中共有3个变量和18个观察数据，3个变量分别代表年龄(age)、尿肌酐(ucr)与是否患病(group)。

1. 线性回归条件的判断

线性：

#基础包绘图
plot(ucr~age,data = UCR) #满足线性回归第一个条件线性
plot(ucr~age,data = UCR,xlab="Age in years",ylab="Urine creatinine(mmol)")

#ggplot绘制散点图
library(ggplot2)  #启动ggplot2软件包
ggplot(data=UCR,aes(x=age,y=ucr))+
  geom_point()+
  stat_smooth(method="lm",se=TRUE) #绘制散点图

散点大致呈一条直线，提示变量“age”和“ucr”存在线性关系，该条件满足

独立性：独立性需要针对实验设计和观察指标进行判断。

正态性（四种方法）：

lm.UCR<-lm(age~ucr,data=UCR) #回归分析
shapiro.test(lm.UCR$residuals)#进行shapiro正态性检验

hist(lm.UCR$residuals)#方法二：通过直方图展示残差的分布判断

#方法三：对于小样本，可以使用qq图展示，是否聚集在一条直线上
qqnorm(lm.UCR$residuals)
qqline(lm.UCR$residuals)

#方法四：残差图：做残差与拟合值的散点图，是否对称
plot(fitted(lm.UCR),lm.UCR$residuals,xlab = "Fitted values",type="n") #绘制散点图，type=n表示不显示散点
text(fitted(lm.UCR),lm.UCR$residuals,labels=rownames(UCR))#标记每个点的名称，text表示加上文本的意思
abline(h=0,col="blue")

残差的方差齐性：

par(mfrow=c(1,3)) #绘制一行3个图
plot(lm.UCR$residuals)
plot(UCR$age,lm.UCR$residuals)
plot(UCR$ucr,lm.UCR$residuals)

预测值和各变量值的残差分布较为均匀，并未出现特殊的分布形式(如漏斗或者扇形)，提示残差的方差齐。

2. 线性模型拟合

lm.UCR=lm(ucr~age,data=UCR);lm.UCR #输出的数据有限
attributes(lm.UCR)#通过函数attributes()查看模型对象的属性

上述就是线性拟合后生成的拟合模型包含的内容，可以通过$提取具体属性的值。比如lm.UCR$residual就能将残差提出出来。

lm.UCR$coefficients#显示的是截距和斜率（回归系数）
lm.UCR$residuals#每一个样本的残差

summary(lm.UCR)#生成主要的线性拟合参数

上述的总结，显示了拟合公式lm是ucr~age的拟合，第二部分是残差的基本信息，第三部分Coefficients部分”列出了截距和自变量的“Estimate (非标准化系数)”、“Std.Error (标准误)”，统计量t值及P值。Multiple R-squared：决定系数。F统计量：对线性回归方程进行统计检验。

回归方差模型假设检验包含两个方面：一是检验回归模型是否成立的模型检验（F检验），二是检验总体回归系数β是否为零的参数检验（t检验）。

决定系数是回归分析中重要的统计量，定义为回归平方和与总平方和之比，记为 $\boldsymbol{R}^{^{2}}$ 。其值的大小反映了自变量对回归效果的贡献，在Y的总变异中回归关系所能解释的百分比。常用作反映拟合优度的指标。当X与Y均为随机变量时，决定系数等于相关系数r的平方和。

该案例中，年龄（age）的斜率为0.149，检验p<0.001，即age的回归系数不等于0。F统计量=60.95，p<0.001，说明回归模型成立。决定系数为0.7921，说明年龄有很大程度可以解释尿肌酐的水平。项在该案例中没有实际意义，p值仅表示常数项与0之间有显著差异。但是 $\mathrm{R^{2}}$ 会受自变量个数的影响，夸大自变量对因变量变异的解释程度，自变量越多， $\mathrm{R^{2}}$ 越大。Adjusted $\mathrm{R^{2}}$ 调整了自变量个数对结果的影响，一般小于 $\mathrm{R^{2}}$ 。Adjusted $\mathrm{R^{2}}$ =0.7791，说明自变量 (年龄)可以解释77.91%的因变量的变异 (尿肌酐)。

如果要具体分析回归模型的假设检验内容：

summary(aov(lm.UCR))

回归决定系数 $\boldsymbol{R}^{^{2}}$ = $\frac{SS_{回归}}{SS_{总}}$

注意：服从双变量正态分布的同样一组资料，同时做了相关分析和回归分析，则相关系数的t检验与回归系数t检验等价

分层线性回归

分层线性回归：在分组变量下，两个变量之间的线性回归关系

本质：对回归方程进行比较的问题

（1）：不同分组中的回归直线是否平行；（2）：若平行，检验其截距是否相等

#假定年龄对尿肌酐在两组中的影响是相同的

mod1 <- lm(ucr~age+group,data = UCR)
summary(mod1)

年龄每增长1岁，尿肌酐就增加0.16156，相同年龄情况下，患病儿童的尿肌酐比正常儿童低0.232356。

#分组绘制散点图
col=ifelse(UCR$group==0,"blue","red")#设置两个组的颜色
pch=ifelse(UCR$group==0,1,19)#设置两个组点的形状
plot(ucr~age,data = UCR,col=col,pch=pch,xlab="Age in years",ylab="Urine Creatinine")
legend("topleft",
       legend = c("Normal children","Diseased children"),
       col=c("blue","red"),
       pch=c(1,19))

分别对两组绘制拟合曲线，分别求解不同模型的斜率和截距。

mod1$coefficients

第一个是截距，第二个分别是age的斜率，前文假设了年龄在两组之间尿肌酐的影响是相同的，所以斜率都是相同的b=0.16156，但截距没有患病组group=0，截距a=1.45，患病组group=1，截距a=1.45-0.233。

mod1$coefficients
b <- mod1$coefficients[2]
a0 <- mod1$coefficients[1]
a1 <- mod1$coefficients[1]+mod1$coefficients[3]
abline(a=a0,b=b,col="blue")
abline(a=a1,b=b,col="red")

在上面的模型里，假定的是儿童年龄对尿肌酐含量具有恒定效应。进行假设检验，尝试不同的斜率来拟合回归直线，然后比较两者斜率的差异是否显著

！！！！建立一个group和age之间存在交互的模型

mod2 <- lm(ucr~age+group+age:group,data = UCR)
#建立了一个age和group存在交互的模型
mod2$coefficients

所以回归方程中y=1.66+0.139*age-0.566*group+0.033*age*group

再次对存在交互模型的斜率和截距进行计算。

健康组中group=0，所以截距a2，斜率b2

a2 <- mod2$coefficients[1]
b2 <- mod2$coefficients[2]

在患病组中group=1，所以截距a3，斜率b3

a3=mod2$coefficients[1]+mod2$coefficients[3]
b3=mod2$coefficients[2]+mod2$coefficients[4]

绘制散点图和拟合曲线：

col=ifelse(UCR$group==0,"blue","red")
pch=ifelse(UCR$group==0,1,19)
plot(ucr~age,data = UCR,col=col,pch=pch,xlab="Age in years",ylab="Urine Creatinine")
legend("topleft",
       legend = c("Normal children","Diseased children"),
       col=c("blue","red"),
       pch=c(1,19))
abline(a=a2,b=b2,col="blue")
abline(a=a3,b=b3,col="red")