GLMM for Poisson ratio of means
glmm_poisson.RdGeneralized linear mixed model for two dependent Poisson outcomes.
Arguments
- data
(list)
A list whose first element is the vector of Poisson values from sample 1 and the second element is the vector of Poisson values from sample 2. Each vector must be sorted by the subject/item index and must be the same sample size. NAs are silently excluded. The default output fromsim_bnb().- test
(String:
"wald";c("wald", "lrt"))
The statistical method used for the test results.test = "wald"performs a Wald test and optionally the Wald confidence intervals.test = "lrt"performs a likelihood ratio test and optionally the profile likelihood confidence intervals (means and ratio). The Wald confidence interval is always used for the limits of the mean of sample 2 and standard deviation of the item (subject) random intercept.- ci_level
(Scalar numeric:
NULL;(0, 1))
IfNULL, confidence intervals are set asNA. If in(0, 1), confidence intervals are calculated at the specified level. Profile likelihood intervals are computationally intensive, so intervals fromtest = "lrt"may be slow.- ...
Optional arguments passed to
glmmTMB::glmmTMB().
Value
A list with the following elements:
| Slot | Subslot | Name | Description |
| 1 | chisq | \(\chi^2\) test statistic for the ratio of means. | |
| 2 | df | Degrees of freedom. | |
| 3 | p | p-value. | |
| 4 | ratio | Estimated ratio of means (sample 2 / sample 1). | |
| 4 | 1 | estimate | Point estimate. |
| 4 | 2 | lower | Confidence interval lower bound. |
| 4 | 3 | upper | Confidence interval upper bound. |
| 5 | mean1 | Estimated mean of sample 1. | |
| 5 | 1 | estimate | Point estimate. |
| 5 | 2 | lower | Confidence interval lower bound. |
| 5 | 3 | upper | Confidence interval upper bound. |
| 6 | mean2 | Estimated mean of sample 2. | |
| 6 | 1 | estimate | Point estimate. |
| 6 | 2 | lower | Wald confidence interval lower bound. |
| 6 | 3 | upper | Wald confidence interval upper bound. |
| 7 | item_sd | Estimated standard deviation of the item (subject) random intercept. | |
| 7 | 1 | estimate | Point estimate. |
| 7 | 2 | lower | Confidence interval lower bound. |
| 7 | 3 | upper | Confidence interval upper bound. |
| 8 | n1 | Sample size of sample 1. | |
| 9 | n2 | Sample size of sample 2. | |
| 10 | method | Method used for the results. | |
| 11 | test | Type of hypothesis test. | |
| 12 | alternative | The alternative hypothesis. | |
| 13 | ci_level | Confidence level of the interval. | |
| 14 | hessian | Information about the Hessian matrix. | |
| 15 | convergence | Information about convergence. |
Details
Uses glmmTMB::glmmTMB() in the form
glmmTMB(
formula = value ~ condition + (1 | item),
data = data,
family = stats::poisson
)to model dependent Poisson outcomes \(X_1 \sim \text{Poisson}(\mu)\) and \(X_2 \sim \text{Poisson}(r \mu)\) where \(\mu\) is the mean of sample 1 and \(r\) is the ratio of the means of sample 2 with respect to sample 1.
The hypotheses for the LRT and Wald test of \(r\) are
$$ \begin{aligned} H_{null} &: log(r) = 0 \\ H_{alt} &: log(r) \neq 0 \end{aligned} $$
where \(r = \frac{\bar{X}_2}{\bar{X}_1}\) is the population ratio of arithmetic means for sample 2 with respect to sample 1 and \(log(r_{null}) = 0\) assumes the population means are identical.
When simulating data from sim_bnb(), the mean is a function of the
item (subject) random effect which in turn is a function of the dispersion
parameter. Thus, glmm_poisson() has biased mean estimates. The bias
increases as the dispersion parameter gets smaller and decreases as the
dispersion parameter gets larger. However, estimates of the ratio and
standard deviation of the random intercept tend to be accurate. In summary,
the Poisson mixed-effects model fit by glmm_poisson() is not recommended
for the BNB data simulated by sim_bnb(). Instead, wald_test_bnb() or
lrt_bnb() should typically be used instead.
References
Hilbe JM (2011). Negative Binomial Regression, 2 edition. Cambridge University Press. ISBN 9780521198158 9780511973420, doi:10.1017/CBO9780511973420 .
Hilbe JM (2014). Modeling count data. Cambridge University Press, New York, NY. ISBN 9781107028333 9781107611252, doi:10.1017/CBO9781139236065 .
Examples
#----------------------------------------------------------------------------
# glmm_poisson() examples
#----------------------------------------------------------------------------
library(depower)
set.seed(1234)
d <- sim_bnb(
n = 40,
mean1 = 10,
ratio = 1.2,
dispersion = 2
)
lrt <- glmm_poisson(d, test = "lrt")
lrt
#> $chisq
#> [1] 9.020431
#>
#> $df
#> [1] 1
#>
#> $p
#> [1] 0.002669784
#>
#> $ratio
#> $ratio$estimate
#> [1] 1.227983
#>
#> $ratio$lower
#> [1] NA
#>
#> $ratio$upper
#> [1] NA
#>
#>
#> $mean1
#> $mean1$estimate
#> [1] 6.959751
#>
#> $mean1$lower
#> [1] NA
#>
#> $mean1$upper
#> [1] NA
#>
#>
#> $mean2
#> $mean2$estimate
#> [1] 8.546456
#>
#> $mean2$lower
#> [1] NA
#>
#> $mean2$upper
#> [1] NA
#>
#>
#> $item_sd
#> $item_sd$estimate
#> [1] 0.8313157
#>
#> $item_sd$lower
#> [1] NA
#>
#> $item_sd$upper
#> [1] NA
#>
#>
#> $n1
#> [1] 40
#>
#> $n2
#> [1] 40
#>
#> $method
#> [1] "GLMM for two dependent Poisson ratio of means"
#>
#> $test
#> [1] "lrt"
#>
#> $alternative
#> [1] "two.sided"
#>
#> $ci_level
#> NULL
#>
#> $hessian
#> [1] "Hessian appears to be positive definite."
#>
#> $convergence
#> [1] "relative convergence (4)"
#>
wald <- glmm_poisson(d, test = "wald", ci_level = 0.95)
wald
#> $chisq
#> [1] 8.973328
#>
#> $df
#> [1] 1
#>
#> $p
#> [1] 0.002739491
#>
#> $ratio
#> $ratio$estimate
#> [1] 1.227983
#>
#> $ratio$lower
#> [1] 1.07358
#>
#> $ratio$upper
#> [1] 1.404592
#>
#>
#> $mean1
#> $mean1$estimate
#> [1] 6.959751
#>
#> $mean1$lower
#> [1] 5.236813
#>
#> $mean1$upper
#> [1] 9.249545
#>
#>
#> $mean2
#> $mean2$estimate
#> [1] 8.546456
#>
#> $mean2$lower
#> [1] 6.451759
#>
#> $mean2$upper
#> [1] 11.32124
#>
#>
#> $item_sd
#> $item_sd$estimate
#> [1] 0.8313157
#>
#> $item_sd$lower
#> [1] 0.6371535
#>
#> $item_sd$upper
#> [1] 1.084646
#>
#>
#> $n1
#> [1] 40
#>
#> $n2
#> [1] 40
#>
#> $method
#> [1] "GLMM for two dependent Poisson ratio of means"
#>
#> $test
#> [1] "wald"
#>
#> $alternative
#> [1] "two.sided"
#>
#> $ci_level
#> [1] 0.95
#>
#> $hessian
#> [1] "Hessian appears to be positive definite."
#>
#> $convergence
#> [1] "relative convergence (4)"
#>
#----------------------------------------------------------------------------
# Compare results to manual calculation of chi-square statistic
#----------------------------------------------------------------------------
# Use the same data, but as a data frame instead of list
set.seed(1234)
d <- sim_bnb(
n = 40,
mean1 = 10,
ratio = 1.2,
dispersion = 2,
return_type = "data.frame"
)
mod_alt <- glmmTMB::glmmTMB(
formula = value ~ condition + (1 | item),
data = d,
family = stats::poisson
)
mod_null <- glmmTMB::glmmTMB(
formula = value ~ 1 + (1 | item),
data = d,
family = stats::poisson
)
lrt_chisq <- as.numeric(-2 * (logLik(mod_null) - logLik(mod_alt)))
lrt_chisq
#> [1] 9.020431
wald_chisq <- summary(mod_alt)$coefficients$cond["condition2", "z value"]^2
wald_chisq
#> [1] 8.973328
anova(mod_null, mod_alt)
#> Data: d
#> Models:
#> mod_null: value ~ 1 + (1 | item), zi=~0, disp=~1
#> mod_alt: value ~ condition + (1 | item), zi=~0, disp=~1
#> Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
#> mod_null 2 497.10 501.86 -246.55 493.10
#> mod_alt 3 490.08 497.22 -242.04 484.08 9.0204 1 0.00267 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1