I have a dataset where parents and staff were asked to rate how different variables would be expected to affect how much time parents were present at the bedside. The original purpose was to create informative priors for a group of parents we were recruiting with measurements of time at bedside. In the interim, we are interested in how parent and staff responses differ.
I fit the following:
stan1 = stan_lmer(acuity ~ role + (1|role), data = all_merged)
stan2 = stan_lmer(acuity ~ role + (1|role), prior = laplace(), data = all_merged)
But there were seven problematic observations in model 2 when using loo. I then specify:
loo(stan2, k_threshold = 0.7)
But get the following error:
Fitting model 1 out of 7 (leaving out observation 44)
Error in model.frame.default(delete.response(Terms), newdata, xlev = orig_levs) :
factor role has new level Social Worker
I’m sure this is because we only have one social worker observation. Is there any way around this or should I just be re-thinking the model instead?
This is a known missing feature in loo and kfold, that is if all observations belonging to a group are left out, the fitted model is missing a parameter for the group needed in the prediction. We have discussed about this with @jonah, and I think he had an idea how to fix this, but it’s again about limited amount of time.
loo is going to be quite difficult for that only social worker, and it’s likely that high k for that observation was also because of that. Could you group that social worker in some other group? Otherwise you would need to code this in Stan yourself, fix this in rstanarm, or wait probably some months so that we have time to fix this. The fix is certainly coming because that is needed also for leave-one-group-out functionality.
Thanks, I thought this might be the case. Things behave better when I organize into larger groups (parents, bedside, consultants, admin). I will keep an eye out for the update though!
@avehtari For stan_(g)lmer models maybe we can do what we do for
posterior_predict when there are new levels. rstanarm always estimates the posterior of an extra group-level parameter (without it entering the likelihood) so we can use it when predicting for new levels. But that’s essentially what we need here too, right?
Ok, this should be doable but will require enough coding that I don’t think I’ll have time to do it before the loo and rstanarm updates we release later this week. Definitely something to get sorted out before the subsequent releases though.
So I think this may be a bit more subtle that we originally thought.
It turns out that I had already implemented what I suggested but it actually doesn’t solve the problem here. The real issue here is not with the
(1|role) term (I’ll call this the Z part of the model to follow the lme4 convention) but with the other
role term (in the X part of the model).
I was thinking about this a bit more and discussing it with @bgoodri and I actually don’t think a model like
acuity ~ role + (1|role) (or more generally,
y ~ f + (1|f), where
f is a grouping variable) makes much sense to estimate. It certainly doesn’t make sense for a frequentist since including
f in the X part of the model assumes that under repeated sampling from the population the levels of
f are fixed while including
(1|f) in the Z part of the model allows for that to not be the case. Seems like a contradiction to have them both in the same model. From a Bayesian perspective I guess it’s not quite as bad, but I’m still not sure it makes much sense. What is the interpretation of such a model?
So, regarding the error from running
k_threshold, this is not because we’re missing the posterior of the group-level parameter (this is handled in the way I mentioned in my previous post) but because we’re missing the posterior of the so-called “fixed” parameter, i.e., the coefficient on the column of X corresponding to the social worker in the original model. And it’s not clear to me what to do in that situation, or if it makes sense to allow leaving out that observation at all.
Thank you jonah, your comments re: model specification are helpful. It’s funny that the mistake came from a model I was trying to fit as
y~ 1 + (f|g) being identified as not making sense without f in the fixed parts, so I think I understand this a little better.
I do have one question though. Really what I am trying to fit is:
y ~ acuity + role
role ~ N(u,sigma)
But I wasn’t sure how to have this make sense in the lmer syntax, i.e. what goes to the right of the ‘|’? Would it just be
y ~ acuity + (1|role)?
I am getting the same error (regarding new factor levels) when using kfold with leave one group out. The model looks like:
y~factor+x+(x|group). I use ‘group’ to generate my folds but recieve an error when ‘factor’ has a new level. I can try to post a reproducible example but seems from the discussion that the developers are aware and will be addressing this issue. Has it been addressed in the dev. version?
I am using