3
$\begingroup$

Gelman and Hill 2007 present a varying using WinBUGS and R:

Simplifying the data description hugely, the data is from 85 counties, with log radon level measured in the basement and the first floor of houses in these counties (multiple houses in each county). y is log radon level, and x is floor (0 for basement, 1 for first floor, so treatment contrasts), county marks the id of the relevant county. This is data from Gelman and Hill's book.

If this were pure R, we could do:

lmer(y~x+(1|county)).

Let's say the research question is: does being in the basement vs 1st floor affect radon emissions? I.e., is the slope b significantly different from 0 (it should be negative, because the higher up you are, the less the radon).

I have WinBUGS code from Gelman and Hill's book (radon.1.bug), and modified code (radon.1a.bug), modified by me because I want an informative prior for the slope for county. I just used some rather large value for b as a prior, and a fairly low sd:

b ~ dnorm (10, .0001)

I first thought I could get away with making this one change; I expected the posterior estimate for b to be skewed towards my prior because the prior has such a low sd (reflects my strong belief that beta=10?). But this had absolutely no impact on the posterior!

Can someone help me understand what I'm doing wrong?

To run the model:

  1. Download the following files from this site:

http://www.stat.columbia.edu/~gelman/arm/examples/radon/

cty.dat, srrs2.dat, radon.1.bug (the WinBUGS model).

Then download and run radon_setup.R to set up the data.

  1. The code I used (Modified slightly from Gelman and Hill):

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

    library(R2WinBUGS)

    ## Specify the location of your WinBUGS directory:

    bugsdir<-"C:/Users/Shravan/Desktop/winbugs14/WinBUGS14"

    radon.1 <- bugs (radon.data, radon.inits, radon.parameters, "radon.1.bug", n.chains=3, n.iter=500,bugs.directory=bugsdir)

    plot (radon.1) radon.1.noburnin <- bugs (radon.data, radon.inits, radon.parameters, "radon.1.bug", n.chains=3, n.iter=500, n.burnin=0, bugs.directory=bugsdir)

    plot (radon.1)

I modified radon.1.bug as follows (call this radon.1a.bug and replace in the above code to see what happens---it seems precisely nothing changes in the estimated slope):

### Start WinBUGS model, call it radon.1a.bug:
model {
  for (i in 1:n){
    y[i] ~ dnorm (y.hat[i], tau.y)
    y.hat[i] <- a[county[i]] + b*x[i]
  }
  #b ~ dnorm (0, .0001)  ## <- Change this to an "informative prior"
  b ~  dnorm (10, .0001)  ## Shouldn't this influence the posterior?
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif (0, 100)

  for (j in 1:J){
    a[j] ~ dnorm (mu.a, tau.a)
  }
  mu.a ~ dnorm (0, .0001)
  tau.a <- pow(sigma.a, -2)
  sigma.a ~ dunif (0, 100)
}
### End WinBUGS code.

I'd be grateful for any help, including pointing me to an example online or a book where an informative prior is used in the context of a normal-normal conjugate prior, in either a single t-test setting or a varying intercepts model of the sort described above.

$\endgroup$
0

1 Answer 1

1
$\begingroup$

Note that WinBUGS' parameterization of the normal distribution is based on precision rather than variance. Therefore, the original prior b ~ dnorm (0, .0001) and your modified prior b ~ dnorm (10, .0001) are essentially the same, since you use a rather large variance (i.e., small precision). In order to make it informative (and thus, influential wrt to posterior), you would need to make the variance small, i.e. the precision large (e.g. b ~ dnorm(0, 4). Since I don't know how you summarize the historical information available for $b$, I can't tell you exactly what variance would be appropriate; for some formalized ways to summarize historical information, I recommend e.g. Neuenschwander B, Capkun-Niggli G, Branson M, Spiegelhalter DJ. Summarizing historical information on controls in clinical trials.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.