## Bayesian regression – part II

In the earlier posting (part I), I demonstrated how one obtains a prior using a portion of the data and then uses this as a prior with the remaining portion of the data to obtain the posterior.

One advantage of using proper priors is that one can implement model selection.

Recall that the output of the MCMCregress function was fit2 — this is the regression model including all four predictors Food, Decor, Service, and East.

Suppose I want to compare this model with the model that removes the Food variable. I just adjust my prior. If I don’t believe this variable should be in the model, my prior for that regression coefficient should be 0 with a tiny variance or a large precision. Currently my prior mean is stored in the variable beta0 and the precision matrix is stored in B0. I change the component of beta0 corresponding to Food to be 0 and then assign a really big value to the (Food, Food) entry in the precision matrix B0. I rerun MCMCregress using this new model and the output is stored in fit2.a.

B0.a=B0; B0.a[2,2]=200000

fit2.a=MCMCregress(Price~Food+Decor+Service+East,

data = italian2, burnin = 0, mcmc = 10000,

thin = 1, verbose = 0, seed = NA, beta.start = NA,

b0 = beta0.a, B0 = B0.a, c0 = c0, d0 = d0,

marginal.likelihood = c(“Chib95”))

To compare the two models (the full one and the one with Food removed), I use the BayesFactor function in MCMCpack. This gives the log marginal likelihood for each model and compares models by Bayes factors. Here is the key output.

BayesFactor(fit2,fit2.a)

The matrix of the natural log Bayes Factors is:

fit2 fit2.a

fit2 0 303

fit2.a -303 0

## Logistic regression exercise

In my Bayesian class, I assigned a problem (exercise 8.3 from BCWR) where one is using logistic regression to model the trajectory of Mike Schmidt’s home run rates.

First, I have written a function to compute the log posterior for a logistic regression model with a vague prior placed on the regression vector. The function together with a short example how it works can be found here.

Second, it seems that there is a problem getting the laplace function to converge if you don’t use a really good starting value. Here is an alternative way of getting the posterior mode.

Suppose y is a vector containing the home run counts and n is a vector with the sample sizes. Let age be the vector of ages and age2 be the vector of ages squared. Define the response to be a matrix containing the vectors of successes and failures.

response=cbind(y,n-y)

Then the mle of beta is found using the glm command with the family=binomial option

fit=glm(response~age+age2,family=binomial)

The posterior mode is the mle:

beta=coef(fit)

The approximation to the posterior variance-covariance matrix is found by

V=vcov(fit)

Now you should be able to use the random walk metropolis algorithm to simulate from the posterior.

## Bayesian regression — part I

To illustrate Bayesian regression, I’ll use an example from Simon Sheather’s new book on linear regression. A restaurant guide collects a number of variables from a group of Italian restaurants in New York City. One is interested in predicting the Price of a meal based on Food (measure of quality), Decor (measure of decor of the restaurant), Service (measure of quality of service), and East (dummy variable indicating if the restaurant is east or west of 5th Avenue).

Here is the dataset in csv format: http://www.stat.tamu.edu/~sheather/book/docs/datasets/nyc.csv

Assuming the usual regression model

, ,

suppose we assign the following prior to : we assume independence with and . Here is the precision matrix, the inverse of the variance-covariance matrix.

How do we obtain values of the parameters ? We’ll use a small part of the data to construct our prior, and then use the remaining part of the data to do inference and eventually do some model selection.

First, I read in the data, select a random sample of rows, and partition the data into the two datasets italian1 and italian2.

italian=read.csv(“http://www.stat.tamu.edu/~sheather/book/docs/datasets/nyc.csv”)

indices=sample(168,size=25)

italian1=italian[indices,]

italian2=italian[-indices,]

y=italian1$Price

X=with(italian1,cbind(1,Food,Decor,Service,East))

fit=blinreg(y,X,5000)

Sigma0=cov(fit$beta)

B0=solve(Sigma0)

c0=(25-5)

d0=sum((y-X%*%beta0)^2)

data = italian2, burnin = 0, mcmc = 10000,

thin = 1, verbose = 0, seed = NA, beta.start = NA,

b0 = beta0, B0 = B0, c0 = c0, d0 = d0,

marginal.likelihood = c(“Chib95”))

beta=fit2[,1:5]

sigma=sqrt(fit2[,6])

lp=c(beta%*%x)

