the data
data(Kidney, package='INLA')
head(Kidney)
scale the response variable
Kidney$years = Kidney$time / 365.25
run the model
library("INLA")
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, int.strategy='grid'),
control.compute = list(config=TRUE),
control.mode = list(theta=c(0,-5,4), restart=TRUE),
verbose=FALSE)
results
knitr::kable(mod$summary.fixed, digits=3)
knitr::kable(mod$summary.hyper, digits=3)
standard deviation
modSd = Pmisc::priorPost(mod)
knitr::kable(modSd$summary, digits=3)
for(D in grep('^sd', names(modSd), value=TRUE)) {
do.call(matplot, modSd[[D]]$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)
xSeq = seq(0,1/12,len=1000)
myDens = Pmisc::sampleDensHaz(mod, xSeq, n=25)
matplot(xSeq*365.25, myDens[,'dens',]/365.25, type='l',
col = '#00000020', lty=1, xlab='days', ylab='dens')
data('meuse', package='sp')
distBreaks = seq(0, 1, len=101)
meuse$distCut = cut(meuse$dist, breaks=distBreaks,
right=FALSE)
meuse$soilFac = factor(meuse$soil, levels=1:3,
labels=c('Calcerous','Non-Calcerous','Red Brick'))
meuse[1:4,c('cadmium','elev','dist','distCut', 'soilFac')]
## cadmium elev dist distCut soilFac
## 1 11.7 7.909 0.00135803 [0,0.01) Calcerous
## 2 8.6 6.983 0.01222430 [0.01,0.02) Calcerous
## 3 6.5 7.800 0.10302900 [0.1,0.11) Calcerous
## 4 2.6 7.655 0.19009400 [0.19,0.2) Non-Calcerous
library("INLA", quietly=TRUE)
mod = inla(
cadmium ~ elev + soilFac +
f(distCut, model='rw2',
prior = 'pc.prec', param = c(0.01, 0.5),
values = levels(distCut), constr=FALSE,
extraconstr = list(
A = model.matrix(~0+cut(0.25, distBreaks)), e=0
)),
data = meuse, family = 'gamma',
control.family = list(hyper=list(
prec=list(prior = 'pc.prec', param = c(0.1, 0.5)))),
verbose=FALSE
)
knitr::kable(mod$summary.fixed, digits=3)
knitr::kable(mod$summary.hyper, digits=3)
standard deviation
modSd = Pmisc::priorPost(mod)
knitr::kable(modSd$summary, digits=3)
for(D in grep('^sd', names(modSd), value=TRUE)) {
do.call(matplot, modSd[[D]]$matplot)
do.call(legend, modSd$legend)
}
mod$summary.random$distCut[is.na(mod$summary.random$distCut$mean),-1] = 0
matplot(distBreaks[-length(distBreaks)],
exp(mod$summary.random$distCut[,paste0(c('0.5','0.025','0.975'), 'quant')]),
ylab = 'relative rate', xlab='distance', type='l', lty=c(1,2,2),
col='black', log='y')
legend("bottomleft", lty=c(1,2), legend=c('median','95% CI'), bty='n')