My background is in CS and ML, so I am quite familiar with hierarchical Bayesian models -- and probabilistic programming. But I do not have much experience with traditional statistical methods. This makes my understanding of the standard formalisms for linear models involving "factor variables" and "contrasts" a bit flawed, I guess.

A (1D, one variable) linear model of $Y_i$ given $X_i$, as I see it, is of the form

$$ Y_i \sim \mathrm{Normal}(\alpha + \beta X_i, \,\sigma) $$

where the $X_i$ are continuous values, to which is applied an affine transformation, and then some noise around it.

However, when people write a linear model in R as

y ~ x

where x really is a categorical variable (say, with values in $\{1, \ldots, k\}$), this, if I understand correctly, designates the model

$$ Y_i \sim Normal([X_i = 1]\beta_1 + \cdots + [X_i = k]\beta_k,\, \sigma) $$

(using Iverson brackets). Or even something equivalent but stranger, by encoding the k values of $X$ in another orthogonal basis, in general:

$$ Y_i \sim Normal(\mathrm{encode}(X_i) \cdot \vec{\beta}, \sigma) $$

Why isn't the formulation

$$ Y_i \sim Normal(\mu_{X_i},\, \sigma) $$

used instead? (Which I called a "mixture" in the title because it's really a GMM where we know the cluster assignments in advance.)

I find this much more natural to describe and interpret. $Y_i$ does not vary with $X_i$ -- it discontinuously switches to a completely unrelated distribution. The domain of the $X_i$ in contrast coding is not the same as above -- they really have to come from a simplex of some sort. But I have never really seen this kind of interpretation -- "what is the treatment effect for a person that is given 80% of condition A and 20% of condition B"? The value of the $X_i$ loses its inherent "topology" if we go that way.

Also, in practical terms, there is even an R package emmeans that converts a linear model fit to exactly this form given through the "cluster means".

So what am I missing? Is it only computational convenience -- reducing the discrete case to the already known regression form -- or is there an interpretational advantage? Why should I ever construct my mental model of a generative process in terms of contrasts?


I think you are confusing mathematical notation with programming language syntax in here. In R ~ is a function, that creates a formula class object. The formula class is a way of defining the model structure in R in terms of target variable(s) on the left hand side of ~, and features on the right hand side. So y ~ x means "A model where we predict y using feature x". Formula has nothing to do with the statistical model underlying the relationship between variables, same formula can be used for many different models.

As about mathematical notation, the regression model is defined as

$$ y = \mathbf{X}\beta + \varepsilon $$

what is equivalent to probabilistic notation

$$ y \sim \mathsf{Normal}(\mu_Y, \sigma) $$

where $\mu_Y = \mathbf{X}\beta$. Here $\sim$ indeed means "is distributed as". As a slight correction to notation you used, it's $\mu_Y$, the mean of $Y$ variable, not $\mu_X$ ($X$ does not even have to be a random variable).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.