有粉丝在后台问我人年发病率的可信区间怎么计算?如下图:
我对人年也不是很熟悉,于是查一下资料,
人年发病率的计算,就是每个进入队列的患者的随访时间求和,就是总的人年数,然后用发生结局的人数/总的人年数*1000就是per1000 person-years,人年的发病率置信区间是通过泊松回归分布进行计算的,计算原理如下图
有了公式计算就容易多了,主要步骤为:
- 先求出每个队列的发病数和人年总数
- 发病数/人年总数×1000即为人年发病率
- 求出相对发病率RR,
- 求出R的对数log(RR)
- 代入最后的公式计算就可以求出置信区间了
本章我们通过R语言来实现人年计算,继续使用我们的乳腺癌数据(公众号回复:乳腺癌可以获得数据),这里我们对数据进行一点小小改变,就是把它改为csv格式(点另存为带分隔符的CSV),为了以后好和stata做出来进行比较。
先导入数据
library(reshape)
library(dplyr)
bc<-read.csv("E:/r/test/Breast cancer survival agec.csv",sep=',',header=TRUE)
head(bc)
names(bc)
我们先来看看数据:
age表示年龄,pathsize表示病理肿瘤大小(厘米),lnpos表示腋窝淋巴结阳性,histgrad表示病理组织学等级,er表示雌激素受体状态,pr表示孕激素受体状态,status结局事件是否死亡,pathscat表示病理肿瘤大小类别(分组变量),ln_yesno表示是否有淋巴结肿大,time是生存时间,后面的agec是我们自己设定的,不用管它。
我们可以看到第一列ID部分有点乱码,我们帮它改一下名字,改不改都可以,我这该死的强迫症
bc<- rename(bc,c(锘縤d= "id"))
改好了以后我们就可以正式分析了,
假设我们想比较有无淋巴结肿大患者(ln_yesno)的人年发病率及可信区间,首先我们要把它分成两个数据,一个淋巴结肿大的和没有肿大的。
bd<-filter(bc,ln_yesno==0)####没有淋巴结肿大的患者数据
转换变量单位
bd$pyears<-bd$time/12
计算患者死亡人数和人年总数
bd$status<-factor(bd$status)
summary(bd$status)
sum(bd$pyears)
也可以用另一种方法计数,结果都是一样的
by_out<-group_by(bd,status)
summarise(by_out,count = n())
同理算出有淋巴结肿大的患者
be<-filter(bc,ln_yesno==1)
be$pyears<-be$time/12
be$status<-factor(be$status)
summary(be$status)
sum(be$pyears)
目前每个队列的发病率和人年总数我们都算出来了,可以进行计算了
a1=42
b1=3668.403
a2=30
b2=1054.636
c1=a1/b1*1000#c为千人年发病率
c2=a2/b2*1000# c为千人年发病率
所以有淋巴结肿大和没有淋巴结肿大的千人年(per 1,000 person-years)发病率为28千人年,11千人年。好的,下面我们来计算相对发病率及可信区间
RR=c1/c2###相对发病率
d=log(RR)
lower=exp(d - 1.96*sqrt((1/a1)+(1/a2)))####上限
upper=exp(d + 1.96*sqrt((1/a1)+(1/a2)))####下限
所以相对发病率为0.40,可信区间为0.25—0.64
也要同学要问有没有简单一点的,公式什么的太复杂,也是有的。popEpi包的rate_ratio函数可以帮我们快速计算
先导入popEp包
library(popEpi)
使用rate_ratio函数的时候,在分别填入两个队列的死亡人数和人年总数就可以了,我们这里X填入无淋巴结肿大组的,Y填入有淋巴结肿大组的
rate_ratio(x = c(42, 3668.403), y = c(30, 1054.636), SE.method = FALSE)
OK,两个结果一样,演示结束。