R/maturity.R

In pbs-assess/gfplot: Data Extraction and Plotting for DFO PBS Groundfish Stocks

#' Fit and plot maturity ogives
#'
#' @param dat Data from [gfdata::get_survey_samples()].
#' @param type Should this be an age or length fit?
#' @param sample_id_re If `TRUE` then the model will include random intercepts
#'   for sample ID.
#' @param year_re If `TRUE` the model will include random intercepts
#'   for year.
#' @param months A numeric vector indicating which months to include when
#'   fitting the maturity ogive. Defaults to all months.
#' @param custom_maturity_at A numeric vector of two threshold codes to define
#'   maturity at with the first being for males and the second for females.
#'   Defaults to `NULL`, which brings in default values from maturity assignment
#'   dataframe included with this package.
#'   `NA` in either position will also retain the default.
#' @param ageing_method_codes A numeric vector of ageing method codes to filter
#'   on. Defaults to `NULL`, which brings in all valid ageing codes.
#'   See [gfdata::get_age_methods()].
#' @param usability_codes An optional vector of usability codes.
#'   All usability codes not in this vector will be omitted.
#'   Set to `NULL` to include all samples.
#' @param link The link function for the binomial GLM (by default, \code{"logit"}).
#' @rdname plot_mat_ogive
#' @export
#' @examples
#' # d <- gfdata::get_survey_samples("pacific ocean perch", ssid = 1)
#' d <- pop_samples
#'
#' m <- fit_mat_ogive(d, type = "age", sample_id_re = FALSE)
#' plot_mat_ogive(m)
#'
#' m <- fit_mat_ogive(d, type = "length", sample_id_re = FALSE)
#' plot_mat_ogive(m)
#' \dontrun{
#' ## with random intercepts for sample ID:
#' m <- fit_mat_ogive(d, type = "length", sample_id_re = TRUE)
#' plot_mat_ogive(m)
#' }
fit_mat_ogive <- function(dat,
                          type = c("age", "length"),
                          sample_id_re = FALSE,
                          year_re = FALSE,
                          months = seq(1, 12),
                          custom_maturity_at = NULL,
                          ageing_method_codes = NULL,
                          usability_codes = c(0, 1, 2, 6),
                          link = c("logit", "probit", "cloglog", "cauchit", "log")) {
  link <- match.arg(link)

  dat <- mutate(dat, month = lubridate::month(trip_start_date))

  dat <- dat %>% filter(maturity_convention_code != 9)

  if ("maturity_convention_maxvalue" %in% names(dat)) {
    dat <- dat %>%
      filter(maturity_code <= maturity_convention_maxvalue)
  }
  if (!is.null(usability_codes)) {
    dat <- filter(dat, .data$usability_code %in% usability_codes)
  }

  type <- match.arg(type)
  dat <- dat[dat$sex %in% c(1, 2), , drop = FALSE]
  dat <- dat[dat$month %in% months, , drop = FALSE]

  if (type == "age" && !is.null(ageing_method_codes)) {
    dat <- filter(dat, ageing_method %in% ageing_method_codes)
  }

  if (nrow(dat) == 0) {
    warning("No data")
    return(NA)
  }

  dat <- dat[!duplicated(dat$specimen_id), , drop = FALSE] # critical!
  dat <- dat %>%
    select(
      species_common_name,
      year,
      age, length, weight,
      maturity_code, sex,
      maturity_convention_code,
      # survey_series_id,
      specimen_id, sample_id, trip_start_date
    )

  mat_df <- maturity_assignment

  dat <- left_join(dat, mat_df, by = c("sex", "maturity_convention_code"))

  if(!is.null(custom_maturity_at)){
    if(!is.na(custom_maturity_at[1])) dat[dat$sex == 1,]$mature_at <- custom_maturity_at[1]
    if(!is.na(custom_maturity_at[2])) dat[dat$sex == 2,]$mature_at <- custom_maturity_at[2]
  }

  dat <- mutate(dat, mature = maturity_code >= mature_at)

  type <- match.arg(type)
  .d <- switch(type,
    age = filter(dat, !is.na(mature), !is.na(age), !is.na(sex)) %>%
      rename(age_or_length = age),
    length = filter(dat, !is.na(mature), !is.na(length), !is.na(sex), !is.na(year)) %>%
      rename(age_or_length = length)
  )
  .d <- mutate(.d, female = ifelse(sex == 2L, 1L, 0L))
  if (nrow(.d) == 0) {
    warning("No data")
    return(NA)
  }

  if (sample_id_re) {
    if (year_re) {
      m <- glmmTMB::glmmTMB(mature ~ age_or_length * female +
        (1 | sample_id) +
        # (1 | survey_series_id) +
        (1 | year),
      data = .d, family = binomial(link)
      )
      b <- glmmTMB::fixef(m)[[1L]]
      ranef <- glmmTMB::ranef(m)
      re <- ranef$cond$year
    } else {
      m <- glmmTMB::glmmTMB(mature ~ age_or_length * female + (1 | sample_id),
        data = .d, family = binomial(link)
      )

      b <- glmmTMB::fixef(m)[[1L]]
    }
  } else {
    if (year_re) {
      m <- glmmTMB::glmmTMB(mature ~ age_or_length * female + (1 | year),
        data = .d, family = binomial(link)
      )
      b <- glmmTMB::fixef(m)[[1L]]
      ranef <- glmmTMB::ranef(m)
      re <- ranef$cond$year
    } else {
      m <- stats::glm(mature ~ age_or_length * female,
        data = .d, family = binomial(link)
      )
      b <- stats::coef(m)
    }
  }

  age_or_length <- seq(min(.d$age_or_length), max(.d$age_or_length),
    length.out = 300L
  )

  if (year_re) {
    # builds nd using all sample ids to make sure full year distribution is represented
    year <- unique(.d$year)
    s_ids <- .d %>%
      select(year, sample_id) %>%
      distinct()

    nd1 <- expand.grid(
      age_or_length = age_or_length,
      sample_id = s_ids$sample_id,
      female = c(0L, 1L), stringsAsFactors = FALSE
    )

    nd <- left_join(nd1, s_ids)

    year_f <- as.character(nd$year)
  } else {
    if (length(unique(.d$sample_id)) > 100L) {
      s_ids <- sample(unique(.d$sample_id), 100L)
    } else {
      s_ids <- unique(.d$sample_id)
    }

    nd <- expand.grid(
      age_or_length = age_or_length, sample_id = s_ids,
      female = c(0L, 1L), stringsAsFactors = FALSE
    )
  }

  if (sample_id_re) {
    nd$glmm_re <- predict(m, newdata = nd, type = "response", se.fit = FALSE)
  }

  if (year_re) {
    year_f <- as.character(nd$year)
    nd$glmm_re <- predict(m, newdata = nd, type = "response", se.fit = FALSE)
    nd$glmm_re2 <- family(m)$linkinv(b[[1L]] + re[year_f, ] + b[[3L]] * nd$female +
        b[[2L]] * nd$age_or_length + b[[4L]] * nd$age_or_length * nd$female)
  }

  nd$glmm_fe <- family(m)$linkinv(b[[1L]] + b[[3L]] * nd$female +
    b[[2L]] * nd$age_or_length + b[[4L]] * nd$age_or_length * nd$female)

  if (year_re) {
    mat_perc <- extract_maturity_perc_re(b, re, m)
  } else {
    mat_perc <- tryCatch(extract_maturity_perc(b, m),
      error = function(e) {message("error: ", e); return(NA)})
  }

  if (link == "logit") {
    se_l50 <- tryCatch({delta_method(~ -(log((1/0.5) - 1) + x1 + x3) / (x2 + x4),
                                     mean = stats::coef(m), cov = stats::vcov(m))}, error = function(e) NA)
  } else if(requireNamespace("numDeriv", quietly = TRUE)) {
    se_l50 <- mat_par_delta_method(model = m, perc = 0.5) %>% as.numeric() %>% sqrt()
  } else {
    se_l50 <- paste("Download the numDeriv package to get the SE of the", type, "of 50% maturity.")
  }

  list(
    data = .d, pred_data = nd, model = m, sample_id_re = sample_id_re,
    year_re = year_re,
    type = type, mat_perc = mat_perc, se_50 = se_l50
  )
}

