# Radon computations for chapter 15 of the book # 1. No-pooling and complete-pooling models # 2. Multilevel model # fit the model using bugs radon.data <- list ("n", "J", "x", "y", "county") radon.inits <- function (){ list (a=rnorm(J), b=rnorm(1), mu.a=rnorm(1), sigma.y=runif(1), sigma.a=runif(1)) } radon.parameters <- c ("a", "b", "mu.a", "sigma.y", "sigma.a") radon.bugs.1 <- bugs (radon.data, radon.inits, radon.parameters, "radon.1.bug", n.chains=3, n.iter=10) windows () plot (radon.bugs.1) attach.bugs (radon.bugs.1) # 3. Multilevel model including uranium as a county-level predictor # fit the model using bugs radon.data <- list ("n", "J", "x", "y", "county", "u") radon.inits <- function (){ list (a=rnorm(J), b=rnorm(1), g.0=rnorm(1), g.1=rnorm(1), sigma.y=runif(1), sigma.a=runif(1)) } radon.parameters <- c ("a", "b", "g.0", "g.1", "sigma.y", "sigma.a") radon.bugs.2 <- bugs (radon.data, radon.inits, radon.parameters, "radon.2.bug", n.chains=3, n.iter=10) windows () plot (radon.bugs.2) # 4. Predictions # these need to be fixed!!! y.squiggle <- rnorm (1, a.hat[26] + b.hat[26], sigma.y) y.squiggle <- rep (NA, n.sims) for (s in 1:n.sims){ y.squiggle <- rnorm (1, a.hat[26] + b.hat[26], sigma.y) } quantile (y.squiggle, c(.25,.5,.75)) unlogged <- exp (y.squiggle) mean (unlogged)