tern.gee coverage - 90.20%

Files
Source

#' Methods for GEE Models
#'
#' Additional methods which can simplify working with the GEE result object.
#' @name gee_methods
#' @returns `VarCorr()` returns the estimated covariance matrix, and
#'   `QIC()` returns the QIC value.
NULL

#' @rdname gee_methods
#' @importFrom nlme VarCorr
#'
#' @param x (`tern_gee`)\cr result of [fit_gee()].
#' @inheritParams nlme::VarCorr
#'
#' @exportS3Method
VarCorr.tern_gee <- function(x, sigma = 1, ...) { # nolint
  dim_mat <- length(x$visit_levels)
  tmp <- id_mat <- diag(dim_mat)
  corest <- x$geese$alpha

  # Start with lower-triangular matrix part.
  lower_mat <- switch(x$corstr,
    unstructured = , # Since this is the same as exchangeable, we can do this.
    exchangeable = {
      tmp[lower.tri(tmp)] <- corest
      tmp
    },
    ar1 = {
      row_col_diff <- row(tmp) - col(tmp)
      tmp[lower.tri(tmp)] <- corest^(row_col_diff[lower.tri(row_col_diff)])
      tmp
    },
    `m-dependent` = {
      row_col_diff <- row(tmp) - col(tmp)
      tmp[lower.tri(tmp)] <- corest[row_col_diff[lower.tri(row_col_diff)]]
      tmp
    }
  )

  # Construct the full symmetric matrix.
  mat <- lower_mat + t(lower_mat) - id_mat
  rownames(mat) <- colnames(mat) <- x$visit_levels
  mat
}

#' @rdname gee_methods
#' @importFrom geepack QIC
#'
#' @param object (`tern_gee`)\cr result of [fit_gee()].
#' @inheritParams geepack::QIC
#'
#' @exportS3Method
QIC.tern_gee <- function(object, ...) { # nolint
  object$qic
}

#' Set Variables to Use in GEE Model
#'
#' @param response (`character`)\cr name of response variable.
#' @param covariates (`character`)\cr vector of names of variables to use as covariates.
#' @param id (`character`)\cr name of variable to use to identify unique IDs.
#' @param arm (`character`)\cr name of arm variable.
#' @param visit (`character`)\cr name of visit variable.
#'
#' @return A list of variables that can be used as the `vars` argument in [fit_gee()].
#' @export
#'
#' @examples
#' vars_gee()
#'
#' vars_gee(
#'   response = "CHG",
#'   covariates = c("SEX", "RACE"),
#'   id = "SUBJID",
#'   arm = "ARMCD",
#'   visit = "AVISITN"
#' )
vars_gee <- function(response = "AVAL",
                     covariates = c(),
                     id = "USUBJID",
                     arm = "ARM",
                     visit = "AVISIT") {
  list(
    response = response,
    covariates = covariates,
    id = id,
    arm = arm,
    visit = visit
  )
}

#' @keywords internal
build_formula <- function(vars) {
  assert_list(vars)
  arm_part <- if (is.null(vars$arm)) NULL else vars$arm
  rhs_formula <- paste(
    c(arm_part, vars$covariates),
    collapse = " + "
  )
  stats::as.formula(paste(
    vars$response,
    "~",
    rhs_formula
  ))
}

#' @keywords internal
build_family <- function(regression) {
  assert_string(regression)

  result_object <- switch(regression,
    logistic = stats::binomial(link = "logit"),
    stop(paste("regression type", regression, "not supported"))
  )

  result_class <- paste0("tern_gee_", regression)

  list(
    object = result_object,
    class = result_class,
    control = geeasy::geelm.control(scale.fix = TRUE)
  )
}

#' @keywords internal
build_cor_details <- function(cor_str, vars, data) {
  assert_string(cor_str)
  assert_list(vars)
  assert_data_frame(data)

  result_str <- switch(cor_str,
    "unstructured" = "unstructured",
    "toeplitz" = "m-dependent",
    "compound symmetry" = "exchangeable",
    "auto-regressive" = "ar1",
    stop(paste("correlation structure", cor_str, "not available"))
  )

  result_mv <- switch(cor_str,
    "unstructured" = 1,
    "toeplitz" = nlevels(data[[vars$visit]]) - 1,
    "compound symmetry" = 1,
    "auto-regressive" = 1
  )

  list(
    str = result_str,
    mv = result_mv
  )
}

#' @keywords internal
order_data <- function(data, vars) {
  assert_data_frame(data)
  assert_list(vars)

  if (is.character(data[[vars$visit]])) {
    message(paste("visit variable", vars$visit, "will be coerced to factor for ordering"))
    message("order is:")
    data[[vars$visit]] <- factor(data[[vars$visit]])
    message(paste(toString(levels(data[[vars$visit]])), "\n"))
  }
  if (is.factor(data[[vars$id]]) || is.character(data[[vars$id]])) {
    data[[vars$id]] <- as.integer(as.factor(data[[vars$id]]))
  }
  assert_numeric(data[[vars$id]])

  right_order <- order(data[[vars$id]], data[[vars$visit]])
  data[right_order, ]
}

