手动计算校正年龄、性别后的标准化死亡率 (SMR)

分析队列人群有无死亡人数超额，通常应用标准人群死亡率来校正，即刻观察到中的实际死亡数（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] 46
E
## [1] 6.74174

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

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

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

alpha = 0.05

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

SMR.lo <- O/E * (1 - 1/9/O - qnorm(1 - alpha/2)/3/sqrt(O))^3
SMR.up <- (O + 1)/E * (1 - 1/9/(O + 1) + qnorm(1 - alpha/2)/3/sqrt(O + 1))^3
SMR.lo
## [1] 4.994967
SMR.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

两者算得一模一样。