- Trauma Data Example Again
- Parameter Estimation with MCSE
- Reference Priors for Binomial Regression
- O-Ring Data Presentations

- Trauma Data Example Again
- Parameter Estimation with MCSE
- Reference Priors for Binomial Regression
- O-Ring Data Presentations

Dr. Turner Osler, a trauma surgeon and former head of the Burn Unit at the University of New Mexico Trauma Center, provided data on survival of patients admitted to the University of New Mexico Trauma Center between 1991 and 1994.

- The predictor variables are injury severity score (ISS), revised trauma score (RTS), AGE, and type of injury (TI). TI is either blunt (\(TI = 0\)), e.g., the result of a car crash, or penetrating (\(TI = 1\)), e.g., gunshot or stab wounds.
- ISS is an overall index of a patient's injuries based on the approximately 1,300 injuries cataloged in the Abbreviated Injury Scale. Takes on values from 0 for a patient with no injuries to 75 for a patient with severe injuries in three or more body areas.
- RTS is an index of physiologic status derived from the Glasgow Coma Scale. It is a weighted average of a number of measurements on an incoming patient such as systolic blood pressure and respiratory rate. Takes on values from 0 for a patient with no vital signs to 7.84 for a patient with normal vital signs.

- Data available here.

# trauma <- read.table("http://www.ics.uci.edu/~wjohnson/BIDA/Ch8/trauma300.txt, # yheader=TRUE) traumadata <- read.table("trauma300.txt", header=TRUE) attach(traumadata) head(traumadata)

## ID death ISS TI RTS AGE ## 1 2979 0 1 1 7.8408 25 ## 2 1167 0 9 0 7.8408 63 ## 3 116 0 29 0 2.9304 32 ## 4 1118 0 10 0 7.8408 61 ## 5 3126 0 2 1 7.8408 43 ## 6 2173 0 9 0 7.8408 27

Dr. Osler proposed a logistic regression model to estimate the probability of a patient's death using an intercept; the predictors ISS, RTS, AGE, and TI; and an interaction between AGE and TI. AGE is viewed as a surrogate for a patient's physiologic reserve, i.e., their ability to withstand trauma.

The model \[ \text{logit}(\theta_i) = \beta_1 + \beta_2 ISS_i + \beta_3 RTS_i + \beta_4 AGE_i + \beta_5 TI_i + \beta_6 (AGE*TI)_i \] is similar to logistic models used by trauma centers throughout the United States.

Introduce \(x_{i6} = (AGE*TI)_i\) and correct for the mean values of all the continuous predictor variables, then the model becomes \[ \text{logit}(\theta_i) = \alpha + \beta_2 x_{i2} + \beta_3 x_{i3} + \beta_4 x_{i4} + \beta_5 x_{i5} + \beta_6 x_{i6}. \]

x.check <- cbind(ISS-mean(ISS), RTS-mean(RTS), AGE-mean(AGE), TI, AGE*TI-mean(AGE*TI)) colnames(x.check) <- c("x2", "x3", "x4", "x5", "x6") head(x.check)

## x2 x3 x4 x5 x6 ## [1,] -13.28 0.553934 -6.41 1 17.633333 ## [2,] -5.28 0.553934 31.59 0 -7.366667 ## [3,] 14.72 -4.356466 0.59 0 -7.366667 ## [4,] -4.28 0.553934 29.59 0 -7.366667 ## [5,] -12.28 0.553934 11.59 1 35.633333 ## [6,] -5.28 0.553934 -4.41 0 -7.366667

model2 <- glm(death ~ x.check, family = binomial(), x = TRUE)

- For our Bayesian analysis we assume the same data model as the frequentist, i.e.Â we continue to use the logit link function.
- We assume the prior distribution of the parameters (the regression coefficients) makes them independent and identically normally distributed with mean 0 and variance \(\sigma^2 = 25\).

- The log unnormalized posterior density calculated by the following R function.

x <- model2$x y <- model2$y lupost <- function(beta, x, y) { eta <- as.numeric(x %*% beta) logp <- ifelse(eta < 0, eta - log1p(exp(eta)), - log1p(exp(- eta))) logq <- ifelse(eta < 0, - log1p(exp(eta)), - eta - log1p(exp(- eta))) logl <- sum(logp[y == 1]) + sum(logq[y == 0]) return(logl - sum(beta^2) / 50) }

- Use
`mcmc`

package and`metrop`

function to obtain posterior simulations. - Following code runs the Metropolis algorithm to simulate the posterior.

library(mcmc) set.seed(1) beta.init <- as.numeric(coefficients(model2)) out <- metrop(lupost, beta.init, 5e4, scale=0.03, x = x, y = y) out$accept

## [1] 0.19804

library(mcmcse)

## mcmcse: Monte Carlo Standard Errors for MCMC ## Version 1.3-2 created on 2017-07-03. ## copyright (c) 2012, James M. Flegal, University of California, Riverside ## John Hughes, University of Colorado, Denver ## Dootika Vats, University of Warwick ## Ning Dai, University of Minnesota ## For citation information, type citation("mcmcse"). ## Type help("mcmcse-package") to get started.

plot(ts(out$batch))

ess(out$batch)

## [1] 228.1051 485.8939 269.3972 614.3843 225.2397 403.7840

multiESS(out$batch)

## [1] 968.6197

- Effective sample sizes are small!
- See highest correlation components 1, 3, and 5, which also have the smallest ESS.
- Try to improve mixing by using a vector for
`scale`

.

out <- metrop(out, scale=c(.1, .03, .1, .03, .1, .03), x = x, y = y) out$accept

## [1] 0.17346

out$time

## user system elapsed ## 5.736 1.297 7.062

plot(ts(out$batch))

ess(out$batch)

## [1] 265.4184 789.0774 650.0411 686.0127 232.3179 322.1590

multiESS(out$batch)

## [1] 1116.557

- Getting better!

acf(out$batch[,1:3])