#' Fit a GEE Model
#'
#' @param vars (`list`)\cr see [vars_gee()].
#' @param data (`data.frame`)\cr input data.
#' @param regression (`string`)\cr choice of regression model.
#' @param cor_struct (`string`)\cr assumed correlation structure.
#'
#' @details The correlation structure can be:
#' * `unstructured`: No constraints are placed on the correlations.
#' * `toeplitz`: Assumes a banded correlation structure, i.e. the correlation
#'    between two time points depends on the distance between the time indices.
#' * `compound symmetry`: Constant correlation between all time points.
#' * `auto-regressive`: Auto-regressive order 1 correlation matrix.
#'
#' @return Object of class `tern_gee` as well as specific to the kind of regression
#'   which was used.
#' @export
#'
#' @examples
#' df <- fev_data
#' df$AVAL <- as.integer(fev_data$FEV1 > 30)
#'
#' fit_gee(vars = vars_gee(arm = "ARMCD"), data = df)
#'
#' fit_gee(vars = vars_gee(arm = "ARMCD"), data = df, cor_struct = "compound symmetry")
fit_gee <- function(vars = vars_gee(),
                    data,
                    regression = c("logistic"),
                    cor_struct = c("unstructured", "toeplitz", "compound symmetry", "auto-regressive")) {
  formula <- build_formula(vars)

  regression <- match.arg(regression)
  family <- build_family(regression)

  data <- order_data(data, vars)
  data[[".id"]] <- data[[vars$id]]
  data[[".waves"]] <- as.integer(data[[vars$visit]])

  cor_struct <- match.arg(cor_struct)
  cor_details <- build_cor_details(cor_struct, vars, data)

  fit <- geeasy::geelm(
    formula = formula,
    id = .id,
    waves = .waves,
    data = data,
    family = family$object,
    corstr = cor_details$str,
    Mv = cor_details$mv,
    control = family$control
  )

  fit$qic <- geepack::QIC(fit)
  fit$visit_levels <- levels(data[[vars$visit]])
  fit$vars <- vars
  fit$data <- data
  assert_factor(data[[vars$arm]])
  fit$ref_level <- levels(data[[vars$arm]])[1L]

  structure(
    fit,
    class = c(family$class, "tern_gee", class(fit))
  )
}

#' Extract Least Square Means from a GEE Model
#'
#' @param object (`tern_gee`)\cr result of [fit_gee()].
#' @param conf_level (`proportion`)\cr confidence level
#' @param weights (`string`)\cr type of weights to be used for the least square means,
#'   see [emmeans::emmeans()] for details.
#' @param specs (`string` or `formula`) specifications passed to [emmeans::emmeans()]
#' @param ... additional arguments for methods
#'
#' @return A `data.frame` with least-square means and contrasts. Additional
#'   classes allow to dispatch downstream methods correctly, too.
#' @export
#'
#' @examples
#' df <- fev_data
#' df$AVAL <- rbinom(n = nrow(df), size = 1, prob = 0.5)
#' fit <- fit_gee(vars = vars_gee(arm = "ARMCD"), data = df)
#'
#' lsmeans(fit)
#'
#' lsmeans(fit, conf_level = 0.90, weights = "equal")
lsmeans <- function(object,
                    conf_level = 0.95,
                    weights = "proportional",
                    specs = object$vars$arm,
                    ...) {
  UseMethod("lsmeans", object)
}

#' @rdname lsmeans
#' @exportS3Method
lsmeans.tern_gee_logistic <- function(object,
                                      conf_level = 0.95,
                                      weights = "proportional",
                                      specs = object$vars$arm,
                                      ...) {
  prop_emm <- emmeans::emmeans(
    object = object,
    specs = specs,
    weights = weights,
    type = "response",
    data = object$data
  )

  prop_df <- cbind(
    data.frame(stats::confint(prop_emm))[, c(object$vars$arm, "prob", "SE", "lower.CL", "upper.CL")],
    n = as.list(prop_emm)$extras[, ".wgt."]
  )
  names(prop_df) <- c(object$vars$arm, "prop_est", "prop_est_se", "prop_lower_cl", "prop_upper_cl", "n")
  ref_level <- levels(object$data[[object$vars$arm]])[1L]

  or_emm <- stats::confint(
    graphics::pairs(prop_emm, reverse = TRUE),
    level = conf_level
  )
  or_df <- as.data.frame(or_emm)
  or_df$comparator <- gsub(pattern = ".+ / (.+)", replacement = "\\1", x = or_df$contrast)
  or_df[[object$vars$arm]] <- gsub(pattern = "(.+) / .+", replacement = "\\1", x = or_df$contrast)
  or_df <- or_df[or_df$comparator == ref_level, ]
  or_df <- rbind(NA, or_df)
  or_df[1L, object$vars$arm] <- ref_level
  or_df <- or_df[, c(object$vars$arm, "odds.ratio", "lower.CL", "upper.CL")]
  or_df <- cbind(or_df, log(or_df[, -1L]), conf_level)
  names(or_df) <- c(
    object$vars$arm,
    "or_est", "or_lower_cl", "or_upper_cl",
    "log_or_est", "log_or_lower_cl", "log_or_upper_cl",
    "conf_level"
  )

  result <- merge(prop_df, or_df, by = object$vars$arm)

  structure(
    result,
    class = c("lsmeans_logistic", class(result))
  )
}

