2

I am conducting a GLMM for a meta-analysis using the beta distribution with the package glmmTMB. My response variable is a vector of correlations (No exact 0 or 1), but Fisher’s transformation fails to normalize the data so this is why I use the beta distribution. Because this is a meta-analysis, I have included sample size as a weights argument. Below is the code:

model <- glmmTMB(Y ~ unit*class + (1|study/pair), weights = weights, family=list(family="beta",link="logit"), data = data)

where unit is a factor with two levels, class is a factor with four levels and pair is a random effect nested within study. Each study measured correlation for two traits, forming a “pair” within each “study”. Now, I am unsure as to the best approach for generating residuals for diagnostic purposes. Using the residuals function;

plot(y = residuals(model), fitted(model, type = "response"))

yields the following graph. I am however unsure if this approach is appropriate, given the general uncertainty about residuals of beta-distributed models.If it is appropriate, I am equally unsure as whether they violate the assumption of homogeneity or not, as there seems to be a downward trend. I assume the two parallel bounds are from the beta distribution, which is fine as far as I understand it. However, I do not know these models enough to be comfortable with diagnosing the residuals they provide.

Residual plot

Below is the plot of observed values against fitted values, which seems good to me (correlation of 0.776).

Observed vs fitted

Alternatively, using DHARMa seems promising, but the simulate.glmmTMB() function does not seem to be implemented for the beta family.

I have also consulted the following papers, suggested here: Similar question on beta residuals, but I have no idea how they could be implemented in R.

Pereira (2017). On quantile residuals in beta regression. Communications in Statistics - Simulation and Computation. doi: doi.org/10.1080/03610918.2017.1381740. See https://arxiv.org/abs/1704.02917 for a pre-print.

Espinheira, Santos, & Cribari-Neto (2017). On nonlinear beta regression residuals. Biometrical Journal, 59. doi: 10.1002/bimj.201600136

Espinheira, Ferrari, Cribari-Neto (2008). On beta regression residuals. Journal of Applied Statistics, 35. doi: 10.1080/02664760701834931

Espnheira, Ferrari, & Cribari-Neto (2008). Influence diagnostics in beta regression. Computational Statistics and Data Analysis, 52. doi: 10.1016/j.csda.2008.02.028

Would anyone be able to suggest the best approach to generate the appropriate residuals for this type of model , for diagnostic purposes, and how to implement it in R? Or whether the ones I currently have are appropriate?

Jean-Michel Matte
  • 21
  • 3
  • if you implement the simulation yourself, you can still use DHARMa, see example here https://github.com/florianhartig/DHARMa/wiki/Adding-new-R-packages-to-DHARMA – Florian Hartig Jun 08 '18 at 15:49
  • I want to implement the code under the section "Quick-and-dirty handcoded simulations" of the link you posted, but the example provided is a glm only and I have a few follow-up questions as to how to implement it with my current model. 1) How would you simulate the variance with the nested study/pair random effect; 2) How would you implement the weight argument? 3) And finally, how do you implement the two shape parameters in the rbeta() function (which is what I assume I would need to use for that distribution)? – Jean-Michel Matte Jun 28 '18 at 17:22
  • 1
    shape parameters have to be extracted from the fitted model, REs in the usual way (normal distr. on linear predictor), or just test conditional on the fitted REs ... weights is tricky, I would just ignore the weights. it's not perfect, but if you ignore weights and calculate conditional on fitted REs, you could just calculate the quantiles of the beta without simulation – Florian Hartig Jun 29 '18 at 16:54

0 Answers0