亚组分析、P交互、P趋势是什么？如何计算呢？

亚组分析、P交互、P趋势是什么？如何计算呢？

（1）亚组分析如何计算？

（2）P交互作用的计算方法？

（3）P趋势如何计算？

（1）亚组分析如何计算？

本篇不说理论，直接以R语言的lung数据集为例进行介绍。

除P交互的外，亚组分析和P趋势的计算以logistic模型为例进行说明：

自变量：sex（性别），age（年龄），ph.ecog（ECOG评分），wt.loss（过去六个月的体重减轻）

因变量：status（是否死亡）

交互作用我们主要关注性别和年龄的交互。

假设：我们主要关注性别和死亡的关系？

首先：由于要研究性别和死亡的关系，因此我们建立的主要分析模型如下：

Logit(p)=β1sex+β2age+β3ph.ecog+β4wt.loss+残差

在这里，β1是我们关注的主要解释参数。age+ph.ecog+wt.loss都是协变量

R语言代码：

#数据处理

library(rms)

library(dplyr)

data<-lung

#生成二分类新变量

data$age1 <- ifelse(data$age > 60,

                        c("1"), c("2"))

data$status<-data$status-1

str(data)#查看数据类型

#转化因子

a<-c("sex","age1")#填入需要转化的变量

data[,a]<-lapply(data[,a],as.factor)#转因子

#建立模型

#logistic模型：sex和status的关系

fit1 <- glm(status~sex+age1+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit1)

#自编函数，计算OR值

formatFit<-function(fit){

  #取P值

  p<-summary(fit)$coefficients[,4]

  #wald值

  wald<-summary(fit)$coefficients[,3]^2

  #B值

  valueB<-coef(fit)

  #OR值

  valueOR<-exp(coef(fit))

  #OR值得95%CI

  confitOR<-exp(confint(fit))

  data.frame(

    B=round(valueB,3),

    Wald=round(wald,3),

    OR_with_CI=paste(round(valueOR,3),"(",

                     round(confitOR[,1],3),"~",round(confitOR[,2],3),")",sep=""),

    P=format.pval(p,digits = 3,eps=0.001)

  )

}

formatFit(fit1)#得出OR值

结果解释：OR=0.304，因此sex=2的死亡风险更低。

1. 亚组分析如何计算

那么，为了进一步探讨，不同年龄组的性别和死亡的关系，我们需要将人群按照年龄组分为60岁以下和60岁以上两组，分别建立两个模型：

60岁以下：Logit(p)=β1sex+β3ph.ecog+β4wt.loss+残差

60岁以上：Logit(p)=β1sex+β3ph.ecog+β4wt.loss+残差

这时候模型中就去掉了age这个变量，然后会得到两个β1，这其实就是亚组分析，也叫分层分析。

R语言代码：

#亚组分析

data1<-subset(data,age>60)

data2<-subset(data,age<=60)

#60岁以上

fit2 <- glm(status~sex+ph.ecog+wt.loss,data=data1,family = "binomial")

#60岁以下

fit3 <- glm(status~sex+ph.ecog+wt.loss,data=data2,family = "binomial")

formatFit(fit2)

formatFit(fit3)

结果解释：OR=0.284和0.309，因此60岁以上和60岁以下，sex=2的死亡风险均更低。

（2）P交互作用的计算方法？

那么OR=0.284和0.309，两者相等吗？这时候就需要计算P交互值了。

P交互值的计算方法有两种，一种是传统方法，一种是似然比检验。

1. 传统方法：在模型中加入主要解释变量和亚组变量的交互项，交互项的P值即为P交互，本例的模型如下：

Logit(p)=β1sex+β2age+β3sex*age+β4ph.ecog+β5wt.loss+残差

β3的P值即为P交互。

1. 似然比检验：步骤有三步

步骤1：建立有交互的模型，即Logit(p)=β1sex+β2age+β3sex*age+β4ph.ecog+β5wt.loss+残差

步骤2：建立无交互的模型，即Logit(p)=β1sex+β2age+β3ph.ecog+β4wt.loss+残差

步骤3：比较这个模型差异，计算的P值即为P交互