delta_method <- function(g, mean, cov) {
  # simplified from msm::deltamethod
  cov <- as.matrix(cov)
  n <- length(mean)
  g <- list(g)
  syms <- paste0("x", seq_len(n))
  for (i in seq_len(n)) assign(syms[i], mean[i])
  gdashmu <- t(sapply(g, function(form) {
    as.numeric(attr(eval(stats::deriv(form, syms)), "gradient"))
  }))
  new.covar <- gdashmu %*% cov %*% t(gdashmu)
  sqrt(diag(new.covar))
}

sd2cv <- function(.sd) {
  sqrt(exp(.sd^2) - 1)
}


#' @param object Output from [fit_mat_ogive()].
#' @param xlab X axis label.
#' @param col A named character vector declaring the colors for F and M.
#' @param title Title for the plot.
#' @param rug Logical indicating whether rug lines should be added.
#' @param rug_n The number of rug lines to sample from the total number of fish.
#' @param x_max Used in determining the right axis limit.
#' @param prediction_type The prediction lines to show. Useful if you only want
#'   to show model fits when you have sufficient data.
#' @param text_label_size Font size for the labels showing the age-at- values
#'   for either age or length on the plot panel
#' @param show_quant_text Logical. If `TRUE`, show the quantile values for each
#' sex on the panel
#' @param french Translate to French?
#'
#' @importFrom stats binomial plogis predict
#' @export
#' @rdname plot_mat_ogive

