Constraining batches of coefficients for finite population variances in two-factor ANOVA

Modeling
1

Hi,
I’m looking for advice on how to specify an ANOVA model, where I can estimate the finite population variances of batches of coefficients. Suppose I have two factors, X and Z, each with four levels. My general strategy is to make three design matrices X_d (the main effects of X), Z_d (the main effects of Z), and I_d (the interactions between X and Z). I generally use ‘effects’ or ‘dummy’ coding to generate these matrices and keep the model identifiable (although it appears not to matter which I use for estimates of finite population variances for coefficient batches). I then fit the following model:

anova_model = """
data{
    int<lower=1> N;
    int<lower=1> K1;
    int<lower=1> K2;
    int<lower=1> K3;
    matrix[N, K1] X_d;
    matrix[N, K2] Z_d;
    matrix[N, K3] I_d;
    vector[N] y;
}
parameters{
    vector[K1] B_x;
    vector[K2] B_z;
    vector[K3] B_i;
    real<lower=0> sd_y;
    real<lower=0> sd_x;
    real<lower=0> sd_z;
    real<lower=0> sd_i;
    real mu;
}
transformed parameters{
    vector[N] yhat;

    yhat = mu + X_d*B_x + Z_d*B_z + I_d*B_i ;
}
model{
    y ~ normal(yhat, sd_y) ;
    B_x ~ normal(0,sd_x);
    B_z ~ normal(0,sd_z);
    B_i ~ normal(0,sd_i);
}
generated quantities{
    vector[N] e_y;
    real s_e;
    real s_x;
    real s_z;
    real s_i;

    e_y = y - yhat;
    s_e = sd(e_y);
    s_x = sd(B_x);
    s_z = sd(B_z);
    s_i = sd(B_i);
}
"""

No problems there, everything typically works as expected. The model converges, the results make sense, etc.

The issue I run into is when one of the factors only has 2 levels, 1 df. For example, if Z only has two levels, using effects or dummy coding results in an N x 1 design matrix that has no finite population variance. I’ve been reading Gelman (2005), Hill and Gelman, and other references. Both sources can do this kind of ANOVA when a factor only has one level (i.e. the latin square example), and both sources also use unconstrained coefficients. I’ve tried making three new design matrices: X_d2, Z_d2, I_d2, which have no intercept and are fully specified (i.e. X_d2 has four columns, one per level, Z_d2 has two columns, one per level, and I_d2 has 8 columns for the interactions). But the model performs poorly. So I modified the code to constrain the effects:

anova_model = """
data{
    int<lower=1> N;
    int<lower=1> K1;
    int<lower=1> K2;
    int<lower=1> K3;
    matrix[N, K1] X_d2;
    matrix[N, K2] Z_d2;
    matrix[N, K3] I_d2;
    vector[N] y;
}
parameters{
    vector[K1] B_x;
    vector[K2] B_z;
    vector[K3] B_i;
    real<lower=0> sd_y;
    real<lower=0> sd_x;
    real<lower=0> sd_z;
    real<lower=0> sd_i;
    real mu;
}
transformed parameters{
    vector[N] yhat;
    vector[K1] B_x_adj;
    vector[K2] B_z_adj;
    vector[K3] B_i_adj;
    real B0;

    B_x_adj = B_x - mean(B_x) ;
    B_z_adj = B_z - mean(B_z);
    B_i_adj = B_i - mean(B_i);
    B0 = mu + mean(B_x) + mean(B_z) + mean(B_i);

    yhat = B0 + X_d2*B_x_adj + Z_d2*B_z_adj + I_d2*B_i_adj ;
}
model{
    y ~ normal(yhat, sd_y) ;
    B_z ~ normal(0,sd_z);
    B_x ~ normal(0,sd_x);
    B_i ~ normal(0,sd_i);
    mu ~ normal(0,1); # the response has been standardized, so this prior makes sense
}
generated quantities{
    vector[N] e_y;
    real s_e;
    real s_x;
    real s_z;
    real s_i;

    e_y = y - yhat;
    s_e = sd(e_y);
    s_x = sd(B_x_adj);
    s_z = sd(B_z_adj);
    s_i = sd(B_i_adj);
}
"""

This model converges much more rapidly and gives results that make sense. I guess I’m wondering if this is appropriate. Should I be using the constrained parameters both for estimation (i.e. yhat) and for finite population variances (s_x, s_z, s_i)? The finite population variances are the same either way, because s_x = sd(B_x_adj) = sd(B_x - mean(B_x)). But I’m not clear if this is the same model referenced in Gelman (2005), and how that model was able to use unconstrained coefficients. I hope I’m on the right track.