therneau
commented
2 years ago There is no function pyears3B, nor is there a call to such a function. There is a pyears3b, however. I would suggest looking at the tests/expected.R, which hammers through an expected survival computation in excruciating detail. The function computes the expected results by hand, and then compares them to the survexp results as a test of survexp. Also tests/expected2.R, and tests/ratetable.R may be helpful.
The parallel step is to look at noweb/code.pdf, which is a long document explaining much of the mathematics behind the survival package. Read the section on expected survival.
StephenKay67
I cannot replicate by manual programming, expected survival results from the survfit function , and have been unable to uncover the underlying program code: on examining the survfit code, the real workhorse is a call (.Call) to the C function "cpyears3B" (at least for the "individual.h" survexp method I am trying to replicate) and there is no uncompiled version called "cpyears3B.c" in the relevant src folder (or able to find googling - likely this function is held in a .c (or .h) file with another name?). I can replicate one year predictions exactly (unsurprising) but not two years for an individual:
library(survival) library(date) df<- data.frame(t = c(1,2) 365.25, # two patient rows which differ in 1or2yrs predict age = c(70,70) 365.25, sex = c("male","male"), year = c(as.date("01Jan1940"),as.date("01Jan1940")) ) official_1y <- survexp(t ~ 1, data=df, scale = 365.25, method="individual.h") # i.e. using default "survexp.us" ratetable
manual_1y <- survexp.us[age = "70", sex="male", year = "1940"] * 365.25 all.equal(official_1y[1], manual_1y,check.attributes = FALSE) # TRUE
2year assuming "period" rate table hence add 1 year to age and year for second year
manual_2y_a <- (survexp.us[age = "70", sex="male", year = "1940"] + survexp.us[age = "71", sex="male", year = "1941"]) * 365.25 all.equal(official_1y[2], manual_2y_a,check.attributes = FALSE) # FALSE
2year assuming "cohort" rate table hence add 1 year to age only and stick with 1940
manual_2y_b <- (survexp.us[age = "70", sex="male", year = "1940"] + survexp.us[age = "71", sex="male", year = "1940"]) * 365.25 all.equal(official_1y[2], manual_2y_b,check.attributes = FALSE) # FALSE but much closer
My interpretation is that survexp.us is a "period" not a "cohort" ratetable hence manual_2y_a should be correct. Clearly though either something is wrong with my reasoning/calculations or the formula underlying survexp is using some approximation formula that I am not aware of. I cannot replicate results outside of R without knowing this formula (possibly something is amiss with my formula/assumptions)? Thanks for any help.