coxph.RdTools to convert a Cox proportional hazard model into Poisson regression
inla.coxph(formula, data, control.hazard = list(), debug = FALSE, tag = "")
inla.rbind.data.frames(...)The formula for the coxph model where the response must be a
inla.surv-object.
All the data used in the formula, as a list.
Control the model for the baseline-hazard; see
?control.hazard.
Print debug-information
An optional tag added to the names of the new variables created
(to make them unique when combined with several calls of inla.coxph.
Note that E..coxph is not included, as its usually merged into one
vector over different expansions.
Data.frames to be rbind-ed, padding with NA.
inla.coxph returns a list of new expanded variables to be
used in the inla-call. Note that element data and
data.list needs to be merged into a list to be passed as the
data argument. See the example for details.
inla.rbind.data.frames returns the rbinded data.frames padded with
NAs. There is a better implementation in dplyr::bind_rows, which is
used if package dplyr is installed.
## How the cbind.data.frames works:
df1 = data.frame(x=1:2, y=2:3, z=3:4)
df2 = data.frame(x=3:4, yy=4:5, zz=5:6)
inla.rbind.data.frames(df1, df2)
#> x y z yy zz
#> 1 1 2 3 NA NA
#> 2 2 3 4 NA NA
#> 3 3 NA NA 4 5
#> 4 4 NA NA 5 6
## Standard example of how to convert a coxph into a Poisson regression
n = 1000
x = runif(n)
lambda = exp(1+x)
y = rexp(n, rate=lambda)
event = rep(1,n)
data = list(y=y, event=event, x=x)
y.surv = inla.surv(y, event)
intercept1 = rep(1, n)
p = inla.coxph(y.surv ~ -1 + intercept1 + x,
list(y.surv = y.surv, x=x, intercept1 = intercept1))
r = inla(p$formula,
family = p$family,
data=c(as.list(p$data), p$data.list),
E = p$E)
summary(r)
#> Time used:
#> Pre = 0.235, Running = 0.223, Post = 0.0235, Total = 0.482
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> intercept1 1.027 0.065 0.901 1.027 1.155 1.027 0
#> x 0.995 0.108 0.782 0.995 1.207 0.995 0
#>
#> Random effects:
#> Name Model
#> baseline.hazard RW1 model
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant 0.975quant
#> Precision for baseline.hazard 19185.28 19814.77 371.38 13125.80 73080.39
#> mode
#> Precision for baseline.hazard 167.19
#>
#> Marginal log-Likelihood: -2738.12
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
#>
## How to use this in a joint model
intercept2 = rep(1, n)
y = 1 + x + rnorm(n, sd=0.1)
df = data.frame(intercept2, x, y)
## new need to cbind the data.frames, and then add the list-part of
## the data
df.joint = c(as.list(inla.rbind.data.frames(p$data, df)), p$data.list)
df.joint$Y = cbind(df.joint$y..coxph, df.joint$y)
## merge the formulas, recall to add '-1' and to use the new joint
## reponse 'Y'
formula = update(p$formula, Y ~ intercept2 -1 + .)
rr = inla(formula,
family = c(p$family, "gaussian"),
data = df.joint,
E = df.joint$E..coxph)
## A check that automatic and manual approach gives the same result
data(Leuk)
## add some random stuff for testing. Note that variables needs to
## be in 'data' as they are expanded
Leuk$off <- runif(nrow(Leuk), min = -0.5, max = 0.5)
Leuk$off.form <- runif(nrow(Leuk), min = -0.5, max = 0.5)
Leuk$w <- runif(nrow(Leuk), min = 0.5, max = 1.0)
formula <- inla.surv(time, cens) ~ sex + age +
wbc + tpi + offset(off.form)
r <- inla(formula, family = "coxph", data = Leuk,
offset = off, weights = w)
cph <- inla.coxph(formula = formula, data = Leuk)
cph.data = c(as.list(cph$data), cph$data.list)
rr <- inla(cph$formula, family = cph$family, data = cph.data,
E = cph$E, offset = off, weights = w)
print(cbind(r$mlik, rr$mlik, r$mlik - rr$mlik))
#> [,1] [,2] [,3]
#> log marginal-likelihood (integration) -1187.672 -2005.293 817.6212
#> log marginal-likelihood (Gaussian) -1188.052 -2005.670 817.6179