4.4 Example: mortality rates following surgery
This is an illustration of what we seens previously but it also will add few important functions to help us analyze output of inla().
Surg## n r hospital
## 1 47 0 A
## 2 148 18 B
## 3 119 8 C
## 4 810 46 D
## 5 211 8 E
## 6 196 13 F
## 7 148 9 G
## 8 215 31 H
## 9 207 14 I
## 10 97 8 J
## 11 256 29 K
## 12 360 24 LHere we want to investigate mortality rate (r)
4.4.1 Model
\[Y_i \sim Binom(n_i, p_i), i = 1, .., 12\]
\[logit(pi) = \alpha + u_i, u_i \sim N(o, \sigma^2) \]
Reminder: logit is the log odds (\(ln(\frac{p_i}{1 -p_i}\))
Default prior for \(\alpha \sim N(0, \frac{1}{\tau})\)
We will PC for \(\u_i\) with
\(P(\sigma > 1) = 0.01\)
Define prior as a list:
prior.prec <- list(prec = list(prior = "pc.prec",
param = c(1, 0.01)))Define formula object:
formula <- r ~ f(hospital, model = "iid", hyper = prior.prec)Calling inla() :
res <- inla(formula,
data = Surg,
family = "binomial", Ntrials = n,
control.predictor = list(compute = TRUE), # marginal densities of LP
control.compute = list(dic = TRUE) # DIC
)4.4.2 resuts!
summary(res)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 0.489, Running = 0.151, Post = 0.0266, Total = 0.666
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -2.547 0.141 -2.842 -2.543 -2.279 -2.535 0
##
## Random effects:
## Name Model
## hospital IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for hospital 11.53 12.72 2.37 8.27 39.20 5.35
##
## Deviance Information Criterion (DIC) ...............: 74.44
## Deviance Information Criterion (DIC, saturated) ....: 24.68
## Effective number of parameters .....................: 8.10
##
## Marginal log-Likelihood: -41.17
## 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)')You can access specific part of the result using $
res$summary.fixed## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -2.547242 0.1412799 -2.842385 -2.542578 -2.279011 -2.534587
## kld
## (Intercept) 1.468744e-06We asked compute = TRUE so we get :
summary.linear.predictor: df with linear predictorssummary.fitted.values: df with fitted values obtained by transforming the linear predictors by the inverse of the link function (if like me you have trouble with log values)marginals.linear.predictor: list with the posterior marginals of the linear predictorsmarginals.fitted.values: list with the posterior marginals of the fitted values obtained by transforming the linear predictors by the inverse of the link function
You can acces them like that:
res$summary.fitted.values## mean sd 0.025quant 0.5quant 0.975quant
## fitted.Predictor.01 0.05544055 0.018915904 0.02257488 0.05432153 0.09593327
## fitted.Predictor.02 0.10175273 0.021235098 0.06719221 0.09933052 0.14989450
## fitted.Predictor.03 0.07108845 0.017212251 0.04194181 0.06954239 0.10956339
## fitted.Predictor.04 0.05952862 0.007868087 0.04529022 0.05912580 0.07607338
## fitted.Predictor.05 0.05288174 0.013028004 0.03006237 0.05204160 0.08060505
## fitted.Predictor.06 0.06957032 0.014574997 0.04438961 0.06842460 0.10155925
## fitted.Predictor.07 0.06725087 0.015681266 0.04027101 0.06601093 0.10175408
## fitted.Predictor.08 0.12138342 0.021672558 0.08401089 0.11969918 0.16844666
## fitted.Predictor.09 0.07025402 0.014379546 0.04539382 0.06912368 0.10180192
## fitted.Predictor.10 0.07838144 0.019434425 0.04644487 0.07624236 0.12287168
## fitted.Predictor.11 0.10124358 0.017156213 0.07200569 0.09974523 0.13894767
## fitted.Predictor.12 0.06873049 0.011630975 0.04819778 0.06795637 0.09377516
## mode
## fitted.Predictor.01 0.05286946
## fitted.Predictor.02 0.09441323
## fitted.Predictor.03 0.06705311
## fitted.Predictor.04 0.05832318
## fitted.Predictor.05 0.05030748
## fitted.Predictor.06 0.06643517
## fitted.Predictor.07 0.06396695
## fitted.Predictor.08 0.11636261
## fitted.Predictor.09 0.06714129
## fitted.Predictor.10 0.07270945
## fitted.Predictor.11 0.09671112
## fitted.Predictor.12 0.06651510Or using r-inla functions to manipulate them:
inla.L.marginal()L is for “Letters”:- e for expectations
- q for quantiles
- s for spline smoothings
- t for transform (if you want varianec instead of precision)
- z for summary statistics
library(ggplot2)
alpha <- res$marginals.fixed[[1]]
ggplot(data.frame(inla.smarginal(alpha)), aes(x, y)) +
geom_line() +
theme_bw()
marg.variance <- inla.tmarginal(function(x) 1/x,
res$marginals.hyperpar$"Precision for hospital")
ggplot(data.frame(inla.smarginal(marg.variance)), aes(x, y)) +
geom_line() +
theme_bw()
list_marginals <- res$marginals.fitted.values
marginals <- data.frame(do.call(rbind, list_marginals))
marginals$hospital <- rep(names(list_marginals),
times = sapply(list_marginals, nrow))
library(ggplot2)
ggplot(marginals, aes(x = x, y = y)) + geom_line() +
facet_wrap(~ hospital) +
labs(x = "", y = "Density") +
geom_vline(xintercept = 0.1, col = "gray") +
theme_bw()