R语言进行人年发病率可信区间计算（per 1,000 person-years）

有粉丝在后台问我人年发病率的可信区间怎么计算？如下图：


我对人年也不是很熟悉，于是查一下资料，

人年发病率的计算，就是每个进入队列的患者的随访时间求和，就是总的人年数，然后用发生结局的人数/总的人年数*1000就是per1000 person-years，人年的发病率置信区间是通过泊松回归分布进行计算的，计算原理如下图

有了公式计算就容易多了，主要步骤为：

1. 先求出每个队列的发病数和人年总数
2. 发病数/人年总数×1000即为人年发病率
3. 求出相对发病率RR，
4. 求出R的对数log(RR)
5. 代入最后的公式计算就可以求出置信区间了

本章我们通过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，两个结果一样，演示结束。