#' Tabulation of a GEE Model
#'
#' Functions to produce tables from a fitted GEE produced with [fit_gee()].
#'
#' @name tabulate_gee
#' @returns The functions have different purposes:
#'   - `as.rtable()` returns either the coefficient table or the covariance matrix as an
#'   `rtables` object.
#'   - `s_lsmeans_logistic()` returns several least square mean statistics from the GEE.
#'   - `a_lsmeans_logistic()` is the formatted analysis function and returns the formatted statistics.
#'   - `summarize_gee_logistic()` is the analyze function and returns the modified `rtables` layout.
NULL

#' @importFrom tern as.rtable
#' @export
tern::as.rtable

#' @exportS3Method
#' @describeIn tabulate_gee Extracts the coefficient table or covariance matrix estimate from a `tern_gee` object.
#' @inheritParams tern::as.rtable
#' @param type (`character`)\cr type of table to extract from `tern_gee` object.
as.rtable.tern_gee <- function(x, # nolint
                               type = c("coef", "cov"),
                               ...) {
  type <- match.arg(type)
  switch(type,
    coef = h_gee_coef(x, ...),
    cov = h_gee_cov(x, ...)
  )
}

#' @keywords internal
h_gee_coef <- function(x, format = "xx.xxxx", conf_level = 0.95, ...) {
  fixed_table <- as.data.frame(stats::coef(summary(x)))
  assert_number(conf_level, lower = 0.001, upper = 0.999)

  fixed_table[["Std. Error"]] <- fixed_table[["Robust S.E."]]
  fixed_table[["z value"]] <- fixed_table[["Robust z"]]
  fixed_table[["Pr(>|z|)"]] <- 2 * stats::pnorm(abs(fixed_table[["z value"]]), lower.tail = FALSE)
  q <- stats::qnorm((1 + conf_level) / 2)
  ci_string <- tern::f_conf_level(conf_level)
  lower_string <- paste("Lower", ci_string)
  upper_string <- paste("Upper", ci_string)
  fixed_table[[lower_string]] <- fixed_table$Estimate - q * fixed_table[["Std. Error"]]
  fixed_table[[paste("Upper", ci_string)]] <- fixed_table$Estimate + q * fixed_table[["Std. Error"]]

  est_se_ci_table <- as.rtable(
    fixed_table[, c("Estimate", "Std. Error", lower_string, upper_string)],
    format = format
  )
  z_table <- as.rtable(fixed_table[, c("z value"), drop = FALSE], format = format)
  pvalue_table <- as.rtable(fixed_table[, "Pr(>|z|)", drop = FALSE], format = "x.xxxx | (<0.0001)")

  cbind_rtables(est_se_ci_table, z_table, pvalue_table)
}

#' @keywords internal
h_gee_cov <- function(x, format = "xx.xxxx") {
  cov_estimate <- VarCorr(x)
  as.rtable(as.data.frame(cov_estimate), format = format)
}

# lsmeans_logistic ----

#' @describeIn tabulate_gee Statistics function which extracts estimates from a
#'   [lsmeans()] data frame based on a logistic GEE model.
#'
#' @param df (`data.frame`)\cr data set resulting from [lsmeans()].
#' @param .in_ref_col (`logical`)\cr `TRUE` when working with the reference level, `FALSE` otherwise.
#'
#' @export
#'
#' @examples
#' library(dplyr)
#'
#' df <- fev_data %>%
#'   mutate(AVAL = as.integer(fev_data$FEV1 > 30))
#' df_counts <- df %>%
#'   select(USUBJID, ARMCD) %>%
#'   unique()
#'
#' lsmeans_df <- lsmeans(fit_gee(vars = vars_gee(arm = "ARMCD"), data = df))
#'
#' s_lsmeans_logistic(lsmeans_df[1, ], .in_ref_col = TRUE)
#'
#' s_lsmeans_logistic(lsmeans_df[2, ], .in_ref_col = FALSE)
s_lsmeans_logistic <- function(df, .in_ref_col) {
  if_not_ref <- function(x) `if`(.in_ref_col, character(), x)
  list(
    n = df$n,
    adj_prop_se = c(df$prop_est, df$prop_est_se), # to be confirmed
    adj_prop_ci = formatters::with_label(c(df$prop_lower_cl, df$prop_upper_cl), f_conf_level(df$conf_level)),
    odds_ratio_est = if_not_ref(df$or_est),
    odds_ratio_ci = formatters::with_label(
      if_not_ref(c(df$or_lower_cl, df$or_upper_cl)),
      f_conf_level(df$conf_level)
    ),
    log_odds_ratio_est = if_not_ref(df$log_or_est),
    log_odds_ratio_ci = formatters::with_label(
      if_not_ref(c(df$log_or_lower_cl, df$log_or_upper_cl)),
      f_conf_level(df$conf_level)
    )
  )
}

