Again, thank you all for the suggestions.
I was going with the skew-normal in the beginning as I thought it fits the data-generating process quite well. The ratings on the scale represent people’s subjective confidence about a decision. In this scale, the scale is only anchored at 0 and 200 without additional ticks/categories.
Thus, I assume that each observation of the DV represents some true confidence of a person + noise. I thought that a gaussian model would capture that idea quite well, but as I knew these confidence indications are normally negatively skewed, I decided to go with the skew normal. Shouldn’t that, at least in principle, be accurate?
@paul.buerkner: Also, I would like to avoid binning/grouping the data if possible. For example, in another experiment the values represent real monetary values that people have to “recover” on the scale rather than giving ratings in the classical sense. Thus, in that case a gaussian model assmuing observations have a true value / mean + error seems to be most plausible to me at least (but then again maybe my intuition is wrong). Moreover, I also use the spread of ratings of the same objects, i.e. their sd, for predicting other variables in a later stage of the analysis. This probably does not make a lot of sense out of context, but it is to say that the solution of fitting an ordinal model on binned data would be something I would like to avoid. That said, if you think it might be the only reason out of this misery or
I tried some additional stuff in the meantime with some interesting but confusing results again:
I fit a truncated model with an upper bound of 200
m_gauss_trunc <- brm(m_conf | trunc(ub = 200) ~ 1
+ (1|ppid)
, data = d_h1, family = "gaussian"
, chains = n_cores-1, cores = n_cores-1
)
in the pp_check the fit looks reasonable, and looic is better than for the previous models (looic = 11290):
However, for some reason I get an intercept in the summary(m_gauss_trunc) output that is around 212.5 with SE of 9.9 :
Group-Level Effects:
~ppid (Number of levels: 57)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 66.55 9.24 50.86 87.09 1151 1.01
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 212.52 9.86 194.05 233.04 622 1.01
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 36.02 1.58 33.17 39.34 6531 1.00
which seems odd as the intercept is > 200, or am I just misinterpreting was the intercept means here? If I look at pp_check(m_gauss_trunc, type="stat") I get this plot:
which does seem to recover the mean/intercept correctly. I do not really understand what is going on here. Does the intercept in this case represent something like: “if we assume the data are from a normal distribtution but have been truncated at 200, the best normal distribution to describe this would have a mean of 208”?
At least I was able to reconstruct that kind of behavior when I simulating a sort-like situation using a truncated normal with the estimates from m_gauss_trunc
set.seed(1234)
d_trunc <- as.data.frame(rtnorm(1000, 212, 36, 0, 200))
names(d_trunc)[1] <- "y"
mtest_trunc <- brm(y | trunc(ub = 200) ~ 1
, data = d_trunc
, chains = n_cores-1, cores = n_cores-1
)
Again, I get an intercept of > 200 but an accurate mean-recovery for the data.
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 212.59 7.74 200.13 229.76 1287 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 37.07 2.74 32.30 43.04 1366 1.00
So I guess, the truncated model does allow for underlying distributions with mean values that are higher than the upper bound of the truncation? In that case, I wonder whether the truncated model would make sense in this case.
Additionally, I try a truncated skew-normal model
m_SN_trunc <- brm(m_conf | trunc(ub = 200) ~ 1
+ (1|ppid)
, data = d_h1, family = "skew_normal"
, chains = n_cores-1, cores = n_cores-1
, control = list(adapt_delta = 0.999)
)
in which case I also get an intercept > 200 but cannot plot the posterior predictive check with pp_check() as I get the following error:
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Error in if (any(lower > upper)) stop("lower>upper integration limits") :
missing value where TRUE/FALSE needed
Again, no idea why this happens there.
Sorry if this is too much information, but I thought it might help if people experience similar problems in the future maybe.
But to give a short summary of the remaining questions:
-
As I see it, if I were to choose a model purely on a-priori theoretical considerations, I would go for the skew-normal model. However, that model shows some influential observations, which is of course not per se problematic but I do not exactly know where these come from compared to the gaussian model . Moreover, it does not recover the mean of the data in the posterior predictions which is at odds with the posterior predictive check of the response distribution looking better compared to the gaussian model and performing better in k-fold CV.
-
If, on the other hand, I were to choose a model based purely on best model-fit (of the models I have tried so far) in terms of looic, I would go with the truncated normal model. However, that model would result in (that is, if my above explanation is sensible) a posterior with a mean/intercept that is higher than the upper bound of the scale, i.e 200, though it does indeed seem to recover the mean in the observed as well as simulated data correctly.
@paul.buerkner Do you think any of these models make (enough) sense to go on using them? Do, in your experience, results change drastically among the different specifications? In my case, the estimates for the predictors of interest were really similar with both the gaussian and skew-normal.
Thanks very much in advance!