hierarchical models

[Ross Boylan]

I estimated a bunch of hierarchical models, and am looking for feedback
about whether the models are appropriate. I have some doubts.

The stan code is below, but the key point is that each batch of
coefficients had a prior with mean 0 and and estimated variance. Call
two such distribution F and G. Let b1, b2, ... be coefficients for
covariants x1, x2, ... (there were about 10) and let b1A, b1B, and b1C
be coefficients for x1 a categorical covariate takes values A, B, or C
respectively (exactly 3 levels). So
b1A ~ F
b1B ~ F
b2A ~ F.
The coefficients without interaction were in a separate batch, so
b1 ~ G
b2 ~ G.

I believe the usual procedure would be
b1A, b1B, b1C ~ F(m1, s1)
b2A, b2B, b2C ~ F(m2, s2)
with no separate b1 and b2 coefficient, and the m's and s's (means and
standard deviations) in turn coming from other distributions.

I think the b1 and b2 coefficients in my formulation are equivalent to
m1 and m2, so the difference there is merely formal.

However, my batching scheme is equivalent to s1=s2. It didn't seem a
good idea to try fitting the standard deviation when each group had only
3 members. It was often challenging to fit the models I did; adding
more parameters seemed pushing my luck.

Finally, I'm a bit concerned that since at the highest hiearchical level
(or do I mean shallowest? the most meta level, the distribution of the
b1, b2, ... coefficients) I simply assumed a mean of 0, the estimates
may be unduly small. They might also be too small if the estimated
variance parameters are too small, e.g., if there were notable
differences in the interaction term for only 1 variable, the batching of
those interactions with all the others would lead to a small sigma.

So, does the procedure I followed seem reasonable? If not, what other
models would be advisable?

Here's a representative model:
data {
int<lower=0> N; // number of respondents (~400)
int<lower=0> J; //number of covariates (~10, not counting clusters)
int<lower=0> NSD; // number of classes of beta for sd (~3 classes)
int<lower=1,upper=NSD> SDClass[J]; // beta[j, k] has sigma_beta[SDClass[j]]
int<lower=0> NCLUSTER; //number of seeds (~200)
int iCluster[N]; // cluster id
matrix[N, J] x; // covariates
int<lower=1,upper=15> nDetail[N]; // outcome
}
parameters {
vector[J] beta_raw; // standard covariate coefficients
vector[NCLUSTER] alpha_raw; // cluster effects
real<lower=0> sigma;
real<lower=0> sigma_beta[NSD];
}
transformed parameters {
vector[J] beta;
vector[NCLUSTER] alpha;
for(ir in 1:J) {
beta[ir] <- beta_raw[ir]*sigma_beta[SDClass[ir]];
}
alpha <- alpha_raw*sigma;
}
model {
vector[N] beta_x;
sigma ~ cauchy(0, 5);
sigma_beta ~ cauchy(0, 5);
beta_raw[1] ~ uniform(-5, 5);
for (j in 2:J) {
beta_raw[j] ~ normal(0, 1);
}
alpha_raw ~ normal(0, 1);
for (n in 1:N) {
beta_x[n] <- (x[n] * beta)+alpha[iCluster[n]];
//truncated poisson can't be vectorized
// and is not available in log form
nDetail[n] ~ poisson(exp(beta_x[n])) T[1, 15];
}
}

Some of the models were more complex, including multinomial model with 5
outcome categories and a model with 2 nested levels of clustering. The
more complex models had considerably more observations, about 2000.

Thanks for any advice.
Ross Boylan

[Bob Carpenter]

It can be hard to fit groups with few members, and 3 is
difficult. It's not impossible, but you'll need reasonable
priors.

Otherwise, the basic idea of grouping parameters in a hierarchical
model is exchangeability (given the inputs for modeling).

As to assuming a mean of 0 at the top of the hierarchy, you
probably want a fairly diffuse standard deviation. And
that should be chosen given what you know about the problem.
Things like whether the predictors standardized.

There's also a major techncial problem for Stan
with the model as formulated.

For the sampler to work effectively, all parameters should
have support at any point that satisfies their constraint.

In the attached model, a variable is declared as

vector[J] beta_raw; // standard covariate coefficients

without constraints and then given a constrained density:

beta_raw[1] ~ uniform(-5, 5);

It would probably make sense to have a more meaningful prior
than uniform.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>