In Pmisc: Various Utilities for knitr and inla

Weibull With Cure Fraction

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="weibullsurv", 
  control.family = list(hyper=list(
      'log alpha' = list(
        prior='loggamma', param=c(1,1)))
    ),
  data=Kidney, 
  control.inla=list(h=0.001, int.strategy='grid'), 
  control.compute = list(config=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)

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')

Gamma

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')]

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
  )

standard deviation

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')