3
$\begingroup$

I am working with a very unbalanced dataset that represents algae response to pollution. This dataset is the reunion of data that came from other studies. Algae are expressed as cell abundance counting, so, the response varies from 0 to more than 20000. Pollution is approached as z-scores because I needed to standardize many different variables, resulting in only one variable to compare. My random effect is related to the kind of measurement performed in the studies to compare treatment and controls: if repeated in time (measured on day1, then on day 2, etc, at the same place), if the measurements were performed in different places (polluted vs non-polluted), or if not specified.

This is what algae abundance looks like. As can be noticed, there are many zeros.

enter image description here

This is how z-scores data look like:

enter image description here

Link to dataset: https://drive.google.com/file/d/1lBrEseqDq4K0pNGp0Gvirn3J8lE-oNDu/view?usp=sharing

I need to measure the effect of pollution on algae considering the random effect, so I'm using this model:

model.abund.phyt<-glmer(response ~ z_scores+ (1|random) , data= dataset, family = "poisson")

I don't know what is going wrong with this analysis but when I run model.abund.phyt this is what I get:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']

Family: poisson ( log )

Formula: response ~ z_scores + (1 | random)

Data: dataset

 AIC      BIC   logLik deviance df.resid 

 Inf      Inf     -Inf      Inf      441 

Random effects:

Groups Name Std.Dev.

random (Intercept) 1

Number of obs: 444, groups: random, 3

Fixed Effects:

(Intercept) z_scores

  2.333        1.021  

optimizer (Nelder_Mead) convergence code: 0 (OK) ; 22820 optimizer warnings; 1 lme4 warnings

When I run the summary of the model, I get several lines of errors like this:

non-integer x = 0.210000

What do I need to do to get the model to run correctly? The numerical variables are recognized as numeric when I run str function on the dataset. Is there a problem with my data or the model specifications?

$\endgroup$
6
  • 3
    $\begingroup$ Unless your outcome values are all non-negative integers, a Poisson model can't be counted on to give sensible results. The warning suggests that you have some non-integer values. Also, 3 groups (in your "random" variable) is rather small for a mixed model. You would probably be better off treating those 3 groups as fixed effects. $\endgroup$
    – EdM
    Jul 29, 2022 at 19:13
  • $\begingroup$ How can I handle these non-integer values? Should I change the family at the model? $\endgroup$ Jul 29, 2022 at 19:15
  • 1
    $\begingroup$ If you don't have integer counts for outcomes then you should use something other than a Poisson model. The proper choice of family depends on the nature of your outcome values. Please edit your question to say more about what that "response" values represent physically, the source of the "z-scores" used as predictors, and what the random groups are (and how many). Please provide that information by editing the question, as comments are easy to overlook and can be deleted. $\endgroup$
    – EdM
    Jul 29, 2022 at 19:25
  • $\begingroup$ I edited the question and added more details. Hope it is clarified. $\endgroup$ Jul 29, 2022 at 19:32
  • $\begingroup$ Can you also explain what you mean by "if repeated in time" and "if in different places"? $\endgroup$
    – dipetkov
    Jul 29, 2022 at 19:54

2 Answers 2

3
$\begingroup$

We use random effects to encode information about structure in the data that implies that observations are not independent. For example:

  • We have measurements from different sites in the same geographic area and we expect that two measurements from the same area are more similar than measurements from two different areas (with the same values for the predictor variables). We can deal with this by adding area random effects.
  • We have repeated measurements from the same site at various points in time and we expect that measurements from the same site at two different times are more similar than measurements from different sites. We can deal with this by adding site random effects.

You don't have this level of information as you seem to be putting together a meta-dataset by combining data collected under different conditions.

It's not meaningful to conceptualize the different conditions as "random effects". As @EdM advises, it would be better that you treat them as fixed effects. This simplifies the model and it'll be easier to deal with other errors, if any.

No matter what model you choose, if you don't have information about the structure of the data (which location were measurements collected from? at what time?), you cannot model correlations between observations appropriately. And if you assume that observations are independent when they are in fact correlated, the inference from your model won't be quite right: p-values too small, confidence intervals too narrow. Be aware of this limitation and don't overinterpret the results.

$\endgroup$
1
  • $\begingroup$ Thank you very much dipetkov and @EdM for the help. I am closing the question now. $\endgroup$ Jul 30, 2022 at 13:50
2
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Thank you very much! Your response completed everything I needed. $\endgroup$ Aug 2, 2022 at 21:52

Your Answer

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

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