## a_lsmeans_logistic ----

#' @describeIn tabulate_gee Formatted Analysis function which can be further customized by calling
#'   [rtables::make_afun()] on it. It is used as `afun` in [rtables::analyze()].
#'
#' @export
a_lsmeans_logistic <- make_afun(
  s_lsmeans_logistic,
  .labels = c(
    adj_prop_se = "Adjusted Mean Proportion (SE)",
    odds_ratio_est = "Odds Ratio",
    log_odds_ratio_est = "Log Odds Ratio"
  ),
  .formats = c(
    n = "xx.",
    adj_prop_se = sprintf_format("%.2f (%.2f)"),
    adj_prop_ci = "(xx.xx, xx.xx)",
    odds_ratio_est = "xx.xx",
    odds_ratio_ci = "(xx.xx, xx.xx)",
    log_odds_ratio_est = "xx.xx",
    log_odds_ratio_ci = "(xx.xx, xx.xx)"
  ),
  .indent_mods = c(
    adj_prop_ci = 1L,
    odds_ratio_ci = 1L,
    log_odds_ratio_ci = 1L
  ),
  .null_ref_cells = FALSE
)

# Note: In production it would be nice to allow an S3 dispatch according to the
# class of the lsmeans input, however for now in the prototype we keep it simple.
# see later then to tern::summarize_variables for how to do that.

#' @describeIn tabulate_gee Analyze function for tabulating least-squares means estimates
#'   from logistic GEE least square mean results.
#'
#' @param lyt (`layout`)\cr input layout where analyses will be added to.
#' @param table_names (`character`)\cr this can be customized in case that the same `vars`
#'   are analyzed multiple times, to avoid warnings from `rtables`.
#' @param .stats (`character`)\cr statistics to select for the table.
#' @param .formats (named `character` or `list`)\cr formats for the statistics.
#' @param .indent_mods (named `integer`)\cr indent modifiers for the labels.
#' @param .labels (named `character`)\cr labels for the statistics (without indent).
#'
#' @export
#'
#' @examples
#' basic_table() %>%
#'   split_cols_by("ARMCD") %>%
#'   add_colcounts() %>%
#'   summarize_gee_logistic(
#'     .in_ref_col = FALSE
#'   ) %>%
#'   build_table(lsmeans_df, alt_counts_df = df_counts)
summarize_gee_logistic <- function(lyt,
                                   ...,
                                   table_names = "lsmeans_logistic_summary",
                                   .stats = NULL,
                                   .formats = NULL,
                                   .indent_mods = NULL,
                                   .labels = NULL) {
  afun <- make_afun(
    a_lsmeans_logistic,
    .stats = .stats,
    .formats = .formats,
    .indent_mods = .indent_mods,
    .labels = .labels
  )
  analyze(
    lyt = lyt,
    vars = "n",
    afun = afun,
    table_names = table_names,
    extra_args = list(...)
  )
}

1		#' Methods for GEE Models
2		#'
3		#' Additional methods which can simplify working with the GEE result object.
4		#' @name gee_methods
5		#' @returns `VarCorr()` returns the estimated covariance matrix, and
6		#' `QIC()` returns the QIC value.
7		NULL
8
9		#' @rdname gee_methods
10		#' @importFrom nlme VarCorr
11		#'
12		#' @param x (`tern_gee`)\cr result of [fit_gee()].
13		#' @inheritParams nlme::VarCorr
14		#'
15		#' @exportS3Method
16		VarCorr.tern_gee <- function(x, sigma = 1, ...) { # nolint
17	6x	dim_mat <- length(x$visit_levels)
18	6x	tmp <- id_mat <- diag(dim_mat)
19	6x	corest <- x$geese$alpha
20
21		# Start with lower-triangular matrix part.
22	6x	lower_mat <- switch(x$corstr,
23	6x	unstructured = , # Since this is the same as exchangeable, we can do this.
24	6x	exchangeable = {
25	4x	tmp[lower.tri(tmp)] <- corest
26	4x	tmp
27		},
28	6x	ar1 = {
29	1x	row_col_diff <- row(tmp) - col(tmp)
30	1x	tmp[lower.tri(tmp)] <- corest^(row_col_diff[lower.tri(row_col_diff)])
31	1x	tmp
32		},
33	6x	`m-dependent` = {
34	1x	row_col_diff <- row(tmp) - col(tmp)
35	1x	tmp[lower.tri(tmp)] <- corest[row_col_diff[lower.tri(row_col_diff)]]
36	1x	tmp
37		}
38		)
39
40		# Construct the full symmetric matrix.
41	6x	mat <- lower_mat + t(lower_mat) - id_mat
42	6x	rownames(mat) <- colnames(mat) <- x$visit_levels
43	6x	mat
44		}
45
46		#' @rdname gee_methods
47		#' @importFrom geepack QIC
48		#'
49		#' @param object (`tern_gee`)\cr result of [fit_gee()].
50		#' @inheritParams geepack::QIC
51		#'
52		#' @exportS3Method
53		QIC.tern_gee <- function(object, ...) { # nolint
54	!	object$qic
55		}

