Generalized linear mixed models are generalized linear models that include random effects varying between clusters or 'higher-level' units of hierarchically structured data. Such models can be estimated using gllamm. The prediction command gllapred can be used to obtain empirical Bayes predictions of the random effects, interpretable as higher-level residuals. Combined with approximate sampling standard deviations, these residuals can be used for identifying unusual higher-level units. However, since the distribution of these predictions is generally not known, we recommend simulating responses from the model using gllasim and comparing 'observed' and simulated residuals. We also discuss different types of level 1 residuals and influence diagnostics.