I have some questions on performing mixed models on multi-rater data when residuals are heteroskedastic. I've found some of the info on Cross Validated confusing and quite technical-- would be very grateful for some pragmatic tips to help with model selection.
My original mixed model formula, written in R for lme4::lmer was:
y ~ x + sex + source + x:sex + x:source + x:sex:source + (1 | ID)
…where “y” is score on a questionnaire (continuous var ranging from 0 to 40), “source” is a within-person binary variable indicating who answered the questionnaire (self- or parent-report) and “x” is my primary variable of interest (continuous var, standardized). “x” varies drastically by sex (it is a hormone measure). The random intercept by “ID” allows the two observations (own + parents) to be considered non-independant measures, nested within each subject ID (n=90). Thus, there should be 180 observations of the outcome (2 per subject).
The distribution of the outcome “y” is positively skewed (though I realise it’s the normality of residuals we’re interested in)…
The model fitted values (x-axis) vs residuals (y-axis) look like this:
Heteroskedasticity can be obseved. Granted, it’s not by a large order of magnitude (max residual = ~2) but presumably its the visible pattern in residuals along predicted values of the outcome that’s the problem.
Say I want to improve the model fit but do not want to log-transform the outcome variable due to the complications it poses to interpretation. As far as I can make out, this leaves me with Generalized Linear Mixed Models, where I can change the underlying distribution & link function.
Distribution: Based on the observed distribution of the outcome (see density plot above), I decided that a gamma distibution could be a better fit than gaussian. Note: Outcome “y” was standardized to have mean 2 and SD 1, such that no 0-value existed.
Link Function: Based on the fact that log-transforming the outcome results in a more normal distribution, I thought that it may be a better link function than the identity link.
I experimented with the 4 permutations of gaussian/gamma distributions & identity/log link function, using the glmer() function instead of lmer() to model 3 of those permutations (gaussian-identity modelled using lmer). Example of glmer syntax:
glmer(y ~ x + sex + source + x:sex + x:source + x:sex:source + (1 | ID) data = d, family = Gamma(link=log), control = glmerControl(optimizer="bobyqa", optCtrl = list(maxfun= 100000)), nAGQ = 20)
My main issue: I don’t understand how the distribution & link function componants interact and how I should make my decision about which combo is best for my data. Should I choose purely based on observation of residuals? If so, (see residual vs fitted & QQ plots below) I would probably choose the gamma distribution with identity link as this gives me the smallest residuals with no heteroskedasticity.
￼ … Or should I base my decision on BIC/AIC values (see below)? If so, it would be the gamma distribution with the log link function that gives me the smallest absolute AIC & BIC (but see how the residuals above look oddly grouped for this distrib/link)…
One final issue (which may or may not be relevant) is that my main variable of interest “x” is strongly correlated with sex (r ~ 0.85), leading to very high correlation estimates between fixed effects (see correlation table below). Is this relevant to the model fit? We could spit the analyses by sex if this were a problem.