Applied 36-724: Bayesian and Computational Statistics Homework 4: Due Friday March 10, 2006
Announcements: In Gelman, please look at Chapter 6 (if you haven't already). Please download and examine the following papers from jstor: Kass, R. E. and Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 430, 773795. Cowles, M. K. and Carlin, B. P. (1996). Markov Chain Monte Carlo Convergence Diagnostics: A Comparative Review. Journal of the American Statistical Association, 91, 883904 . Kass, R. E., Carlin, B. P., Gelman, A., Neal, R. M. (1998). Markov Chain Monte Carlo in Practice: A Roundtable Discussion. The American Statistician, 52, 93100. (1998). In the remaining lectures for this part of the course I will discuss model fitting, model checking and model comparison, using examples from Parts IV and V of Gelman (and other sources). Problems: 1. Refer again to problem #12, pp. 99-100 of Gelman (the data is in Table 2.2, p. 69). priors on and , the joint model for the parameters and data is p(y, , |t) =
i
If we put independent
Pois(yi | + ti )p()p()
If you were happy with your informative prior from this problem for HW #3, you may use it. Otherwise try something simple like p() = p() = with appropriate choices for , 2 , , and 2 . (a) Give expressions for the complete conditionals for and . If there is a meaningful way to simplify your expressions, please do so. (b) Construct an MCMC algorithm for sampling from the posterior p(, |y) for your model. Chances are you will have to use Metropolis steps for both and . Use the "normal random walk" proposal distributions, i.e., qm (|(m-1) ) = qm (|(m-1) ) = N(|(m-1) , c ) N(|(m-1) , c ) N(| , 2 ) 2 N(| , )
where c and c are tuning parameters that you want to tune so that the acceptance rate is between 25% and 50% (and, really, closer to 50%; see e.g. pp 306-307 of Gelman). Code your algorithm in R. Notes: You may find you need to work with the logs of the complete conditionals in order to avoid numerical problems when multiplying together the terms of the likelihood and prior. Similarly, you may find you need to work with the log of the acceptance ratio for each parameter. (c) Generate an appropriate sample from the stationary distribution of your MCMC routine (i.e. throw away ^ an initial burn-in segment). Check, using time series plots, autocorrelation plots, R, etc., that your chain has reached stationarity, and decide whether your chain has to be thinned or block-averaged to reduce auto-correlation between draws. 1
(d) Produce a final sample from your MCMC chain, and produce appropriate numerical and graphical summaries of the p(, |y) using your sample. Compare with your results from HW #3. 2. The file "gss-data.r" defines a data frame gss.data containing the responses of 3121 respondents in the 1972, 1973, and 1974 General Social Survey. The responses are cross-classified by year of study (3 levels) and response (0=no, 1=yes) to each of three questions. The questions all begin: "Please tell me whether or not you think it should be possible for a pregnant woman to obtain a legal abortion if..." and end with: y1 : "...she is married and does not want any more children." y2 : "...the family has a very low income and cannot afford any more children." y3 : "...she is not married and does not want to marry the man." This data set is discussed in Tanner's (1996) book "Tools for Statistical Inference" and was considered by Goodman (1974, Biometrika) and by Haberman (1979, Analysis of Qualitative Data, Vol 2). (a) Use par(mfrow=c(2,2)) and barplot() to make barplots of the frequencies of the eight response patterns in each of the three years 1972, 1973, 1974. Because of the preponderance of "all yes's" and "all no's" responses, it is plausible to think that there are basically two classes of people responding to the survey, a class C1 of people who are disinclined to accept abortion under any of the three listed conditions, and a class C2 who would be inclined to accept abortion under all of them. This suggests a model for one respondent of the form
W = J
P[Yi1 yi1 , Yi2 = yi2 , Yi3 = yi3 , | p, ] =
w=1
w
j=1
pwi jj [1 - pw j ]1-yi j
y
()
where w = P[i Cw ], yi j = 1 if respondent i responds "yes" to question j, and yi j = 0 for "no", and pi j = P[yi j = 1|i Cw ], and where there are W = 2 classes of respondents and J = 3 questions. We wish to estimate the parameters w (w = 1, 2), and pw j (w = 1, 2, j = 1, 2, 3), for the N = 1055 respondents in the 1972 portion of the survey. Assuming independence between respondents, explicitly write out the likelihood g(Y | p, ), for the observed data matrix Y = [yi j : i = 1, . . . , N, j = 1, 2, 3]. (b) The likelihood is greatly simplified by use of data augmentation. Consider the latent variables (a.k.a. missing data, or augmentation data) zi2 = Show that if we assume zi2 yi j
iid
0, 1,
if i C1 if i C2
,
zi1 = 1 - zi2
indep
Bin(1, 2) Bin(1, p1 j), if zi2 = 0 Bin(1, p2 j), if zi2 = 1
then the distribution of the observed data is as in (). Show, therefore, that we may take the complete data likelihood to be ziw yi j 1-yi j pw j [1 - pw j ] w f (Y, Z|p, ) = w j i n w c = pwwjj [1 - pw j ]nw -cw j w
w j
where nw =
i ziw ,
and cw j =
i ziw yi j .
2
(c) The model in part (c) is equivalent to the hierarchical model yi j |z, p zi2 | pw j 2
indep
Bernoulli(pzi2+1 , j) Bernoulli(2) (some priors) (some prior)
iid
Choose reasonable priors for the p's and for 2 (but not also 1 --why?). Estimate this model in BUGS, for the 1972 data only. Notes: This is an example of a "latent class" model, also called a "discrete mixture" or "clustering" model (the clases C1 and C2 represent latent classes, or clusters, of response patterns). As with many such models, there is a relabeling ambiguity: if we interchange the class indices w = 1 and w = 2 throughout, the fit to the data is the same; this can lead to odd bimodalities in posterior distributions, for example. A common way to resolve the ambiguity is to arrange that p2 j > p1 j for one particular question j. It may be worthwhile to try to run BUGS from R in this problem, following the examples in lecture and on Gelman's website. If you cannot get BUGS to work (with or without R), this is still a fairly reasonable problem to do in R, as in problem 1 (with judicious choice of priors, everything is pure Gibbs steps). Generate an appropriate sample from the stationary distribution of your MCMC routine (i.e. throw away ^ an initial burn-in segment). Check, using time series plots, autocorrelation plots, R, etc., that your chain has reached stationarity, and decide whether your chain has to be thinned or block-averaged to reduce auto-correlation between draws. Interesting things to look at would be: 2 The pw j 's Perhaps a few "interesting" ziw 's, or P(zi2 = 1|yi1 , yi2 , yi3 ) for a few interesting cases (how many such probabilities are there?). Plot estimated probabilities of the eight response patterns in each latent class, and compare with the barplots you made in part (a). 3. One hypothesis about the different years of the survey is that the proportions of respondents in C1 and C2 change from year to year (but the response behavior in each class stays the same). Letting t be the year, this suggests the complete data likelihood zitw 1974 y f (Y, Z, |p, ) = pwi jj [1 - pw j ]1-yi j tw w j i t=1972 1974 W ntw J cw j (nw -cw j ) = pw j (1 - ptw j ) tw
t=1972 w=1 j=1
now = (tw , t = 1972, 1973, 1974, w = 1, 2); zitw = 1 if respondent i is from year t and is in class Cw (and 0 otherwise); ntw = i zitw , nw = t i zitw ; and cw j = t i zitw yi j .
(a) Fit this model in BUGS or R (whatever you used in problem 2 should work here) and summarize the posterior distributions of the 's. Compare with the posterior distribution of 2 in problem 2. (b) A way to test the hypothesis about the model is to compare the two models tw freely varying as t = 1972, 1973, 1974, w = 1, 2 tw w as t = 1972, 1973, 1974, w = 1, 2 (i.e., tw is allowed to vary in w but not in t). Carry out this comparison, using Bayes Factors.
3