GLM for NB ratio of means

Generalized linear model for two independent negative binomial outcomes.

Usage

glm_nb(data, equal_dispersion = FALSE, test = "wald", ci_level = NULL, ...)

Arguments

data: (list)
A list whose first element is the vector of negative binomial values from group 1 and the second element is the vector of negative binomial values from group 2. NAs are silently excluded. The default output from sim_nb().
equal_dispersion: (Scalar logical: FALSE)
If TRUE, the model is fit assuming both groups have the same population dispersion parameter. If FALSE (default), the model is fit assuming different dispersions.
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. Wald confidence intervals are always returned for the limits of the mean of sample 2 and, when equal_dispersion=FALSE, for the limits of the dispersion of sample 2.
ci_level: (Scalar numeric: NULL; (0, 1))
If NULL, confidence intervals are set as NA. If in (0, 1), confidence intervals are calculated at the specified level. Profile likelihood intervals are computationally intensive, so intervals from test = "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 (group 2 / group 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 group 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 group 2.
6	1	`estimate`	Point estimate.
6	2	`lower`	Wald confidence interval lower bound.
6	3	`upper`	Wald confidence interval upper bound.
7		`dispersion1`	Estimated dispersion of group 1.
7	1	`estimate`	Point estimate.
7	2	`lower`	Confidence interval lower bound.
7	3	`upper`	Confidence interval upper bound.
8		`dispersion2`	Estimated dispersion of group 2.
8	1	`estimate`	Point estimate.
8	2	`lower`	Wald confidence interval lower bound.
8	3	`upper`	Wald confidence interval upper bound.
9		`n1`	Sample size of group 1.
10		`n2`	Sample size of group 2.
11		`method`	Method used for the results.
12		`test`	Type of hypothesis test.
13		`alternative`	The alternative hypothesis.
14		`equal_dispersion`	Whether or not equal dispersions were assumed.
15		`ci_level`	Confidence level of the intervals.
16		`hessian`	Information about the Hessian matrix.
17		`convergence`	Information about convergence.

Details

Uses glmmTMB::glmmTMB() in the form

glmmTMB(
  formula = value ~ condition,
  data = data,
  dispformula = ~ condition,
  family = nbinom2
)

to model independent negative binomial outcomes $X_1 \sim \text{NB}(\mu, \theta_1)$ and $X_2 \sim \text{NB}(r\mu, \theta_2)$ where $\mu$ is the mean of group 1, $r$ is the ratio of the means of group 2 with respect to group 1, $\theta_1$ is the dispersion parameter of group 1, and $\theta_2$ is the dispersion parameter of group 2.

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 group 2 with respect to group 1 and $log(r_{null}) = 0$ assumes the population means are identical.

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

#----------------------------------------------------------------------------
# glm_nb() examples
#----------------------------------------------------------------------------
library(depower)

set.seed(1234)
d <- sim_nb(
  n1 = 60,
  n2 = 40,
  mean1 = 10,
  ratio = 1.5,
  dispersion1 = 2,
  dispersion2 = 8
)

lrt <- glm_nb(d, equal_dispersion = FALSE, test = "lrt", ci_level = 0.95)
lrt
#> $chisq
#> [1] 10.02464
#> 
#> $df
#> [1] 1
#> 
#> $p
#> [1] 0.001544595
#> 
#> $ratio
#> $ratio$estimate
#> [1] 1.542934
#> 
#> $ratio$lower
#> [1] 1.189824
#> 
#> $ratio$upper
#> [1] 1.982785
#> 
#> 
#> $mean1
#> $mean1$estimate
#> [1] 9.31667
#> 
#> $mean1$lower
#> [1] 7.493765
#> 
#> $mean1$upper
#> [1] 11.72399
#> 
#> 
#> $mean2
#> $mean2$estimate
#> [1] 14.37501
#> 
#> $mean2$lower
#> [1] 12.70001
#> 
#> $mean2$upper
#> [1] 16.27092
#> 
#> 
#> $dispersion1
#> $dispersion1$estimate
#> [1] 1.545422
#> 
#> $dispersion1$lower
#> [1] 1.005316
#> 
#> $dispersion1$upper
#> [1] 2.3586
#> 
#> 
#> $dispersion2
#> $dispersion2$estimate
#> [1] 11.07998
#> 
#> $dispersion2$lower
#> [1] 4.980081
#> 
#> $dispersion2$upper
#> [1] 24.65142
#> 
#> 
#> $n1
#> [1] 60
#> 
#> $n2
#> [1] 40
#> 
#> $method
#> [1] "GLM for independent negative binomial ratio of means"
#> 
#> $test
#> [1] "lrt"
#> 
#> $alternative
#> [1] "two.sided"
#> 
#> $equal_dispersion
#> [1] FALSE
#> 
#> $ci_level
#> [1] 0.95
#> 
#> $hessian
#> [1] "Hessian appears to be positive definite."
#> 
#> $convergence
#> [1] "relative convergence (4)"
#> 

