comparison of methods for computing the mean of the logit-normal

logit_normal.R

## Comparison of bias-correction methods

## n.b. both "logit-normal" and "logistic-normal" are used
##  for logit(Y) ~ Gaussian(mu, sigma)
library(logitnorm)
library(emmeans)

##
eta <- -2
ranef_sd <- 1

res <- numeric(0) ## grow the vector, to hell with efficiency

## brute force 'Monte Carlo'
set.seed(101)
res[["mc1"]] <- mean(plogis(rnorm(10000, eta, ranef_sd))) ## 0.1556

## numerical integration
res[["numint"]] <- momentsLogitnorm(mu=eta,sigma=ranef_sd)[["mean"]]

## delta method à la emmeans
## (finite difference of link$mu.eta())
## https://github.com/rvlenth/emmeans/blob/b4145ad81e2d5bb3fedc82b1e5db68e1305aeb53/R/summary.R#L914-L925
logit_link <- make.link("logit")
dfun <- function(eta, sd) {
    link_ba <- emmeans:::.make.bias.adj.link(logit_link, sd)
    link_ba$linkinv(eta)
}
res[["delta_emmeans"]] <- dfun(eta, ranef_sd) ## 0.1591842

## home-grown delta method
d2fun <- deriv(D(expression(1/(1+exp(-x))), "x"), "x", function.arg = "x")
res[["delta_mine"]] <- plogis(eta) + ranef_sd^2/2*attr(d2fun(eta), "gradient")

res
##          mc1        numint delta_emmeans    delta_mine 
##    0.1556848     0.1554625     0.1591842     0.1591842

## compare delta method to num integration

library(tidyverse); theme_set(theme_bw())
library(cowplot)
library(colorspace)

res <- (crossing(eta = seq(-4, 4, length = 41),
          sd = 10^(seq(-1, 0.5, length = 41)))
    |> pmap_dfr(~ tibble(eta = .x, sd = .y,
                         delta = dfun(.x, .y),
                         numint = momentsLogitnorm(mu=.x,
                                                   sigma=.y)[["mean"]]))
    |> mutate(delta_dev = (delta-numint)/numint)
)


brkvec1 <- c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 10)
brkvec <- c(-rev(brkvec1), 0, brkvec1)

gg0 <- ggplot(res, aes(eta, sd)) +
    scale_x_continuous(expand = c(0,0)) +
    scale_y_log10(expand = c(0,0))

pg1 <- plot_grid(nrow = 1,
          rel_widths = c(0.3, 0.4, 0.4),
          gg0 + geom_contour_filled(aes(z=numint)) + labs(title = "numint") +
         theme(legend.position = "none"),
         gg0 + geom_contour_filled(aes(z=delta)) + labs(title = "delta"),
         gg0 + geom_contour_filled(aes(z=delta_dev), breaks = brkvec) +
         labs(title = "relative diff") +
         scale_fill_discrete_diverging()
)

print(pg1)
save_plot(pg1, file = "delta.png", base_width = 12)