plot_mat_ogive <- function(object,
                           col = c("M" = "grey50", "F" = "black"),
                           xlab = if (object$type[[1]] == "age") "Age (years)" else "Length (cm)",
                           title =
                             if (object$type[[1]] == "age") "Age at maturity" else "Length at maturity",
                           rug = TRUE, rug_n = 1500, x_max = 1.75,
                           prediction_type = c("all", "male", "female", "none"),
                           text_label_size = 3,
                           show_quant_text = TRUE,
                           french = FALSE) {
  if (object$sample_id_re) {
    b <- glmmTMB::fixef(object$model)[[1L]]
  } else {
    if (object$year_re) {
      b <- glmmTMB::fixef(object$model)[[1L]]
    } else {
      b <- stats::coef(object$model)
    }
  }

  nd_re <- object$pred_data
  nd_fe <- object$pred_data |>
    dplyr::distinct(age_or_length, female, .keep_all = TRUE) # Remove duplicate predictions (predictions were done for each sample_id)
  # nd_fe <- filter(nd_re, year == nd_re$year[[1L]]) # fake; all same
  nd_fe$glmm_re <- NULL # also may not exist if no random effects

  prediction_type <- match.arg(prediction_type)
  if (prediction_type == "male") {
    nd_fe <- filter(nd_fe, female == 0L)
  }
  if (prediction_type == "female") {
    nd_fe <- filter(nd_fe, female == 1L)
  }

  # if (prediction_type == "none") {
  #   nd_fe <- filter(nd_fe, !female %in% c(0L, 1L)) # i.e. nothing
  # }

  m_perc <- data.frame(
    p0.5 = binomial_perc(a = b[[1]], b = b[[2]], perc = 0.5, linkinv = family(object$model)$linkinv)
  )
  m_perc$p0.95 <- binomial_perc(a = b[[1]], b = b[[2]], perc = 0.95, linkinv = family(object$model)$linkinv)
  m_perc$p0.05 <- binomial_perc(a = b[[1]], b = b[[2]], perc = 0.05, linkinv = family(object$model)$linkinv)

  f_perc <- data.frame(
    p0.5 = binomial_perc(a = b[[1]] + b[[3]], b = b[[2]] + b[[4]], perc = 0.5, linkinv = family(object$model)$linkinv)
  )
  f_perc$p0.95 <- binomial_perc(a = b[[1]] + b[[3]], b = b[[2]] + b[[4]], perc = 0.95, linkinv = family(object$model)$linkinv)
  f_perc$p0.05 <- binomial_perc(a = b[[1]] + b[[3]], b = b[[2]] + b[[4]], perc = 0.05, linkinv = family(object$model)$linkinv)

  labs_f <- tibble(
    p = c("05", "50", "95"),
    value = c(f_perc$p0.05, f_perc$p0.5, f_perc$p0.95),
    x = 0.75 * max(nd_re$age_or_length), # re-calculated below
    y = seq(0.75, 0.6, length.out = 3L),
    sex = "F"
  )

  labs_m <- tibble(
    p = c("05", "50", "95"),
    value = c(m_perc$p0.05, m_perc$p0.5, m_perc$p0.95),
    x = 0.75 * max(nd_re$age_or_length), # re-calculated below
    y = seq(0.4, 0.25, length.out = 3),
    sex = "M"
  )
  labs <- bind_rows(labs_m, labs_f)

  if (prediction_type == "male") {
    labs <- filter(labs, sex == "M")
  }
  if (prediction_type == "female") {
    labs <- filter(labs, sex == "F")
  }

  nd_fe <- mutate(nd_fe, sex = ifelse(female == 1L, "F", "M"))
  nd_re <- mutate(nd_re, sex = ifelse(female == 1L, "F", "M"))
  object$data <- mutate(object$data, sex = ifelse(female == 1L, "F", "M"))

  if (object$type[[1]] == "age") {
    labs <- mutate(labs,
      label =
        paste0(
          sex, " ", p, " = ",
          sprintf("%.1f", round(value, 1L)), en2fr("y", translate = french)
        )
    )
  } else {
    labs <- mutate(labs,
      label =
        paste0(sex, " ", p, " = ", sprintf("%.1f", round(value, 1L)), "cm")
    )
  }
  max_x <- min(c(max(labs$value) * x_max, max(nd_fe$age_or_length)))

  if (object$type[[1]] == "age") {
    labs <- mutate(labs, x = max_x * 0.7) # actual x position calculation
  } else {
    labs <- mutate(labs, x = max_x * 0.05) # actual x position calculation
  }

  nd_fe$sex <- factor(nd_fe$sex, levels = c("F", "M"))
  nd_re$sex <- factor(nd_re$sex, levels = c("F", "M"))
  labs$sex <- factor(labs$sex, levels = c("F", "M"))
  labs$sex <- factor(labs$sex, levels = c("F", "M"))

  if (isTRUE(french)) {
    labs$label <- gsub("\\.", ",", labs$label)
  }
  g <- ggplot(nd_fe, aes_string("age_or_length", "glmm_fe", colour = "sex", lty = "sex"))

  if ("glmm_re" %in% names(nd_re)) {
    if (object$sample_id_re) {
      n_re <- length(unique(nd_re$sample_id))/5
      n_re2 <- ifelse(n_re < 15, 15, n_re)
      g <- g + geom_line(
        data = nd_re,
        aes_string("age_or_length", "glmm_re",
          group = "paste(sample_id, sex)",
          colour = "sex", lty = "sex"
        ), inherit.aes = FALSE, alpha = 1/n_re2,
        show.legend = FALSE
      )
    } else {
      if (object$year_re) {
        n_re <- length(unique(nd_re$year))
        g <- g + geom_line(
          data = nd_re,
          aes_string("age_or_length", "glmm_re",
            group = "paste(year, sex)",
            colour = "sex", lty = "sex"
          ), inherit.aes = FALSE, alpha = 2/n_re,
          show.legend = FALSE
        )
      }
    }
  }

  if (prediction_type != "none") {
    g <- g + geom_line(size = 1.0) #lwd = 1, alpha = 0.8
    g <- g + geom_vline(
      data = filter(labs, p == "50"),
      aes_string(xintercept = "value", colour = "sex", lty = "sex"), lwd = 0.8,
      alpha = 0.8,
      show.legend = FALSE
    )
    if(show_quant_text){
      g <- g + geom_text(
        data = labs, aes_string(
          x = "x", y = "y", label = "label"
        ),
        hjust = 0, show.legend = FALSE, size = text_label_size
      )
    }
  }

  g <- g + scale_colour_manual(values = col,
    breaks = c("F", "M"), drop = FALSE) +
    ggplot2::scale_linetype_discrete(breaks = c('F', 'M'), drop = FALSE) +
    xlab(en2fr(xlab, french)) + ylab(en2fr("Probability mature", french)) +
    theme_pbs() +
    coord_cartesian(
      expand = FALSE, ylim = c(-0.005, 1.005),
      xlim = c(0, max_x)
    ) +
    labs(colour = en2fr("Sex", french), lty = en2fr("Sex", french)) +
    ggplot2::ggtitle(en2fr(title, french))

  if (rug) {
    if (nrow(object$data) > rug_n) {
      temp <- object$data[sample(seq_len(nrow(object$data)), rug_n), , drop = FALSE]
    } else {
      temp <- object$data
    }
    position <- if (object$type == "age") "jitter" else "identity"
    g <- g + ggplot2::geom_rug(
      data = filter(temp, mature == 0L),
      sides = "b", position = position, alpha = 0.5, lty = 1,
      aes_string(x = "age_or_length", y = "as.numeric(mature)", colour = "sex", lty = "sex"), show.legend = FALSE
    )
    g <- g + ggplot2::geom_rug(
      data = filter(temp, mature == 1L),
      sides = "t", position = position, alpha = 0.5, lty = 1,
      aes_string(x = "age_or_length", y = "as.numeric(mature)", colour = "sex", lty = "sex"), show.legend = FALSE
    )
  }

  g
}