wald <- glm_nb(d, equal_dispersion = FALSE, test = "wald", ci_level = 0.95)
wald
#> $chisq
#> [1] 11.3516
#> 
#> $df
#> [1] 1
#> 
#> $p
#> [1] 0.0007538306
#> 
#> $ratio
#> $ratio$estimate
#> [1] 1.542934
#> 
#> $ratio$lower
#> [1] 1.198893
#> 
#> $ratio$upper
#> [1] 1.985703
#> 
#> 
#> $mean1
#> $mean1$estimate
#> [1] 9.31667
#> 
#> $mean1$lower
#> [1] 7.478497
#> 
#> $mean1$upper
#> [1] 11.60666
#> 
#> 
#> $mean2
#> $mean2$estimate
#> [1] 14.37501
#> 
#> $mean2$lower
#> [1] 12.70001
#> 
#> $mean2$upper
#> [1] 16.27092
#> 
#> 
#> $dispersion1
#> $dispersion1$estimate
#> [1] 1.545422
#> 
#> $dispersion1$lower
#> [1] 1.010286
#> 
#> $dispersion1$upper
#> [1] 2.364013
#> 
#> 
#> $dispersion2
#> $dispersion2$estimate
#> [1] 11.07998
#> 
#> $dispersion2$lower
#> [1] 4.980081
#> 
#> $dispersion2$upper
#> [1] 24.65142
#> 
#> 
#> $n1
#> [1] 60
#> 
#> $n2
#> [1] 40
#> 
#> $method
#> [1] "GLM for independent negative binomial ratio of means"
#> 
#> $test
#> [1] "wald"
#> 
#> $alternative
#> [1] "two.sided"
#> 
#> $equal_dispersion
#> [1] FALSE
#> 
#> $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)
df <- sim_nb(
  n1 = 60,
  n2 = 40,
  mean1 = 10,
  ratio = 1.5,
  dispersion1 = 2,
  dispersion2 = 8,
  return_type = "data.frame"
)

mod_alt <- glmmTMB::glmmTMB(
  formula = value ~ condition,
  data = df,
  dispformula = ~ condition,
  family = glmmTMB::nbinom2
)
mod_null <- glmmTMB::glmmTMB(
  formula = value ~ 1,
  data = df,
  dispformula = ~ condition,
  family = glmmTMB::nbinom2
)

lrt_chisq <- as.numeric(-2 * (logLik(mod_null) - logLik(mod_alt)))
lrt_chisq
#> [1] 10.02464
wald_chisq <- summary(mod_alt)$coefficients$cond["condition2", "z value"]^2
wald_chisq
#> [1] 11.3516

anova(mod_null, mod_alt)
#> Data: df
#> Models:
#> mod_null: value ~ 1, zi=~0, disp=~condition
#> mod_alt: value ~ condition, zi=~0, disp=~condition
#>          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
#> mod_null  3 657.69 665.50 -325.84   651.69                            
#> mod_alt   4 649.66 660.08 -320.83   641.66 10.025      1   0.001545 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#----------------------------------------------------------------------------
# Compare results to wald_test_nb()
#----------------------------------------------------------------------------
wald2 <- wald_test_nb(d, equal_dispersion = FALSE, ci_level = 0.95)
all.equal(wald$chisq, wald2$chisq, tolerance = 0.01)
#> [1] TRUE