I'm using poisson regressions to analyze count data. I have two groups of patients in a clinical trial, and I'm comparing numbers of brain lesions that can be detected on their MRIs at 3 different assessment points. Here is an example of what the frequency tables look like: 

enter image description here

I was first looking at lesion frequency over time as a function of group, while having several other covariates in the model. Here is an example of a model I was using:  m1 <- glmer(lesions ~ Groupvisit + Agegender + DiseaseDuration + Site + (visit|id), family=poisson, data=visit1_2_3)  I am interested in both the Group*visit interaction and the main effect of Group. Because the model was not converging, I started simplifying it, and as a sanity check, I tried running separate poisson regression glms on each visit's data to see if I would get similar results for the Group effect as I get form the mixed effects models. The models produce discrepant findings with the simple glms yielding smaller errors and significant effects of Group. Below are examples of the two types of models using the same data (i.e. fitting the mixed effects model to just a single visits's data). 

enter image description here ____________________________________________ 

As you can see, the results are quite different between the two models. I understand that glmer fits the models using maximum likelihood and glm does it using least squares (IWLS), with glmer tending to produce larger error estimates. However not being a quantitative expert, I'm not able to draw any meaningful conclusions from that for the interpretation of my data. I am trying to decide which findings I should report, and I was hoping to get some well-informed opinion or advice on which type of model should be more accurate or valid in my case.  Thank you! 


1 Answer 1


In terms of reporting, your original model that includes time (visit) is much closer to what you want.

Don't bother with the individual visit models. When you specify the random effect of the intercept in your glmer using "(1|id)" you allow the intercept to shift for each individual at a lower level of analysis (your random effects) (to see the random effect of each individual use lme4::coeff(model)). It is akin to setting up a repeated measures analysis without a repeated measure because you are only looking at one visit, but that is not exactly what is happening here. If you are specifically interested in your random effect trajectories of individual patients, you would need some time variable for that, as you mentioned above. In this case much of the difference between individuals would be captured in the random effect of the intercept and we can see that in your output. Look at your intercept estimate. Ignore the p-values for a moment. The only big difference in the estimates between the models is your Intercept estimate. The others are pretty negligible depending on your scale. You also have more error in the glmer because it is using maximum likelihood estimation (MLE) so the line wont fit as well as a method that minimizes error in that data set alone (glm). The MLE model's estimates should be better at out of sample prediction. The only other thing that changed is the significance because of the change in the error.

Your model structure is interesting if you want to look at frequency or change across time. One would assume you would include an interaction with visit and treat visit as a continuous variable (depending on a few factors). Are you running inferential stats on these to get p-values or building a predictive model? What is your total number of data points and how many participants do you have? Could you also explain your variables a bit more? Are they categorical or continuous, how many levels if categorical, do you have any missing data from patients that only came for one or two visits, etc... It will be easier to help. I would assume you would want something more like this structure:

m1 <- glmer(lesions ~ visit*Age*gender*DiseaseDuration + 
           (visit|id) + (1|Site), family=poisson, data=visit1_2_3) 

I would also be careful as DiseaseDuration may have essential collinearity with visit and result in variance inflation that can cause issues with parameter estimates. But it is hard to know without more details about the variables.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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