How to account for incrementation in a log-linear model

Question

I want to perform a mixed regression analysis with random intercept and uncorrelated random slope after multiple imputation.

The dependent variable is continuous, namely a duration as number of days (‘days’).

A lot of observations have the value 0 as you can the following table()-output depicts:

   0 0.25  0.5 0.75    1    2    3    4    5    6    7    8    9   10   11   12   14   16   18   19   20   24   29   35   38   56 
 220    1    8    2  199  128   35    8    6    2    5    7    4    5    4    3    1    1    1    1    1    1    1    1    1    1

When I use a linear model

with(imp, lmer(days ~ binary_predictor + ... + 
                  (1 | group) + (0 + binary_predictor | group)))

I get for 20 of 20 imputations the message “singular fit”. Using qqp(days, "lnorm") shows that the probability distributions fits quite well after logarithmic transformation. However, due to the “0”-values I have to increment days by 1 for regression analysis.

with(imp, lmer(log(days + 1) ~ binary_predictor + ... + 
                  (1 | group) + (0 + binary_predictor | group)))

The resulting estimate for binary_predictoris is 0.595 with exp(0.596) = 1.813031. Unfortunately, the estimate differs if I increment days by 2. Its then exp(0.408) = 1.503807.

Is there a way to account for the incrementation of days in the regression model? Or would you suggest another regression model instead?

I tried mixed Cox regression but the assumption of proportional hazards is violated.

UPDATE #1 after the answer of Dimitris Rizopoulos

Thanks for you answer.

After I found this website, I checked the fit of the probability distributions.

days.1 <- days + 1
days.int <- ceiling(days)

From the resulting plots I figured that lnorm would be the best probability distributions to use.

I also calculated a Poisson model

with(imp, glmer(days.int ~ binary_predictor + ... + 
     (1 | group) + (0 + binary_predictor | group)), family=poisson ))

which lead to this result:

                Estimate Std.Error   t.value        df   P(>|t|)       RIV       FMI 
(Intercept)        0.167     0.190     0.879  5911.367     0.380     0.060     0.057 
binary_predictor     1.037     0.178     5.814 74393.785     0.000     0.016     0.016

Unfortunately, I do not know how to interpret the estimate of the Poisson model. Is it in days or some kind of ratio?

exp(1.037) -> 2.820742

The result of the log-linear model

with(imp, lmer(log(days + 1) ~ binary_predictor + ... + 
      (1 | group) + (0 + binary_predictor | group)))

is quite diferent

                  Estimate   Std.Error     t.value          df     P(>|t|)         RIV         FMI 
(Intercept)          0.494       0.180       2.745  551034.661       0.006       0.006       0.006 
binary_predictor       0.595       0.158       3.763 4528779.275       0.000       0.002       0.002 

with

exp(0.595) -> 1.813031

Answer 1

Your outcome variable seems to be discrete and not continuous as is the assumption behind the linear mixed model. Sometimes when you have many discrete categories (e.g., more than 10), assuming a normal distribution for this outcome (better put for the error terms of the model) may still provide a reasonable fit. However, in your case and because you have too many zeros this could be problematic.

You could consider fitting a model that is more appropriate for the nature of your outcome variable. For example a Poisson or Negative Binomial mixed model, treating your outcome as a count or collapse some categories and treat as an ordinal variable.

Both can be done with, for example, the GLMMadaptive package in R. For count data, check examples here, and for ordinal data here.