#' @importFrom stats binomial plogis predict
#' @export
#' @rdname plot_mat_ogive

plot_mat_annual_ogives <- function(object,
                                   xlab = if (object$type[[1]] == "age") "Age (years)" else "Length (cm)",
                                   title =
                                     if (object$type[[1]] == "age") "Age at maturity" else "Length at maturity",
                                   rug = TRUE, rug_n = 1500, x_max = 1.75,
                                   prediction_type = c("all", "male", "female", "none"),
                                   french = FALSE) {
  if (object$year_re) {
    b <- glmmTMB::fixef(object$model)[[1L]]
    ranef <- glmmTMB::ranef(object$model)
    re <- ranef$cond$year
  } else {
    stop("Year is not a random effect in this model; use plot_mat_ogive instead.")
  }

  nd_re <- object$pred_data
  nd_fe <- object$pred_data
  # nd_fe <- filter(nd_re, year == nd_re$year[[1L]]) # fake; all same
  nd_fe$glmm_re <- NULL # also may not exist if no random effects

  prediction_type <- match.arg(prediction_type)
  if (prediction_type == "male") {
    nd_fe <- filter(nd_fe, female == 0L)
  }
  if (prediction_type == "female") {
    nd_fe <- filter(nd_fe, female == 1L)
  }

  labs_year <- list()
  for (i in (unique(as.character(nd_re$year)))) {
    m_perc <- data.frame(
      p0.5 = binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.5, linkinv = family(object$model)$linkinv)
    )
    m_perc$p0.95 <- binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.95, linkinv = family(object$model)$linkinv)
    m_perc$p0.05 <- binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.05, linkinv = family(object$model)$linkinv)

    f_perc <- data.frame(
      p0.5 = binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.5, linkinv = family(object$model)$linkinv)
    )
    f_perc$p0.95 <- binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.95, linkinv = family(object$model)$linkinv)
    f_perc$p0.05 <- binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.05, linkinv = family(object$model)$linkinv)

    labs_f <- tibble(
      p = c("05", "50", "95"),
      value = c(f_perc$p0.05, f_perc$p0.5, f_perc$p0.95),
      x = 0.75 * max(nd_re$age_or_length), # re-calculated below
      y = seq(0.75, 0.6, length.out = 3L),
      sex = "F",
      year = i
    )

    labs_m <- tibble(
      p = c("05", "50", "95"),
      value = c(m_perc$p0.05, m_perc$p0.5, m_perc$p0.95),
      x = 0.75 * max(nd_re$age_or_length), # re-calculated below
      y = seq(0.4, 0.25, length.out = 3),
      sex = "M",
      year = i
    )
    labs_year[[i]] <- bind_rows(labs_m, labs_f)
  }
  labs <- do.call(rbind, labs_year)

  if (prediction_type == "male") {
    labs <- filter(labs, sex == "M")
  }
  if (prediction_type == "female") {
    labs <- filter(labs, sex == "F")
  }

  nd_fe <- mutate(nd_fe, sex = ifelse(female == 1L, "F", "M"))
  nd_re <- mutate(nd_re, sex = ifelse(female == 1L, "F", "M"))
  object$data <- mutate(object$data, sex = ifelse(female == 1L, "F", "M"))

  if (object$type[[1]] == "age") {
    labs <- mutate(labs,
      label =
        paste0(
          sex, " ", p, " = ",
          sprintf("%.1f", round(value, 1L)), en2fr("y", translate = french)
        )
    )
  } else {
    labs <- mutate(labs,
      label =
        paste0(sex, " ", p, " = ", sprintf("%.1f", round(value, 1L)), "cm")
    )
  }
  max_x <- min(c(max(labs$value) * x_max, max(nd_fe$age_or_length)))

  if (object$type[[1]] == "age") {
    labs <- mutate(labs, x = max_x * 0.7) # actual x position calculation
  } else {
    labs <- mutate(labs, x = max_x * 0.05) # actual x position calculation
  }

  nd_fe$sex <- factor(nd_fe$sex, levels = c("F", "M"))
  nd_re$sex <- factor(nd_re$sex, levels = c("F", "M"))
  nd_re$year <- as.factor(nd_re$year)

  g <- ggplot(nd_re, aes_string("age_or_length", "glmm_re2", colour = "year"))
  g <- g + geom_vline(
    data = filter(labs, p == "50"),
    aes_string(xintercept = "value", colour = "year"), lwd = 0.8,
    alpha = 0.3, show.legend = FALSE
  )
  g <- g + geom_line(size = 2, alpha = 0.5)
  g <- g + facet_wrap(~sex, nrow = 2)
  g <- g + scale_colour_viridis_d() +
    labs(colour = "Year") +
    coord_cartesian(
      expand = FALSE, ylim = c(-0.005, 1.005),
      xlim = c(0, max_x)
    ) + gfplot::theme_pbs()

  if (rug) {
    if (nrow(object$data) > rug_n) {
      temp <- object$data[sample(seq_len(nrow(object$data)), rug_n), , drop = FALSE]
    } else {
      temp <- object$data
    }
    position <- if (object$type == "age") "jitter" else "identity"
    g <- g + ggplot2::geom_rug(
      data = filter(temp, mature == 0L),
      sides = "b", position = position, alpha = 0.5, lty = 1, lwd = 2,
      aes_string(x = "age_or_length", y = "as.numeric(mature)", colour = "as.character(year)")
    )
    g <- g + ggplot2::geom_rug(
      data = filter(temp, mature == 1L),
      sides = "t", position = position, alpha = 0.5, lty = 1, lwd = 2,
      aes_string(x = "age_or_length", y = "as.numeric(mature)", colour = "as.character(year)")
    )
  }
  g <- g + xlab(xlab) + ylab("Probability mature") +
    ggplot2::ggtitle(title)
  g
}

