## Stat 295, 3/26/09

Please note: The homework for next week has been posted in the immediately preceding blog entry. It is due on Thursday.

Comments on the homework: Some groups have been having difficulty writing down the likelihood. The likelihood for observation j is $\lambda^{y_j} \exp(-\lambda)/y_j!$ Some were writing the first term over all data as $\lambda^{\vec{y}}$. What does that mean? There is no such thing. They should have written $\lambda^{\sum{y_j}}$, which is the product of the individual likelihoods for each observation, so that the likelihood to use should have been $\lambda^{\sum{y_j}}\exp(-n\lambda)$, where n is the number of observations.

Other problem: People were programming the Gibbs sampler…They were not updating the value of the variables sent to sample1(). Need betas[i]=beta=z\$beta; that middle assignment is critical. Without it, you are always sending the same position to the sampler, and it won’t move around as it should.

There was a perfect paper by undergrads in the course; the grad students need to match this!

We were talking about the linear regression model from last time, which was, you recall:

$y_i=\beta_0+\beta_1 x_{1j}+\beta_2 x_{2j}$

So we looked at the plots. We fit the least squares model using R function lm(), which gives the frequentist result. We looked at confidence intervals and point estimates.

For our Bayesian analysis we’ll put independent priors as per the blog last time. They were very broad priors that are almost flat. Jeff has learned that it is possible to put flat priors on the betas (but BUGS didn’t let you do that before so that’s why the priors we talked about on Tuesday are proper).

We analyzed the R code and the WinBUGS code.

Then we ran it from R (consult the R code). We looked at the output of the bugs() function in the arm package, which we called linreg.sim (see R code for this); We want Rhat close to 1 to verify good mixing. We looked at the traces, which looked good.

You can define a function and have BUGS evaluate that function over the chain and return a sample. So for example, the usual Rsquare is given by

$\begin{array}{rclr} R^2 & = & 1-\frac{\mbox{SSE}}{\mbox{SST}} \\ & = & 1-\frac{\sum{(y_j-\hat{y}_j)^2}}{\sum{(y_j-\bar{y})^2}}\end{array}$

Adjusted $R^2$: Divide the fraction’s numerator by $(n-p)$ and its denominator by $(n-1)$

In the Bayesian case define

$R_B^2=1-\frac{\sigma^2}{s_y^2}$ where $s_y^2 = \frac{1}{n-1}\sum{(y_j-\bar{y})^2}$

(closest to the adjusted $R^2$).

This is in the R code sent to WinBUGS.

typical.y (see the WinBUGS file) computes a posterior distribution on the estimated mean time it takes to service, using mean(cases[]), mean(distance[]).

We talked about other possibly useful features, like the log file, tracing, and the listserv.

Jeff remarked that the best way to learn BUGS is to look at the models that they have as examples.

Next example:

ANOVA models. One-way ANOVA is a useful linear model.

Our example: Interested in whether different exercise programs affect time to walking in infants.

Four different programs, called Active, Passive, No Exercise, Control Group.

Six infants in each group (randomized). Outcome is time in months when infants first walk.

Model:

$\begin{array}{rcl}y_{ij}&=&\alpha_i+\epsilon_{ij}\\ \alpha_i & \sim & N(\mu_{\alpha}, \sigma_{\alpha})\\ \epsilon_{ij} & \sim & N(0,\sigma) \end{array}$

This is a random effects model and also a hierarchical model.

$i=\mbox{group}, j=\mbox{infant}$ within the group.

$\alpha_i$ denotes the mean time to walking in group i. Without the second of these three lines it’s a one-way fixed effects ANOVA model. With it it is a random effects ANOVA model. It’s natural for a Bayesian to put a distribution on $\alpha_i$ so it’s natural (and simple) for her to consider a random effects model.