the data
data(Kidney, package='INLA')
head(Kidney)
## time event age sex disease dis1 dis2 dis3 ID
## 1 8 1 28 0 1 1 0 0 1
## 2 16 1 28 0 1 1 0 0 1
## 3 23 1 48 1 2 0 1 0 2
## 4 13 0 48 1 2 0 1 0 2
## 5 22 1 32 0 1 1 0 0 3
## 6 28 1 32 0 1 1 0 0 3
scale the response variable
Kidney$years = Kidney$time / 365.25
run the model
library("INLA")
## Loading required package: sp
## Loading required package: Matrix
## This is INLA 0.0-1485844051, dated 2017-01-31 (09:14:12+0300).
## See www.r-inla.org/contact-us for how to get help.
##
## Attaching package: 'INLA'
## The following object is masked _by_ '.GlobalEnv':
##
## Kidney
mod = inla(
inla.surv(years, event) ~ age * sex +
f(ID, model="iid", prior='pc.prec', param=c(0.5, 0.1)),
family="weibullcure",
control.family = list(hyper=list(
'log alpha' = list(
prior='loggamma', param=c(1,1)),
'logit probability' = list(
prior='logitbeta', param=c(1,1))
)),
data=Kidney,
control.inla=list(h=0.001),
control.mode = list(theta=c(0,-5,4), restart=TRUE),
verbose=TRUE)
results
knitr::kable(mod$summary.fixed, digits=3)
mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | kld | |
---|---|---|---|---|---|---|---|
(Intercept) | 3.654 | 1.078 | 1.529 | 3.651 | 5.800 | 3.644 | 0 |
age | -0.027 | 0.023 | -0.071 | -0.028 | 0.019 | -0.029 | 0 |
sex | -3.418 | 1.245 | -5.918 | -3.404 | -0.994 | -3.379 | 0 |
age:sex | 0.044 | 0.027 | -0.009 | 0.044 | 0.097 | 0.044 | 0 |
knitr::kable(mod$summary.hyper, digits=3)
mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | |
---|---|---|---|---|---|---|
alpha parameter for ps | 1.130 | 0.149 | 0.842 | 1.130 | 1.424 | 1.141 |
p parameter for ps | 0.015 | 0.014 | 0.001 | 0.011 | 0.052 | 0.002 |
Precision for ID | 4.649 | 4.943 | 0.922 | 3.154 | 17.338 | 1.803 |
standard deviation
modSd = Pmisc::priorPost(mod)
knitr::kable(modSd$summary, digits=2)
mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | |
---|---|---|---|---|---|---|
alpha for weibullcure | 1.13 | 0.15 | 0.84 | 1.13 | 1.42 | 1.14 |
prob for weibullcure | 0.01 | 0.01 | 0.00 | 0.01 | 0.05 | 0.00 |
sd for ID | 0.58 | 0.21 | 0.24 | 0.56 | 1.04 | 0.74 |
do.call(matplot, modSd$sd$matplot)
do.call(legend, modSd$legend)
do.call(matplot, modSd$alpha$matplot)
do.call(legend, modSd$legend)
do.call(matplot, modSd$prob$matplot)
do.call(legend, modSd$legend)