1		#' Set Variables to Use in GEE Model
2		#'
3		#' @param response (`character`)\cr name of response variable.
4		#' @param covariates (`character`)\cr vector of names of variables to use as covariates.
5		#' @param id (`character`)\cr name of variable to use to identify unique IDs.
6		#' @param arm (`character`)\cr name of arm variable.
7		#' @param visit (`character`)\cr name of visit variable.
8		#'
9		#' @return A list of variables that can be used as the `vars` argument in [fit_gee()].
10		#' @export
11		#'
12		#' @examples
13		#' vars_gee()
14		#'
15		#' vars_gee(
16		#' response = "CHG",
17		#' covariates = c("SEX", "RACE"),
18		#' id = "SUBJID",
19		#' arm = "ARMCD",
20		#' visit = "AVISITN"
21		#' )
22		vars_gee <- function(response = "AVAL",
23		covariates = c(),
24		id = "USUBJID",
25		arm = "ARM",
26		visit = "AVISIT") {
27	3x	list(
28	3x	response = response,
29	3x	covariates = covariates,
30	3x	id = id,
31	3x	arm = arm,
32	3x	visit = visit
33		)
34		}
35
36		#' @keywords internal
37		build_formula <- function(vars) {
38	14x	assert_list(vars)
39	14x	arm_part <- if (is.null(vars$arm)) NULL else vars$arm
40	14x	rhs_formula <- paste(
41	14x	c(arm_part, vars$covariates),
42	14x	collapse = " + "
43		)
44	14x	stats::as.formula(paste(
45	14x	vars$response,
46		"~",
47	14x	rhs_formula
48		))
49		}
50
51		#' @keywords internal
52		build_family <- function(regression) {
53	14x	assert_string(regression)
54
55	14x	result_object <- switch(regression,
56	14x	logistic = stats::binomial(link = "logit"),
57	14x	stop(paste("regression type", regression, "not supported"))
58		)
59
60	13x	result_class <- paste0("tern_gee_", regression)
61
62	13x	list(
63	13x	object = result_object,
64	13x	class = result_class,
65	13x	control = geeasy::geelm.control(scale.fix = TRUE)
66		)
67		}
68
69		#' @keywords internal
70		build_cor_details <- function(cor_str, vars, data) {
71	14x	assert_string(cor_str)
72	14x	assert_list(vars)
73	14x	assert_data_frame(data)
74
75	14x	result_str <- switch(cor_str,
76	14x	"unstructured" = "unstructured",
77	14x	"toeplitz" = "m-dependent",
78	14x	"compound symmetry" = "exchangeable",
79	14x	"auto-regressive" = "ar1",
80	14x	stop(paste("correlation structure", cor_str, "not available"))
81		)
82
83	13x	result_mv <- switch(cor_str,
84	13x	"unstructured" = 1,
85	13x	"toeplitz" = nlevels(data[[vars$visit]]) - 1,
86	13x	"compound symmetry" = 1,
87	13x	"auto-regressive" = 1
88		)
89
90	13x	list(
91	13x	str = result_str,
92	13x	mv = result_mv
93		)
94		}
95
96		#' @keywords internal
97		order_data <- function(data, vars) {
98	14x	assert_data_frame(data)
99	14x	assert_list(vars)
100
101	14x	if (is.character(data[[vars$visit]])) {
102	1x	message(paste("visit variable", vars$visit, "will be coerced to factor for ordering"))
103	1x	message("order is:")
104	1x	data[[vars$visit]] <- factor(data[[vars$visit]])
105	1x	message(paste(toString(levels(data[[vars$visit]])), "\n"))
106		}
107	14x	if (is.factor(data[[vars$id]]) \|\| is.character(data[[vars$id]])) {
108	14x	data[[vars$id]] <- as.integer(as.factor(data[[vars$id]]))
109		}
110	14x	assert_numeric(data[[vars$id]])
111
112	14x	right_order <- order(data[[vars$id]], data[[vars$visit]])
113	14x	data[right_order, ]
114		}
115
116		#' Fit a GEE Model
117		#'
118		#' @param vars (`list`)\cr see [vars_gee()].
119		#' @param data (`data.frame`)\cr input data.
120		#' @param regression (`string`)\cr choice of regression model.
121		#' @param cor_struct (`string`)\cr assumed correlation structure.
122		#'
123		#' @details The correlation structure can be:
124		#' * `unstructured`: No constraints are placed on the correlations.
125		#' * `toeplitz`: Assumes a banded correlation structure, i.e. the correlation
126		#' between two time points depends on the distance between the time indices.
127		#' * `compound symmetry`: Constant correlation between all time points.
128		#' * `auto-regressive`: Auto-regressive order 1 correlation matrix.
129		#'
130		#' @return Object of class `tern_gee` as well as specific to the kind of regression
131		#' which was used.
132		#' @export
133		#'
134		#' @examples
135		#' df <- fev_data
136		#' df$AVAL <- as.integer(fev_data$FEV1 > 30)
137		#'
138		#' fit_gee(vars = vars_gee(arm = "ARMCD"), data = df)
139		#'
140		#' fit_gee(vars = vars_gee(arm = "ARMCD"), data = df, cor_struct = "compound symmetry")
141		fit_gee <- function(vars = vars_gee(),
142		data,
143		regression = c("logistic"),
144		cor_struct = c("unstructured", "toeplitz", "compound symmetry", "auto-regressive")) {
145	12x	formula <- build_formula(vars)
146
147	12x	regression <- match.arg(regression)
148	12x	family <- build_family(regression)
149
150	12x	data <- order_data(data, vars)
151	12x	data[[".id"]] <- data[[vars$id]]
152	12x	data[[".waves"]] <- as.integer(data[[vars$visit]])
153
154	12x	cor_struct <- match.arg(cor_struct)
155	12x	cor_details <- build_cor_details(cor_struct, vars, data)
156
157	12x	fit <- geeasy::geelm(
158	12x	formula = formula,
159	12x	id = .id,
160	12x	waves = .waves,
161	12x	data = data,
162	12x	family = family$object,
163	12x	corstr = cor_details$str,
164	12x	Mv = cor_details$mv,
165	12x	control = family$control
166		)
167
168	12x	fit$qic <- geepack::QIC(fit)
169	12x	fit$visit_levels <- levels(data[[vars$visit]])
170	12x	fit$vars <- vars
171	12x	fit$data <- data
172	12x	assert_factor(data[[vars$arm]])
173	12x	fit$ref_level <- levels(data[[vars$arm]])[1L]
174
175	12x	structure(
176	12x	fit,
177	12x	class = c(family$class, "tern_gee", class(fit))
178		)
179		}

