i don't have any references for this technique tho

gistfile1.md

library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#> 
#> Attaching package: 'gamlss.data'
#> The following object is masked from 'package:datasets':
#> 
#>     sleep
#> Loading required package: gamlss.dist
#> Loading required package: MASS
#> Loading required package: nlme
#> Loading required package: parallel
#>  **********   GAMLSS Version 5.4-12  **********
#> For more on GAMLSS look at https://www.gamlss.com/
#> Type gamlssNews() to see new features/changes/bug fixes.
library(marginaleffects)
library(tidyverse)

#Helper function
inv_logit <- function(p) exp(p) / (1 + exp(p))

set.seed(123)

num <- 10^2

N = 28

d <- tibble(dafi = c(rZIBB(num, mu =.8, sigma = .1, nu = .13, bd = N), 
                     rZIBB(num, mu =.5, sigma = .9, nu = .20, bd = N)), 
            arm = rep(c("control", "intervention"), each = num))

# No intercepts so group diffs are easy to compute
fit <- gamlss(
  cbind(dafi,N-dafi) ~ -1 + arm,
  sigma.formula = ~ -1 + arm, 
  nu.formula = ~ -1 + arm, 
  data = d,
  family = ZIBB(
    mu.link = "logit", 
    sigma.link = "log",
    nu.link = "logit"
  ), 
  control = gamlss.control(trace = FALSE)
)

mu <- coefAll(fit) |> unlist()
sigma <- vcov(fit)

# where the SEs come from
summary(fit)[, 2]
#>      armcontrol armintervention      armcontrol armintervention      armcontrol 
#>      0.09782647      0.24909966      0.20658828      0.24045186      0.30773091 
#> armintervention 
#>      0.87344668
sqrt(diag(sigma))
#>      armcontrol armintervention      armcontrol armintervention      armcontrol 
#>      0.09782647      0.24909966      0.20658828      0.24045186      0.30773091 
#> armintervention 
#>      0.87344668

draws <- MASS::mvrnorm(100000, mu, sigma) |> as.data.frame()
epred_control <- plogis(draws$mu.armcontrol) * N * (1 - plogis(draws$nu.armcontrol))
epred_intervention <- plogis(draws$mu.armintervention) * N * (1 - plogis(draws$nu.armintervention))

summary(epred_control)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   14.09   18.96   19.55   19.49   20.08   22.28
summary(epred_intervention)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   2.735  10.095  10.864  10.818  11.604  15.885
summary(epred_control - epred_intervention)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   1.496   7.716   8.654   8.670   9.603  17.315
quantile(epred_control - epred_intervention, c(.025, .975))
#>     2.5%    97.5% 
#>  5.87795 11.54691


expression <- "exp(mu.armcontrol) / (1 + exp(mu.armcontrol)) * N * (1 - exp(nu.armcontrol) / (1 + exp(nu.armcontrol))) - exp(mu.armintervention) / (1 + exp(mu.armintervention)) * N * (1 - exp(nu.armintervention) / (1 + exp(nu.armintervention)))"
car::deltaMethod(
  mu, expression,
  vcov(fit),
  func = "diff"
)
#>      Estimate     SE  2.5 % 97.5 %
#> diff   8.2950 1.3004 5.7463 10.844

^{Created on 2023-05-31 with reprex v2.0.2}