鸿鹄网站建设,做网站开封,门店库存管理软件,中国最权威的网站排名分析队列人群有无死亡人数超额#xff0c;通常应用标准人群死亡率来校正#xff0c;即刻观察到中的实际死亡数#xff08;D#xff09;与定一个标准的死亡人数#xff08;E#xff09;#xff0c;D与E之比称为死亡比#xff08;standarized Mortality ratio#xff0c…分析队列人群有无死亡人数超额通常应用标准人群死亡率来校正即刻观察到中的实际死亡数D与定一个标准的死亡人数ED与E之比称为死亡比standarized Mortality ratioSMR. 标准化死亡率 (SMR) 是观察到的病例与预期病例的比率。 今天我们介绍一下怎么通过手动计算校正年龄、性别后的标准化死亡率 (SMR)。 我们先导入R包和数据也不能算完全手动还是需要用到survival包。
library(survival)
bc-read.csv(E:/r/test/smr1.csv,sep,,headerTRUE)
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为46E为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两者算得一模一样。
