1
$\begingroup$

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.

Normal distribution

Log normmal distribution

Poisson distribution

Negative binomial distribution

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
$\endgroup$
  • $\begingroup$ One technique that can sometime be used is to replace zero values with very small values, such as 1.0E-10, which will allow logs to be used. As a practical matter, values such as 1.0E-10 days are effectively zero. $\endgroup$ – James Phillips Mar 24 at 11:26
  • $\begingroup$ From the Poisson model the coefficients have an intepretation for the log expected number of day. Hence, if you expontiate them you get an intepretation for the expected number of days. With regard to the fit of the distribution to your data, check if you need to account for over-dispersion and/or extra zeros. Examples are in the links I provided in my answer. $\endgroup$ – Dimitris Rizopoulos Mar 24 at 18:08
1
$\begingroup$

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.

$\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.