5

If $X$ is a categorical variable, and I am interested in the posterior distributions of $\beta_1$, where $\beta_1$ is a vector of coefficients, one for each level of X, are these equivalent models?

Model 1: $$Y \sim ( \beta_0 + \beta_1X_1,\sigma^2)$$ $$\beta_1 \sim N(0,\tau)$$

Model 2: $$Y \sim(\beta_1X_1,\sigma^2)$$ $$\beta_1 \sim N(\beta_0,\tau)$$

In both models, $\beta_0$ has an informed prior, $\tau$ has a weak prior: $$\beta_0\sim \text{Gamma}(6,3)$$ $$\tau \sim \text{Gamma}(0.1,0.1)$$

If these are not equivalent, why not?

I ask because although I thought these would be the same, but when implemented in JAGS, the MCMC chain for $\beta_0$ from the first model has much higher autocorrelation than the second.

I (originally) put my R code here SO.test.tgz, but it is no longer available. I apologize ...

David LeBauer
  • 7,060
  • 6
  • 44
  • 89

4 Answers4

6

Updated Response: You still don't have a full specification for model #2. However, I can sort of guess what you mean -- correct me if I'm wrong. The trouble is that the statements $Y = \beta_1 X$ & $Y = \beta_0 + \beta_1 X$ are not probabilistic.

[ Aside: In a mathematical sense, you're defining a set of linear equations. If you have only one datapoint then you can solve $Y = \beta_1 X$ for the only value of $\beta_1$ which satisfies this equation. If you have more than one value of $Y$ and $X_1$ then this represents an overconstrained system -- without a solution. ]

My hunch is that you mean the following: $$Y \sim N( \beta_1 X_1, \sigma^2)$$ for some value of $\sigma^2$ which is assumed to be known -- (or perhaps you also want to estimate it). This places your random effects model within the context of a standard linear mixed effects model.

If we then assume that $Y$ is distributed as I have stated, then:

Answer: No, they are not equal. In the first model, $\beta_0$ is essentially an intercept (Imagine multiplying $\beta_0$ by $X_0$ where $X_0$ is always equal to 1). In the second model, $\beta_0$ represents a random offset from $\beta_1$. Note that this model (#2) wouldn't make sense to someone who doesn't do Bayes... It would be identical to running a linear regression where you have two predictors that are perfectly multi-collinear, but since you've made distributional assumptions within a Bayesian model, you can estimate it. That said, I'm not sure you should.

Note: You can run the first model without the added distributional assumption (which I included above more for model #2). However, I've never seen this sort of specification. I think it would be identical to stating that $(Y - \beta_1 X_1) \sim Gamma(6,3)$. In other words, an error term which is distributed as a gamma random variable. My hunch is that you instead want a random intercept model (with normally distributed errors). If that be the case, use model #1, don't use model #2. Yes $\beta_0$ and $\beta_1$ will exhibit autocorrelation until you standardize the levels of your categorical variable -- make sure that the summation across all individuals is equal to zero (which you can do by subtracting the mean of $X_1$ from every single observation of $X_1$).

M. Tibbits
  • 1,693
  • 14
  • 25
4

The only similarity in the two models is the general type of models they belong to, otherwise they are not similar in general as pointed out by M. Tibbits.

Both these models belong the class of hierarchical models with varying slope (cf Gelman and Hill 2006 for detailed treatment)

The answer for "why not" are many and one of them is pointed out by M. Tibbits. A few more are:

  • In the model 2 the slope on an average is 6/3, where as it is 0 for model 1.(This description can be made much better if we had data but it is "almost" accurate given the limited description)

  • You would expect the data in model 1 to be distributed randomly along a approximately horizontal line where as in model 2 you would expect data to be randomly distributed along a line of slope approx. beta_0 and zero intercept.

These answers can be best varified if you similated the data based on your hierarchical setting and look at the results of the analysis.

Thanks,

S.

suncoolsu
  • 6,202
  • 30
  • 46
2

You already have good answers and have accepted one, but I'm not sure anyone's put it plainly enough for even me to intuit. At the core, your two models are:

[1] $Y=\beta_0 + \beta_1X_1 + \epsilon$

[2] $Y=(\beta_0 + \beta_1)X_1 + \epsilon$ = $\beta_0X_1 + \beta_1X_1 + \epsilon$

Yes, the $\beta$'s (and $\epsilon$) have distributions in your example, but in my mind that only smears the value out and doesn't change the form of the formulae. In that light, they're obviously different: your second model essentially fixes the intercept $\beta_0=0$ and then makes the slope $\beta_1$ a bit more complicated.

David LeBauer
  • 7,060
  • 6
  • 44
  • 89
Wayne
  • 19,981
  • 4
  • 50
  • 99
1

If $x_1=0$, then $Y$~$(\beta_0,\sigma^2)$ in model 1 and $Y$~$(0,\sigma^2)$ in model 2.

If $x_1=1$, then $Y$~$(\beta_0+\beta_1,\sigma^2)$ in model 1 and $Y$~$(\beta_1,\sigma^2)$ in model 2.

Look for example at the first line: is $\beta_0$ a random variable or is a zero constant?

glassy
  • 1,012
  • 2
  • 10
  • 16