Description Usage Arguments Details Note Author(s) References Examples
This function tranform survival data into a format compatible with
the glm()
function for fitting an auxiliary Poisson model,
providing the parameter estimates of the associated proportional hazard model.
1 2 3 4 5 |
data |
a data frame with columns:
|
all.breaks |
the breakpoints between time intervals |
interval.width |
the width of the time intervals on which the risks will be assumed consant,
in case of intervals of the same length.
This parameter is ignored if |
nInts |
the number of intervals containing the same expected number of events
(used only if |
factors |
a vector of characters, containing the names of the factors to be kept in the transformed data set |
compress |
a logical, indicating whether the record with the same factor profile should be summarized into one record, i.e. whether the data should be expressed in a short form |
x |
The fitted Poisson model on the poissonized data |
type |
the type of plot, either 'haz' for the hazard function or 'Surv', for the survival curve |
add |
should the plot added to the active device? |
xscale |
scaling factor for the time (x) axis |
by |
covariate for which a different curve per level has to be plotted |
col, ... |
other graphical parameters |
If interval.width
is not null, the study period is diveded into
equal-length intervals of length interval.width
.
Otherwise, nInts
intervals are used, and the location of their bounds
is computed based on the empirical quantiles of the survival function.
This code is hugely inspired by original code made publicly available by Stephanie Kovalchik [web link]
Federico Rotolo [aut, cre], Xavier Paoletti [ctr], Marc Buyse [ctr], Tomasz Burzykowski [ctr], Stefan Michiels [ctr]
Whitehead, J. Fitting Cox's regression model to survival data using GLIM. J Roy Stat Soc C Appl Stat 1980; 29(3):268-275. http://www.jstor.org/stable/2346901.
Crowther MJ, Riley RD, Staessen JA, Wang J, Gueyffier F, Lambert PC. Individual patient data meta-analysis of survival data using Poisson regression models. BMC Medical Research Methodology 2012; 12:34. doi: 10.1186/1471-2288-12-34.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 | ################################################################################
# Example 1 - KIDNEY data #
################################################################################
library(survival)
data(kidney)
kidney <- kidney[1:(nrow(kidney)/2)*2,]
head(kidney)
par(mfrow=c(1, 3))
for (int in c(50, 20, 10)) {
head(wdata1 <- poissonize(kidney, interval.width = int,
factors = c('disease'), compress = FALSE))
head(wdata2 <- poissonize(kidney, interval.width = int,
factors = c('disease'), compress = TRUE))
fitcox <- (coxph(Surv(time, status) ~ disease, data = kidney))
fitpoi1 <- glm(event ~ -1 + interval + disease + offset(log(time)),
data = wdata1, family = 'poisson')
fitpoi2 <- glm(m ~ -1 + interval + offset(log(Rt)) + disease,
data = wdata2, family = 'poisson')
cox.base <- basehaz(fitcox, centered = FALSE)
plot(stepfun(cox.base$time[-nrow(cox.base)], exp(-cox.base$hazard)),
ylim = 0:1, xlim = c(0, max(cox.base$time)),
do.points = FALSE, verticals = FALSE, xaxs = 'i',
main = paste0('KIDNEY data set\nInterval width = ', int),
xlab = 'Time', ylab = 'Survival probability')
plotsson(fitpoi1, 'Surv', add = TRUE, col = 2, lty = 2)
plotsson(fitpoi2, 'Surv', add = TRUE, col = 3, lty = 3)
legend('topright', col = 1:3, lty = 1:3,
legend = c('Breslow (Cox)', 'Poisson',
'Poisson (compressed dataset)'))
}
print(cbind(Cox = coef(fitcox),
Poisson = rev(rev(coef(fitpoi1))[1:3]),
Poisson_Compressed = rev(rev(coef(fitpoi2))[1:3])), digits = 2)
################################################################################
# Example 2 - COLON data #
################################################################################
library(survival)
data(colon)
head(wdata1 <- poissonize(subset(colon, etype == 1), interval.width = 365.25,
factors=c('surg', 'sex', 'age'), compress = FALSE))
head(wdata2 <- poissonize(subset(colon, etype == 1), interval.width = 365.25,
factors=c('surg', 'sex', 'age'), compress = TRUE))
fitcox <- coxph(Surv(time, status) ~ surg + sex + age,
data = subset(colon, etype == 1))
system.time({
fitpoi1 <- glm(event ~ -1 + interval + surg + sex + age + offset(log(time)),
data = wdata1, fam = 'poisson')
})
system.time({
fitpoi2 <- glm(m ~ -1 + interval + offset(log(Rt)) + surg + sex + age,
data = wdata2, family = 'poisson')
})
{
cox.base <- basehaz(fitcox, centered = FALSE)
par(mfrow = c(1, 1))
plot(stepfun(cox.base$time[-nrow(cox.base)], exp(-cox.base$hazard)),
ylim = 0:1, xlim = c(0, max(cox.base$time)),
do.points = FALSE, verticals = FALSE, xaxs = 'i',
main = 'COLON data set', xlab = 'Time', ylab = 'Survival probability')
plotsson(fitpoi1, 'Surv', add = TRUE, col = 2, lty = 2)
plotsson(fitpoi2, 'Surv', add = TRUE, col = 3, lty = 3)
legend('topright', col = 1:3, lty = 1:3,
legend = c('Cox', 'Poisson', 'Poisson (compressed dataset)'))
}
print(cbind(Cox = coef(fitcox),
Poisson = rev(rev(coef(fitpoi1))[1:3]),
Poisson_Compressed = rev(rev(coef(fitpoi2))[1:3])), digits = 2)
################################################################################
# Example 3 - LUNG data #
################################################################################
library(survival)
data(lung)
lung$status <- lung$status - 1
lung$id <- 1:nrow(lung)
head(wdata1 <- poissonize(lung, interval.width = 365.25/12,
factors = c('pat.karno', 'sex', 'age'),
compress = FALSE))
head(wdata2 <- poissonize(lung, interval.width = 365.25/12,
factors = c('pat.karno', 'sex', 'age'),
compress = TRUE))
fitcox <- coxph(Surv(time, status) ~ pat.karno + sex + age, data = lung)
system.time({
fitpoi1 <- glm(event ~ -1 + interval + pat.karno + sex + age +
offset(log(time)),
data = wdata1, family = 'poisson')
})
system.time({
fitpoi2 <- glm(m ~ -1 + interval + pat.karno + sex + age + offset(log(Rt)),
data = wdata2, family = 'poisson')
})
{
cox.base <- basehaz(fitcox, centered = FALSE)
plot(stepfun(cox.base$time[-nrow(cox.base)], exp(-cox.base$hazard)),
ylim = 0:1, xlim = c(0, max(cox.base$time)),
do.points = FALSE, verticals = FALSE, xaxs = 'i',
main = 'LUNG data set', xlab = 'Time', ylab = 'Survival probability')
plotsson(fitpoi1, 'Surv', add = TRUE, col = 2, lty = 2)
plotsson(fitpoi2, 'Surv', add = TRUE, col = 3, lty = 3)
legend('topright', col = 1:3, lty = 1:3,
legend = c('Cox', 'Poisson', 'Poisson (compressed dataset)'))
}
print(cbind(Cox = coef(fitcox),
Poisson = rev(rev(coef(fitpoi1))[1:3]),
Poisson_Compressed = rev(rev(coef(fitpoi2))[1:3])), digits = 2)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.