1		#' Extract Least Square Means from a GEE Model
2		#'
3		#' @param object (`tern_gee`)\cr result of [fit_gee()].
4		#' @param conf_level (`proportion`)\cr confidence level
5		#' @param weights (`string`)\cr type of weights to be used for the least square means,
6		#' see [emmeans::emmeans()] for details.
7		#' @param specs (`string` or `formula`) specifications passed to [emmeans::emmeans()]
8		#' @param ... additional arguments for methods
9		#'
10		#' @return A `data.frame` with least-square means and contrasts. Additional
11		#' classes allow to dispatch downstream methods correctly, too.
12		#' @export
13		#'
14		#' @examples
15		#' df <- fev_data
16		#' df$AVAL <- rbinom(n = nrow(df), size = 1, prob = 0.5)
17		#' fit <- fit_gee(vars = vars_gee(arm = "ARMCD"), data = df)
18		#'
19		#' lsmeans(fit)
20		#'
21		#' lsmeans(fit, conf_level = 0.90, weights = "equal")
22		lsmeans <- function(object,
23		conf_level = 0.95,
24		weights = "proportional",
25		specs = object$vars$arm,
26		...) {
27	5x	UseMethod("lsmeans", object)
28		}
29
30		#' @rdname lsmeans
31		#' @exportS3Method
32		lsmeans.tern_gee_logistic <- function(object,
33		conf_level = 0.95,
34		weights = "proportional",
35		specs = object$vars$arm,
36		...) {
37	5x	prop_emm <- emmeans::emmeans(
38	5x	object = object,
39	5x	specs = specs,
40	5x	weights = weights,
41	5x	type = "response",
42	5x	data = object$data
43		)
44
45	5x	prop_df <- cbind(
46	5x	data.frame(stats::confint(prop_emm))[, c(object$vars$arm, "prob", "SE", "lower.CL", "upper.CL")],
47	5x	n = as.list(prop_emm)$extras[, ".wgt."]
48		)
49	5x	names(prop_df) <- c(object$vars$arm, "prop_est", "prop_est_se", "prop_lower_cl", "prop_upper_cl", "n")
50	5x	ref_level <- levels(object$data[[object$vars$arm]])[1L]
51
52	5x	or_emm <- stats::confint(
53	5x	graphics::pairs(prop_emm, reverse = TRUE),
54	5x	level = conf_level
55		)
56	5x	or_df <- as.data.frame(or_emm)
57	5x	or_df$comparator <- gsub(pattern = ".+ / (.+)", replacement = "\\1", x = or_df$contrast)
58	5x	or_df[[object$vars$arm]] <- gsub(pattern = "(.+) / .+", replacement = "\\1", x = or_df$contrast)
59	5x	or_df <- or_df[or_df$comparator == ref_level, ]
60	5x	or_df <- rbind(NA, or_df)
61	5x	or_df[1L, object$vars$arm] <- ref_level
62	5x	or_df <- or_df[, c(object$vars$arm, "odds.ratio", "lower.CL", "upper.CL")]
63	5x	or_df <- cbind(or_df, log(or_df[, -1L]), conf_level)
64	5x	names(or_df) <- c(
65	5x	object$vars$arm,
66	5x	"or_est", "or_lower_cl", "or_upper_cl",
67	5x	"log_or_est", "log_or_lower_cl", "log_or_upper_cl",
68	5x	"conf_level"
69		)
70
71	5x	result <- merge(prop_df, or_df, by = object$vars$arm)
72
73	5x	structure(
74	5x	result,
75	5x	class = c("lsmeans_logistic", class(result))
76		)
77		}

