I am attempting to build a model for the following dataset:

Level 1 Observations (Product-Level): 89000
Level 2 Observations ("BU_SBU" Department-Level): 135

Unfortunately I cannot share a sample of my data, since it is confidential.

The dependent variable in the model is a percentage (Delivery Reliability, 0-100%). Fixed effects include roughly 20 variables at level 1 and 5 variables at level 2. The only random effects are the intercepts at level 2. Having run the regression, I have a number of questions regarding the violation of model assumptions which I cannot answer myself:

  1. Constant variance of residuals: The graphic shows that there appears to be an upper- and lower-bound of the residuals. My guess is that this is due to the limitation of the dependent variable. But do the upper- and lower-bounds shown in the graphic actually indicate a violation of model assumptions? I have also run a GLMER model with a binomial(logit)-link but this did not resolve the issue. The diagnostic plots look almost identical in all three cases.

  2. Distribution of residuals: Is there a way to compute confidence intervals for residual QQ-plots of LMER models? And is it possible to compute heteroskedasticity-robust standard errors via the lme4-package?

  3. Normal distribution of level-2 intercepts: The level-2 intercepts do not appear to be Normally distributed. Is this an issue and if so, how can I resolve it?

I would greatly appreciate, if someone could help me at least with some of these questions. I am currently stuck and was not able to find any resources that provide answers. I am also grateful for recommendations to helpful literature.

enter image description here

enter image description here

enter image description here

Dropbox to Diagnostic Plots


It's hard to give a more concrete answer without actually looking at the data, but if your dependent variable is bounded between 0 and 1 I'd try a beta regression mixed model. You can try to fit it with the glmmADMB package. Here's a short example of a beta regression model with a random intercept:

data(GasolineYield)   #response variable yield is between 0 and 1
mod1<-glmmadmb(yield ~ temp + (1|batch), family="beta", data=GasolineYield)

This should help with the residuals and also produce predictions bounded within the interval of your response variable.

  • $\begingroup$ Hi Aghila, thank you for your answer. For some reason I cannot get the glmmADMB package to work. Are there other packages that allow beta regressions for hierarchical models? And, more importantly, wouldn't a linear hierarchical model with heteroskedasticity robust standard errors also work fine? And, if so, is it possible to run a linear hierarchical regression with robust standard errors in R? $\endgroup$ – Christian Mar 26 '15 at 13:10
  • $\begingroup$ If you have problems installing the package try with this: install.packages("glmmADMB", repos=c("htt p://glmmadmb.r-forge.r-project.org/repos", getOption("repos")),type="source") $\endgroup$ – Aghila Mar 26 '15 at 13:24
  • $\begingroup$ The installation worked fine. However, I get the error " The function maximizer failed (couldn't find STD file) Troubleshooting steps include (1) run with 'save.dir' set and inspect output files; (2) change run parameters: see '?admbControl' In addition: Warning message: running command 'C:\Windows\system32\cmd.exe /c "C:/Program Files/R/R-2.14.1/library/glmmADMB/bin/windows32/glmmadmb.exe" -maxfn 500 -maxph 5 -noinit -shess' had status 1 " when I include more than 5000 level-1 observations in the regression. I could not find any resources only that specify how to deal with this.Can you help? $\endgroup$ – Christian Mar 26 '15 at 15:07
  • $\begingroup$ @Christian Stackoverflow is better suited to providing troubleshooting for this type of issue. You might get better answers there. $\endgroup$ – Sycorax Apr 9 '15 at 15:43

Based on some further research I have conducted, I have derived the following options for dealing with the heteroskedasticity in this setting:

  1. Using heteroskedasticity-robust standard errors. As far as I know these are also not available in R for HLM models. However, SPSS and SAS are able to provide them.
  2. Leaving the model as it is, as HLM models with maximum likelihood estimators are relatively robust towards heteroskedasticity (see http://www.sciencedirect.com/science/article/pii/S016794730600185X and http://www.sciencedirect.com/science/article/pii/S0167947303001816).

The simulation studies also show that the other two violations of model assumptions mentioned in the questions do not seem to a big issue.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.