extract_maturity_perc <- function(object, model) {

  m.p0.5 <- binomial_perc(a = object[[1]], b = object[[2]], perc = 0.5, linkinv = family(model)$linkinv)
  m.p0.95 <- binomial_perc(a = object[[1]], b = object[[2]], perc = 0.95, linkinv = family(model)$linkinv)
  m.p0.05 <- binomial_perc(a = object[[1]], b = object[[2]], perc = 0.05, linkinv = family(model)$linkinv)

  f.p0.5 <- binomial_perc(
    a = object[[1]] + object[[3]],
    b = object[[2]] + object[[4]], perc = 0.5,
    linkinv =family(model)$linkinv
  )

  f.p0.95 <- binomial_perc(
    a = object[[1]] + object[[3]],
    b = object[[2]] + object[[4]], perc = 0.95,
    linkinv =family(model)$linkinv
  )

  f.p0.05 <- binomial_perc(
    a = object[[1]] + object[[3]],
    b = object[[2]] + object[[4]], perc = 0.05,
    linkinv =family(model)$linkinv
  )
  list(
    m.p0.5 = m.p0.5, m.p0.95 = m.p0.95, m.p0.05 = m.p0.05, f.p0.5 = f.p0.5,
    f.p0.95 = f.p0.95, f.p0.05 = f.p0.05
  )
}