1		#' Tabulation of a GEE Model
2		#'
3		#' Functions to produce tables from a fitted GEE produced with [fit_gee()].
4		#'
5		#' @name tabulate_gee
6		#' @returns The functions have different purposes:
7		#' - `as.rtable()` returns either the coefficient table or the covariance matrix as an
8		#' `rtables` object.
9		#' - `s_lsmeans_logistic()` returns several least square mean statistics from the GEE.
10		#' - `a_lsmeans_logistic()` is the formatted analysis function and returns the formatted statistics.
11		#' - `summarize_gee_logistic()` is the analyze function and returns the modified `rtables` layout.
12		NULL
13
14		#' @importFrom tern as.rtable
15		#' @export
16		tern::as.rtable
17
18		#' @exportS3Method
19		#' @describeIn tabulate_gee Extracts the coefficient table or covariance matrix estimate from a `tern_gee` object.
20		#' @inheritParams tern::as.rtable
21		#' @param type (`character`)\cr type of table to extract from `tern_gee` object.
22		as.rtable.tern_gee <- function(x, # nolint
23		type = c("coef", "cov"),
24		...) {
25	1x	type <- match.arg(type)
26	1x	switch(type,
27	!	coef = h_gee_coef(x, ...),
28	1x	cov = h_gee_cov(x, ...)
29		)
30		}
31
32		#' @keywords internal
33		h_gee_coef <- function(x, format = "xx.xxxx", conf_level = 0.95, ...) {
34	!	fixed_table <- as.data.frame(stats::coef(summary(x)))
35	!	assert_number(conf_level, lower = 0.001, upper = 0.999)
36
37	!	fixed_table[["Std. Error"]] <- fixed_table[["Robust S.E."]]
38	!	fixed_table[["z value"]] <- fixed_table[["Robust z"]]
39	!	fixed_table[["Pr(>\|z\|)"]] <- 2 * stats::pnorm(abs(fixed_table[["z value"]]), lower.tail = FALSE)
40	!	q <- stats::qnorm((1 + conf_level) / 2)
41	!	ci_string <- tern::f_conf_level(conf_level)
42	!	lower_string <- paste("Lower", ci_string)
43	!	upper_string <- paste("Upper", ci_string)
44	!	fixed_table[[lower_string]] <- fixed_table$Estimate - q * fixed_table[["Std. Error"]]
45	!	fixed_table[[paste("Upper", ci_string)]] <- fixed_table$Estimate + q * fixed_table[["Std. Error"]]
46
47	!	est_se_ci_table <- as.rtable(
48	!	fixed_table[, c("Estimate", "Std. Error", lower_string, upper_string)],
49	!	format = format
50		)
51	!	z_table <- as.rtable(fixed_table[, c("z value"), drop = FALSE], format = format)
52	!	pvalue_table <- as.rtable(fixed_table[, "Pr(>\|z\|)", drop = FALSE], format = "x.xxxx \| (<0.0001)")
53
54	!	cbind_rtables(est_se_ci_table, z_table, pvalue_table)
55		}
56
57		#' @keywords internal
58		h_gee_cov <- function(x, format = "xx.xxxx") {
59	2x	cov_estimate <- VarCorr(x)
60	2x	as.rtable(as.data.frame(cov_estimate), format = format)
61		}
62
63		# lsmeans_logistic ----
64
65		#' @describeIn tabulate_gee Statistics function which extracts estimates from a
66		#' [lsmeans()] data frame based on a logistic GEE model.
67		#'
68		#' @param df (`data.frame`)\cr data set resulting from [lsmeans()].
69		#' @param .in_ref_col (`logical`)\cr `TRUE` when working with the reference level, `FALSE` otherwise.
70		#'
71		#' @export
72		#'
73		#' @examples
74		#' library(dplyr)
75		#'
76		#' df <- fev_data %>%
77		#' mutate(AVAL = as.integer(fev_data$FEV1 > 30))
78		#' df_counts <- df %>%
79		#' select(USUBJID, ARMCD) %>%
80		#' unique()
81		#'
82		#' lsmeans_df <- lsmeans(fit_gee(vars = vars_gee(arm = "ARMCD"), data = df))
83		#'
84		#' s_lsmeans_logistic(lsmeans_df[1, ], .in_ref_col = TRUE)
85		#'
86		#' s_lsmeans_logistic(lsmeans_df[2, ], .in_ref_col = FALSE)
87		s_lsmeans_logistic <- function(df, .in_ref_col) {
88	2x	if_not_ref <- function(x) `if`(.in_ref_col, character(), x)
89	2x	list(
90	2x	n = df$n,
91	2x	adj_prop_se = c(df$prop_est, df$prop_est_se), # to be confirmed
92	2x	adj_prop_ci = formatters::with_label(c(df$prop_lower_cl, df$prop_upper_cl), f_conf_level(df$conf_level)),
93	2x	odds_ratio_est = if_not_ref(df$or_est),
94	2x	odds_ratio_ci = formatters::with_label(
95	2x	if_not_ref(c(df$or_lower_cl, df$or_upper_cl)),
96	2x	f_conf_level(df$conf_level)
97		),
98	2x	log_odds_ratio_est = if_not_ref(df$log_or_est),
99	2x	log_odds_ratio_ci = formatters::with_label(
100	2x	if_not_ref(c(df$log_or_lower_cl, df$log_or_upper_cl)),
101	2x	f_conf_level(df$conf_level)
102		)
103		)
104		}
105
106		## a_lsmeans_logistic ----
107
108		#' @describeIn tabulate_gee Formatted Analysis function which can be further customized by calling
109		#' [rtables::make_afun()] on it. It is used as `afun` in [rtables::analyze()].
110		#'
111		#' @export
112		a_lsmeans_logistic <- make_afun(
113		s_lsmeans_logistic,
114		.labels = c(
115		adj_prop_se = "Adjusted Mean Proportion (SE)",
116		odds_ratio_est = "Odds Ratio",
117		log_odds_ratio_est = "Log Odds Ratio"
118		),
119		.formats = c(
120		n = "xx.",
121		adj_prop_se = sprintf_format("%.2f (%.2f)"),
122		adj_prop_ci = "(xx.xx, xx.xx)",
123		odds_ratio_est = "xx.xx",
124		odds_ratio_ci = "(xx.xx, xx.xx)",
125		log_odds_ratio_est = "xx.xx",
126		log_odds_ratio_ci = "(xx.xx, xx.xx)"
127		),
128		.indent_mods = c(
129		adj_prop_ci = 1L,
130		odds_ratio_ci = 1L,
131		log_odds_ratio_ci = 1L
132		),
133		.null_ref_cells = FALSE
134		)
135
136		# Note: In production it would be nice to allow an S3 dispatch according to the
137		# class of the lsmeans input, however for now in the prototype we keep it simple.
138		# see later then to tern::summarize_variables for how to do that.
139
140		#' @describeIn tabulate_gee Analyze function for tabulating least-squares means estimates
141		#' from logistic GEE least square mean results.
142		#'
143		#' @param lyt (`layout`)\cr input layout where analyses will be added to.
144		#' @param table_names (`character`)\cr this can be customized in case that the same `vars`
145		#' are analyzed multiple times, to avoid warnings from `rtables`.
146		#' @param .stats (`character`)\cr statistics to select for the table.
147		#' @param .formats (named `character` or `list`)\cr formats for the statistics.
148		#' @param .indent_mods (named `integer`)\cr indent modifiers for the labels.
149		#' @param .labels (named `character`)\cr labels for the statistics (without indent).
150		#'
151		#' @export
152		#'
153		#' @examples
154		#' basic_table() %>%
155		#' split_cols_by("ARMCD") %>%
156		#' add_colcounts() %>%
157		#' summarize_gee_logistic(
158		#' .in_ref_col = FALSE
159		#' ) %>%
160		#' build_table(lsmeans_df, alt_counts_df = df_counts)
161		summarize_gee_logistic <- function(lyt,
162		...,
163		table_names = "lsmeans_logistic_summary",
164		.stats = NULL,
165		.formats = NULL,
166		.indent_mods = NULL,
167		.labels = NULL) {
168	2x	afun <- make_afun(
169	2x	a_lsmeans_logistic,
170	2x	.stats = .stats,
171	2x	.formats = .formats,
172	2x	.indent_mods = .indent_mods,
173	2x	.labels = .labels
174		)
175	2x	analyze(
176	2x	lyt = lyt,
177	2x	vars = "n",
178	2x	afun = afun,
179	2x	table_names = table_names,
180	2x	extra_args = list(...)
181		)
182		}