I want to analyze the impact of the rain on smoking probability. I observed people in two cities on the streets and marked the following parameters: city, gender of the person, duration of observing person, how many times (s)he smokes, was it rain during observation or not. To deal with different length of observations I made the new parameter (smokes per hour (smkHour = smokes/duration)) and used the following code in R:

Model <- lmer (smkHour ~ rain + (1|city/gender))

However, there are many zeros in the smkHour. BoxCox transformation is not an option here and I found that it was possible to analyze first zeros versus non-zeros and a binomial model and after analyse non-zeros with log-transformation. Therefore, the first question:

Is it correct to analyze this data in that way?

If so, probability of zeros depends on duration of the observation (probability of smoke depends on the time of observing person). Thus, the second question:

How to take into account in the first model the duration of the observation?

Is it possible to include it to the error effect or not? Could it be like that?

smkHour2 - binomial zero/non-zero transformation of smkHour

smkHourNZ - only non-zeros of smkHour

model1 <- glmer (smkHour2 ~ rain + (duration|city/gender), family=binomial)
model2 <- lmer (log(smkHourNZ) ~ rain + (1|city/gender)

The first model doesn't respect the constraint that rates must be non-negative. It also assumes a Gaussian conditional distribution for smoke rate, given the regressors. But, your calculated rates were obtained from counts, and will be quantized. Quantization will be especially pronounced because counts are low, and a Gaussian distribution isn't a good approximation in this regime. Finally, the first model treats the calculated rates as ground truth, when in fact they're estimates. And, these estimates have greater uncertainty for subjects that were observed for shorter periods of time. Neglecting this uncertainty means short observations will have undue influence on estimated model parameters.

To address these concerns, an alternative approch is to model counts directly. Let $x_i$ be a vector containing regressors for the $i$th subject, who smoked $y_i$ times over an interval of length $t_i$. Suppose we postulate a linear relationship between the regressors and $r_i$, the smoking rate for this subject:

$$r_i = \beta \cdot x_i$$

where $\beta$ is a vector of coefficients (common to all subjects). Multiplying the smoking rate by the length of the observation interval gives the expected number of times the subject smoked during this interval:

$$\lambda_i = r_i t_i$$

The expected smoke count is then used to produce a conditional distribution over the number of smokes in the observed interval. For example, a Poisson distribution with mean $\lambda_i$. This can be written as:

$$p(y_i \mid x_i, t_i, \beta) = \text{Poiss}\big( y_i \mid (\beta \cdot x_i) t_i \big)$$

Assuming smoke counts are conditionally independent given the regressors, the log likelihood for all $n$ subjects is:

$$\mathcal{L}(\beta) = \sum_i \log p(y_i \mid x_i, t_i, \beta)$$

Maximum likelihood estimates of the coefficients $\beta$ can be obtained by maximizing $\mathcal{L}(\beta)$.

This is equivalent to performing Poisson regression using the observed smoke counts as targets and a transformed set of regressors. The transformed regressors are obtained by multiplying the original regressors for each subject by the time that subject was observed. The equivalence follows from the fact that $(\beta \cdot x_i) t_i = \beta \cdot (t_i x_i)$. This trick won't work for nonlinear models.

If the Poisson distribution is too restrictive, you could consider zero-inflated Poisson or negative binomial models.

  • $\begingroup$ Thank you for your answer. I think it would be better to use zero-inflated negative binomial model. I tried glmmTMB package from R and it works. $\endgroup$ – Ivan Oct 18 at 11:18

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.