R语言泊松回归并分层人年发病率统计分析

在既往文章中，我们已经介绍了R语言计算人年及可信区间的计算。但是计算的是总的人年发病率的比较情况，假如我们想知道分层发病率的情况呢？拿既往乳腺癌的数据为例子，我们已经知道了有淋巴结肿大和没有淋巴结肿大患者总的生存率的比较，但是如果我们想了解在每个年龄段有淋巴结肿大和没有淋巴结肿大患者生存率有无区别？如下图

我们以R语言survival包演示泊松回归年龄分层发病率统计，继续使用我们的乳腺癌数据（公众号回复：乳腺癌可以获得数据），首先把数据和包导入

library(reshape)
library(survival)
library(dplyr)
bc<-read.csv("E:/r/test/Breast cancer survival agec.csv",sep=',',header=TRUE)
bc<- rename(bc,c(锘縤d= "id")) 

我们先来看看数据：
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。

有些变量是分类变量，我们先把它转换成因子，把时间除以12，让它的单位变成年

bc$pyears<-bc$time/12
bc$ln_yesno<-factor(bc$ln_yesno)

因为我们想知道不同年龄的生存率，所以我们需要对年龄进行分割把它变成0-30岁，30-40岁，40-50岁，50-60岁，60-70岁，70-100岁这么多组

bc$age1<-cut(bc$age,breaks=c(0,30,40,50,60,70,100),include.lowest=T,
             labels = c(1,2,3,4,5,6))
bc$age1<-factor(bc$age1)

建立泊松回归

fit1 <- glm(status ~ln_yesno+age1+offset(log(pyears)),
            data=bc, family=poisson)
summary(fit1)

提取系数转换

coef(fit1)
exp(coef(fit1))

解读和其他回归是一样的，可以看出有无淋巴结肿大对生存率有影响，随着年龄增大对生存率影响不大。如果年龄的P是＜0.05有意义的，从OR表明60-70岁组是生存率最低的。
泊松回归也可以作图

par(mfrow=c(1,4))
plot(fit1)

接着我们可以使用survival包的pyears函数进行人年统计

fit2 <- pyears(Surv(time/12, status) ~cut(age, c(0,30,40,50,60,70,100))+
                 ln_yesno, data =bc, scale = 1)
summary(fit2)

这样我们的人年统计的表就出来啦，继续统计分层人年发病率

fit2$event/fit2$pyears

这样的话发病率也出来啦，想求1000人年发病率乘以1000就可以了

fit2$event/fit2$pyears*1000

同理看出60-70岁组是生存率最低的。
我们可以把这部分数据提取出来作图表示，这样直观一点
先提取数据，并对它进行格式转换

bb1<-fit2$pyears
bb1<-as.data.frame(bb1)

给它增加行名并转成宽格式

bb1<- tibble::rownames_to_column(bb1, "age")
be<-melt(bb1,id=c("age"),
         measure.vars = (c("0","1")),
         variable.name = "lnyesno",
         value.name = "value") ##ID为固定不变的变量，measure.vars为需要整合的变量，variable.name 为新变量名字

转换好格式以后就可以作图了

library(ggplot2)
ggplot(be, aes(be$age, be$value,group=variable)) +  
  geom_point()+geom_line()

美化一下

ggplot(be, aes(be$age, be$value,fill=variable,group=variable)) +  
  geom_point(shape=21,size=4,col="black")+geom_line(linetype=1,size=1)

0-30岁这组因为没有发病所以显示为0，并不是这组生存率最低，我们可以把这个数据去掉

be[1,3]<-NA
ggplot(be, aes(be$age, be$value,fill=variable,group=variable)) +  
  geom_point(shape=21,size=4,col="black")+geom_line(linetype=1,size=1)

OK，表明随着年龄增大，生存率降低，60-70岁这个年龄最低。