3
$\begingroup$

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")

model validation plots - real data 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")

validation plots on simulated data 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?

$\endgroup$
1
$\begingroup$

The idea to simulate from the fitted model to see what kind of residuals would be expected is a good one. However, there are better ways to compare the simulated to your observed residuals than eyeballing.

Those are implemented in the DHARMa R package, which uses simulations from the fitted model to calculate standardized residuals. An explanation of the approach and examples are provided in the vignette of the package.

$\endgroup$
0
$\begingroup$

An article that discusses the idea of simulating and comparing is:

Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne, D.F and Wickham, H. (2009) Statistical Inference for exploratory data analysis and model diagnostics Phil. Trans. R. Soc. A 2009 367, 4361-4383 doi: 10.1098/rsta.2009.0120

$\endgroup$

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.