Package 'depower'

Description

Generalized linear model for two independent negative binomial outcomes.

Usage

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

Arguments

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.

Value

A list with the following elements:

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.

See Also

glmmTMB::glmmTMB()

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

wald <- glm_nb(d, equal_dispersion = FALSE, 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_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 = d,
  dispformula = ~ condition,
  family = glmmTMB::nbinom2,
)
mod_null <- glmmTMB::glmmTMB(
  formula = value ~ 1,
  data = d,
  dispformula = ~ condition,
  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)

Description

Generalized linear mixed model for bivariate negative binomial outcomes.

Usage

glmm_bnb(data, test = "wald", ci_level = NULL, ...)

Arguments

Details

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.

Value

A list with the following elements:

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.

See Also

glmmTMB::glmmTMB()

Examples

#----------------------------------------------------------------------------
# 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)

Description

Generalized linear mixed model for two dependent Poisson outcomes.

Usage

glmm_poisson(data, test = "wald", ci_level = NULL, ...)

Arguments

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.

Value

A list with the following elements:

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.

See Also

glmmTMB::glmmTMB()

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

wald <- glmm_poisson(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,
  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
wald_chisq <- summary(mod_alt)$coefficients$cond["condition2", "z value"]^2
wald_chisq

anova(mod_null, mod_alt)

Description

Likelihood ratio test for the ratio of means from bivariate negative binomial outcomes.

Usage

lrt_bnb(data, ratio_null = 1, ...)

Arguments

Details

This function is primarily designed for speed in simulation. Missing values are silently excluded.

Suppose $X_1 \mid G = g \sim \text{Poisson}(\mu g)$ and $X_2 \mid G = g \sim \text{Poisson}(r \mu g)$ where $G \sim \text{Gamma}(\theta, \theta^{-1})$ is the random item (subject) effect. Then $X_1, X_2 \sim \text{BNB}(\mu, r, \theta)$ is the joint distribution where $X_1$ and $X_2$ are dependent (though conditionally independent), $X_1$ is the count outcome for sample 1 of the items (subjects), $X_2$ is the count outcome for sample 2 of the items (subjects), $\mu$ is the conditional mean of sample 1, $r$ is the ratio of the conditional means of sample 2 with respect to sample 1, and $\theta$ is the gamma distribution shape parameter which controls the dispersion and the correlation between sample 1 and 2.

The hypotheses for the LRT of $r$ are

$\begin{aligned} H_{null} &: r = r_{null} \\ H_{alt} &: r \neq r_{null} \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 $r_{null}$ is a constant for the assumed null population ratio of means (typically $r_{null} = 1$).

The LRT statistic is

$\begin{aligned} \lambda &= -2 \ln \frac{\text{sup}_{\Theta_{null}} L(r, \mu, \theta)}{\text{sup}_{\Theta} L(r, \mu, \theta)} \\ &= -2 \left[ \ln \text{sup}_{\Theta_{null}} L(r, \mu, \theta) - \ln \text{sup}_{\Theta} L(r, \mu, \theta) \right] \\ &= -2(l(r_{null}, \tilde{\mu}, \tilde{\theta}) - l(\hat{r}, \hat{\mu}, \hat{\theta})) \end{aligned}$

Under $H_{null}$, the LRT test statistic is asymptotically distributed as $\chi^2_1$. The approximate level $\alpha$ test rejects $H_{null}$ if $\lambda \geq \chi^2_1(1 - \alpha)$. Note that the asymptotic critical value is known to underestimate the exact critical value. Hence, the nominal significance level may not be achieved for small sample sizes (possibly $n \leq 10$ or $n \leq 50$).

Value

A list with the following elements:

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

Examples

#----------------------------------------------------------------------------
# lrt_bnb() examples
#----------------------------------------------------------------------------
library(depower)

set.seed(1234)
sim_bnb(
  n = 40,
  mean1 = 10,
  ratio = 1.2,
  dispersion = 2
) |>
  lrt_bnb()

Description

Likelihood ratio test for the ratio of means from two independent negative binomial outcomes.

Usage

lrt_nb(data, equal_dispersion = FALSE, ratio_null = 1, ...)

Arguments

Details

This function is primarily designed for speed in simulation. Missing values are silently excluded.

Suppose $X_1 \sim NB(\mu, \theta_1)$ and $X_2 \sim NB(r\mu, \theta_2)$ where $X_1$ and $X_2$ are independent, $X_1$ is the count outcome for items in group 1, $X_2$ is the count outcome for items in group 2, $\mu$ is the arithmetic mean count in group 1, $r$ is the ratio of arithmetic means for 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 of $r$ are

$\begin{aligned} H_{null} &: r = r_{null} \\ H_{alt} &: r \neq r_{null} \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 $r_{null}$ is a constant for the assumed null population ratio of means (typically $r_{null} = 1$).

The LRT statistic is

$\begin{aligned} \lambda &= -2 \ln \frac{\text{sup}_{\Theta_{null}} L(r, \mu, \theta_1, \theta_2)}{\text{sup}_{\Theta} L(r, \mu, \theta_1, \theta_2)} \\ &= -2 \left[ \ln \text{sup}_{\Theta_{null}} L(r, \mu, \theta_1, \theta_2) - \ln \text{sup}_{\Theta} L(r, \mu, \theta_1, \theta_2) \right] \\ &= -2(l(r_{null}, \tilde{\mu}, \tilde{\theta}_1, \tilde{\theta}_2) - l(\hat{r}, \hat{\mu}, \hat{\theta}_1, \hat{\theta}_2)) \end{aligned}$

Under $H_{null}$, the LRT test statistic is asymptotically distributed as $\chi^2_1$. The approximate level $\alpha$ test rejects $H_{null}$ if $\lambda \geq \chi^2_1(1 - \alpha)$. Note that the asymptotic critical value is known to underestimate the exact critical value. Hence, the nominal significance level may not be achieved for small sample sizes (possibly $n \leq 10$ or $n \leq 50$).

Value

A list with the following elements:

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

Examples

#----------------------------------------------------------------------------
# lrt_nb() examples
#----------------------------------------------------------------------------
library(depower)

set.seed(1234)
sim_nb(
  n1 = 60,
  n2 = 40,
  mean1 = 10,
  ratio = 1.5,
  dispersion1 = 2,
  dispersion2 = 8
) |>
  lrt_nb()

Description

Maximum likelihood estimates (MLE) for bivariate negative binomial outcomes.

Usage

mle_bnb_null(data, ratio_null = 1, method = "nlm_constrained", ...)

mle_bnb_alt(data, method = "nlm_constrained", ...)

Arguments

Details

These functions are primarily designed for speed in simulation. Missing values are silently excluded.

Suppose $X_1 \mid G = g \sim \text{Poisson}(\mu g)$ and $X_2 \mid G = g \sim \text{Poisson}(r \mu g)$ where $G \sim \text{Gamma}(\theta, \theta^{-1})$ is the random item (subject) effect. Then $X_1, X_2 \sim \text{BNB}(\mu, r, \theta)$ is the joint distribution where $X_1$ and $X_2$ are dependent (though conditionally independent), $X_1$ is the count outcome for sample 1 of the items (subjects), $X_2$ is the count outcome for sample 2 of the items (subjects), $\mu$ is the conditional mean of sample 1, $r$ is the ratio of the conditional means of sample 2 with respect to sample 1, and $\theta$ is the gamma distribution shape parameter which controls the dispersion and the correlation between sample 1 and 2.

The MLEs of $r$ and $\mu$ are $\hat{r} = \frac{\bar{x}_2}{\bar{x}_1}$ and $\hat{\mu} = \bar{x}_1$. The MLE of $\theta$ is found by maximizing the profile log-likelihood $l(\hat{r}, \hat{\mu}, \theta)$ with respect to $\theta$. When $r = r_{null}$ is known, the MLE of $\mu$ is $\tilde{\mu} = \frac{\bar{x}_1 + \bar{x}_2}{1 + r_{null}}$ and $\tilde{\theta}$ is obtained by maximizing the profile log-likelihood $l(r_{null}, \tilde{\mu}, \theta)$ with respect to $\theta$.

The backend method for numerical optimization is controlled by argument method which refers to stats::nlm(), stats::nlminb(), or stats::optim(). If you would like to see warnings from the optimizer, include argument warnings = TRUE.

Value

• For mle_bnb_alt, a list with the following elements:

  Slot	Name	Description
  1	mean1	MLE for mean of sample 1.
  2	mean2	MLE for mean of sample 2.
  3	ratio	MLE for ratio of means.
  4	dispersion	MLE for BNB dispersion.
  5	nll	Minimum of negative log-likelihood.
  6	nparams	Number of estimated parameters.
  7	n1	Sample size of sample 1.
  8	n2	Sample size of sample 2.
  9	method	Method used for the results.
  10	mle_method	Method used for optimization.
  11	mle_code	Integer indicating why the optimization process terminated.
  12	mle_message	Additional information from the optimizer.

• For mle_bnb_null, a list with the following elements:

  Slot	Name	Description
  1	mean1	MLE for mean of sample 1.
  2	mean2	MLE for mean of sample 2.
  3	ratio_null	Population ratio of means assumed for null hypothesis. mean2 = mean1 * ratio_null.
  4	dispersion	MLE for BNB dispersion.
  5	nll	Minimum of negative log-likelihood.
  6	nparams	Number of estimated parameters.
  7	n1	Sample size of sample 1.
  8	n2	Sample size of sample 2.
  9	method	Method used for the results.
  10	mle_method	Method used for optimization.
  11	mle_code	Integer indicating why the optimization process terminated.
  12	mle_message	Additional information from the optimizer.

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

See Also

sim_bnb(), nll_bnb

Examples

#----------------------------------------------------------------------------
# mle_bnb() examples
#----------------------------------------------------------------------------
library(depower)

set.seed(1234)
d <- sim_bnb(
  n = 40,
  mean1 = 10,
  ratio = 1.2,
  dispersion = 2
)

mle_alt <- d |>
  mle_bnb_alt()

mle_null <- d |>
  mle_bnb_null()

mle_alt
mle_null

Description

Maximum likelihood estimates (MLE) for two independent negative binomial outcomes.

Usage

mle_nb_null(
  data,
  equal_dispersion = FALSE,
  ratio_null = 1,
  method = "nlm_constrained",
  ...
)

mle_nb_alt(data, equal_dispersion = FALSE, method = "nlm_constrained", ...)

Arguments

Details

These functions are primarily designed for speed in simulation. Missing values are silently excluded.

Suppose $X_1 \sim \text{NB}(\mu, \theta_1)$ and $X_2 \sim \text{NB}(r\mu, \theta_2)$, where $X_1$ and $X_2$ are independent, $X_1$ is the count outcome for items in group 1, $X_2$ is the count outcome for items in group 2, $\mu$ is the arithmetic mean count in group 1, $r$ is the ratio of arithmetic means for 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 MLEs of $r$ and $\mu$ are $\hat{r} = \frac{\bar{x}_2}{\bar{x}_1}$ and $\hat{\mu} = \bar{x}_1$. The MLEs of $\theta_1$ and $\theta_2$ are found by maximizing the profile log-likelihood $l(\hat{r}, \hat{\mu}, \theta_1, \theta_2)$ with respect to $\theta_1$ and $\theta_2$. When $r = r_{null}$ is known, the MLE of $\mu$ is $\tilde{\mu} = \frac{n_1 \bar{x}_1 + n_2 \bar{x}_2}{n_1 + n_2}$ and $\tilde{\theta}_1$ and $\tilde{\theta}_2$ are obtained by maximizing the profile log-likelihood $l(r_{null}, \tilde{\mu}, \theta_1, \theta_2)$.

The backend method for numerical optimization is controlled by argument method which refers to stats::nlm(), stats::nlminb(), or stats::optim(). If you would like to see warnings from the optimizer, include argument warnings = TRUE.

Value

• For mle_nb_alt(), a list with the following elements:

  Slot	Name	Description
  1	mean1	MLE for mean of group 1.
  2	mean2	MLE for mean of group 2.
  3	ratio	MLE for ratio of means.
  4	dispersion1	MLE for dispersion of group 1.
  5	dispersion2	MLE for dispersion of group 2.
  6	equal_dispersion	Were equal dispersions assumed.
  7	n1	Sample size of group 1.
  8	n2	Sample size of group 2.
  9	nll	Minimum of negative log-likelihood.
  10	nparams	Number of estimated parameters.
  11	method	Method used for the results.
  12	mle_method	Method used for optimization.
  13	mle_code	Integer indicating why the optimization process terminated.
  14	mle_message	Additional information from the optimizer.

• For mle_nb_null(), a list with the following elements:

  Slot	Name	Description
  1	mean1	MLE for mean of group 1.
  2	mean2	MLE for mean of group 2.
  3	ratio_null	Population ratio of means assumed for null hypothesis. mean2 = mean1 * ratio_null.
  4	dispersion1	MLE for dispersion of group 1.
  5	dispersion2	MLE for dispersion of group 2.
  6	equal_dispersion	Were equal dispersions assumed.
  7	n1	Sample size of group 1.
  8	n2	Sample size of group 2.
  9	nll	Minimum of negative log-likelihood.
  10	nparams	Number of estimated parameters.
  11	method	Method used for the results.
  12	mle_method	Method used for optimization.
  13	mle_code	Integer indicating why the optimization process terminated.
  14	mle_message	Additional information from the optimizer.

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

See Also

sim_nb(), nll_nb

Examples

#----------------------------------------------------------------------------
# mle_nb() examples
#----------------------------------------------------------------------------
library(depower)

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

mle_alt <- d |>
  mle_nb_alt()

mle_null <- d |>
  mle_nb_null()

mle_alt
mle_null

Description

The negative log-likelihood for bivariate negative binomial outcomes.

Usage

nll_bnb_null(param, value1, value2, ratio_null)

nll_bnb_alt(param, value1, value2)

Arguments

Details

These functions are primarily designed for speed in simulation. Arguments are not checked.

Suppose $X_1 \mid G = g \sim \text{Poisson}(\mu g)$ and $X_2 \mid G = g \sim \text{Poisson}(r \mu g)$ where $G \sim \text{Gamma}(\theta, \theta^{-1})$ is the random item (subject) effect. Then $X_1, X_2 \sim \text{BNB}(\mu, r, \theta)$ is the joint distribution where $X_1$ and $X_2$ are dependent (though conditionally independent), $X_1$ is the count outcome for sample 1 of the items (subjects), $X_2$ is the count outcome for sample 2 of the items (subjects), $\mu$ is the conditional mean of sample 1, $r$ is the ratio of the conditional means of sample 2 with respect to sample 1, and $\theta$ is the gamma distribution shape parameter which controls the dispersion and the correlation between sample 1 and 2.

The likelihood is

$\begin{aligned} L(r, \mu, \theta \mid X_1, X_2) = & \left( \frac{\theta^{\theta}}{\Gamma(\theta)} \right)^{n} \times \\ & \frac{\mu^{\sum{x_{1i}} + \sum{x_{2i}}}}{\prod_{i=1}^{n} x_{1i}!} \frac{r^{\sum{x_{2i}}}}{\prod_{i=1}^{n} x_{2i}!} \times \\ & \frac{\prod_{i = 1}^{n} \Gamma(x_{1i} + x_{2i} + \theta)}{(\mu + r \mu + \theta)^{\sum x_{1i} + x_{2i} + \theta}} \end{aligned}$

and the parameter space is $\Theta = \left\{ (r, \mu, \theta) : r, \mu, \theta > 0 \right\}$. The log-likelihood is

$\begin{aligned} l(r, \mu, \theta) = \ & n \left[ \theta \ln \theta - \ln \Gamma(\theta) \right] + \\ & n (\bar{x}_1 + \bar{x}_2) \ln(\mu) + n \bar{x}_2 \ln r + \\ & \sum_{i=1}^{n}{\ln \Gamma(x_{1i} + x_{2i} + \theta)} - \\ & n (\bar{x}_1 + \bar{x}_2 + \theta) \ln(\mu + r\mu + \theta) - \\ & \sum_{i = 1}^{n}{\ln x_{1i}!} - \sum_{i = 1}^{n}{\ln x_{2i}!} \end{aligned}$

Value

Scalar numeric negative log-likelihood.

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

See Also

sim_nb(), stats::nlminb(), stats::nlm(), stats::optim()

Examples

#----------------------------------------------------------------------------
# nll_bnb*() examples
#----------------------------------------------------------------------------
library(depower)

set.seed(1234)
d <- sim_bnb(
  n = 40,
  mean1 = 10,
  ratio = 1.2,
  dispersion = 2
)

nll_bnb_alt(
  param = c(mean1 = 10, mean2 = 12, dispersion = 2),
  value1 = d[[1L]],
  value2 = d[[2L]]
)

nll_bnb_null(
  param = c(mean = 10, dispersion = 2),
  value1 = d[[1L]],
  value2 = d[[2L]],
  ratio_null = 1
)

Description

The negative log-likelihood for two independent samples of negative binomial distributions.

Usage

nll_nb_null(param, value1, value2, equal_dispersion, ratio_null)

nll_nb_alt(param, value1, value2, equal_dispersion)

Arguments

Details

These functions are primarily designed for speed in simulation. Arguments are not checked.

Suppose $X_1 \sim \text{NB}(\mu, \theta_1)$ and $X_2 \sim \text{NB}(r\mu, \theta_2)$ where $X_1$ and $X_2$ are independent, $X_1$ is the count outcome for items in group 1, $X_2$ is the count outcome for items in group 2, $\mu$ is the arithmetic mean count in group 1, $r$ is the ratio of arithmetic means for 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.

Unequal dispersion parameters

When the dispersion parameters are not equal, the likelihood is

$\begin{aligned} L(r, \mu, \theta_1, \theta_2 \mid X_1, X_2) = & \left( \frac{\theta_1^{\theta_1}}{\Gamma(\theta_1)} \right)^{n_1} \frac{\mu^{\sum{x_{1i}}}}{(\mu + \theta_1)^{\sum{x_{1i} + n_1 \theta_1}}} \times \\ & \left( \frac{\theta_2^{\theta_2}}{\Gamma(\theta_2)} \right)^{n_2} \frac{(r \mu)^{\sum{x_{2j}}}}{(r \mu + \theta_2)^{\sum{x_{2j} + n_2 \theta_2}}} \times \\ & \prod_{i = 1}^{n_1}{\frac{\Gamma(x_{1i} + \theta_1)}{x_{1i}!}} \prod_{j = 1}^{n_2}{\frac{\Gamma(x_{2j} + \theta_2)}{x_{2j}!}} \end{aligned}$

and the parameter space is $\Theta = \left\{ (r, \mu, \theta_1, \theta_2) : r, \mu, \theta_1, \theta_2 > 0 \right\}$. The log-likelihood is

$\begin{aligned} l(r, \mu, \theta_1, \theta_2) = \ &n_1 \left[ \theta_1 \ln \theta_1 - \ln \Gamma(\theta_1) \right] + \\ &n_2 \left[ \theta_2 \ln \theta_2 - \ln \Gamma(\theta_2) \right] + \\ &(n_1 \bar{x}_1 + n_2 \bar{x}_2) \ln(\mu) - n_1 (\bar{x}_1 + \theta_1) \ln(\mu + \theta_1) + \\ &n_2 \bar{x}_2 \ln(r) - n_2 (\bar{x}_2 + \theta_2) \ln(r \mu + \theta_2) + \\ &\sum_{i = 1}^{n_1}{\left( \ln \Gamma(y_{1i} + \theta_1) - \ln(y_{1i}!) \right)} + \\ &\sum_{j = 1}^{n_2}{\left( \ln \Gamma(y_{2j} + \theta_2) - \ln(y_{2j}!) \right)} \end{aligned}$

Equal dispersion parameters

When the dispersion parameters are equal, the likelihood is

$\begin{aligned} L(r, \mu, \theta \mid X_1, X_2) = & \left( \frac{\theta^{\theta}}{\Gamma(\theta)} \right)^{n_1 + n_2} \times \\ & \frac{\mu^{\sum{x_{1i}}}}{(\mu + \theta)^{\sum{x_{1i} + n_1 \theta}}} \frac{(r \mu)^{\sum{x_{2j}}}}{(r \mu + \theta)^{\sum{x_{2j} + n_2 \theta}}} \times \\ & \prod_{i = 1}^{n_1}{\frac{\Gamma(x_{1i} + \theta)}{x_{1i}!}} \prod_{j = 1}^{n_2}{\frac{\Gamma(x_{2j} + \theta)}{x_{2j}!}} \end{aligned}$

and the parameter space is $\Theta = \left\{ (r, \mu, \theta) : r, \mu, \theta > 0 \right\}$. The log-likelihood is

$\begin{aligned} l(r, \mu, \theta) = \ &(n_1 + n_2) \left[ \theta \ln \theta - \ln \Gamma(\theta) \right] + \\ &(n_1 \bar{x}_1 + n_2 \bar{x}_2) \ln(\mu) - n_1 (\bar{x}_1 + \theta) \ln(\mu + \theta) + \\ &n_2 \bar{x}_2 \ln(r) - n_2 (\bar{x}_2 + \theta) \ln(r \mu + \theta) + \\ &\sum_{i = 1}^{n_1}{\left( \ln \Gamma(y_{1i} + \theta) - \ln(y_{1i}!) \right)} + \\ &\sum_{j = 1}^{n_2}{\left( \ln \Gamma(y_{2j} + \theta) - \ln(y_{2j}!) \right)} \end{aligned}$

Value

Scalar numeric negative log-likelihood.

References

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

See Also

sim_nb(), stats::nlminb(), stats::nlm(), stats::optim()

Examples

#----------------------------------------------------------------------------
# nll_nb_*() examples
#----------------------------------------------------------------------------
library(depower)

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

nll_nb_alt(
  param = c(mean1 = 10, mean2 = 15, dispersion = 2, dispersion2 = 8),
  value1 = d[[1L]],
  value2 = d[[2L]],
  equal_dispersion = FALSE
)

nll_nb_null(
  param = c(mean = 10, dispersion = 2, dispersion2 = 8),
  value1 = d[[1L]],
  value2 = d[[2L]],
  equal_dispersion = FALSE,
  ratio_null = 1
)

Description

An automatic plot method for objects returned by power().

Usage

## S3 method for class 'depower'
plot(
  x,
  x_axis = NULL,
  y_axis = NULL,
  color = NULL,
  facet_row = NULL,
  facet_col = NULL,
  hline = NULL,
  caption = TRUE,
  caption_width = 70L,
  ...
)

Arguments

Details

If you are limited by the output from plot.depower(), keep in mind that the object returned by power() is a standard data frame. This allows you to easily plot all results with standard plotting functions. In addition, because plot.depower() uses ggplot2, you can modify the plot as you normally would. For example:

set.seed(1234)
sim_log_lognormal(
  n1 = c(10, 15),
  n2 = c(10, 15),
  ratio = c(1.3, 1.5),
  cv1 = c(0.3),
  cv2 = c(0.3, 0.5),
  nsims = 1000
) |>
  power(alpha = 0.05) |>
  plot(hline = 0.8, caption_width = 60) +
  ggplot2::theme_bw() +
  ggplot2::theme(plot.caption = ggplot2::element_text(hjust = 0)) +
  ggplot2::labs(title = "Power for the ratio of geometric means")

Value

A ggplot2::ggplot() object.

See Also

power()

Examples

#----------------------------------------------------------------------------
# plot() examples
#----------------------------------------------------------------------------
library(depower)

# Power for independent two-sample t-test
set.seed(1234)
sim_log_lognormal(
  n1 = c(10, 15),
  n2 = c(10, 15),
  ratio = c(1.3, 1.5),
  cv1 = c(0.3),
  cv2 = c(0.3, 0.5),
  nsims = 500
) |>
  power(alpha = 0.05) |>
  plot()

# Power for dependent two-sample t-test
set.seed(1234)
sim_log_lognormal(
  n1 = c(10, 15),
  n2 = c(10, 15),
  ratio = c(1.3, 1.5),
  cv1 = c(0.3, 0.5),
  cv2 = c(0.3, 0.5),
  cor = c(0.3),
  nsims = 500
) |>
  power(alpha = 0.01) |>
  plot()

# Power for two-sample independent AND two-sample dependent t-test
set.seed(1234)
sim_log_lognormal(
  n1 = c(10, 15),
  n2 = c(10, 15),
  ratio = c(1.3, 1.5),
  cv1 = c(0.3),
  cv2 = c(0.3),
  cor = c(0, 0.3, 0.6),
  nsims = 500
) |>
  power(alpha = c(0.05, 0.01)) |>
  plot(facet_row = "cor", color = "test")

# Power for one-sample t-test
set.seed(1234)
sim_log_lognormal(
  n1 = c(10, 15),
  ratio = c(1.2, 1.4),
  cv1 = c(0.3, 0.5),
  nsims = 500
) |>
  power(alpha = c(0.05, 0.01)) |>
  plot()


# Power for independent two-sample NB test
set.seed(1234)
sim_nb(
  n1 = c(10, 15),
  mean1 = 10,
  ratio = c(1.8, 2),
  dispersion1 = 10,
  dispersion2 = 3,
  nsims = 100
) |>
  power(alpha = 0.01) |>
  plot()

# Power for BNB test
set.seed(1234)
sim_bnb(
  n = c(10, 12),
  mean1 = 10,
  ratio = c(1.3, 1.5),
  dispersion = 5,
  nsims = 100
) |>
  power(alpha = 0.01) |>
  plot()

Description

A method to calculate power for objects returned by sim_log_lognormal(), sim_nb(), and sim_bnb().

Usage

power(data, ..., alpha = 0.05, list_column = FALSE, ncores = 1L)

Arguments

Details

Power is calculated as the proportion of hypothesis tests which result in a p-value less than or equal to alpha. e.g.

sum(p <= alpha, na.rm = TRUE) / nsims

Power is defined as the expected probability of rejecting the null hypothesis for a chosen value of the unknown effect. In a multiple comparisons scenario, power is defined as the marginal power, which is the expected power of the test for each individual null hypothesis assumed to be false.

Other forms of power under the multiple comparisons scenario include disjunctive or conjunctive power. Disjunctive power is defined as the expected probability of correctly rejecting one or more null hypotheses. Conjunctive power is defined as the expected probability of correctly rejecting all null hypotheses. In the simplest case, and where all hypotheses are independent, if the marginal power is defined as $\pi$ and $m$ is the number of null hypotheses assumed to be false, then disjunctive power may be calculated as $1 - (1 - \pi)^m$ and conjunctive power may be calculated as $\pi^m$. Disjunctive power tends to decrease with increasingly correlated hypotheses and conjunctive power tends to increase with increasingly correlated hypotheses.

Argument ...

... are the name-value pairs for the functions used to perform the tests. If not named, the functions coerced to character will be used for the name-value pairs. Typical in non-standard evaluation, ... accepts bare functions and converts them to a list of expressions. Each element in this list will be validated as a call and then evaluated on the simulated data. A base::call() is simply an unevaluated function. Below are some examples of specifying ... in power().

# Examples of specifying ... in power()
data <- sim_nb(
  n1 = 10,
  mean1 = 10,
  ratio = c(1.6, 2),
  dispersion1 = 2,
  dispersion2 = 2,
  nsims = 200
)

# ... is empty, so a an appropriate default function will be provided
power(data)

# This is equivalent to leaving ... empty
power(data, "NB Wald test" = wald_test_nb())

# If not named, "wald_test_nb()" will be used to label the function
power(data, wald_test_nb())

# You can specify any parameters in the call. The data argument
# will automatically be inserted or overwritten.
data |>
  power("NB Wald test" = wald_test_nb(equal_dispersion=TRUE, link="log"))

# Multiple functions may be used.
data |>
  power(
    wald_test_nb(link='log'),
    wald_test_nb(link='sqrt'),
    wald_test_nb(link='squared'),
    wald_test_nb(link='identity')
  )

# Just like functions in a pipe, the parentheses are required.
# This will error because wald_test_nb is missing parentheses.
try(power(data, wald_test_nb))

In most cases*, any user created test function may be utilized in ... if the following conditions are satisfied:

1. The function contains argument data which is defined as a list with the first and second elements for simulated data.

2. The return object is a list with element p for the p-value of the hypothesis test.

Validate with test cases beforehand.

*Simulated data of class log_lognormal_mixed_two_sample has both independent and dependent data. To ensure the appropriate test function is used, power.log_lognormal_mixed_two_sample() allows only t_test_welch() and t_test_paired() in .... Each will be evaluated on the simulated data according to column data$cor. If one or both of these functions are not included in ..., the corresponding default function will be used automatically. If any other test function is included, an error will be returned.

Argument alpha

$\alpha$ is known as the type I assertion probability and is defined as the expected probability of rejecting a null hypothesis when it was actually true. $\alpha$ is compared with the p-value and used as the decision boundary for rejecting or not rejecting the null hypothesis.

The family-wise error rate is the expected probability of making one or more type I assertions among a family of hypotheses. Using Bonferroni's method, $\alpha$ is chosen for the family of hypotheses then divided by the number of tests performed ($m$). Each individual hypothesis is tested at $\frac{\alpha}{m}$. For example, if you plan to conduct 30 hypothesis tests and want to control the family-wise error rate to no greater than $\alpha = 0.05$, you would set alpha = 0.05/30.

The per-family error rate is the expected number of type I assertions among a family of hypotheses. If you calculate power for the scenario where you perform 1,000 hypotheses and want to control the per-family error rate to no greater than 10 type I assertions, you would choose alpha = 10/1000. This implicitly assumes all 1,000 hypotheses are truly null. Alternatively, if you assume 800 of these hypotheses are truly null and 200 are not, alpha = 10/1000 would control the per-family error rate to no greater than 8 type I assertions. If it is acceptable to keep the per-family error rate as 10, setting alpha = 10/800 would provide greater marginal power than the previous scenario.

These two methods assume that the distribution of p-values for the truly null hypotheses are uniform(0,1), but remain valid under various other testing scenarios (such as dependent tests). Other multiple comparison methods, such as FDR control, are common in practice but don't directly fit into this power simulation framework.

Column nsims

The final number of valid simulations per unique set of simulation parameters may be less than the original number requested. This may occur when the test results in a missing p-value. For wald_test_bnb(), pathological MLE estimates, generally from small sample sizes and very small dispersions, may result in a negative estimated standard deviation of the ratio. Thus the test statistic and p-value would not be calculated. Note that simulated data from sim_nb() and sim_bnb() may also reduce nsims during the data simulation phase.

The nsims column in the return data frame is the effective number of simulations for power results.

Value

A data frame with the following columns appended to the data object:

For power(list_column = FALSE), columns data, and result are excluded from the data frame.

References

Yu L, Fernandez S, Brock G (2017). “Power analysis for RNA-Seq differential expression studies.” BMC Bioinformatics, 18(1), 234. ISSN 1471-2105, doi:10.1186/s12859-017-1648-2.

Yu L, Fernandez S, Brock G (2020). “Power analysis for RNA-Seq differential expression studies using generalized linear mixed effects models.” BMC Bioinformatics, 21(1), 198. ISSN 1471-2105, doi:10.1186/s12859-020-3541-7.

Rettiganti M, Nagaraja HN (2012). “Power Analyses for Negative Binomial Models with Application to Multiple Sclerosis Clinical Trials.” Journal of Biopharmaceutical Statistics, 22(2), 237–259. ISSN 1054-3406, 1520-5711, doi:10.1080/10543406.2010.528105.

Aban IB, Cutter GR, Mavinga N (2009). “Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data.” Computational Statistics & Data Analysis, 53(3), 820–833. ISSN 01679473, doi:10.1016/j.csda.2008.07.034.

Julious SA (2004). “Sample sizes for clinical trials with Normal data.” Statistics in Medicine, 23(12), 1921–1986. doi:10.1002/sim.1783.

Vickerstaff V, Omar RZ, Ambler G (2019). “Methods to adjust for multiple comparisons in the analysis and sample size calculation of randomised controlled trials with multiple primary outcomes.” BMC Medical Research Methodology, 19(1), 129. ISSN 1471-2288, doi:10.1186/s12874-019-0754-4.

See Also

plot.depower()

Examples

#----------------------------------------------------------------------------
# power() examples
#----------------------------------------------------------------------------
library(depower)

# Power for independent two-sample t-Test
set.seed(1234)
data <- sim_log_lognormal(
  n1 = 20,
  n2 = 20,
  ratio = c(1.2, 1.4),
  cv1 = 0.4,
  cv2 = 0.4,
  cor = 0,
  nsims = 1000
)

# Welch's t-test is used by default
power(data)

# But you can specify anything else that is needed
power(
  data = data,
  "Welch's t-Test" = t_test_welch(alternative = "greater"),
  alpha = 0.01
)

# Power for dependent two-sample t-Test
set.seed(1234)
sim_log_lognormal(
  n1 = 20,
  n2 = 20,
  ratio = c(1.2, 1.4),
  cv1 = 0.4,
  cv2 = 0.4,
  cor = 0.5,
  nsims = 1000
) |>
  power()

# Power for mixed-type two-sample t-Test
set.seed(1234)
sim_log_lognormal(
  n1 = 20,
  n2 = 20,
  ratio = c(1.2, 1.4),
  cv1 = 0.4,
  cv2 = 0.4,
  cor = c(0, 0.5),
  nsims = 1000
) |>
  power()

# Power for one-sample t-Test
set.seed(1234)
sim_log_lognormal(
  n1 = 20,
  ratio = c(1.2, 1.4),
  cv1 = 0.4,
  nsims = 1000
) |>
  power()

# Power for independent two-sample NB test
set.seed(1234)
sim_nb(
  n1 = 10,
  mean1 = 10,
  ratio = c(1.6, 2),
  dispersion1 = 2,
  dispersion2 = 2,
  nsims = 200
) |>
  power()

# Power for BNB test
set.seed(1234)
sim_bnb(
  n = 10,
  mean1 = 10,
  ratio = c(1.4, 1.6),
  dispersion = 10,
  nsims = 200
) |>
  power()

Description

Simulate data from the bivariate negative binomial (BNB) distribution. The BNB distribution is used to simulate count data where the event counts are jointly dependent (correlated). For independent data, see sim_nb().

Usage

sim_bnb(
  n,
  mean1,
  mean2,
  ratio,
  dispersion,
  nsims = 1L,
  return_type = "list",
  max_zeros = 0.99,
  ncores = 1L
)

Arguments

Details

The negative binomial distribution may be defined using a gamma-Poisson mixture distribution. In this case, the Poisson parameter $\lambda$ is a random variable with gamma distribution. Equivalence between different parameterizations are demonstrated below:

# Define constants and their relationships
n <- 10000
dispersion <- 8
mu <- 4
p <- dispersion / (dispersion + mu)
q <- mu / (mu + dispersion)
variance <- mu + (mu^2 / dispersion)
rate <- p / (1 - p)
scale <- (1 - p) / p

# alternative formula for mu
mu_alt <- (dispersion * (1 - p)) / p
stopifnot(isTRUE(all.equal(mu, mu_alt)))

set.seed(20240321)

# Using built-in rnbinom with dispersion and mean
w <- rnbinom(n = n, size = dispersion, mu = mu)

# Using gamma-Poisson mixture with gamma rate parameter
x <- rpois(
  n = n,
  lambda = rgamma(n = n, shape = dispersion, rate = rate)
)

# Using gamma-Poisson mixture with gamma scale parameter
y <- rpois(
  n = n,
  lambda = rgamma(n = n, shape = dispersion, scale = scale)
)

# Using gamma-Poisson mixture with multiplicative mean and
# gamma scale parameter
z <- rpois(
  n = n,
  lambda = mu * rgamma(n = n, shape = dispersion, scale = 1/dispersion)
)

# Compare CDFs
par(mar=c(4,4,1,1))
plot(
  x = sort(w),
  y = (1:n)/n,
  xlim = range(c(w,x,y,z)),
  ylim = c(0,1),
  col = 'green',
  lwd = 4,
  type = 'l',
  main = 'CDF'
)
lines(x = sort(x), y = (1:n)/n, col = 'red', lwd = 2)
lines(x = sort(y), y = (1:n)/n, col = 'yellow', lwd = 1.5)
lines(x = sort(z), y = (1:n)/n, col = 'black')

The BNB distribution is implemented by compounding two conditionally independent Poisson random variables $X_1 \mid G = g \sim \text{Poisson}(\mu g)$ and $X_2 \mid G = g \sim \text{Poisson}(r \mu g)$ with a gamma random variable $G \sim \text{Gamma}(\theta, \theta^{-1})$. The probability mass function for the joint distribution of $X_1, X_2$ is

$P(X_1 = x_1, X_2 = x_2) = \frac{\Gamma(x_1 + x_2 + \theta)}{(\mu + r \mu + \theta)^{x_1 + x_2 + \theta}} \frac{\mu^{x_1}}{x_1!} \frac{(r \mu)^{x_2}}{x_2!} \frac{\theta^{\theta}}{\Gamma(\theta)}$

where $x_1,x_2 \in \mathbb{N}^{\geq 0}$ are specific values of the count outcomes, $\theta \in \mathbb{R}^{> 0}$ is the dispersion parameter which controls the dispersion and level of correlation between the two samples (otherwise known as the shape parameter of the gamma distribution), $\mu \in \mathbb{R}^{> 0}$ is the mean parameter, and $r = \frac{\mu_2}{\mu_1} \in \mathbb{R}^{> 0}$ is the ratio parameter representing the multiplicative change in the mean of the second sample relative to the first sample. $G$ denotes the random subject effect and the gamma distribution scale parameter is assumed to be the inverse of the dispersion parameter ($\theta^{-1}$) for identifiability.