分析队列人群有无死亡人数超额，通常应用标准人群死亡率来校正，即刻观察到中的实际死亡数（D）与定一个标准的死亡人数（E），D与E之比称为死亡比（standarized Mortality ratio，SMR）. 标准化死亡率 (SMR) 是观察到的病例与预期病例的比率。

今天我们介绍一下怎么通过手动计算校正年龄、性别后的标准化死亡率 (SMR)。 我们先导入R包和数据，也不能算完全手动，还是需要用到survival包。

library(survival)bc<-read.csv("E:/r/test/smr1.csv",sep=',',header=TRUE)head(bc,6)## sex age entry_date status futime## 1 male 11213.45 1985/2/26 1 2750.988## 2 female 6681.62 1994/3/1 0 1981.679## 3 male 10411.76 1992/6/6 0 2979.385## 4 female 10665.13 1985/9/22 0 2576.967## 5 male 19065.91 1986/11/2 0 4993.524## 6 male 10154.70 1993/4/23 0 1821.558

这是一个很简单的数据，sex为性别，age是年龄，entry_date为诊断也就是进入这个队列的时间，status为结局变量，futime为生存时间。（公众号回复：SMR1可以获得该数据） 先处理一下生存时间

bc$entry_date<-as.Date(bc$entry_date)

我们首先要算出它的haz， 每人年的平均人口死亡率（d/(pyrs)，其中d是死亡人数，pyrs是人年），可以通过survival包的survexp函数计算。sex这个指标不能少，如果全是女的也要设置.

bc$risk1 <- -log(survexp(futime ~ 1, data = bc, rmap = list(year = entry_date, age = age, sex = sex), cohort = FALSE, conditional = TRUE))

survexp函数中选项还有一个费率表，就是用来调整生成的预测概率的，可以理解为校正功能。 survival函数自带了3个survexp.us, survexp.usr和survexp.mn. survexp.usr就是1940年至2014年按年龄、性别和种族划分的美国人口。这里有黑人和白人的数据集，我这里只取白人的

survexp.uswhite <- survexp.usr[,,"white",]head(survexp.uswhite ,6)## Rate table with dimension(s): age sex year ## , , year = 1940## ## sex## age male female## 0 1.350207e-04 1.057536e-04## 1 1.336591e-05 1.185314e-05## 2 7.264935e-06 6.029907e-06## 3 5.206865e-06 4.411492e-06## 4 4.192119e-06 3.506694e-06## 5 3.780843e-06 3.013293e-06## ## , , year = 1941## ## sex## age male female## 0 1.300527e-04 1.017030e-04## 1 1.260761e-05 1.118577e-05## 2 6.919267e-06 5.733727e-06## 3 4.976544e-06 4.208639e-06## 4 4.019425e-06 3.345001e-06## 5 3.624615e-06 2.879024e-06## ## , , year = 1942##

如果你想算某个地方的生存率，可以使用当地人口普查数据，这里我们把费率表加进去。

bc$risk2 <- -log(survexp(futime ~ 1, data = bc, rmap = list(year = entry_date, age = age, sex = sex), cohort = FALSE, ratetable = survexp.uswhite , conditional = TRUE))head(bc,6)## sex age entry_date status futime risk1 risk2## 1 male 11213.45 1985-02-26 1 2750.988 0.017979337 0.015404330## 2 female 6681.62 1994-03-01 0 1981.679 0.002648331 0.002482978## 3 male 10411.76 1992-06-06 0 2979.385 0.015381638 0.013551117## 4 female 10665.13 1985-09-22 0 2576.967 0.006089193 0.005058961## 5 male 19065.91 1986-11-02 0 4993.524 0.186826509 0.175753223## 6 male 10154.70 1993-04-23 0 1821.558 0.008510061 0.007464535

可以看到，加入费率表的概率和不加的是不一样的。算出了预测概率后我们就可以进一步计算了。

O <- sum(bc$status)E <- sum(bc$risk2)O;## [1] 46E## [1] 6.74174

可以得到O为46，E为6.74.SMR等于O/E

SMR <- O/ESMR## [1] 6.823164

接下来计算可信区间，先设置一下alpha

alpha = 0.05

接下来计算可信区间，公式是固定的，直接放进去就可以了

SMR.lo <- O/E * (1 - 1/9/O - qnorm(1 - alpha/2)/3/sqrt(O))^3SMR.up <- (O + 1)/E * (1 - 1/9/(O + 1) + qnorm(1 - alpha/2)/3/sqrt(O + 1))^3SMR.lo## [1] 4.994967SMR.up## [1] 9.101355

这样全部结果就计算出来啦，计算结果我们使用survexp.fr包来验证一下

library(survexp.fr)attach(bc)bc$entry_date<-as.Date(bc$entry_date)SMR(futime, status, age, sex, entry_date,ratetable =survexp.uswhite)## $O## [1] 46## ## $E## [1] 6.74174## ## $SMR.classic## $SMR.classic$SMR## [1] 6.823164## ## $SMR.classic$SMR.lo## [1] 4.994967## ## $SMR.classic$SMR.up## [1] 9.101355## ## $SMR.classic$p.value## [1] 0## ## ## $SMR.poisson## $SMR.poisson$SMR## [1] 6.823164## ## $SMR.poisson$SMR.lo## [1] 5.110733## ## $SMR.poisson$SMR.up## [1] 9.109373## ## $SMR.poisson$p.value## [1] 8.903473e-39

两者算得一模一样。