This follows on to the helpful answer by @dipetkov (+1) about random effects, which you rightly accepted.
Even if you treat the measurement-type groups
* as fixed effects, there are additional issues with your model and data that need your attention.
First, you need to look at your data to see if they make sense. For example, of your 444 data rows only 8 (all in the "locais_dif" group
) have response
values greater than 2500. These 8 contain exact duplicates of 4 sets of values (rows 5 to 8 of your data are identical to rows 11 to 14). Overall, 93 of your rows (62 with non-0 response
values) exactly duplicate other rows. The extreme values and duplicates might be correct, but it's important to start with validated data before you model.
After you ensure that your data are clean, look at them directly. With ggplot2
you could do something like
ggplot(data=subset(algae[!duplicated(algae[,]),],
subset=response<2500),
mapping=aes(y=response,x=z_scores,color=group)) +
geom_point() + geom_smooth()
where for now I removed the duplicated rows and the 8 extreme outliers to see how responses depend jointly on z_scores
and your measurement group
. (This gives a lot of warnings probably due to the large number of 0 response
values, but there is a plot.)
Second, if you just include a "main effect" for each of your 3 measurement-type groups
, you're assuming that they have different baseline response
values but that the slope of the relationship between response
and z_scores
is the same for all 3 groups. That isn't always a good assumption, and it seems to be violated in your data (as my suggested plot above will show). (That's also the underlying assumption in using random intercepts in a mixed model.) You should include an interaction between z_scores
and the groups
. (If a mixed model were appropriate, that would mean adding a random slope to the model.)
Third, your model assumes a strictly linear association between response
and z_scores
. Such a simple assumption seldom holds, and it doesn't seem to hold well in your data. It makes sense to try flexible fitting of continuous predictors, for example with regression splines.
Fourth, your "response" values in your data aren't count values. The non-zero values all seem to have fractional parts. The Poisson distribution is for non-negative integer counts. You can still get a model to fit, but the likelihood calculations needed to get things like AIC can't be done properly. That's why you got the warnings that originally attracted your attention. If the response
values were based on calculations from integer counts and you want to use a Poisson model, you should go back to model the original counts instead.
Fifth, if you specify family = "quasipoisson"
instead, you can get a fit without any warnings. The method used to fit doesn't use maximum likelihood and thus also doesn't give AIC. (I think it's the method that your Poisson model used when it found non-integer response values.) The quasi-Poisson model estimates a "dispersion parameter": how much greater the variance is versus the Poisson assumption that the variance equals the mean (dispersion parameter of 1). When I tried that on your data with an interaction between group
and linear z_scores
, the estimated dispersion parameter was over 2700! Even after removing the 8 apparent outliers the dispersion parameter was 383. That brings into question whether you should be trying to force your data into a (quasi) Poisson framework. A model that makes no assumptions about the variance of response
values, like ordinal logistic regression, could be a better choice. This page discusses other ways to approach non-negative outcome data with a lot of zeros.
Finally, be careful about how you determined the z_scores
you use as predictor values. If you determined those based on your knowledge of the subject matter independently of the response
values in these data, that's fine. But if you used these response
data in some way to determine the z_scores
then your estimates of "statistical significance" will be invalid and your model might not extend well to new data.
*To remove ambiguity, I've renamed your "random" column describing the type of measurement to be "group" instead.