# Using simulated data to check when patterns in GLMM residual plots are acceptable

I have run the following Poisson GLMM:

m1 <- glmer(kereru ~ logNBT*PatchPC5yr + Urban.YN + DaysSinceNov1 + MinSinceSunrise
+ Rain + Noise + (1|PatchID) + (1|Observer) + (1|ID),
data=bird_Rescaled, family='poisson')


And have made the following validation plots to examine whether model assumptions appear to be met:

par(mfrow=c(3,3))
PearsonResids <- resid(m1, type="pearson")
plot(PearsonResids~predict(m1, type = "response"), main = " Resids vs fitted values")
plot(PearsonResids~bird_Rescaled$logNBT, main = " Resids vs Xvar") plot(PearsonResids~bird_Rescaled$PatchPC5yr, main = " Resids vs Xvar")
plot(PearsonResids~bird_Rescaled$Urban.YN, main = " Resids vs Xvar") plot(PearsonResids~bird_Rescaled$DaysSinceNov1, main = " Resids vs Xvar")
plot(PearsonResids~bird_Rescaled$MinSinceSunrise, main = " Resids vs Xvar") plot(PearsonResids~bird_Rescaled$Rain, main = " Resids vs Xvar")
plot(PearsonResids~bird_Rescaled$Noise, main = " Resids vs Xvar")  There appear to be quite a few patterns in the data. One of these patterns – the tendency for large positive residuals at low fitted values – does seem to indicate misspecification. However, I wondered if some of the other patterns are not violations of model assumptions per se, but due to the unbalanced nature of the data. For example, in the plot of residuals versus logNBT, maybe the band of points running from roughly the centre of the plot to the lower right corner just represents the large number of zeros in the response variable? In the same plot, maybe the ‘fanning out’ of points in the plot of residuals versus logNBT just occurs because there are few data points for low values of logNBT, and many data points for high values? Similarly, maybe the apparently differing variances suggested by each of the box and whisker plots results from the fact that some of my factor levels had far fewer observations than others? So, I’m wondering if a useful approach for discriminating ‘real’ versus ‘imagined’ violation of model assumptions is to (1) simulate a matching dataset (in terms of sample size for the level of each predictor, estimated influence of each predictor, etc) for which I know that the assumptions of my model have been met, (2) repeat my model and model validation plots on these simulated data, then (3) compare patterns in these plots with the patterns in the plots from the real data (and if necessary to re-simulate the data many times, to see how often I get plots with more extreme patterning than that observed for my real data). I made the following plots under this approach: # simulating data: Xvars <- subset(bird_Rescaled, select=c(logNBT, PatchPC5yr, Urban.YN, DaysSinceNov1, MinSinceSunrise, Rain, Noise, Observer, PatchID, ID)) #predictor vars mu <- predict(m1, type='response', re.form=NULL) #re.form=NULL includes random effects in predictions yObs <- rpois(lambda=mu, n=length(mu)) #simulating response var: Poisson distributed, with mean conditional on values of Xvars myDat <- cbind(yObs,Xvars) # Modelling simulated data, then checking patterns in validation plots: simMod <- glmer(yObs ~ logNBT*PatchPC5yr + Urban.YN + DaysSinceNov1 + MinSinceSunrise + Rain + Noise + (1|PatchID) + (1|Observer) + (1|ID), data=myDat, family='poisson', glmerControl(optimizer="bobyqa")) # model fails to converge otherwise par(mfrow=c(3,3)) PearsonResids <- resid(simMod, type="pearson") plot(PearsonResids~predict(simMod, type = "response"), main = " Resids vs fitted values") plot(PearsonResids~myDat$logNBT, main = " Resids vs Xvar")
plot(PearsonResids~myDat$PatchPC5yr, main = " Resids vs Xvar") plot(PearsonResids~myDat$Urban.YN, main = " Resids vs Xvar")
plot(PearsonResids~myDat$DaysSinceNov1, main = " Resids vs Xvar") plot(PearsonResids~myDat$MinSinceSunrise, main = " Resids vs Xvar")
plot(PearsonResids~myDat$Rain, main = " Resids vs Xvar") plot(PearsonResids~myDat$Noise, main = " Resids vs Xvar")


Based on this exercise it looks like many of the ‘patterns’ in the original plots do not indicate violation of model assumptions, because they also appear in the plots for which I know that assumptions are valid. However, I’m sceptical of my results! In the very least I would have expected that my simulated plots would be symmetrical around zero and would not be full of outliers. So, my question is: is the approach that I’ve used valid (is the approach valid in general, and did I simulate the data correctly), and if so, why all the asymmetry and outliers in the plots?