glmm_bnb | R Documentation |
Generalized linear mixed model for bivariate negative binomial outcomes.
glmm_bnb(data, test = "wald", ci_level = NULL, ...)
data |
(list) |
test |
(String: |
ci_level |
(Scalar numeric: |
... |
Optional arguments passed to |
Uses glmmTMB::glmmTMB()
in the form
glmmTMB( formula = value ~ condition + (1 | item), data = data, dispformula = ~ 1, family = nbinom2 )
to model dependent negative binomial outcomes
X_1, X_2 \sim \text{BNB}(\mu, r, \theta)
where \mu
is the mean of
sample 1, r
is the ratio of the means of sample 2 with respect to
sample 1, and \theta
is the dispersion parameter.
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_bnb()
has biased mean and dispersion 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. The p-value
for glmm_bnb()
is generally overconservative compared to glmm_poisson()
,
wald_test_bnb()
and lrt_bnb()
. In summary, the negative binomial
mixed-effects model fit by glmm_bnb()
is not recommended for the BNB data
simulated by sim_bnb()
. Instead, wald_test_bnb()
or lrt_bnb()
should
typically be used instead.
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 | Confidence interval lower bound. |
6 | 3 | upper | Confidence interval upper bound. |
7 | dispersion | Estimated dispersion. | |
7 | 1 | estimate | Point estimate. |
7 | 2 | lower | Confidence interval lower bound. |
7 | 3 | upper | Confidence interval upper bound. |
8 | item_sd | Estimated standard deviation of the item (subject) random intercept. | |
8 | 1 | estimate | Point estimate. |
8 | 2 | lower | Confidence interval lower bound. |
8 | 3 | upper | Confidence interval upper bound. |
9 | n1 | Sample size of sample 1. | |
10 | n2 | Sample size of sample 2. | |
11 | method | Method used for the results. | |
12 | test | Type of hypothesis test. | |
13 | alternative | The alternative hypothesis. | |
14 | ci_level | Confidence level of the interval. | |
15 | hessian | Information about the Hessian matrix. | |
16 | convergence | Information about convergence. |
hilbe_2011depower
\insertRefhilbe_2014depower
glmmTMB::glmmTMB()
#----------------------------------------------------------------------------
# glmm_bnb() examples
#----------------------------------------------------------------------------
library(depower)
set.seed(1234)
d <- sim_bnb(
n = 40,
mean1 = 10,
ratio = 1.2,
dispersion = 2
)
lrt <- glmm_bnb(d, test = "lrt")
lrt
wald <- glmm_bnb(d, test = "wald", ci_level = 0.95)
wald
#----------------------------------------------------------------------------
# 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,
dispformula = ~ 1,
family = glmmTMB::nbinom2,
)
mod_null <- glmmTMB::glmmTMB(
formula = value ~ 1 + (1 | item),
data = d,
dispformula = ~ 1,
family = glmmTMB::nbinom2,
)
lrt_chisq <- as.numeric(-2 * (logLik(mod_null) - logLik(mod_alt)))
lrt_chisq
wald_chisq <- summary(mod_alt)$coefficients$cond["condition2", "z value"]^2
wald_chisq
anova(mod_null, mod_alt)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.