STAT 295 2/12/09

Happy 200th birthday, Charles Darwin!

Today we discussed the problem (#3) regarding the non-independent data. If you can't write f(x,y|θ)=k(x,y)f_y(y|θ)f_x(x|θ), then you cannot treat the sequential problem as if the y variable is independent of x (that is, the calculation given in the analogous problem in problem set #2 does not work). I need to do a bit of research on the k(x,y) piece of this.

Another example: Bioassay experiment. First: transforming parameters to (-∞,∞) is often a good idea. In this case we used a log transform. Similarly, the probabilities are transformed from (0,1) to (-∞,∞) using the logit transformation:

q=logit(p)=log(p/(1-p))

with q=(β₀+β₁x)=F(q) where p=F(q)=(1+exp(-q))^-1

log(p/(1-p))=logit(p)=β₀+β₁x

β₁ is interpreted as the difference in log odds for a 1-unit change in x (the dose).

Logistic regression is very important in statistics. Working statisticians will see it over and over again.

Note that F'(v)=F(v)(1-F(v))

We have animal data from the book, and wish to make inference on the βs.

The likelihood is binomial in p_i (i is the group) and then the product over all groups of animals. This is because we have two situations (the animal lived or it died) with probability p (given by the inverse logit) that it lived and (1-p) that it died. So the likelihood is the product of p's and (1-p)'s over all the animals in the study.

We'll use a sampling approach on a discretized version of β₀, β₁. But need to get an idea of how far to go. Go for the mode (MLE since priors are flat). Then estimate the standard errors, go out several standard errors and use that.

We did this in R using the R function glm (generalized linear model).

We recommend looking at functions in Albert's LearnBayes package to see how he did it. You'll learn new things about R, and have a better idea about what Albert is doing.

We looked at the plot (Figure 4.3 of the book) and saw that we've got to go pretty far out (farther than Jeff guessed) in order to get a decent sample.

How to sample from the distribution. The contour plot gives us an idea. Chop up region into a lot of little squares. Have something that is proportional to the density. Height of this surface is sort of like the probability that theta takes that value. Can use the sample() function to do this (sampling with replacement). We did this and the sample in Fig. 4.4 is the result. Looks pretty good.

We want to compute LD-50, the dose at which 50% of the animals die. Thus

Pr(y|x)=1/(1+exp(-(β₀+β₁x)))=0.5

1=exp(-(β₀+β₁x_LD50))

0=-β₀-β₁x_LD50, so x=-β₀/β₁

Now look at the Bayesian analysis. From the sample, divide the samples of β₀ by the samples of β₁, look at the 95% credible interval for the sample. We did this, looked at the histogram and computed the 95% credible interval. (The result is in the log dose).

We also used the density() function to draw the posterior density for the marginal of LD-50. This uses a nonparametric kernel estimator to draw a smooth line. β₁ is the log odds ratio for a unit change in dose. To get the odds ratio, exponentiate.

(The nonparametric kernel estimator works this way: At each point of the sample, we draw a peaked, narrow bell-shaped curve with unit area. We then add the curves for the whole sample and divide by the number of samples to get an approximation to the probability density curve.)

Sheila asked about hypothesis testing. We can get one-sided credible intervals (analogous to one-sided confidence intervals) quite simply. For example if s$y contains the samples of some quantity β, can compute probability that β>1 as mean((s$y)>0).

I pointed out that there is no Bayesian analog to frequentist two-sided hypothesis testing. We'll get to this later, when we discuss hypothesis testing.