extract_maturity_perc_re <- function(betas, random_intercepts, model) {
  b <- betas
  re <- random_intercepts
  out <- list()
  for (i in rownames(re)) {
    m.p0.5 <- binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.5, linkinv = family(model)$linkinv)
    m.p0.95 <- binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.95, linkinv = family(model)$linkinv)
    m.p0.05 <- binomial_perc(a = b[[1]] + re[i, ], b = b[[2]], perc = 0.05, linkinv = family(model)$linkinv)

    f.p0.5 <- binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.5, linkinv = family(model)$linkinv)
    f.p0.95 <- binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.95, linkinv = family(model)$linkinv)
    f.p0.05 <- binomial_perc(a = b[[1]] + b[[3]] + re[i, ], b = b[[2]] + b[[4]], perc = 0.05, linkinv = family(model)$linkinv)
    out[[i]] <- list(
      m.p0.5 = m.p0.5, m.p0.95 = m.p0.95, m.p0.05 = m.p0.05,
      f.p0.5 = f.p0.5, f.p0.95 = f.p0.95, f.p0.05 = f.p0.05
    )
  }
  m.mean.p0.5 <- binomial_perc(a = b[[1]], b = b[[2]], perc = 0.5, linkinv = family(model)$linkinv)
  f.mean.p0.5 <- binomial_perc(
    a = b[[1]] + b[[3]],
    b = b[[2]] + b[[4]], perc = 0.5,
    linkinv = family(model)$linkinv
  )
  # nrow(re) + 1
  out[["mean"]] <- list(m.mean.p0.5 = m.mean.p0.5, f.mean.p0.5 = f.mean.p0.5)
  out
}

binomial_perc <- function(x, a, b, perc = 0.5, linkinv, ...) {
  f <- function(x) linkinv(a + b * x) - perc
  uniroot(f, interval = c(-100, 200))$root
}

mat_par_delta_method <- function(model, perc = 0.5) {
  stopifnot(requireNamespace("numDeriv", quietly = FALSE))
  f <- function(x) binomial_perc(a = x[1] + x[3], b = x[2] + x[4], perc = perc, linkinv = family(model)$linkinv)
  gradient <- numDeriv::grad(f, x = stats::coef(model))
  gradient %*% stats::vcov(model) %*% gradient
}