2
$\begingroup$

I have a dataset of 2,000 subjects, each providing daily social contact data for one week (14,000 observations = 2000 subjects * 7 days). I aim to describe the number of contacts using variables such as the subject's age, household size, day of the week, and potentially other factors (e.g., region of residence or occupation).

The outcome variable, the daily number of contacts per capita, is overdispersed, so I wanted to try a generalized linear model assuming negative binomial distribution. However, the data includes repeated measures for the same subject (i.e., 7 observations per subject), which violates the independence assumption. So I tried a generalized linear mixed models (glmer.nb in lme4 package in R) but I never got the optimizing algorithm to converge despite adjusting various control parameters. I am looking for alternative solutions to the issues without resorting to the generalized linear mixed modeling.

Here is my question: Would the following approach be reasonable? I sample from the original data such that only one observation per subject exists by uniformly randomly selecting a day of the week for each subject. This reduces the dataset from 14,000 to 2,000 observations. I then perform negative binomial regression on this subset and repeat this process multiple times (e.g., 20,000 times).

As an exercise, I performed a negative binomial regression on 14000 observations (ignoring the violation of independence assumptions of repeated measures) and also ran negative binomial regression on 4,000 samples of 2,000 observations. The estimates (95% confidence interval) were similar across both approaches though the standard errors were wider for the sampled subsets (below). confidence intervals for the estimates_2000 were taken as the central 95% percentile from 4,000 simulations.

var estimates_14000 estimates_2000
1 (Intercept) 3.05 (2.78 - 3.34) 3.04 (2.51 - 3.66)
2 age_grp5-9 1.67 (1.53 - 1.82) 1.66 (1.41 - 1.96)
3 age_grp10-14 1.54 (1.41 - 1.68) 1.54 (1.34 - 1.79)
4 age_grp15-19 1.45 (1.33 - 1.59) 1.45 (1.25 - 1.69)
5 age_grp20-29 0.60 (0.55 - 0.65) 0.60 (0.52 - 0.69)
6 age_grp30-39 0.69 (0.63 - 0.74) 0.69 (0.60 - 0.79)
7 age_grp40-49 0.75 (0.69 - 0.81) 0.75 (0.66 - 0.86)
8 age_grp50-59 0.82 (0.75 - 0.88) 0.82 (0.72 - 0.94)
9 age_grp60-69 1.05 (0.97 - 1.13) 1.04 (0.92 - 1.21)
10 age_grp70-79 1.11 (1.02 - 1.21) 1.11 (0.96 - 1.29)
11 sexM 1.00 (0.98 - 1.02) 1.00 (0.96 - 1.04)
12 dayofweekMonday 1.47 (1.41 - 1.54) 1.47 (1.28 - 1.68)
13 dayofweekTuesday 1.46 (1.39 - 1.52) 1.46 (1.26 - 1.67)
14 dayofweekWednesday 1.74 (1.67 - 1.82) 1.75 (1.52 - 2.00)
15 dayofweekThursday 1.58 (1.51 - 1.66) 1.58 (1.38 - 1.82)
16 dayofweekFriday 1.55 (1.48 - 1.63) 1.56 (1.35 - 1.78)
17 dayofweekSaturday 1.13 (1.08 - 1.19) 1.14 (0.98 - 1.31)
18 hhsize_grp2 1.15 (1.10 - 1.20) 1.15 (1.05 - 1.26)
19 hhsize_grp3 1.27 (1.22 - 1.33) 1.27 (1.16 - 1.39)
20 hhsize_grp4 1.47 (1.40 - 1.53) 1.47 (1.34 - 1.60)
21 hhsize_grp5+ 1.66 (1.57 - 1.74) 1.65 (1.48 - 1.85)
Improve this question
$\endgroup$
3
  • $\begingroup$ This seems reasonable to me, although I don't base that on anything more than intuition. But, if you have age in years, I would use that. Categorizing a continuous variable is rarely a good idea. $\endgroup$
    – Peter Flom
    Commented Jun 5 at 10:42
  • $\begingroup$ I think something that may be helpful is to show how you fit the model to determine how it may have failed to converge (along with important information like how much random effect variance was estimable). Otherwise, one could potentially use a Bayesian NB model with priors which allow the model to converge. $\endgroup$
    – Shawn Hemelstrand
    Commented Jun 5 at 11:22
  • $\begingroup$ Thank you for shraing your opinions and suggestions! I am trying a Bayesian and it seems to work well. $\endgroup$
    – Jong-Hoon Kim
    Commented Jun 6 at 5:09

1 Answer 1

Reset to default
1
$\begingroup$

If it's just overdispersed between people, you could use a Poisson mixed effects model with a (normally-distributed) random subject effect on the intercept (on the log-rate scale). For a single observation that behaves very similarly to Negative Binomial (log-normal dist. can approximate gamma dist. well).

Alternatively, if it's convergence issues with maximum likelhood estimation, you can consider doing a Bayesian model, which e.g. the brms package makes easy (it has family=negbinomial() option).

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you your comments! I didn't realize that I could have used poisson regression, which in fact easily solves a algorithm convergence issue. I am also testing a brms package and it seems to work well, too. $\endgroup$
    – Jong-Hoon Kim
    Commented Jun 6 at 5:10
  • $\begingroup$ Just to point out that only a Poisson with a random effect (to exactly get a negative binomial, if would have to be a gamma-distributed random effect on the Poisson-rate, or to do something very similar a normally distributed random effect on the intercept on the log-rate scale of a Poisson regression) captures overdispersion. $\endgroup$
    – Björn
    Commented Jun 7 at 12:58
  • 1
    $\begingroup$ Thank you for the clarification! I used the following codes, glmer(contacts ~ age_grp + sex + dayofweek + hhsize_grp + (1 |id), weights=wt, data=dat, family=poisson(link = "log")). I think this corresponds to the latter case of the two cases you listed. $\endgroup$
    – Jong-Hoon Kim
    Commented Jun 7 at 20:10

Your Answer

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

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