格式化生命表以用于生存分析
Formating life-tables to use in survival analysis
我正在尝试使用 R 中的 'relsurv' 包来比较队列的生存与国民生活 tables。下面的代码使用 relsurv 中的示例显示了我的问题,但更改了 life-table 数据。我刚刚在下面的 life-table 数据中使用了两年和两个年龄,实际数据要大得多但给出了相同的错误。错误是 'invalid ratetable argument' 但我已经按照示例 life-tables 'slopop' 和 'survexp.us'.
对其进行了格式化
library(survival)
library(relsurv)
data(rdata) # example data from relsurv
raw = read.table(header=T, stringsAsFactors = F, sep=' ', text='
Year Age sex qx
1980 30 1 0.00189
1980 31 1 0.00188
1981 30 1 0.00191
1981 31 1 0.00191
1980 30 2 0.00077
1980 31 2 0.00078
1981 30 2 0.00076
1981 31 2 0.00074
')
ages = c(30,40) # in years
years = c(1980, 1990)
rtab = array(data=NA, dim=c(length(ages), 2, length(years))) # set up blank array: ages, sexes, years
for (y in unique(raw$Year)){
for (s in 1:2){
rtab[ , s, y-min(years)+1] = -1 * log(1-subset(raw, Year==y&sex==s)$qx) / 365.24 # probability of death in next year, transformed to hazard (see ratetables help)
}
}
attributes(rtab)$dimnames[[1]] = as.character(ages)
attributes(rtab)$dimnames[[2]] = c('male','female')
attributes(rtab)$dimnames[[3]] = as.character(years)
attributes(rtab)$dimid <- c("age", "sex", 'year')
attributes(rtab)$dim <- c(length(ages), 2, length(years))
attributes(rtab)$factor = c(0,0,1)
attributes(rtab)$type = c(2,1,4)
attributes(rtab)$cutpoints[[1]] = ages*365.24 # must be in days
attributes(rtab)$cutpoints[[2]] = NULL
attributes(rtab)$cutpoints[[3]] = as.date(paste("1Jan", years, sep='')) # must be date
attributes(rtab)$class = "ratetable"
# example from relsurv
rsmul(Surv(time,cens) ~ sex+as.factor(agegr)+
ratetable(age=age*365.24, sex=sex, year=year),
data=rdata, ratetable=rtab, int=1)
尝试使用 relsurv 包中的转换函数重新格式化数据。那应该会给你一个兼容的数据集。
此致,
乔什
要补充三点:
您应该设置 attributes(rtab)$factor = c(0,1,0)
,因为性别(第二维)是一个因素(即不随时间改变)。
检查某个值是否有效的好方法 table 是使用 is.ratetable()
函数。 is.ratetable(rtab, verbose = TRUE)
甚至会 return 一条消息说明错误。
不使用verbose
先检查is.ratetable
的结果,因为它会谎言[=56] =]关于有效率tables.
此评论的其余部分是关于这个谎言的。
如果没有给出 type
属性,is.ratetable
将使用 factor
属性计算;你可以通过打印函数来看到这一点。但是,这样做似乎不正确。它使用 type <- 1 * (fac == 1) + 2 * (fac == 0) + 4 * (fac > 0)
,其中 fac
是 attributes(rtab)$factor
.
但是下一节检查 type
属性(如果提供)说唯一有效的值是 1
、2
、3
和 4
。从上面的代码中不可能得到1
。
例如,让我们检查一下 relsurv
包中提供的 slopop
速率 table。
library(relsurv)
data(slopop)
is.ratetable(slopop)
# [1] TRUE
is.ratetable(slopop, verbose = TRUE)
# [1] "wrong length for cutpoints 3"
我认为这是您的费率 table 被挂起的地方。
我正在尝试使用 R 中的 'relsurv' 包来比较队列的生存与国民生活 tables。下面的代码使用 relsurv 中的示例显示了我的问题,但更改了 life-table 数据。我刚刚在下面的 life-table 数据中使用了两年和两个年龄,实际数据要大得多但给出了相同的错误。错误是 'invalid ratetable argument' 但我已经按照示例 life-tables 'slopop' 和 'survexp.us'.
对其进行了格式化library(survival)
library(relsurv)
data(rdata) # example data from relsurv
raw = read.table(header=T, stringsAsFactors = F, sep=' ', text='
Year Age sex qx
1980 30 1 0.00189
1980 31 1 0.00188
1981 30 1 0.00191
1981 31 1 0.00191
1980 30 2 0.00077
1980 31 2 0.00078
1981 30 2 0.00076
1981 31 2 0.00074
')
ages = c(30,40) # in years
years = c(1980, 1990)
rtab = array(data=NA, dim=c(length(ages), 2, length(years))) # set up blank array: ages, sexes, years
for (y in unique(raw$Year)){
for (s in 1:2){
rtab[ , s, y-min(years)+1] = -1 * log(1-subset(raw, Year==y&sex==s)$qx) / 365.24 # probability of death in next year, transformed to hazard (see ratetables help)
}
}
attributes(rtab)$dimnames[[1]] = as.character(ages)
attributes(rtab)$dimnames[[2]] = c('male','female')
attributes(rtab)$dimnames[[3]] = as.character(years)
attributes(rtab)$dimid <- c("age", "sex", 'year')
attributes(rtab)$dim <- c(length(ages), 2, length(years))
attributes(rtab)$factor = c(0,0,1)
attributes(rtab)$type = c(2,1,4)
attributes(rtab)$cutpoints[[1]] = ages*365.24 # must be in days
attributes(rtab)$cutpoints[[2]] = NULL
attributes(rtab)$cutpoints[[3]] = as.date(paste("1Jan", years, sep='')) # must be date
attributes(rtab)$class = "ratetable"
# example from relsurv
rsmul(Surv(time,cens) ~ sex+as.factor(agegr)+
ratetable(age=age*365.24, sex=sex, year=year),
data=rdata, ratetable=rtab, int=1)
尝试使用 relsurv 包中的转换函数重新格式化数据。那应该会给你一个兼容的数据集。
此致, 乔什
要补充三点:
您应该设置
attributes(rtab)$factor = c(0,1,0)
,因为性别(第二维)是一个因素(即不随时间改变)。检查某个值是否有效的好方法 table 是使用
is.ratetable()
函数。is.ratetable(rtab, verbose = TRUE)
甚至会 return 一条消息说明错误。不使用
verbose
先检查is.ratetable
的结果,因为它会谎言[=56] =]关于有效率tables.
此评论的其余部分是关于这个谎言的。
如果没有给出 type
属性,is.ratetable
将使用 factor
属性计算;你可以通过打印函数来看到这一点。但是,这样做似乎不正确。它使用 type <- 1 * (fac == 1) + 2 * (fac == 0) + 4 * (fac > 0)
,其中 fac
是 attributes(rtab)$factor
.
但是下一节检查 type
属性(如果提供)说唯一有效的值是 1
、2
、3
和 4
。从上面的代码中不可能得到1
。
例如,让我们检查一下 relsurv
包中提供的 slopop
速率 table。
library(relsurv)
data(slopop)
is.ratetable(slopop)
# [1] TRUE
is.ratetable(slopop, verbose = TRUE)
# [1] "wrong length for cutpoints 3"
我认为这是您的费率 table 被挂起的地方。