若P交互<0.05，则说明OR=0.284和0.309在统计学上有差异，否则无差异。

下面分别给出线性模型、logistic模型、cox模型传统方法和似然比检验计算交互P的代码：

#模型1：线性模型

#有交互模型

fit1 <- glm(meal.cal~sex*age1+sex+age1+ph.ecog+wt.loss,data=data,family = "gaussian")

broom::tidy(fit1)

#无交互模型

fit2 <- glm(meal.cal~sex+age1+ph.ecog+wt.loss,data=data,family = "gaussian")

broom::tidy(fit2)

#似然比检验

anova(fit1,fit2,test="Chisq")

#模型2：logistic模型

# 有交互模型

fit1 <- glm(status~sex*age1+sex+age1+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit1)

#无交互模型

fit2 <- glm(status~sex+age1+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit2)

#似然比检验

anova(fit1,fit2,test="Chisq")

#模型3：Cox模型

# 有交互模型

fit1 <- coxph(Surv(time,status)~sex*age1+sex+age1+ph.ecog+wt.loss,data=data)

broom::tidy(fit1)

#无交互模型

fit2 <- coxph(Surv(time,status)~sex+age1+ph.ecog+wt.loss,data=data)

broom::tidy(fit2)

#似然比检验

anova(fit1,fit2,test="Chisq")

（3）P趋势如何计算？

由于P趋势的计算要求自变量至少为三种，因此，这里研究目的改为年龄和死亡的关系。

首先将年龄划分为<50，50~60，60~70，70及以上四组。

然后首先将年龄设置为哑变量，建立以下模型：

Logit(p)=β1age12+β2age13+β3age14+β4sex++β5ph.ecog+β6wt.loss+残差

我们会得到β1、β2、β3三个参数估计值，然后进一步得出OR1、OR2、OR3.

R语言代码：

###########P趋势的计算

#年龄分组

data <- within(data,{

  age1[age<50]=1

  age1[age>=50 & age < 60] = 2

  age1[age>=60 & age < 70] = 3

  age1[age>= 70] = 4

})

#变因子

data$age1<-as.factor(data$age1)

#设参照

data$age1<-relevel(data$age1, ref="1")

#建立模型

fit1 <- glm(status~age1+sex+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit1)

formatFit(fit1)

OR分别为2.055、1.454、2.433，这三个值在统计学有没有趋势性变化呢？这时候就需要计算P趋势值了。

目前有两种方法计算P趋势，以本例为例子进行说明

（1）将age变量编码为1、2、3、4，作为等级变量纳入模型，计算出的P值就是P趋势，P趋势的方向是正向还是反向，可以看估计的β值或OR值

（2）先将age变量编码为1、2、3、4，然后计算每一组的中位数，本例计算的中位数为。。。。，因此将原来的1、2、3、4替换成。。。。/，再作为等级变量纳入模型，计算出的P值就是P趋势，P趋势的方向是正向还是反向，可以看估计的β值或OR值

R语言代码：

############方法1：age1变为等级变量

#age1变为数值

data$age1<-as.numeric(data$age1)

#建立模型

fit1 <- glm(status~age1+sex+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit1)

formatFit(fit1)

############方法2：age1各组的中位数变为等级变量

library(dplyr)

# median of groups

summarise(group_by(data,age1),

          min(age),

          max(age),

          mean(age),

          sd(age),

          median(age),

          quantile(age,0.25),

          quantile(age,0.75),

          )

#age的中位数:44，56，64，7

data$age2[data$age1 == 1]<- 44

data$age2[data$age1 == 2]<- 56

data$age2[data$age1 == 3]<- 64

data$age2[data$age1 == 4]<- 73

#建立模型

fit1 <- glm(status~age2+sex+ph.ecog+wt.loss,data=data,family = "binomial")

broom::tidy(fit1)

formatFit(fit1)

结果解释：P趋势值>0.05，因此不同年龄的死亡风险无趋势性变化。

这里可以发现两种做法结果不太一致，目前主流的做法基本是第2种做法，因此，小编建议大家参考第2种做法，当然第1种做法也没错。