Saturday, February 27, 2016

Bayesian Analysis: Smear and Compare

This blog post is co-written with Julia Haaf (@JuliaHaaf).

Suppose Theory A makes a prediction about main effects and Theory B makes a prediction about interactions.  Can we compare the theories with data?  This question was posed on the Facebook Psychological Methods Group by Rickard Carlsson.

Uli Schimmak (@R_Index) put the discussion on general terms with this setup:
Should patients take Drug X to increase their life expectancy?
Theory A = All patients benefit equally (unlikely to be true).
Theory B1 = Women benefit, but men's LE is not affected.
Theory B2 = Women benefit, but men's LE decreases (cross-over).
Theory B3 = Women benefit, and men benefit too, but less.

We are going to assume that each of these statements represents a theoretically interesting position of constraint. The goal is to use data to state the relative evidence for or against these positions.  This question is pretty hard from a frequentist perspective as there are difficult order constraints to be considered.  Fortunately, it is relatively simple from a Bayesian model-comparison perspective.

Model Specifications

The first and perhaps most important steps is representing these verbal positions as competing statistical models, or performing model specification.  Model specification is a bit of an art, and here is our approach:

Let \(Y_{ijk}\) be life expectancy where \(i\) denotes gender (\(i=w\) for women; \(i=m\) for men), where \(j\) denotes drug status (\(j=0\) for placebo; \(j=1\) for treatment), and where \(k\) denote the replicate as there are several people per cell. We can start with the standard setup:
\[
Y_{ijk} \sim \mbox{Normal}(\mu_{ij},\sigma^2).
\]

The next step is building in meaningful constraints on true cell means \(\mu_{ij}\).  The standard approach is to think in terms of the grand mean, main effects, and interactions.  We think in this case and for these positions, the standard approach is not as suitable as the following two-cornerstone approach:

\(\mu_{w0} = \alpha\)
\(\mu_{w1} = \alpha + \beta\)
\(\mu_{m0} = \gamma \)
\(\mu_{m1} = \gamma+\delta \)

With this parameterization, all the models can be expressed as various constraints on the relationship between \(\beta\), the effect for women, and \(\delta\), the effect for men.

Model A: The constraint in Theory A is instantiated by setting \(\beta=\delta\). We place a half normal on this single parameter, the equal effect of the drug on men and women. See the Figure below.

Model B1: The constraint in Theory B1 instantiated by setting \(\beta>0\) and \(\delta=0\).  A half normal on \(\beta\) will do.

Model B2: The constraint in Theory B2 is instantiated by setting \(\beta>0\) and \(\delta<0\).  We used independent half normals here.

Model B3: The constrain in Theory B3 is that \(0<\delta<\beta\). This model is also shown, and it is similar to the preceding one; the difference is in the form of the constraints.



Of course, there are other models which might be useful including the null, models with no commitment to benefit for women, or models that do not assume women benefit more than men. Adding them presents no additional difficulties in the Bayesian approach.

Analysis: Smear & Compare

If we are willing to make fine specifications, as above, then it is straight forward to derive predictions for data of a set sample size. These predictions are shown as a function of the sample effect for men and women, that is, the change in lifespan between treatment and placebo for each gender.   These effects are denoted as as \(\hat{\beta}\) and \(\hat{\delta}\), respectively.  Here are the predictions:



Notice how these predictions are smeared versions of the model.  That is what sample noise does!

With these predictions, we are ready to observe data. Suppose we observe that the treatment extends women's lives by 2 years and men's lives by one year.  We have now included this observed value as a red dot in the below figures.


As can be seen, the observation is best predicted by Model B3.

The Bayes factor is the relative comparisons of these predictions.  We can ask how much B3 beats the other models.  Here it is:

B3 beats A by 3.7-to-1
B3 beats B1 by 18.2-to-1
B3 beats B2 by 212-to-1


2 comments:

David Kellen said...

Hi Jeff and Julia, nice post.

Quick question: How exactly did you impose the order constraints on B3? I ask because the way one does this can influence the prior predictive and the Bayes Factor.

Jeff Rouder said...

David, We don't have data here. So we just convolved the B3 model with white noise to derive the predictions. Then we renormalized. Here is the code.

#Model

normBiTb.1=function(theta1,theta2,Sigma) ifelse(theta1>theta2,2*dtmvnorm(cbind(theta1,theta2),c(0,0),Sigma,lower=c(0,0)),0)

normBiTb=function(theta1,theta2,Sigma)
{
temp=outer(theta1,theta2,normBiTb.1,Sigma)
tot=sum(temp)
return(temp/tot)
}


sd0=1.5
theta1=seq(-4,4,.05)
theta2=seq(-4,4,.05)
priorWedge=normBiTb(theta1,theta2,Sigma=matrix(ncol=2,c(sd0^2,0,0,sd0^2)))

#Here are the predictions:

kern = convKernel(sigma = 8,k="gaussian")

nrmlz=function(mat)
{
tot=sum(mat)
mat/tot
}

predWedge=nrmlz(applyFilter(priorWedge,kern))


Best,
Jeff