Min. 1st Qu. Median Mean 3rd Qu. Max.

43.59 46.44 47.00 46.99 47.53 50.17

summary(ys)

Min. 1st Qu. Median Mean 3rd Qu. Max.

25.10 43.12 46.90 46.93 50.75 67.82

## Are birth weights normal?

A few weeks ago, I asked my class for an example of data that would be normally distributed and someone mentioned weights of people.

Are weights, more specifically birthweights of newborns normally distributed?

Peter Dunn in a JSE article here analyzes statistics for 44 babies born in a single day in an Australian hospital. The article has a link to the dataset babyboom.dat that contains the data.

If we construct a histogram of the birthweights, we get the following display.

We see some left skewness in this distribution — does this indicate that our normality assumption is false?

One can check the normality assumption by means of the posterior predictive distribution. Here’s an outline of the approach. (There are special functions in the LearnBayes package that make this process simple to implement.)

1. First we fit the normal() model using the traditional noninformative prior proportional to . A simulated sample from the posterior is found using the function normpostsim.

fit=normpostsim(birth.weight,10000)

2. Next we think of an appropriate testing function. Since the birth weights don’t appear symmetric, a natural checking function is the sample skewness. We write a short function “skew” to compute this statistic.

sum((y-mean(y))^3)/sd(y)^3/(length(y)-1)

## A “robustness” illustration of Bayes factors

Here is another illustration of Bayes factors. We are worried about the possibility of outliers and so we assume the observations are distributed from a t distribution with location , scale , and degrees of freedom . We have some prior beliefs about the location and scale: we assume independent with and .

We observe the following data:

112.6,86.0,99.1,106.3,101.5,80.5,99.6,108.6,95.9)

{

mu=theta[1]

sigma=exp(theta[2])

y=stuff$y

df=stuff$df

loglike=sum(dt((y-mu)/sigma, df, log=TRUE)-log(sigma))

logprior=dnorm(mu,100,10,log=TRUE)+dnorm(log(sigma),0,2, log=TRUE)

loglike+logprior

}

## A Bayesian test of the equality of proportions

One problem with a p-value is that it is hard to interpret if you don’t know the sample size. Andy Gelman talks about this problem on his blog.

Here I’ll contrast the message communicated by a p-value and a Bayesian test by looking at the following three tables:

Table A: [32 18; 18 32]

Table B: [272 228; 228,272]

Table C: [2570 2430; 2430,2570]

[1] 0.9436246

[1] -0.3035413

[1] -1.411366

## Bayesian testing of a point null hypothesis

Here is today’s example illustrating a Bayesian test of a point null hypothesis. In Agresti and Franklin’s intro stats book, there is an example wondering if Americans work, on average, 40 hours per week. If represents the mean hours worked per week for all Americans, we are interested in testing against the alternative .

In the example in Agresti and Franklin’s text, a sample of 868 workers were selected and one observes that the workers worked an average of 39.11 hours and one assumes the population standard deviation is equal to . The test statistic is

and one computes a p-value of 0.061. There seems to be moderately significant evidence that the mean working hours for all workers is different from 40.

Here’s an outline of a Bayesian test (I give more details in class):

1. We assign probabilities to the hypotheses and . Here we assume each hypothesis is equally likely.

2. Under the hypothesis , the prior for is a point mass at the hypotheized value of 40.

3. Under the alternative hypothesis, we’ll place a N(40, ) prior on . This reflects the belief that, under the alternative, values of near 40 are more likely that values far from 40.

Under these assumptions, we give (in class) the posterior probability of the hypothesis . This probability is programmed in the function mnormt.twosided in the LearnBayes package. The inputs to this function are (1) the value to be tested, (2) the prior probability of , (3) a vector of values of and (4) a data vector consisting of . The outputs are the Bayes factor in support of and the posterior probability of .

Since it is not clear how to specify the prior standard deviation — we compute the posterior probability for a range of values of . I write a short function to output the posterior probability as a function of and then I use the curve function to plot this probability for values of from 0.01 to 2.

What we see is that the posterior probability of exceeds 0.35 for all values of . This suggests that the p-value overstates the evidence against the null hypothesis.

> data=c(39.11, 868, 14)

> prob.H=function(tau)

+ mnormt.twosided(40,.5,tau,data)$post

> curve(prob.H(x),from=.01,to=2,xlab=”TAU”,ylab=”P(H)”)