Comparing generalized linear mixed models (varying the distribution & link function)

Question

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.

Thank you.

Answer 1

With the models you have shown, I agree that the Gamma w/ identity link looks most promising. However, I'd like to suggest another option for you to consider.

If you would like to preserve the linear interpretation and deal with non-normality (i.e., heteroskedascitiy) in the residuals, you could use so-called robust estimators of the variance-covariance matrix. This matrix is used to calculate standard errors of your fixed effect coefficients. In short, these methods are designed to deal with the nuisance of non-normality in the residuals without explicitly modeling it. See this very informative article by Pek et al. for an introduction to this topic.

The only challenge is that these are hard to get in mixed models in R. However, since you have a simple random intercept model without random slopes, you can use the plm() package designed for panel data modeling to get the appropriate standard errors. See code below for "HC1" standard errors (you can also get "HC2", "HC3", etc. versions). Here I use the sleepstudy data included in lme4:

library(lme4)
library(lmtest)
library(plm)
plm_random <- plm(Reaction ~ 1 + Days , data = sleepstudy, 
                  effect = "individual", index = "Subject",
                  model = "random") # note that index is your grouping factor
coeftest(plm_random, vcov.=function(x) vcovHC(x, method="arellano", type="HC1"))

If you want to add random slopes, you cannot use plm and will have to look into other options. This paper on the package robustlmm is helpful and lays out some of the options.

Regarding the issue of sex, you could estimate the models separately by sex, however you should consider whether this will adversely impact your sample size.