2016-07-01 24 views
0

我試圖在R中使用'relsurv'軟件包來比較隊列與國家生命表的生存期。下面的代碼使用relsurv中的示例顯示了我的問題,但更改了生命表數據。我剛剛在下面的生命表數據中使用了兩年和兩年的時間,實際數據要大得多,但也會出現相同的錯誤。該錯誤是'無效的速率表參數',但我已經按照示例生命表'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) 

回答

0

嘗試使用relsurv包中的transrate函數重新格式化數據。這應該給你一個兼容的數據集。

問候, 喬希

0

三件事補充:

  1. 您應該設置attributes(rtab)$factor = c(0,1,0),因爲性(第二維)是一個因素(即不隨時間變化)。

  2. 檢查是否有效的費率表的一個好方法是使用is.ratetable()函數。 is.ratetable(rtab, verbose = TRUE)甚至會返回一條消息,指出錯誤。

  3. 檢查的is.ratetable結果,而無需使用verbose第一,因爲它會謊言有關有效的費率表。

該評論的其餘部分是關於這個謊言。

如果type屬性沒有給出,is.ratetable將使用factor屬性計算它;你可以通過打印該功能來看到這一點。但是,它似乎做錯了。它使用type <- 1 * (fac == 1) + 2 * (fac == 0) + 4 * (fac > 0),其中facattributes(rtab)$factor

但接下來的部分,它檢查是否商提供的它type屬性,說,唯一有效的值是1234。從上面的代碼中獲取1是不可能的。

例如,我們來看一下relsurv包裝提供的slopop率表。

library(relsurv) 
data(slopop) 
is.ratetable(slopop) 
# [1] TRUE 

is.ratetable(slopop, verbose = TRUE) 
# [1] "wrong length for cutpoints 3" 

我認爲這是您的費率表掛起的位置。