From b66d6d8831ec2f0d9965eb0001dcb1506a056374 Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 28 Apr 2024 15:45:06 +0000 Subject: [PATCH 01/12] add in all draft functions --- .Rbuildignore | 3 + .gitignore | 3 + .../package_structure/all_draft_functions.R | 221 +++++++++++ .../{design.Rmd => design_lml.Rmd} | 0 design/package_structure/design_xgc.R | 355 ++++++++++++++++++ design/package_structure/design_xgc.Rmd | 232 ++++++++++++ .../test_all_draft_functions.R | 42 +++ 7 files changed, 856 insertions(+) create mode 100644 design/package_structure/all_draft_functions.R rename design/package_structure/{design.Rmd => design_lml.Rmd} (100%) create mode 100644 design/package_structure/design_xgc.R create mode 100644 design/package_structure/design_xgc.Rmd create mode 100644 design/package_structure/test_all_draft_functions.R diff --git a/.Rbuildignore b/.Rbuildignore index c0f0cf0..8cd5db0 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -17,3 +17,6 @@ coverage.* init.sh workflows.md images +^RobinCar2\.Rcheck$ +^RobinCar2.*\.tar\.gz$ +^RobinCar2.*\.tgz$ diff --git a/.gitignore b/.gitignore index 7065c0b..d381fc6 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,6 @@ coverage.* node_modules/ package-lock.json package.json +RobinCar2.Rcheck/ +RobinCar2*.tar.gz +RobinCar2*.tgz diff --git a/design/package_structure/all_draft_functions.R b/design/package_structure/all_draft_functions.R new file mode 100644 index 0000000..cd20e07 --- /dev/null +++ b/design/package_structure/all_draft_functions.R @@ -0,0 +1,221 @@ + +get_eff_measure <- function(fit) { + if (identical(class(fit), "lm")) { + out_eff <- "diff" + } else if ("glm" %in% class(fit)) { + use_family <- fit$family$family + use_link <- fit$family$link + if (use_family == "binomial" & use_link == "logit") { + out_eff <- "odds ratio" + } else if (use_family == "binomial" & use_link == "identity") { + out_eff <- "diff" + } else if (use_family %in% c("binomial", "poisson") & use_link == "log") { + out_eff <- "risk ratio" + } else { + stop("pity, specified family and link function in glm() are not supported, may be in a future version...") + } + } else { + stop("fit object should be either from stat::lm() or stat::glm()") + } + out_eff +} + +### g-computation and its VCOV ------------------------------------------------- + +# func: prepare data +get_countfact_pred <- function(fit, trt_var_name){ + treatments <- fit$xlevels[[trt_var_name]] + use_dat <- fit$model + + dat_counterfact <- sapply(treatments, function(t) { + tmp <- use_dat + tmp$y_observed <- use_dat[, 1] + tmp$trt_observed <- use_dat[[trt_var_name]] + tmp$trt_counterfact <- t + tmp + }, simplify = FALSE) + + out_counterfact <- sapply(treatments, function(t) { + outtmp <- tmp <- dat_counterfact[[t]] + tmp[[trt_var_name]] <- tmp$trt_counterfact + outtmp$mu_t <- predict(fit, newdata = tmp, type = "response") + outtmp$theta_t <- mean(outtmp$mu_t) + outtmp$resid_mu_t <- outtmp[, 1] - outtmp$mu_t + outtmp$resid_theta_t <- outtmp[, 1] - outtmp$mu_t + outtmp + }, simplify = FALSE) + + out_counterfact +} + +# func: calculate V +calculate_V <- function(t, s, counterfact_res) { + if (t == s) { # : calculate v_tt + # v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + + # 2*cov(Y_j,mu_t_j | j \in I_t) - + # var(mu_t) + calres <- counterfact_res[[t]] + + n <- nrow(calres) + It_ind <- (calres$trt_observed == t) + pi_t <- sum(It_ind) / n + + var_y_mut_It <- var(calres$resid_mu_t[It_ind]) + cov_y_mut_It <- cov(calres$y_observed[It_ind], calres$mu_t[It_ind]) + var_mut <- var(calres$mu_t) + + v_out <- var_y_mut_It / pi_t + 2 * cov_y_mut_It - var_mut + } else { # : calculate v_ts + # v_ts = cov(Y_j,mu_s_j | j \in I_t) + + # cov(Y_j,mu_t_j | j \in I_s) - + # cov(mu_t,mu_s) + + calres_t <- counterfact_res[[t]] + calres_s <- counterfact_res[[s]] + + n <- nrow(calres_t) + It_ind <- (calres_s$trt_observed == t) + Is_ind <- (calres_t$trt_observed == s) + pi_t <- sum(It_ind) / n + pi_s <- sum(Is_ind) / n + + cov_y_mus_It <- cov(calres_s$y_observed[It_ind], calres_s$mu_t[It_ind]) + cov_y_mut_Is <- cov(calres_t$y_observed[Is_ind], calres_t$mu_t[Is_ind]) + cov_mut_mus <- cov(calres_s$mu_t, calres_t$mu_t) + + v_out <- cov_y_mus_It + cov_y_mut_Is - cov_mut_mus + } + v_out +} + +vectorized_calculate_V <- Vectorize(calculate_V, c("t", "s")) + +# func: calculate n^(-1)f'Vf +op <- function(x, square = FALSE) (x * (1 - x))^ifelse(square, 2, 1) +odds <- function(x) (x / (1 - x)) + +calculate_fVf <- function(t, s, theta, V, n, eff_measure) { + t_ind <- which(names(theta) == t) + s_ind <- which(names(theta) == s) + + # depneding on eff_measure (1.diff; 2.risk ratio; 3.odds ratio) + # 1. v_tt - 2*v_ts + v_ss + # 2. v_tt/(theta_t^2) - 2*v_ts/(theta_t*theta_s) + v_ss/(theta_s^2) + # 3. v_tt/(op(theta_t)^2) - 2*v_ts/(op(theta_t)*op(theta_s)) + v_ss/(op(theta_s)^2) + if (eff_measure == "diff") { + outtmp <- V[t_ind, t_ind] + V[s_ind, s_ind] - 2 * V[t_ind, s_ind] + } else if (eff_measure == "risk ratio") { + outtmp <- V[t_ind, t_ind] / (theta[t_ind]^2) + + V[s_ind, s_ind] / (theta[s_ind]^2) - + 2 * V[t_ind, s_ind] / (theta[t_ind] * theta[s_ind]) + } else if (eff_measure == "odds ratio") { + outtmp <- V[t_ind, t_ind] / op(theta[t_ind], square = T) + + V[s_ind, s_ind] / op(theta[s_ind], square = T) - + 2 * V[t_ind, s_ind] / (op(theta[t_ind]) * op(theta[s_ind])) + } + outtmp / sqrt(n) +} + +vectorized_calculate_fVf <- Vectorize(calculate_fVf, c("t", "s")) + +# func: wrapper function given fit and pred_counterfact +calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ + treatments <- fit$xlevels[[trt_var_name]] + + theta_t <- sapply(counterfact_pred, function(dat_t) unique(dat_t$theta_t)) + + V <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_V", + counterfact_res = counterfact_pred + ) + rownames(V) <- colnames(V) <- treatments + + n_fVf <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_fVf", + theta = theta_t, + V = V, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit) + ) + rownames(n_fVf) <- colnames(n_fVf) <- treatments + + out <- list( theta = theta_t, + V = V, + f_vcov = n_fVf, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit), + trt_var_name = trt_var_name, + fit = fit ) + class(out) <- c(class(out),"vcov_rc") + out +} + +# func: reporting function +calculate_f <- function(theta_t, theta_s, eff_measure) { + switch(eff_measure, + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) + ) +} + +report <- function(fit, trt_var_name, theta, f_vcov, digits = 3) { + eff_measure <- get_eff_measure(fit) + n <- nrow(fit$model) + treatments <- fit$xlevels[[trt_var_name]] + n.arm <- summary(fit$model[[trt_var_name]]) + n.arm.perc <- round(n.arm * 100 / n, 1) |> format(nsmall = 1) + + out_cases <- expand.grid(seq_along(theta), seq_along(theta)) + out_cases <- out_cases[out_cases[, 1] != out_cases[, 2], , drop = FALSE] + t1 <- out_cases[, 1] + t0 <- out_cases[, 2] + + outtmp <- data.frame( + Total = n, + Arm_1 = treatments[t1], + N_Arm_1 = paste0(n.arm[t1], "(", n.arm.perc[t1], "%)"), + Arm_0 = treatments[t0], + N_Arm_0 = paste0(n.arm[t0], "(", n.arm.perc[t0], "%)"), + Eff_measure = eff_measure, + Estimate = calculate_f(theta[t1], theta[t0], eff_measure) + ) + outtmp[["S.E."]] <- sqrt(f_vcov[cbind(t1, t0)]) + outtmp[["95%CI_lower"]] <- outtmp$Estimate + qnorm(0.025) * outtmp[["S.E."]] + outtmp[["95%CI_upper"]] <- outtmp$Estimate + qnorm(0.975) * outtmp[["S.E."]] + + if (eff_measure == "diff") { + outtmp["Estimate"] <- round(outtmp["Estimate"], digits) |> format(nsmall = digits) + outtmp["S.E."] <- round(outtmp["S.E."], digits) |> format(nsmall = digits) + outtmp["95%CI_lower"] <- round(outtmp["95%CI_lower"], digits) |> format(nsmall = digits) + outtmp["95%CI_upper"] <- round(outtmp["95%CI_upper"], digits) |> format(nsmall = digits) + } else { + outtmp["S.E."] <- ((exp(outtmp["S.E."]) - 1) * exp(2 * outtmp["Estimate"] + outtmp["S.E."]^2)) |> + round(digits) |> + format(nsmall = digits) + outtmp["Estimate"] <- exp(outtmp["Estimate"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_lower"] <- exp(outtmp["95%CI_lower"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_upper"] <- exp(outtmp["95%CI_upper"]) |> + round(digits) |> + format(nsmall = digits) + } + outtmp +} + +report_fvcov <- function(result, digits = 3){ + report( + fit = result$fit, + trt_var_name = result$trt_var_name, + theta = result$theta, + f_vcov = result$f_vcov, + digits = digits + ) +} diff --git a/design/package_structure/design.Rmd b/design/package_structure/design_lml.Rmd similarity index 100% rename from design/package_structure/design.Rmd rename to design/package_structure/design_lml.Rmd diff --git a/design/package_structure/design_xgc.R b/design/package_structure/design_xgc.R new file mode 100644 index 0000000..1463fe1 --- /dev/null +++ b/design/package_structure/design_xgc.R @@ -0,0 +1,355 @@ +# Design note +# Author: X. Gregory Chen + +# [Refenrece] +# Ye, T., Bannick, M., Yi, Y. and Shao, J., 2023. Robust variance estimation for +# covariate-adjusted unconditional treatment effect in randomized clinical trials +# with binary outcomes. Statistical Theory and Related Fields, 7(2), pp.159-163. +# DOI: https://doi.org/10.1080/24754269.2023.2205802 + +# Notation in this script tries to follow the above reference +# (_hat used to denote an estimate is dropped for brevity, but shout be clear in the context) +# +# Y is a continuous or binary random variate +# +# E(Y|A,X) = g^(-1)( beta_A'A + beta_X'X ), +# where A is a dummy vector for treatment 1,..,or K +# X is a covariate vector to be adjusted +# beta is regression coef vector, in particular beta_A relates to treatment effect +# g(.) is the link function, depending on distribution of Y and interested eff measure +# +# mu_t = g^(-1)( beta_A'a_t + beta_X'X ), an estimated/predicted E(Y|A=t,X) per subject +# where t \in {1,...,K} +# a_t is a vector of length K with t^th element being 1 and 0 otherwise +# +# theta_t = mean(mu_t) +# +# f(theta;t,s) = +# 1. theta_t - theta_s (diff, risk diff) +# 2. theta_t / theta_s (risk ratio) +# 3. log(odds(theta_t) / odds(theta_s)) (odds ratio) +# where odds(p) = p/(1-p) +# t,s are treatment index taking distinct values from the set {1,...,K} +# +# asympt.cov of estimated f(theta;t,s) = vcov_f_core/sqrt(n) +# and asympt.cov of theta vcov_f_core = +# 1. v_tt - 2*v_ts + v_ss +# 2. v_tt/(theta_t^2) - 2*v_ts/(theta_t*theta_s) + v_ss/(theta_s^2) +# 3. v_tt/(op(theta_t)^2) - 2*v_ts/(op(theta_t)*op(theta_s)) + v_ss/(op(theta_s)^2) +# where n is the sample size +# op(p) = p*(1-p) +# +# let pi_t be the sample proportion of subjects received treatment t +# I_t = a subset of all subjects {1,...,n} that received treatment t +# +# v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + +# 2*cov(Y_j,mu_t_j | j \in I_t) - +# var(mu_t) +# +# v_ts = cov(Y_j,mu_s_j | j \in I_t) + +# cov(Y_j,mu_t_j | j \in I_s) - +# cov(mu_t,mu_s) + + +### prepare data and 'fit' ----------------------------------------------------- +# quote from "design_lml.Rmd" +n <- 200 +block_rand <- function(block = c(0, 0, 1, 1)) { + function(n) { + r <- lapply(seq_len(ceiling(n / length(block))), function(i) { + sample(block) + }) + unlist(r)[1:n] + } +} +by_strata <- function(f, strata) { + ret <- rep(NA, length(strata)) + for (x in split(seq_len(length(strata)), strata)) { + ret[x] <- f(length(x)) + } + return(ret) +} +sim_data <- function(n) { + cov1 <- rnorm(n) + z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) + z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) + permute_block <- c(0, 0, 1, 1) + trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) + trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) + df <- data.frame(trt, z1, z2, cov1) + x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) + coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) + theta <- x_mat %*% coef + y <- plogis(theta) + y_bin <- as.integer(runif(n) < y) + + df$y <- y_bin + df +} + +d <- sim_data(500) +fit <- glm(y ~ trt:z1 + trt, family = binomial(), data = d) + +### Actual start-point, with an input of +### an R object 'fit', from stat::lm() or stat::glm() + +### general helper -------------------------------------------------------------- + +get_eff_measure <- function(fit) { + if (identical(class(fit), "lm")) { + out_eff <- "diff" + } else if ("glm" %in% class(fit)) { + use_family <- fit$family$family + use_link <- fit$family$link + if (use_family == "binomial" & use_link == "logit") { + out_eff <- "odds ratio" + } else if (use_family == "binomial" & use_link == "identity") { + out_eff <- "diff" + } else if (use_family %in% c("binomial", "poisson") & use_link == "log") { + out_eff <- "risk ratio" + } else { + stop("pity, specified family and link function in glm() are not supported, may be in a future version...") + } + } else { + stop("fit object should be either from stat::lm() or stat::glm()") + } + out_eff +} + +### g-computation and its VCOV ------------------------------------------------- + +# func: prepare data +get_countfact_pred <- function(fit, trt_var_name){ + treatments <- fit$xlevels[[trt_var_name]] + use_dat <- fit$model + + dat_counterfact <- sapply(treatments, function(t) { + tmp <- use_dat + tmp$y_observed <- use_dat[, 1] + tmp$trt_observed <- use_dat[[trt_var_name]] + tmp$trt_counterfact <- t + tmp + }, simplify = FALSE) + + out_counterfact <- sapply(treatments, function(t) { + outtmp <- tmp <- dat_counterfact[[t]] + tmp[[trt_var_name]] <- tmp$trt_counterfact + outtmp$mu_t <- predict(fit, newdata = tmp, type = "response") + outtmp$theta_t <- mean(outtmp$mu_t) + outtmp$resid_mu_t <- outtmp[, 1] - outtmp$mu_t + outtmp$resid_theta_t <- outtmp[, 1] - outtmp$mu_t + outtmp + }, simplify = FALSE) + + out_counterfact +} + +# Step 1: get g-computation theta_t +pred_counterfact <- get_countfact_pred(fit,"trt") +theta_t <- sapply(pred_counterfact, function(dat_t) unique(dat_t$theta_t)) + +# func: calculate V +calculate_V <- function(t, s, counterfact_res) { + if (t == s) { # : calculate v_tt + # v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + + # 2*cov(Y_j,mu_t_j | j \in I_t) - + # var(mu_t) + calres <- counterfact_res[[t]] + + n <- nrow(calres) + It_ind <- (calres$trt_observed == t) + pi_t <- sum(It_ind) / n + + var_y_mut_It <- var(calres$resid_mu_t[It_ind]) + cov_y_mut_It <- cov(calres$y_observed[It_ind], calres$mu_t[It_ind]) + var_mut <- var(calres$mu_t) + + v_out <- var_y_mut_It / pi_t + 2 * cov_y_mut_It - var_mut + } else { # : calculate v_ts + # v_ts = cov(Y_j,mu_s_j | j \in I_t) + + # cov(Y_j,mu_t_j | j \in I_s) - + # cov(mu_t,mu_s) + + calres_t <- counterfact_res[[t]] + calres_s <- counterfact_res[[s]] + + n <- nrow(calres_t) + It_ind <- (calres_s$trt_observed == t) + Is_ind <- (calres_t$trt_observed == s) + pi_t <- sum(It_ind) / n + pi_s <- sum(Is_ind) / n + + cov_y_mus_It <- cov(calres_s$y_observed[It_ind], calres_s$mu_t[It_ind]) + cov_y_mut_Is <- cov(calres_t$y_observed[Is_ind], calres_t$mu_t[Is_ind]) + cov_mut_mus <- cov(calres_s$mu_t, calres_t$mu_t) + + v_out <- cov_y_mus_It + cov_y_mut_Is - cov_mut_mus + } + v_out +} + +vectorized_calculate_V <- Vectorize(calculate_V, c("t", "s")) + +# Substep (a) in Step 2: get V +treatments <- fit$xlevels[["trt"]] +V <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_V", + counterfact_res = pred_counterfact +) +rownames(V) <- colnames(V) <- treatments + +# func: calculate n^(-1)f'Vf +op <- function(x, square = FALSE) (x * (1 - x))^ifelse(square, 2, 1) +odds <- function(x) (x / (1 - x)) + +calculate_fVf <- function(t, s, theta, V, n, eff_measure) { + t_ind <- which(names(theta) == t) + s_ind <- which(names(theta) == s) + + # depneding on eff_measure (1.diff; 2.risk ratio; 3.odds ratio) + # 1. v_tt - 2*v_ts + v_ss + # 2. v_tt/(theta_t^2) - 2*v_ts/(theta_t*theta_s) + v_ss/(theta_s^2) + # 3. v_tt/(op(theta_t)^2) - 2*v_ts/(op(theta_t)*op(theta_s)) + v_ss/(op(theta_s)^2) + if (eff_measure == "diff") { + outtmp <- V[t_ind, t_ind] + V[s_ind, s_ind] - 2 * V[t_ind, s_ind] + } else if (eff_measure == "risk ratio") { + outtmp <- V[t_ind, t_ind] / (theta[t_ind]^2) + + V[s_ind, s_ind] / (theta[s_ind]^2) - + 2 * V[t_ind, s_ind] / (theta[t_ind] * theta[s_ind]) + } else if (eff_measure == "odds ratio") { + outtmp <- V[t_ind, t_ind] / op(theta[t_ind], square = T) + + V[s_ind, s_ind] / op(theta[s_ind], square = T) - + 2 * V[t_ind, s_ind] / (op(theta[t_ind]) * op(theta[s_ind])) + } + outtmp / sqrt(n) +} + +vectorized_calculate_fVf <- Vectorize(calculate_fVf, c("t", "s")) + +# Substep (b) in Step 2: get V +n_fVf <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_vcov_f", + theta = theta_t, + V = V, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit) +) +rownames(n_fVf) <- colnames(n_fVf) <- treatments + +# func: wrapper function given fit and pred_counterfact +calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ + treatments <- fit$xlevels[[trt_var_name]] + + theta_t <- sapply(counterfact_pred, function(dat_t) unique(dat_t$theta_t)) + + V <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_V", + counterfact_res = counterfact_pred + ) + rownames(V) <- colnames(V) <- treatments + + n_fVf <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_fVf", + theta = theta_t, + V = V, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit) + ) + rownames(n_fVf) <- colnames(n_fVf) <- treatments + + out <- list( theta = theta_t, + V = V, + f_vcov = n_fVf, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit), + trt_var_name = trt_var_name, + fit = fit ) + class(out) <- c(class(out),"vcov_rc") + out +} + +# func: reporting function +calculate_f <- function(theta_t, theta_s, eff_measure) { + switch(eff_measure, + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) + ) +} + +report <- function(fit, trt_var_name, theta, f_vcov, digits = 3) { + eff_measure <- get_eff_measure(fit) + n <- nrow(fit$model) + treatments <- fit$xlevels[[trt_var_name]] + n.arm <- summary(fit$model[[trt_var_name]]) + n.arm.perc <- round(n.arm * 100 / n, 1) |> format(nsmall = 1) + + out_cases <- expand.grid(seq_along(theta), seq_along(theta)) + out_cases <- out_cases[out_cases[, 1] != out_cases[, 2], , drop = FALSE] + t1 <- out_cases[, 1] + t0 <- out_cases[, 2] + + outtmp <- data.frame( + Total = n, + Arm_1 = treatments[t1], + N_Arm_1 = paste0(n.arm[t1], "(", n.arm.perc[t1], "%)"), + Arm_0 = treatments[t0], + N_Arm_0 = paste0(n.arm[t0], "(", n.arm.perc[t0], "%)"), + Eff_measure = eff_measure, + Estimate = calculate_f(theta[t1], theta[t0], eff_measure) + ) + outtmp[["S.E."]] <- sqrt(f_vcov[cbind(t1, t0)]) + outtmp[["95%CI_lower"]] <- outtmp$Estimate + qnorm(0.025) * outtmp[["S.E."]] + outtmp[["95%CI_upper"]] <- outtmp$Estimate + qnorm(0.975) * outtmp[["S.E."]] + + if (eff_measure == "diff") { + outtmp["Estimate"] <- round(outtmp["Estimate"], digits) |> format(nsmall = digits) + outtmp["S.E."] <- round(outtmp["S.E."], digits) |> format(nsmall = digits) + outtmp["95%CI_lower"] <- round(outtmp["95%CI_lower"], digits) |> format(nsmall = digits) + outtmp["95%CI_upper"] <- round(outtmp["95%CI_upper"], digits) |> format(nsmall = digits) + } else { + outtmp["S.E."] <- ((exp(outtmp["S.E."]) - 1) * exp(2 * outtmp["Estimate"] + outtmp["S.E."]^2)) |> + round(digits) |> + format(nsmall = digits) + outtmp["Estimate"] <- exp(outtmp["Estimate"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_lower"] <- exp(outtmp["95%CI_lower"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_upper"] <- exp(outtmp["95%CI_upper"]) |> + round(digits) |> + format(nsmall = digits) + } + outtmp +} + +report_fvcov <- function(result, digits = 3){ + report( + fit = result$fit, + trt_var_name = result$trt_var_name, + theta = result$theta, + f_vcov = result$f_vcov, + digits = digits + ) + +} + +# example: report +report( + fit = fit, + trt_var_name = "trt", + theta = theta_t, + f_vcov = n_fVf, + digits = 3 +) + + diff --git a/design/package_structure/design_xgc.Rmd b/design/package_structure/design_xgc.Rmd new file mode 100644 index 0000000..34fbb52 --- /dev/null +++ b/design/package_structure/design_xgc.Rmd @@ -0,0 +1,232 @@ +--- +title: "Test on RobinCar2" +author: Gregory Chen +--- + +## Create dummy data + +Before we start our model, first we need to prepare our data. + +In the following part we create a simple simulation data. +z1 comes from a bernouli distribution of p = 0.5. +z2 comes from a multinomial distribution of p = (0.4, 0.3, 0.3) +cov1 comes from a normal distribution. +trt is based on stratified block randomization, block size is 6, stratified by z1 and z2. +y is the response, the formula is + +exp(y) = 0.2 + 0.5 * (trt = "trt1") + 0.7 * (trt = "trt2") + 1.5 * (z1 = "b") - 1 * (z2 = y) - 0.5 * (z2 = z) + 0.2 * cov + +```{r} +n <- 200 +block_rand <- function(block = c(0, 0, 1, 1)) { + function(n) { + r <- lapply(seq_len(ceiling(n / length(block))), function(i) { + sample(block) + }) + unlist(r)[1:n] + } +} +by_strata <- function(f, strata) { + ret <- rep(NA, length(strata)) + for (x in split(seq_len(length(strata)), strata)) { + ret[x] <- f(length(x)) + } + return(ret) +} +sim_data <- function(n) { + cov1 <- rnorm(n) + z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) + z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) + permute_block <- c(0, 0, 1, 1) + trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) + trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) + df <- data.frame(trt, z1, z2, cov1) + x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) + coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) + theta <- x_mat %*% coef + y <- plogis(theta) + y_bin <- as.integer(runif(n) < y) + + df$y <- y_bin + df +} + +d <- sim_data(500) +``` + +## Construct the model + +In this part, we will create the model based on the data provided. +Let's begin with binary outcome: + +```{r} +fit <- glm(y ~ trt:z1 + trt, family = binomial(), data = d) +``` + +## Obtain the predictions + +```{r} +predict_cf <- function(fit, trt, data, ...) { + UseMethod("predict_cf") +} + +predict_cf.glm <- function(fit, trt, data, ...) { + lvls <- fit$xlevels[[trt]] + # data has to be a data.frame + df <- rbind(fit$model, fit$model, fit$model) + df[[trt]] <- factor(rep(lvls, each = nrow(fit$model)), levels = lvls) + preds <- predict(fit, type = "response", newdata = df) + matrix(preds, ncol = length(lvls), dimnames = list(row.names(fit$model), lvls)) +} +``` + +with this `predict_cf` function we already obtain the counterfactual response + +```{r} +pred <- predict_cf(fit, "trt") +pred +``` + +## Obtain the biasness + +```{r} +bias <- function(fit, trt, strat, data) { + UseMethod("bias") +} +bias.glm <- function(fit, trt, strat, data) { + trt_var <- fit$model[[trt]] + if (length(strat) != 0) { + strat_var <- fit$data[, strat] + } else { + strat_var <- rep(0L, length(fit$y)) + } + residuals <- fit$y - fit$fitted.values + d <- matrix(NA_real_, nrow = length(fit$y), ncol = length(fit$xlevels[[trt]])) + id_strat <- split(seq_len(length(residuals)), strat_var) + for (i in id_strat) { + df <- vapply(split(residuals[i], trt_var[i]), function(xx) mean(xx), FUN.VALUE = 0) + d[i, ] <- matrix(df, nrow = length(i), ncol = length(df), byrow = TRUE) + } + d +} +``` + +## Ensure the unbiasness + +to ensure the unbiasness, the predict_cf need additional argument + +```{r} +predict_cf.glm <- function(fit, trt, data, ensure_unbias = TRUE, strata = NULL, ...) { + lvls <- fit$xlevels[[trt]] + # data has to be a data.frame + df <- rbind(fit$model, fit$model, fit$model) + df[[trt]] <- factor(rep(lvls, each = nrow(fit$model)), levels = lvls) + preds <- predict(fit, type = "response", newdata = df) + ret <- matrix(preds, ncol = length(lvls), dimnames = list(row.names(fit$model), lvls)) + if (ensure_unbias) { + ret <- ret - bias(fit, trt, strata) + } + ret +} +``` + +```{r} +pred <- predict_cf(fit, "trt") +``` + + +## Robust variance + +```{r} +adjust_pi <- function(pi_t) { + diag(pi_t) - pi_t %*% t(pi_t) +} + +get_erb <- function(resi, strata, trt, pit, randomization) { + if (length(strata) == 0) { + return(0) + } + if (randomization %in% c("simple", "pocock-simon")) { + return(0) + } + # Calculate Omega Z under simple + omegaz_sr <- adjust_pi(pit) + idx <- split(seq_len(length(resi)), cbind(trt, strata)) + resi_per_strata <- vapply(idx, function(ii) mean(resi[ii]), FUN.VALUE = 0) + # Calculate strata levels and proportions for + # the outer expectation + + strata_props <- vapply(idx, length, FUN.VALUE = 0L) + strata_props <- strata_props / sum(strata_props) + # Estimate R(B) by first getting the conditional expectation + # vector for a particular strata (vector contains + # all treatment groups), then dividing by the pi_t + + rb_z <- resi_per_strata / as.numeric(pit) + # Compute the R(B)[Omega_{SR} - Omega_{Z_i}]R(B) | Z_i + # for each Z_i + strata_levels <- length(resi_per_strata) + n_trt <- length(pit) + ind <- matrix(seq_len(strata_levels), byrow = TRUE, ncol = n_trt) + rb_z_sum <- lapply( + seq_len(nrow(ind)), + function(x) rb_z[ind[x, ]] %*% t(rb_z[ind[x, ]]) * sum(strata_props[ind[x, ]]) + ) + rb_z_sum <- Reduce(`+`, rb_z_sum) + rb_z_sum * omegaz_sr +} + +vcov_robin <- function(fit, ...) { + UseMethod("vcov_robin") +} + +are <- function(objs, cls) { + all(vapply(objs, is, class2 = cls, FUN.VALUE = TRUE)) +} + +vcov_robin.glm <- function(fit, trt, strat = NULL, preds = predict_cf(fit, trt, strat), sr_decompose = TRUE, randomization = "simple", ...) { + raw_data <- fit$data + resi <- residuals(fit, type = "response") + est <- colMeans(preds) + var_preds <- var(preds) + pit <- as.numeric(table(raw_data[[trt]]) / nrow(raw_data)) + trt_lvls <- fit$xlevels[[trt]] + + idx <- split(seq_len(nrow(raw_data)), raw_data[[trt]]) + cov_ymu <- vapply(idx, function(is) stats::cov(fit$y[is], preds[is, ]), FUN.VALUE = rep(0, ncol(preds))) + + vcov_sr <- vapply(idx, function(is) stats::var(fit$y[is]), FUN.VALUE = 0) / pit + if (sr_decompose) { + vcov_sr <- vcov_sr + (diag(var_preds) - 2 * diag(cov_ymu)) / pit + } + v <- diag(vcov_sr) + cov_ymu + t(cov_ymu) - var_preds + if (!randomization %in% c("simple", "pocock-simon") && length(strat) > 0) { + v <- v - get_erb(resi, raw_data[strat], raw_data[[trt]], pit, randomization) + } + ret <- v / length(resi) + dimnames(ret) <- list(trt_lvls, trt_lvls) + return(ret) +} +vcov_robin(fit, "trt", "z1") +``` + + +## robincar_glm wrapper + +```{r} +robincar_glm2 <- function(data, formula, trt, strata, car_scheme, ...) { + fit <- glm(formula = formula, data = data, ...) + pred <- predict_cf(fit, trt, ensure_unbias = TRUE) + var <- vcov_robin(fit, trt = trt, strat = strata, randomization = car_scheme, preds = pred) + list( + pred = colMeans(pred), + var = var + ) +} +robincar_glm(d, formula = y ~ z1:trt + trt, "trt", "z1", "coin") +``` + +potentially, there can be other arguments just like `RobinCar::robincar_glm`. +The formula can be inferred from the choices of "ANOVA", "ANCOVA", "ANHECOVA", etc. if not provided. + +All these functions can be exported for direct usage, and also the wrapper is provided. diff --git a/design/package_structure/test_all_draft_functions.R b/design/package_structure/test_all_draft_functions.R new file mode 100644 index 0000000..3fa1760 --- /dev/null +++ b/design/package_structure/test_all_draft_functions.R @@ -0,0 +1,42 @@ +n <- 200 +block_rand <- function(block = c(0, 0, 1, 1)) { + function(n) { + r <- lapply(seq_len(ceiling(n / length(block))), function(i) { + sample(block) + }) + unlist(r)[1:n] + } +} +by_strata <- function(f, strata) { + ret <- rep(NA, length(strata)) + for (x in split(seq_len(length(strata)), strata)) { + ret[x] <- f(length(x)) + } + return(ret) +} +sim_data <- function(n) { + cov1 <- rnorm(n) + z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) + z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) + permute_block <- c(0, 0, 1, 1) + trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) + trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) + df <- data.frame(trt, z1, z2, cov1) + x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) + coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) + theta <- x_mat %*% coef + y <- plogis(theta) + y_bin <- as.integer(runif(n) < y) + + df$y <- y_bin + df +} + +d <- sim_data(500) +fit <- glm(y ~ trt:z1 + trt, family = binomial(link="logit"), data = d) + +## start of the pkg +source(file.path(getwd(),"design","package_structure","all_draft_functions.R")) +pred_counterfact <- get_countfact_pred(fit,"trt") +fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact) +report_fvcov(fit.fvcov,3) From 5543640be8492b0e4a46c3f04f9c34bc28fcedc6 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Mon, 6 May 2024 06:42:41 +0000 Subject: [PATCH 02/12] [skip style] [skip vbump] Restyle files --- .../package_structure/all_draft_functions.R | 31 ++++++++++--------- design/package_structure/design_xgc.R | 29 +++++++++-------- .../test_all_draft_functions.R | 8 ++--- 3 files changed, 34 insertions(+), 34 deletions(-) diff --git a/design/package_structure/all_draft_functions.R b/design/package_structure/all_draft_functions.R index cd20e07..c900512 100644 --- a/design/package_structure/all_draft_functions.R +++ b/design/package_structure/all_draft_functions.R @@ -1,4 +1,3 @@ - get_eff_measure <- function(fit) { if (identical(class(fit), "lm")) { out_eff <- "diff" @@ -23,7 +22,7 @@ get_eff_measure <- function(fit) { ### g-computation and its VCOV ------------------------------------------------- # func: prepare data -get_countfact_pred <- function(fit, trt_var_name){ +get_countfact_pred <- function(fit, trt_var_name) { treatments <- fit$xlevels[[trt_var_name]] use_dat <- fit$model @@ -119,7 +118,7 @@ calculate_fVf <- function(t, s, theta, V, n, eff_measure) { vectorized_calculate_fVf <- Vectorize(calculate_fVf, c("t", "s")) # func: wrapper function given fit and pred_counterfact -calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ +calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred) { treatments <- fit$xlevels[[trt_var_name]] theta_t <- sapply(counterfact_pred, function(dat_t) unique(dat_t$theta_t)) @@ -143,23 +142,25 @@ calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ ) rownames(n_fVf) <- colnames(n_fVf) <- treatments - out <- list( theta = theta_t, - V = V, - f_vcov = n_fVf, - n = nrow(fit$model), - eff_measure = get_eff_measure(fit), - trt_var_name = trt_var_name, - fit = fit ) - class(out) <- c(class(out),"vcov_rc") + out <- list( + theta = theta_t, + V = V, + f_vcov = n_fVf, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit), + trt_var_name = trt_var_name, + fit = fit + ) + class(out) <- c(class(out), "vcov_rc") out } # func: reporting function calculate_f <- function(theta_t, theta_s, eff_measure) { switch(eff_measure, - "diff" = theta_t - theta_s, - "risk ratio" = log(theta_t / theta_s), - "odds ratio" = log(odds(theta_t) / odds(theta_s)) + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) ) } @@ -210,7 +211,7 @@ report <- function(fit, trt_var_name, theta, f_vcov, digits = 3) { outtmp } -report_fvcov <- function(result, digits = 3){ +report_fvcov <- function(result, digits = 3) { report( fit = result$fit, trt_var_name = result$trt_var_name, diff --git a/design/package_structure/design_xgc.R b/design/package_structure/design_xgc.R index 1463fe1..f70c6e1 100644 --- a/design/package_structure/design_xgc.R +++ b/design/package_structure/design_xgc.R @@ -119,7 +119,7 @@ get_eff_measure <- function(fit) { ### g-computation and its VCOV ------------------------------------------------- # func: prepare data -get_countfact_pred <- function(fit, trt_var_name){ +get_countfact_pred <- function(fit, trt_var_name) { treatments <- fit$xlevels[[trt_var_name]] use_dat <- fit$model @@ -145,7 +145,7 @@ get_countfact_pred <- function(fit, trt_var_name){ } # Step 1: get g-computation theta_t -pred_counterfact <- get_countfact_pred(fit,"trt") +pred_counterfact <- get_countfact_pred(fit, "trt") theta_t <- sapply(pred_counterfact, function(dat_t) unique(dat_t$theta_t)) # func: calculate V @@ -241,7 +241,7 @@ n_fVf <- outer( rownames(n_fVf) <- colnames(n_fVf) <- treatments # func: wrapper function given fit and pred_counterfact -calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ +calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred) { treatments <- fit$xlevels[[trt_var_name]] theta_t <- sapply(counterfact_pred, function(dat_t) unique(dat_t$theta_t)) @@ -265,14 +265,16 @@ calculate_f_vcov <- function(fit, trt_var_name, counterfact_pred){ ) rownames(n_fVf) <- colnames(n_fVf) <- treatments - out <- list( theta = theta_t, - V = V, - f_vcov = n_fVf, - n = nrow(fit$model), - eff_measure = get_eff_measure(fit), - trt_var_name = trt_var_name, - fit = fit ) - class(out) <- c(class(out),"vcov_rc") + out <- list( + theta = theta_t, + V = V, + f_vcov = n_fVf, + n = nrow(fit$model), + eff_measure = get_eff_measure(fit), + trt_var_name = trt_var_name, + fit = fit + ) + class(out) <- c(class(out), "vcov_rc") out } @@ -332,7 +334,7 @@ report <- function(fit, trt_var_name, theta, f_vcov, digits = 3) { outtmp } -report_fvcov <- function(result, digits = 3){ +report_fvcov <- function(result, digits = 3) { report( fit = result$fit, trt_var_name = result$trt_var_name, @@ -340,7 +342,6 @@ report_fvcov <- function(result, digits = 3){ f_vcov = result$f_vcov, digits = digits ) - } # example: report @@ -351,5 +352,3 @@ report( f_vcov = n_fVf, digits = 3 ) - - diff --git a/design/package_structure/test_all_draft_functions.R b/design/package_structure/test_all_draft_functions.R index 3fa1760..73edf68 100644 --- a/design/package_structure/test_all_draft_functions.R +++ b/design/package_structure/test_all_draft_functions.R @@ -33,10 +33,10 @@ sim_data <- function(n) { } d <- sim_data(500) -fit <- glm(y ~ trt:z1 + trt, family = binomial(link="logit"), data = d) +fit <- glm(y ~ trt:z1 + trt, family = binomial(link = "logit"), data = d) ## start of the pkg -source(file.path(getwd(),"design","package_structure","all_draft_functions.R")) -pred_counterfact <- get_countfact_pred(fit,"trt") +source(file.path(getwd(), "design", "package_structure", "all_draft_functions.R")) +pred_counterfact <- get_countfact_pred(fit, "trt") fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact) -report_fvcov(fit.fvcov,3) +report_fvcov(fit.fvcov, 3) From 8c157fe1a305725da586e4a3fbf9532605782263 Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 19 May 2024 13:38:42 +0000 Subject: [PATCH 03/12] clean up design folder for merging --- .../{design_lml.Rmd => design.Rmd} | 0 .../{design_xgc.R => design_note_xgc_v1.R} | 2 +- design/package_structure/design_xgc.Rmd | 232 ------------------ ...ft_functions.R => developer_note_xgc_v1.R} | 52 +++- .../test_all_draft_functions.R | 42 ---- 5 files changed, 51 insertions(+), 277 deletions(-) rename design/package_structure/{design_lml.Rmd => design.Rmd} (100%) rename design/package_structure/{design_xgc.R => design_note_xgc_v1.R} (99%) delete mode 100644 design/package_structure/design_xgc.Rmd rename design/package_structure/{all_draft_functions.R => developer_note_xgc_v1.R} (82%) delete mode 100644 design/package_structure/test_all_draft_functions.R diff --git a/design/package_structure/design_lml.Rmd b/design/package_structure/design.Rmd similarity index 100% rename from design/package_structure/design_lml.Rmd rename to design/package_structure/design.Rmd diff --git a/design/package_structure/design_xgc.R b/design/package_structure/design_note_xgc_v1.R similarity index 99% rename from design/package_structure/design_xgc.R rename to design/package_structure/design_note_xgc_v1.R index f70c6e1..5e5d998 100644 --- a/design/package_structure/design_xgc.R +++ b/design/package_structure/design_note_xgc_v1.R @@ -1,4 +1,4 @@ -# Design note +# Design note for developing v1 # Author: X. Gregory Chen # [Refenrece] diff --git a/design/package_structure/design_xgc.Rmd b/design/package_structure/design_xgc.Rmd deleted file mode 100644 index 34fbb52..0000000 --- a/design/package_structure/design_xgc.Rmd +++ /dev/null @@ -1,232 +0,0 @@ ---- -title: "Test on RobinCar2" -author: Gregory Chen ---- - -## Create dummy data - -Before we start our model, first we need to prepare our data. - -In the following part we create a simple simulation data. -z1 comes from a bernouli distribution of p = 0.5. -z2 comes from a multinomial distribution of p = (0.4, 0.3, 0.3) -cov1 comes from a normal distribution. -trt is based on stratified block randomization, block size is 6, stratified by z1 and z2. -y is the response, the formula is - -exp(y) = 0.2 + 0.5 * (trt = "trt1") + 0.7 * (trt = "trt2") + 1.5 * (z1 = "b") - 1 * (z2 = y) - 0.5 * (z2 = z) + 0.2 * cov - -```{r} -n <- 200 -block_rand <- function(block = c(0, 0, 1, 1)) { - function(n) { - r <- lapply(seq_len(ceiling(n / length(block))), function(i) { - sample(block) - }) - unlist(r)[1:n] - } -} -by_strata <- function(f, strata) { - ret <- rep(NA, length(strata)) - for (x in split(seq_len(length(strata)), strata)) { - ret[x] <- f(length(x)) - } - return(ret) -} -sim_data <- function(n) { - cov1 <- rnorm(n) - z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) - z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) - permute_block <- c(0, 0, 1, 1) - trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) - trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) - df <- data.frame(trt, z1, z2, cov1) - x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) - coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) - theta <- x_mat %*% coef - y <- plogis(theta) - y_bin <- as.integer(runif(n) < y) - - df$y <- y_bin - df -} - -d <- sim_data(500) -``` - -## Construct the model - -In this part, we will create the model based on the data provided. -Let's begin with binary outcome: - -```{r} -fit <- glm(y ~ trt:z1 + trt, family = binomial(), data = d) -``` - -## Obtain the predictions - -```{r} -predict_cf <- function(fit, trt, data, ...) { - UseMethod("predict_cf") -} - -predict_cf.glm <- function(fit, trt, data, ...) { - lvls <- fit$xlevels[[trt]] - # data has to be a data.frame - df <- rbind(fit$model, fit$model, fit$model) - df[[trt]] <- factor(rep(lvls, each = nrow(fit$model)), levels = lvls) - preds <- predict(fit, type = "response", newdata = df) - matrix(preds, ncol = length(lvls), dimnames = list(row.names(fit$model), lvls)) -} -``` - -with this `predict_cf` function we already obtain the counterfactual response - -```{r} -pred <- predict_cf(fit, "trt") -pred -``` - -## Obtain the biasness - -```{r} -bias <- function(fit, trt, strat, data) { - UseMethod("bias") -} -bias.glm <- function(fit, trt, strat, data) { - trt_var <- fit$model[[trt]] - if (length(strat) != 0) { - strat_var <- fit$data[, strat] - } else { - strat_var <- rep(0L, length(fit$y)) - } - residuals <- fit$y - fit$fitted.values - d <- matrix(NA_real_, nrow = length(fit$y), ncol = length(fit$xlevels[[trt]])) - id_strat <- split(seq_len(length(residuals)), strat_var) - for (i in id_strat) { - df <- vapply(split(residuals[i], trt_var[i]), function(xx) mean(xx), FUN.VALUE = 0) - d[i, ] <- matrix(df, nrow = length(i), ncol = length(df), byrow = TRUE) - } - d -} -``` - -## Ensure the unbiasness - -to ensure the unbiasness, the predict_cf need additional argument - -```{r} -predict_cf.glm <- function(fit, trt, data, ensure_unbias = TRUE, strata = NULL, ...) { - lvls <- fit$xlevels[[trt]] - # data has to be a data.frame - df <- rbind(fit$model, fit$model, fit$model) - df[[trt]] <- factor(rep(lvls, each = nrow(fit$model)), levels = lvls) - preds <- predict(fit, type = "response", newdata = df) - ret <- matrix(preds, ncol = length(lvls), dimnames = list(row.names(fit$model), lvls)) - if (ensure_unbias) { - ret <- ret - bias(fit, trt, strata) - } - ret -} -``` - -```{r} -pred <- predict_cf(fit, "trt") -``` - - -## Robust variance - -```{r} -adjust_pi <- function(pi_t) { - diag(pi_t) - pi_t %*% t(pi_t) -} - -get_erb <- function(resi, strata, trt, pit, randomization) { - if (length(strata) == 0) { - return(0) - } - if (randomization %in% c("simple", "pocock-simon")) { - return(0) - } - # Calculate Omega Z under simple - omegaz_sr <- adjust_pi(pit) - idx <- split(seq_len(length(resi)), cbind(trt, strata)) - resi_per_strata <- vapply(idx, function(ii) mean(resi[ii]), FUN.VALUE = 0) - # Calculate strata levels and proportions for - # the outer expectation - - strata_props <- vapply(idx, length, FUN.VALUE = 0L) - strata_props <- strata_props / sum(strata_props) - # Estimate R(B) by first getting the conditional expectation - # vector for a particular strata (vector contains - # all treatment groups), then dividing by the pi_t - - rb_z <- resi_per_strata / as.numeric(pit) - # Compute the R(B)[Omega_{SR} - Omega_{Z_i}]R(B) | Z_i - # for each Z_i - strata_levels <- length(resi_per_strata) - n_trt <- length(pit) - ind <- matrix(seq_len(strata_levels), byrow = TRUE, ncol = n_trt) - rb_z_sum <- lapply( - seq_len(nrow(ind)), - function(x) rb_z[ind[x, ]] %*% t(rb_z[ind[x, ]]) * sum(strata_props[ind[x, ]]) - ) - rb_z_sum <- Reduce(`+`, rb_z_sum) - rb_z_sum * omegaz_sr -} - -vcov_robin <- function(fit, ...) { - UseMethod("vcov_robin") -} - -are <- function(objs, cls) { - all(vapply(objs, is, class2 = cls, FUN.VALUE = TRUE)) -} - -vcov_robin.glm <- function(fit, trt, strat = NULL, preds = predict_cf(fit, trt, strat), sr_decompose = TRUE, randomization = "simple", ...) { - raw_data <- fit$data - resi <- residuals(fit, type = "response") - est <- colMeans(preds) - var_preds <- var(preds) - pit <- as.numeric(table(raw_data[[trt]]) / nrow(raw_data)) - trt_lvls <- fit$xlevels[[trt]] - - idx <- split(seq_len(nrow(raw_data)), raw_data[[trt]]) - cov_ymu <- vapply(idx, function(is) stats::cov(fit$y[is], preds[is, ]), FUN.VALUE = rep(0, ncol(preds))) - - vcov_sr <- vapply(idx, function(is) stats::var(fit$y[is]), FUN.VALUE = 0) / pit - if (sr_decompose) { - vcov_sr <- vcov_sr + (diag(var_preds) - 2 * diag(cov_ymu)) / pit - } - v <- diag(vcov_sr) + cov_ymu + t(cov_ymu) - var_preds - if (!randomization %in% c("simple", "pocock-simon") && length(strat) > 0) { - v <- v - get_erb(resi, raw_data[strat], raw_data[[trt]], pit, randomization) - } - ret <- v / length(resi) - dimnames(ret) <- list(trt_lvls, trt_lvls) - return(ret) -} -vcov_robin(fit, "trt", "z1") -``` - - -## robincar_glm wrapper - -```{r} -robincar_glm2 <- function(data, formula, trt, strata, car_scheme, ...) { - fit <- glm(formula = formula, data = data, ...) - pred <- predict_cf(fit, trt, ensure_unbias = TRUE) - var <- vcov_robin(fit, trt = trt, strat = strata, randomization = car_scheme, preds = pred) - list( - pred = colMeans(pred), - var = var - ) -} -robincar_glm(d, formula = y ~ z1:trt + trt, "trt", "z1", "coin") -``` - -potentially, there can be other arguments just like `RobinCar::robincar_glm`. -The formula can be inferred from the choices of "ANOVA", "ANCOVA", "ANHECOVA", etc. if not provided. - -All these functions can be exported for direct usage, and also the wrapper is provided. diff --git a/design/package_structure/all_draft_functions.R b/design/package_structure/developer_note_xgc_v1.R similarity index 82% rename from design/package_structure/all_draft_functions.R rename to design/package_structure/developer_note_xgc_v1.R index c900512..0bdbda5 100644 --- a/design/package_structure/all_draft_functions.R +++ b/design/package_structure/developer_note_xgc_v1.R @@ -1,3 +1,8 @@ +# Developer note for v1 +# Author: X. Gregory Chen + +#### all draft functions #### + get_eff_measure <- function(fit) { if (identical(class(fit), "lm")) { out_eff <- "diff" @@ -19,8 +24,6 @@ get_eff_measure <- function(fit) { out_eff } -### g-computation and its VCOV ------------------------------------------------- - # func: prepare data get_countfact_pred <- function(fit, trt_var_name) { treatments <- fit$xlevels[[trt_var_name]] @@ -220,3 +223,48 @@ report_fvcov <- function(result, digits = 3) { digits = digits ) } + +#### Example run using draft functions above #### +n <- 200 +block_rand <- function(block = c(0, 0, 1, 1)) { + function(n) { + r <- lapply(seq_len(ceiling(n / length(block))), function(i) { + sample(block) + }) + unlist(r)[1:n] + } +} +by_strata <- function(f, strata) { + ret <- rep(NA, length(strata)) + for (x in split(seq_len(length(strata)), strata)) { + ret[x] <- f(length(x)) + } + return(ret) +} +sim_data <- function(n) { + cov1 <- rnorm(n) + z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) + z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) + permute_block <- c(0, 0, 1, 1) + trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) + trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) + df <- data.frame(trt, z1, z2, cov1) + x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) + coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) + theta <- x_mat %*% coef + y <- plogis(theta) + y_bin <- as.integer(runif(n) < y) + + df$y <- y_bin + df +} + +d <- sim_data(500) +fit <- glm(y ~ trt:z1 + trt, family = binomial(link = "logit"), data = d) + +## start of the pkg +source(file.path(getwd(), "design", "package_structure", "all_draft_functions.R")) +pred_counterfact <- get_countfact_pred(fit, "trt") +fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact) +report_fvcov(fit.fvcov, 3) + diff --git a/design/package_structure/test_all_draft_functions.R b/design/package_structure/test_all_draft_functions.R deleted file mode 100644 index 73edf68..0000000 --- a/design/package_structure/test_all_draft_functions.R +++ /dev/null @@ -1,42 +0,0 @@ -n <- 200 -block_rand <- function(block = c(0, 0, 1, 1)) { - function(n) { - r <- lapply(seq_len(ceiling(n / length(block))), function(i) { - sample(block) - }) - unlist(r)[1:n] - } -} -by_strata <- function(f, strata) { - ret <- rep(NA, length(strata)) - for (x in split(seq_len(length(strata)), strata)) { - ret[x] <- f(length(x)) - } - return(ret) -} -sim_data <- function(n) { - cov1 <- rnorm(n) - z1 <- sample(size = n, factor(c("a", "b")), c(0.5, 0.5), replace = TRUE) - z2 <- sample(size = n, factor(c("x", "y", "z")), c(0.4, 0.3, 0.3), replace = TRUE) - permute_block <- c(0, 0, 1, 1) - trt <- by_strata(block_rand(c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo")), interaction(z1, z2)) - trt <- factor(trt, levels = c("pbo", "trt1", "trt2")) - df <- data.frame(trt, z1, z2, cov1) - x_mat <- model.matrix(~ trt + z1 + z2 + cov1, data = df) - coef <- c(0.2, 0.5, 0.7, 1.5, -1, -0.5, 0.2) - theta <- x_mat %*% coef - y <- plogis(theta) - y_bin <- as.integer(runif(n) < y) - - df$y <- y_bin - df -} - -d <- sim_data(500) -fit <- glm(y ~ trt:z1 + trt, family = binomial(link = "logit"), data = d) - -## start of the pkg -source(file.path(getwd(), "design", "package_structure", "all_draft_functions.R")) -pred_counterfact <- get_countfact_pred(fit, "trt") -fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact) -report_fvcov(fit.fvcov, 3) From 4a582533bdfd339059187fc82bf8502872208666 Mon Sep 17 00:00:00 2001 From: Gregory Chen <65814078+hoppanda@users.noreply.github.com> Date: Sun, 19 May 2024 15:55:00 +0200 Subject: [PATCH 04/12] Update branch 9 with main (#23) * update data and fix design (#17) * update data and fix design * update treatment to factor * fix styler and lint * 7 counter factual prediction (#14) * add function prototype * add basic tests on functionality * add readme * add lm support * update tests * update prediction * simplify * fix typo * fix check issues * update vignettes skeleton * use dummy data in test * update print of cf predictions * fix checks * fix check issues * update bias args * [skip style] [skip vbump] Restyle files * add treatment to counter factual prediction * fix check * 18 sandwich variance (#19) * add vcov * fix check issues * correct issue --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> * 20 add another developers design note (#22) * add in all draft functions * [skip style] [skip vbump] Restyle files * clean up design folder for merging * [skip style] [skip vbump] Restyle files * Update .gitignore --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Liming <36079400+clarkliming@users.noreply.github.com> --------- Co-authored-by: Liming <36079400+clarkliming@users.noreply.github.com> Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> --- .Rbuildignore | 5 +- DESCRIPTION | 11 +- NAMESPACE | 15 +- NEWS.md | 11 +- R/RobinCar2-package.R | 11 + R/bias.R | 24 + R/data.R | 10 +- R/hello.R | 15 - R/predict_couterfactual.R | 77 +++ R/prediction_cf.R | 13 + R/utils.R | 37 ++ R/variance.R | 12 + README.Rmd | 46 ++ data-raw/create_dummy.R | 12 +- data/dummy_data.rda | Bin 8072 -> 11593 bytes design/package_structure/design.Rmd | 12 +- inst/WORDLIST | 2 + man/RobinCar2-package.Rd | 22 + man/bias.Rd | 21 + man/dummy_data.Rd | 10 +- man/h_get_vars.Rd | 20 + man/hello.Rd | 20 - man/predict_counterfactual.Rd | 24 + man/prediction_cf_methods.Rd | 21 + tests/testthat/_snaps/bias.md | 608 ++++++++++++++++++ .../testthat/_snaps/predict_counterfactual.md | 33 + tests/testthat/helper-fit.R | 8 + tests/testthat/test-bias.R | 5 + tests/testthat/test-hello.R | 5 - tests/testthat/test-predict_counterfactual.R | 11 + tests/testthat/test-utils.R | 34 + vignettes/hello.Rmd | 12 - vignettes/intro.Rmd | 20 + 33 files changed, 1098 insertions(+), 89 deletions(-) create mode 100644 R/RobinCar2-package.R create mode 100644 R/bias.R delete mode 100644 R/hello.R create mode 100644 R/predict_couterfactual.R create mode 100644 R/prediction_cf.R create mode 100644 R/utils.R create mode 100644 R/variance.R create mode 100644 README.Rmd create mode 100644 man/RobinCar2-package.Rd create mode 100644 man/bias.Rd create mode 100644 man/h_get_vars.Rd delete mode 100644 man/hello.Rd create mode 100644 man/predict_counterfactual.Rd create mode 100644 man/prediction_cf_methods.Rd create mode 100644 tests/testthat/_snaps/bias.md create mode 100644 tests/testthat/_snaps/predict_counterfactual.md create mode 100644 tests/testthat/helper-fit.R create mode 100644 tests/testthat/test-bias.R delete mode 100644 tests/testthat/test-hello.R create mode 100644 tests/testthat/test-predict_counterfactual.R create mode 100644 tests/testthat/test-utils.R delete mode 100644 vignettes/hello.Rmd create mode 100644 vignettes/intro.Rmd diff --git a/.Rbuildignore b/.Rbuildignore index 8cd5db0..5df0095 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -17,6 +17,5 @@ coverage.* init.sh workflows.md images -^RobinCar2\.Rcheck$ -^RobinCar2.*\.tar\.gz$ -^RobinCar2.*\.tgz$ +design +data-raw diff --git a/DESCRIPTION b/DESCRIPTION index 63327a5..11a5aa1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Type: Package Package: RobinCar2 Title: R Package Template Version: 0.1.0.9000 -Date: 2023-06-20 +Date: 2024-05-08 Authors@R: person("openpharma", , , "openpharma@example.com", role = c("aut", "cre")) Description: R package template with GitHub Actions workflows included. @@ -12,13 +12,16 @@ BugReports: https://github.com/openpharma/RobinCar2/issues Depends: R (>= 3.6) Imports: - stringr + checkmate, + sandwich, + stats Suggests: knitr, - testthat (>= 2.0) + rmarkdown, + testthat (>= 3.0) VignetteBuilder: knitr -biocViews: +Config/testthat/edition: 3 Encoding: UTF-8 Language: en-US LazyData: true diff --git a/NAMESPACE b/NAMESPACE index 01843dc..55f3734 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,3 +1,16 @@ # Generated by roxygen2: do not edit by hand -export(hello) +S3method(predict_counterfactual,glm) +S3method(predict_counterfactual,lm) +S3method(print,prediction_cf) +S3method(vcovHC,prediction_cf) +export(bias) +export(predict_counterfactual) +import(checkmate) +importFrom(sandwich,vcovHC) +importFrom(stats,coefficients) +importFrom(stats,fitted) +importFrom(stats,model.matrix) +importFrom(stats,model.response) +importFrom(stats,predict) +importFrom(stats,residuals) diff --git a/NEWS.md b/NEWS.md index f770305..0cc6e5e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,12 +2,5 @@ ### New features -* Add an initializer script. - -### Enhancements - -* Documentation on how to use the initialize a package. - -### Bug fixes - -* None. \ No newline at end of file +* Add unbiased counter-factual prediction. +* Add robust sandwich variance for the marginal treatment effect. diff --git a/R/RobinCar2-package.R b/R/RobinCar2-package.R new file mode 100644 index 0000000..c47f10d --- /dev/null +++ b/R/RobinCar2-package.R @@ -0,0 +1,11 @@ +#' `RobinCar2` Package +#' +#' `RobinCar2` implements unbiased prediction and robust inference of variance of a fit in R. +#' +#' @aliases RobinCar2-package +"_PACKAGE" + +#' @import checkmate +#' @importFrom stats predict residuals fitted model.response model.matrix coefficients +#' @importFrom sandwich vcovHC +NULL diff --git a/R/bias.R b/R/bias.R new file mode 100644 index 0000000..f6b3c4c --- /dev/null +++ b/R/bias.R @@ -0,0 +1,24 @@ +#' Prediction Bias +#' @description Obtain prediction bias within each stratum. +#' @param residual (`numeric`) residuals. +#' @param treatment (`factor`) treatment. +#' @param group_idx (`character`) stratum index. +#' +#' @return Numeric matrix of bias in each stratum. +#' @export +bias <- function(residual, treatment, group_idx) { + assert_numeric(residual) + assert_factor(treatment, len = length(residual)) + assert_list(group_idx, types = "integer") + grp <- unlist(group_idx) + assert_integer(grp, lower = 1L, upper = max(grp)) + + trt_lvls <- levels(treatment) + + d <- matrix(NA_real_, nrow = length(grp), ncol = length(trt_lvls)) + for (i in group_idx) { + mval <- vapply(split(residual[i], treatment[i]), mean, FUN.VALUE = 0) + d[i, ] <- matrix(mval, nrow = length(i), ncol = length(mval), byrow = TRUE) + } + d +} diff --git a/R/data.R b/R/data.R index 9cc292f..a3aef96 100644 --- a/R/data.R +++ b/R/data.R @@ -2,13 +2,13 @@ #' #' This dataset contains the dummy trial data with permute block randomization. #' -#' @format A data frame with 400 rows and 3 columns: +#' @format A data frame with 600 rows and 7 columns: #' \describe{ #' \item{id}{The ID of the patients.} -#' \item{treatment}{The treatment assignment.} -#' \item{s1}{The first stratification variable.} -#' \item{s2}{The second stratification variable.} -#' \item{covar}{The covariate.} +#' \item{treatment}{The treatment assignment, "pbo", "trt1" and "trt2"} +#' \item{s1}{The first stratification variable, "a" and "b".} +#' \item{s2}{The second stratification variable, "c" and "d".} +#' \item{covar}{The covariate following normal distribution.} #' \item{y}{The continuous response.} #' \item{y_b}{The binary response.} #' } diff --git a/R/hello.R b/R/hello.R deleted file mode 100644 index ae4284d..0000000 --- a/R/hello.R +++ /dev/null @@ -1,15 +0,0 @@ -#' Personal greeting -#' -#' @description Greet a person and appropriately capitalize their name. -#' -#' @param name Your name (character string; e.g. "john doe"). -#' -#' @return A character string, capitalized to title case. -#' @export -#' -#' @examples -#' hello("james bond") -hello <- function(name = "your name") { - name <- stringr::str_to_title(name) - print(paste("Hello,", name)) -} diff --git a/R/predict_couterfactual.R b/R/predict_couterfactual.R new file mode 100644 index 0000000..445013d --- /dev/null +++ b/R/predict_couterfactual.R @@ -0,0 +1,77 @@ +#' Counterfactual Prediction +#' @description Obtain counterfactual prediction of a fit. +#' +#' @param fit fitted object. +#' @param treatment (`string` or `formula`) treatment variable in string, or a formula of form +#' treatment ~ strata(s). +#' @param data (`data.frame`) raw dataset. +#' @param unbiased (`flag`) indicator of whether to remove potential bias of the prediction. +#' +#' @return Numeric matrix of counter factual prediction. +#' +#' @export +predict_counterfactual <- function(fit, treatment, data, unbiased) { + UseMethod("predict_counterfactual") +} + +#' @export +predict_counterfactual.lm <- function(fit, treatment, data, unbiased = TRUE) { + treatment <- h_get_vars(treatment) + assert_data_frame(data) + assert_subset(unlist(treatment), colnames(data)) + formula <- formula(fit) + assert_subset(treatment$treatment, all.vars(formula[[3]])) + assert( + test_character(data[[treatment$treatment]]), + test_factor(data[[treatment$treatment]]) + ) + data[[treatment$treatment]] <- as.factor(data[[treatment$treatment]]) + assert_flag(unbiased) + + trt_lvls <- levels(data[[treatment$treatment]]) + n_lvls <- length(trt_lvls) + + df <- lapply( + data, + function(i) { + rep(i, times = n_lvls) + } + ) + + df[[treatment$treatment]] <- rep(trt_lvls, each = nrow(data)) + + mm <- model.matrix(fit, data = df) + pred_linear <- mm %*% coefficients(fit) + preds <- family(fit)$linkinv(pred_linear) + + ret <- matrix(preds, ncol = n_lvls, dimnames = list(row.names(data), trt_lvls)) + y <- model.response(fit$model) + residual <- y - fitted(fit) + + strata <- data[treatment$strata] + if (ncol(strata) == 0L) { + strata <- integer(nrow(strata)) + } + group_idx <- split(seq_len(nrow(data)), strata) + + if (unbiased) { + ret <- ret - bias(residual, data[[treatment$treatment]], group_idx) + } + structure( + .Data = colMeans(ret), + residual = residual, + predictions = ret, + predictions_linear = pred_linear, + response = y, + fit = fit, + model_matrix = mm, + treatment = data[[treatment$treatment]], + group_idx = group_idx, + class = "prediction_cf" + ) +} + +#' @export +predict_counterfactual.glm <- function(fit, treatment, data = fit$data, unbiased = TRUE) { + predict_counterfactual.lm(fit = fit, data = data, treatment = treatment, unbiased = unbiased) +} diff --git a/R/prediction_cf.R b/R/prediction_cf.R new file mode 100644 index 0000000..d917d78 --- /dev/null +++ b/R/prediction_cf.R @@ -0,0 +1,13 @@ +#' S3 Methods for `prediction_cf` +#' @param x (`prediction_cf`)\cr the obtained counter-factual prediction object. +#' @name prediction_cf_methods +NULL + +#' @describeIn prediction_cf_methods prints the prediction_cf object. +#' @exportS3Method +#' @keywords internal +print.prediction_cf <- function(x, ...) { + cat("counter-factual prediction\n\n") + print(x[seq_along(x)]) + cat("\n") +} diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..891c827 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,37 @@ +#' Extract Variable Names +#' +#' @param treatment (`string` or `formula`) string name of the treatment, or a formula. +#' +#' @details Extract the formula elements, including `treatment` and `strata`. +#' +#' @return A list of two elements, `treatmetn` and `strata`. +h_get_vars <- function(treatment) { + if (test_string(treatment)) { + ret <- list( + treatment = treatment, + strata = character(0) + ) + } else if (test_formula(treatment)) { + if (!identical(length(treatment), 3L)) { + stop("treatment formula must be of type treatment ~ strata") + } + if (!is.name(treatment[[2]])) { + stop("left hand side of the treatment formula should be a single name!") + } + treatvar <- as.character(treatment[[2]]) + strata <- setdiff(all.vars(treatment[[3]]), ".") + ret <- list( + treatment = treatvar, + strata = strata + ) + } else { + stop("treatment must be a length 1 character or a formula of form treatment ~ strata") + } + ret +} + +block_sum <- function(x, n) { + assert_matrix(x) + nr <- nrow(x) / n + matrix(colSums(matrix(x, nrow = n)), nrow = nr) +} diff --git a/R/variance.R b/R/variance.R new file mode 100644 index 0000000..be5f633 --- /dev/null +++ b/R/variance.R @@ -0,0 +1,12 @@ +#' @exportS3Method +vcovHC.prediction_cf <- function(x, type = "HC3",...) { + fit <- attr(x, "fit") + vc <- vcovHC(fit, type = type) + mm <- attr(x, "model_matrix") + n <- length(md) / length(names(x)) + md <- family(fit)$mu.eta(attr(x, "predictions_linear")) / n + z <- block_sum(as.vector(md) * mm, n) + ret <- z %*% vc %*% t(z) + dimnames(ret) <- list(names(x), names(x)) + ret +} diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..283168e --- /dev/null +++ b/README.Rmd @@ -0,0 +1,46 @@ +--- +output: github_document +--- + + + +# RobinCar2 + +```{r, echo = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "README-" +) +``` + + +[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) +[![CRAN status](https://www.r-pkg.org/badges/version-last-release/RobinCar2)](https://www.r-pkg.org/badges/version-last-release/RobinCar2) +[![CRAN monthly downloads](https://cranlogs.r-pkg.org/badges/RobinCar2)](https://cranlogs.r-pkg.org/badges/RobinCar2) +[![CRAN total downloads](https://cranlogs.r-pkg.org/badges/grand-total/RobinCar2)](https://cranlogs.r-pkg.org/badges/grand-total/RobinCar2) +[![Code Coverage](https://raw.githubusercontent.com/openpharma/RobinCar2/_xml_coverage_reports/data/main/badge.svg)](https://raw.githubusercontent.com/openpharma/RobinCar2/_xml_coverage_reports/data/main/coverage.xml) + +\ + +RobinCar2 is a package to provide robust inference of covariate adjusted analysis under stratified randomization. + +## Installation + +### Development + +You can install the current development version from `github` with: + +```{r gh-installation, eval = FALSE} +if (!require("remotes")) { + install.packages("remotes") +} + +remotes::install_github( + "openpharma/RobinCar2" +) +``` + +## Citing `RobinCar2` + +To cite `RobinCar2` please see [here](https://openpharma.github.io/RobinCar2/main/authors.html#citation). diff --git a/data-raw/create_dummy.R b/data-raw/create_dummy.R index 78415c7..59b600a 100644 --- a/data-raw/create_dummy.R +++ b/data-raw/create_dummy.R @@ -4,8 +4,8 @@ library(dplyr) #' Set seed to ensure reproducibility set.seed(20240408) -n <- 400 -block <- c(0L, 0L, 1L, 1L) +n <- 600 +block <- c("trt1", "trt1", "trt2", "trt2", "pbo", "pbo") dummy_data <- tibble( s1 = sample(c("a", "b"), replace = TRUE, size = n), @@ -17,13 +17,15 @@ dummy_data <- tibble( mutate(treatment = unlist(replicate(ceiling(n() / length(block)), sample(block)))[seq_len(n())]) %>% ungroup() %>% mutate( - y = covar * 0.2 + 0.4 * (s1 == "a") - 0.1 * (s2 == "c") + 0.6 * treatment + rnorm(n), + y = covar * 0.2 + 0.4 * (s1 == "a") - 0.1 * (s2 == "c") + + 0.6 * (treatment == "trt1") + 0.8 * (treatment == "trt2") + rnorm(n), y_b = ifelse(y > 0.6, 1L, 0L) ) %>% mutate( s1 = factor(s1), - s2 = factor(s2) + s2 = factor(s2), + treatment = factor(treatment, levels = c("pbo", "trt1", "trt2")) ) %>% select(id, treatment, s1, s2, covar, y, y_b) -usethis::use_data(dummy_data) +usethis::use_data(dummy_data, overwrite = TRUE) diff --git a/data/dummy_data.rda b/data/dummy_data.rda index e344aa2158be748d95f08011c20f3686a8938cfc..a563a25d372ca6dce1533a69d5cd301b50f0818a 100644 GIT binary patch literal 11593 zcmaiaWo#TwtnH=^Hniy`%?29ghMAd}nWS%$(fsN$*PU-#eB( z(y{FEFKgz=f|_O=3<4TtN*a*|cXI$-OyAG{A9e`!0RQ3jvfH(n7Xa0*mO5cPRlmGm z-&|#FK!gx&sihE!U2M9f+JK-@R{&yCS%`_YxoBewl?_xm`7aShf*%0_B7hjEHXEEPT7}NRmI+%>LQ2HW&Av2r(chOiPnx8-+ z8btXm2U6BIWXBQxtWiy{PB zxHFn3D&1w6@G$3OFyAu|>XVU?Nh2i6Ft0-iQMpNz0q4q^WuW^3qv>&EvuU%^X|rix ztc8e0l?-#p(yjzcH4}dFz%9}NV zf)MjFVgmrk|Ct#wtlxhgDi8VV;w+9lUfwLgA7K*<9t0l=F0U0IQV|w079~?x4+|R< zwvZe`SCBU>pvud~K(nN(UuTe%uYe^QLK7vT6s9t$RTc=#E3fTms8lxm77&Q8Y+#Ln zp&V8)sK23U&}wbfOjD&AHuzsb3~)g~r6E;3MiiQgunIa1nqfqFX#q93pu9Zb7kCb> zLX~wrjH%TW&0qvLH4o9gi+HsbUM8s1L3*tN`~ZV-P_!iApPBsc4$#hkz>TtV_xw=;{Rw z0e%q?jLpR~d8P3GwJCxoWu>|%qRNsaEL{G}$*iCtf}%WvA#ATd3S3Y+B#H#Pp_vdi zVWdJ)PZ3aB5Joj9Qi+IRgr=ZBC|p^FCV+;4hOU5Co}}`fD2#wAfF`Q6FQRk^mVzQ& zF%YO=SP)TVRvHKeORXPP&H#fEMvWmr6={O4fM_H-M}s{-Ki{Y#iYyD1$-NV1-DntPls4l=qLTrkH z2w)J;fi1KP#hpS)GuVx|ZuF5OM@^7F&c6_WBqU?zVbhh?Ht=He`0xAkR8QYZ*a-1UW?lULQ=Tj9fw-kKWzsCLd^n#YV}=<=j2!+YnJs6T+@M64N0)=LRPy zjwo!zH2(R#%psSBS#J*h&Z>bQbyfubs{{v0#wCyOY2}sl_I&$E+Z-`=u_bz8YYy2&stgISb%lo5m>fTPty->x@KUHm&P9F(a#Zh_2h zf=5wi2w39BmyF3rl%KqqX}gXnQh`M*@%py5-;q^6b{4JX3f2D~!XKQJ?yV$Dgee8NnT)pU4Gj&^rMyTsyjzfM;ecrrf za2e7j?1*|}!#h3?D=gGUAGjKs@vbt{39Q*VQ11jtiK!Ad6dXtS#PhJ4^lVs?$2@^oW{|Jsh3B* zhF90^5~2f(bJ8?=&_83UH-&@)DeQK-QjSutURAhu->Otu_{t&U{6aI6g>$=jU&!B7 zbj^>2yGYI)jcs-V4FNUt~}#!KFJkH`++(~ z+r)fnV`$c0c70@bKa!F)f__%5Gys6C{8s)jRu2}$*Oe^CD;SaJu{ba~U)Cv3$bQz6DF9`J%gJ8@>%w*c-HE)0&(vXZc0?+M2zK^-^`@Q_GO(Q4;=h zJ(kMG(7D~`5ts{pIB>X8-cL@FRs#Lor7lk8%6U%RZDL6Mi6S$nkn|-RJOCOX!EoZ` z7CL6L(v;vsdzyHAlMQ$-PW;-0MC$0$Rn^RfurHct*UBGIfMp-hsyu4N)ZzV*_kOwI zsBs3$(TC3ZmfjDI7=@9Qdi`7y+;i;_(_JTt8`z+AzsWJlor#8?-mcJ}Qf{2U_VM_e z-pCk9V(9Sj_!f*N*AkN*e%f=3D;1K@=#@U`+6fuuw~!NrhX16uo9IuO-flZDi^myJ zppO_G9?5stoDr;z9o{GS0)u|0V&0yR2p^99G^Whun3VkLg*rxpZ-1EyKB=QNbc&=v z`3Y<<#@bRKCT;FY=-Df(hiDGpO3q3ae84TMbl{|W$DXxG!igeA?+7J79%{`+q!?hWsJ4wc4FDkUX zoHX~EZN$z;6G~XdnDO%Rl4Z-rcz1+gAP>?%W*g%w?AmH!5KojCc6Z`uGJE~K8 zf|X9sVWrYX`SP7z*Ul@-yG}e?9CcM8zs%3K4XOw1XxEi>7_;vRT1hV@z=~a`&^cxx zzdts-%xp|D{^ASytb&uuA%pw+c-rTZ=W<|fAQBe7m$<=Skdwe>D)EU~18%Q2H4ew#Yk0+PDnlW$* zV1vCnixI4wR$sR@#UPXUIcfPXlJz3ZzD>9JHEXA>U^R?KUthy4_(>Zj$@FhphU@pUe}1mi-gGpqniS<$M^BC~~d}B(O$x zFFZl4O8AWks$Jxd`&Rszb0jepD7byabI{uM$m{tDKgcm&!9C8-bID>GN;3qUkjpBZ z6r0DzcWtAvfXoZNBVIX*XFoAtJL`mNkaJfqh& z-z$Z9BH=l8H(F!+NQ?EynPV$KJCz!xWU@Lr!MAB|@KSM0nF6c!U`Lgk2?kJ}MjWY~ zek%#Lc500=)&3t#OjDmHB`1n3ey!x;j7BFa`YKtXEMFLn9$Af}_z#e%uYF%Pj1`?=qYWhD0Shybzg)+QD)zRb5%nTt1Fh2t||g>QF1-nbQ^fh0Dv>p;BZ(b zfU-n>qYO;6dG?8lWQVYx#VkCKxPM!$WvbQ4Pn@WN71< zxp+Den@}VTV{1X!-G3l7>n{6fQDpy&x#15q^+`*8TwHvUKyfC1NzBv_$c!{SI6~nV zZ`cC?g|07~Kw;NJ!vuaS-{|*uhq>L=u1uFNj(RA7hHBOZ%N}*1+_>_HYe@!)EnBIHDURv%uWXudwxTZ$ zuSC{P!p>(}D>6J-p$QOt!ei>b|z=67{rb_L9zNL?Sx`wg(qatLcJ`sUtVv2K@^NmwzQ4cRqpd&rCQ<^%|d*VYSH&^e1)}z zq=J!wcMMfCYf0vQ@ky$pG|Uz6p24zt%KTR{{0LT8%?dVR9+4j@a(1kGYwemkd0stO zs?)YpLDA&Y4@fj{o^xee)&?&apG)Vp?g-^ng@qkDYJs~YpK$y$c4 zNrXr1L)H67r;=LpkK>#BDt~2zav;8QO?67>v+d&K5Px#jtJ$c!O-_tJMMjy=tZg8c z#+HA@)qN)OL4+08x9Pcd!&^pWw%XCh$v!n*C9xTv5FsY~@Ed1)U2UmoVpbi;xZNj2 z;w~!CzEsMtFAD^?-jhtC+Qk2Oj(X3(%!roB_%c;9I?hNj+jotP*(Ekz;ZaV&+raN+C^Kh_@rOENF`qI_X4J|ScQJF=0z8m)v3PoH) zZ)sgSuIq~iPcOHNc(!10X#3mJ<)p_apV`>nbD9809xaopxbt0(`@z`U)Uo0fy?7|= zB4;jRWREV^fT`SgbYXiiJ%xNJ!{+L`s1Vv{FuER5C&oe+x8J=|56|fE?w2x@td$wt z`qtW8LS#@Yh+0<&SuK5T!2Sn?=Y>j*s@XSyeD^>^2ao&2Ff)yr2>cU^mJDIpuL#$0 z#9?f<_d#{}e}+b1b4p@v_O&BzMO3|oj$-rKUyPYH4dLYr{}XYrzAxRxIdUjWr(Ljl zliq`9bG3w>)cZoPdsXy~F~;U({=8C1XqEgk{>VY?EK_X8@@ZjzkD$!-<7w+2Cz7i# zapGmK^|$MAytulIhIurN{#EdoA-}RVedPn>5c`3t@oG(e@^Go<4Ytfpb^uPMgD~%A zqki(zZP1Q&koU@VaUSek%X?R<%aIbE7Fd}^r8BXdrh^Q}|oltjeM@>ltJfKzx=Hg?_K z(C@zZTJ2;#U8tQ3#xND2sT5Zn{ps;u5Mv;Z;Uw*&wlO;9|2fExmx-8_?IjJL&kSur z_ymg|HLDM)2Y5cN`0&Tj@ZGF$Z*a~7I{qP z$We!%r7{KUPv#sTEfGDH4qfbsc631)qvmiUf=P6R_|Gxjj0XYF{70vDC~Kpe*(E1r^Bg_Wno$RiB&amFQeSZQ%ya?1VppJz&YyJlGA>ree`h!yv# z1%1&n7rcm!BexB2mncF`P4;JxYf>V;x+`^sGl${e+jq9V779YYWq#E}fRlEp34?9h zq|OzM1nIFu?_80IdlJZXPGu61;lC-R(3_%KUhXY~v^gf+S;JCog)$XfV~hRTH!OWM zrdIDASH_Z|`SGE$kVQm^nw(L+C7%B~-sD9qbIBCxDAVnCN{Z+UnZf1`VEWo{n;|D8 zq%Uf)wu~luQA5MMdFcJDrG1DYKa!yBgztVfJjI2LRZ@8ypqRM6GD<%l77Btg!n>eO zFq}XOEmhxQmE!f$$6SRk=9i~ q-@mmt3MV~YvrqNu+=we*gMu3ghaQWi_}!w#RW zR)_95ri5C}S+Acx$E7j9mH{Ic;V1cSmqk*(bm3a(2izK^LbQxALJ{#&dpH`RBl|fs zrreFUCn~KNBQnwe>1xk&aL+s$6Rh3AFJRKkHxE(V$~P!?ROskr9HGHOiy8D}f=Yvi zyWJ;bK)Hw zaRw_!ciMS*4n@mrq7DW>T@UqWrep-wo7==!>EQauh}7Gg4Q1P+#SLLFZ{2blX~n2S zP;Wad$xAXl-NG;p=NGH#gs>CCE|`nzCB-*Xc=?PkvDqs2 zxNo{n4;zJi%_Q*|-j&j`4b9i-%Ar%mf@6Or#-J315OCk)SD|$ea+naa-K8=6w0|J) zU88;0v-KhcMrDgvDcK|d!t%P)q;x< z{0KfvRZq{yw3|;xFY`ykQiY!)=SNz){O*7sKWp<$&v|a6M+Ng~js8AsJEb|U#m&+z z;zGkcK3#5^xkvMpM1%V z9(b()J6ABjkXS!!AmGLGoB;fX;xK>)m73o%32zzHJBqjzapyv0l9#2K;LIRbwi6g@ z&1}Ej7Ij|ygmLW~;HHc?M(EAiWPX*JgYt<9kzQg(y_lLTn^7C)8)4(SM^PRj>J;@d zaiJ$@8RELhT5%wYL|?o$NfT)iD5d5lv^?2o8 zJA!#ERwg-YSUcIKd$}}2P5Az`MwoRdbqbMlVM#tUbw@IQimj~2(^;s0;VZH4$_(4U z59S*NbJ@<&lz*AT?uIoB>5;~C$NmcNQxyl))EbW!Vldv+N}WcpskD=Fdhy{B-<{)q zgR5{L@HBx#F8Sze^Pu2jSwrH5gT*ZJ&7;>B)l}YjmGu4A6Phv>i$C!t(TLdj`~rt4 zR_`0_yFLDITipwgvHpQMDyr1g-oEWYbhKn1zS09;4Xzu*z;aWvheWu0WRhp-mED@f ztf3kp@tF=ShTz#`)2sNH&YvRdU<`$bPDY~pbJaMs4>*;;r zmT0%yv^FWupQfNRa@xV^$YGKG`reCDj z**WYH9d;b^P4Y^}snp1(eJleU$XhRUtyyOyg9+FsPVU|dRih#Xoag|XsJ{5|d{qsPCs zqQcOyE=@Fh$tA!={7TO8)e};-VS(9fJ~e98%G*>J2ugouwfS@e*}I7DH{}?*;&>{A z*ZoEl`l^lNSRQEuQj_G>5fHNZwzK@x=O1%`FR@2aL?=?GZ})OL z4h13Xs$>%N2BO1>$Jp{%(#y~KBRN#hbHK=~3H>%;BGi)lnOhiE&M#e=9TPt6Nl}fN$xv|`MH4T3*(gk*&ki}6`J%iQEF;2xy4m1-+T-Pcxugk zYCmj+!Q*|Ul1Hs%K#D%hekQ**CRI}M-VvR6!NW{|t=V*z)diwMo%s{rl9_?1ULS7C zcalX+xmc4xi_}RZ?&EXK0fP2}wo&Knskwy+MTj8SO2bWO$MlCIam!Z%`ply^f`)4A z5|N9?R)q1-JzC&6B+yVOHx&V$ru^#>II`@Ge_#E!of?1#aO0+9Fhwv*ZwRY7aOv5NDMS zt34fK3SpW0oL5zlm%m-i9a}`em-j50vj~p!ll=plvfM3p>{zh)<KwZUMOAl$7AiI7pEn$?@2J&Z zr1n5k4oYR>8JQ-qLG`Xl7!3rsMJvF|u#Vtc5h^RxSQqh5>Vgrmrq-C}zkgryWaPn$Og8a}BM53*vnY7)cc{9_>D?vVYg6WvpEbs4P>r8FF z(42&-YslbM@YX`Q36`j+d_y3`^VNBvLt<=tQc0*0G# ze;3Q6H~<9Rsb*2A%dm!5gccpr+r?7(ypk_vCpx?THk%W{1Di&Eg3bdX4@O{x-&wC% z%C`_zMiMjgwWwa~+U)kWB7&)%#C_3~e06oGaxA>M)T3t-Dxlq+@?jqAT=0c6Vx=Or zMt_E#j<(UAeZZ=4$_Ueic7ntN*-H3U%@D7+ zmjqK>{EyQKy7T@H46NB39~uE5LitsM1Tr!&({p{9eWqlO_MGF(y!OBD%n608wCB>y zH|=0kn4g)n8H+_kXvrIDM|UNgs5VNk5^+qmi)zTIEvmI=SO0eHSb^Rm@RZpz%Hzof z^QX!mIWHuOVI1mwRsvnftGW0+9mbiG?-c@-0>^4RHo-T-{mkXx@tgV7b4RzR14NNV za{n}dvU^u9?Dq8U6$p4y2#KYfi>;2HKXU9KR^0)qADhormP$NhIxP>?NqQc;hNb%h zzcXR;H4v}qB9upM{=D>hnwgxS20>2~$V6;8(xBXsl`*KOG5@Zh+4Y)9Q+Vw{-C>20 zsdM|@p2Z7w?*OhI?M)~KZ}uYBO%&*}#MYQUhVOS@2+xQ)E38{iP2fJZ6%);KF~TM6 zub&5(9J^DM65&_)_K2$wet8A0U~#7-QeiO$uRnH~7(XFSrQo!5qbOFAr%vuvx5LOC zDsd=u>*xF2stLUJGV5YTZAJ7l5Tf(2J2!L;Q!@WMG}%}9^ZiRt&i7W<^&!5cRC6BI zr7+J!7;OpqUmiSXy38ORE6-)t>Uv`jMa97rS}MmpkTrB)CWG!?8Ta#Gc7W_mmuEU2 zA$ep2q~pUF=IL?Z!Z6L)uhc^8X$_ckGyj&cM)_vS60TaQYBKY8^N-W|X!TfaZ6&+A zRKl~y2ivkw{T4{nbQJaFhOJpv3~ny^jc0u!o$$ew~`8muT-PQEw@h{HUM?6 z_jjn)^pJ%r&Naas^;(bb9VZ2j7=2iW06v{%Ej#%ZoCg&U!$F8j`bA62!fS(H}`S9*!I!J5QvVkH2NS}J67T{2@(mGgw`im0DQ>*&#_>CNXi`PnZd{J140GR*}!oWECLx?#mZwpRHj+Qt5-kd#an@0vNpJ@xu{hnr@->b^9! z3Xcd^vSa;|&SZC7m6YFwTQ8}hCE42~zcaGTR$uwP=#%gpJP^5jAqxZ(0!~zEO9+9R~)tFWNt==NU3wUrUkIn@Zy11?QZey;b>Ne!I!&$g_Q{2^7 z>jOM!L({^8wFZ|Ua~ zc)%;5;eH>@Hnlg6=`>%UP$tU|4ZTWzGKOw`9pc(_u8&q&apcrA@(rTrIoa6P&$e2F zsfH!BD}ufF#a6mLGjS{m8roRaP!l6kO=N%UWq$}B(pF6Ga~JpB?srLap#TCI;%WRk z_wdU0-;$E3G4q64_(VIcIRZzY)R&})lDrqWgS?PD=e(Ub`!x~^duVSFa*tsg5 zrFs7$zKNUeX+S3WSEl}iYTmlwM9|{0Fh9O@Pf>bPuFreN(K0%zvGmu42`_^rtt5XK zawlOlfcv?7JzTQM=|YaUveK)O%>stHDw@{)d+cToHuLa){C^V{Owv%D3*75(Zf%sO zZ`<5*#LiD!XE>Jb;4X~;zGZGyDurUMvIq&uwYQnSFZMeG*8h(F#riVv5f~b_Z%E!O zu}Pzk!|LLP9w^W-b;j|$>LQ^2OL$>_W^C}taeQ8E=uv-#)aYfj*{UZiA8?Py741&< z+QyF2bg$|$73)pE%4PW$x+n7E6~4so1&Y3#7u!7mFwiEw1om$Kjr&Xjjoy!a3le<& z+aVg`{;y;Icz4y3OHPJfsg;y>R4mdea%PcA3hKknC|_bCH}FF$07@0^$&*y~#_JN3 zVzNf&n`gccvxVFaA_yk@;9oBlxE3Mk`>7J)sI=yB`)U(<*#>XMcK?UjUs>~1h3Gc_ ziBv2Tg%Ou5%YNMnwT;DUiAC2-IRIL9iff)@Pr-#iP-F;z)`T$bzh6?|H(vP=s0`?@ zsioIT=4PYA#`;)wJ27Dj4p|~USA}t52fFq99CxP3YO-@+cNt$^>O1?2e&+Z3LsT7J z?D`;5x3tH$2TnZI>!GtU&264yY%yY%PM}()UPnyg#&sj#vvypv{YtE1_ses3@o(3q zG$)}49s}qNn47Vbc$GP8XS)Np0@!saIs-wl->Q-OOQC?RkvUzb&oW5EPlN%S!V6aX z22x>ak$QoiWLIyEu1w&y=C*{~&TDaUip}ft_V&-yv!k=lb zjCxjOGLikol__45|yHbZg27xqp%>4p>sy?K%BWjtY4}{b8RHqVnIu}IBjW^H)#Pw zf0@kVaM{DHS0;nudC#X;;|CmF{ibD4jofpyydGt z!|D&&7@*kT&&soGwdZ+)@5VoO1ZUjsUh+{&_S|-S*BDktHssVbXm&xuMvFPsN zwB>tc@j_@GClILZuG8j5_-MnsW2d$-3xhi~&)2hRQtic~qln+7cbMeRSGUfi|4P69 zTWeA0AWH0&00JXZxo`WWQWgdCT%QdubxYZjPc~}%IxS}bR|ITfgSZ0f+b0$=? z5RBko_O$!X1}IVHifIWzG}I?WJ=Gdk@0)O0a4;REN5_hi)NmXn0^l=wkPI)G3-$Ow z22{LvQjIS?_WX=!?G2_-w)H{MXkkauNR8I9MRuRFUEUOCkYQ8wKfNcmX&%2B(sPR3 zCV)L(78FB^rSUI2ds`ckTD;jb38NARRNTG^9YIrv1Sv0pTGI%^ej&rHl2tk2?hjU? zYu6)0DcDyZ)TFXJ^<)F+{b`LHoI+Sb#?=?kPLvBQh2TrI6SolxtfSc*p#KxS*QvY> z&pXRt7Yt6!?Bs647zIF}=b52PN_0fO=xBrvHAFJo*vwAotW!*dW9CRFq|omT6LXKC zJvxH1>+vGld5_Vn=Xt=wt7N4hwcWfGzMfOre1y8*3qP_7`Re!I!6bl+`bX3FJZX4+ zwh68d2o6nGG3kOYn@-&T=b#JX>ck(Gs-OUFGukz^Ui{;3OOT{enl35Bg?1OG*W;SW z_K6bPnDsjn)`+G1 z*{7C&b@P(RE@r*_v@H+G!)_5DJ{VSFqMb~z`h-l*QD+4TD!-nGJb_ujf~&3QxSk~T znxDRSHd0P||4L}DNA*JHO8DA(=={6SJ|m|TFQI={yO7qqvyr5b*z_lNdAa9hqWjB6 zismLK!cH)U@xm2mpO=xIfa%ru4vDk y9$7jxWqpSfpcelA+E45U2I=n!^8pUl5UXhZlvSLL%IctqmGhpq@=q!;H3McySwv10g-O$Zl3S+7yLdm zJF_o#cJ}qk8QY2R$r>^17>)a-r6AyEe0%#}Y!JU8AfEYdZXg8S)u_32h*x(sP_=OB zw>lW-{WYy0Kb>01RM%genrJ9JUhJesKzKwL3{70H$x+DGuv4cg7iNY)%hGNn<>5Lu z(WOP#5jt=gsLEEv!a3U}ODa%Z5Uw-s3k3T?7lp$Z;mm$z6}#`o!Yv`b-M=RYYb;-@Yo&Z=zk z+g)dC08x{hD4cL`gd|iB5TO6_F=3rvj8X_)s#ZQ+JztPmg2FBb7?=HyUPF>DNl@3{PUTR8reag)(q6dYN!6DU zpngL1qH>S}kgVup{80>Dl-2?J32A_xa0CUH6|Q6}4q4Jxz`i1Y&Y+xp>klSp{^kV*BL8@<7lvygMG=E)djK3DoXjY1^a~~VK|5#!1phGKZr(gc2c(zHN>GN5#mb3F2@>&rSu*X$q%X8(0$+1_V7*;FN zE9j#~6^0Q5h;?0|MuSfhs7{4n)7@bMJ+0RMFrmRHioBokTkDgZ3=3o&D& zl8I2r`eF)WlJTuj0J;EF1@l46(Q3;ft@z1XP96COQUs>SxtxmJE>!Q{?arj^7gSc; z8ahDmigG%f8~4rLiXb><@WJ~KLDGGQv#?*w)8+QW(0N@P}8jum!2pnDk-a| zY%rM;mm^6}<_SF$pRZd6?&{%6iu<%#;TNVO^y`A+{0<$ z3N|#%Gz^<$RaJe>l9cLAn?-OI6O8Bk=O%&!Lp`;Xv?d`6W3D^!u&~0gSC@hRE4|^F|c$(k)z3fmQ*Fws=v-c&wn{O01*2wr$WE7OY~^b(dp+R zDzO4OA@J|JvjTwwWv*qQ<`H4G3J68e9yX#qM&~6X#jl2O?U1+m{3>7&OOC)26NZ6g z$F2#6ljwWH;E&bShT>i1HVg7&!kWn*&&;jWA1vJSG(Im>_1|BsARsW_%y88{h%5qI z|FRDuhB`?sT?1qFEaQAJj({fSw8*Tjw3tMQ!2nFNskUAh+4HrZr@H-Jup5Sl0BzKi zBr8IiV~2iP9P`;9Ji+6x!M6HnCPcgfcf^G8RK*8L2tFwqOE7W=pF40&lnix2FDix% z0Fwxz7=)W$sYLT9GSvYtkT9R?OI9bu+mY3Dj}7_~zy}0%rdDGlsW9>m3shUL9rTma z)hn$el(;@#a5K);@D@W9?*R9r!uPG;)&=KKTGj!zI#N&V99Hdmvzz|CF+>QJN|;BW zo*6ja;*T5H=oy_lP$3mN(4FkhZ*+mxnRjJ(*AHC?%mEs(Mc%kzln^8tJa?gTf>Yiv zc?}(+(liq^A~Z+_pM}1gAR?kKXH{h2UY9vQzW6%1Fc!}(yTrUw^H(^(r7UMp5~xr; zs~G4FnnKQ+4|7&ee}UM{7*mp~8i!*-UQ8Ta7<~Q_9-QB&BCv``;zX-N?V%LF5!HP997elc2l{`oEIw}{bVwC2rr zu}euA`uTe$d5ILJ8QV-nCyElEi)KtuKV@4>i{yA@!fj6vXiU)oKhgMxk&zDo^OQ9( zhc%1gnYU;8Xb0iNS`)A9qWRQs^kV(#p1Y%REw)H;dUs`R!p&acm|OgKhOpHLphqk~ zTQV>Zi*IuYxpca;+ngE`a3mb@OG;=r6|!fd)EXCl6;$R7p8}F0V-`Fhf?Tal{(fCX z`R7iZttA&GN2>gDoT4xQgG&j)tQuX!Lj$;J5xnVYn(CO~UongNE@g1J=P&t`9DSu@M9+mv4D zYICEBOD(pMxEC!+*)c0Q(3ZEeGdIw(n19Kr&qOcii1<--`=47wyR;~2H`T`yg$|@a zBI*O$h{wsd0~k;_Lv=mR+67gj$=*uiyj=e?@oWi}<*q^NEhIrsBckV%JR)!3kBQNT z-8QYY4^J3Hq^DL>zgr#o7@E$lQ}2$edM zOdA!7L(s$5%$pwby|l6ccc<(#ExS|T1Zv3F=&6u6j_QW*2RbV-+gkrgvy_adD~(>@UM<@ zalFhBrZvwd^Wk}sl6ZD1)mw4@+&|#P86tGO8YyT&V6PwG&*8N;%2lqAj3WA z+$57xhlIbQ&U!1Xsk#HS!$6AhDL>ga+5F0`XIo;zf`AFZN<#n3Fj;)lFUgz7Q<%(6-Iq4sq)QU)eH5IrO%Y#8VnKyC%6cyX>%= zzrJL$Vb;P_vi*E-+TS%Iu5w!K#qZ*d`h7e1&k%kz9(-lroBVk&jjSQC;#DgI9ZkcZ;y!fZ? z#OM#9Ziy5oo^v~PbseMbq_DY>*C9t!ijCb+nACJ&cz^PEgZgX2-=pf1N4wU((fAil zjzV-Wb|J4Czm`>Nc|Pv$IbTEFHZ2IWEnv^pTA@%rXKMYmvJ%^oD=qk-YFS`x)Q1za4_IpR_;X!eDHCY;Nol}%b*k!AMVj7MBF^J zHYYqF;Mlmu*`8N+=n>Oph}qP9d`PWFwumJZKz{$soKv%U_+WH>BAJzS!k9Z+U*`%t=jVc~qoRGL}vkBQk);WwB2w-$fpFu|ZihHqd< z;+G$cj!kfURu)DxKYxaVW>#f4sD_hg2Dr7?DDum2toz0~W#OvDs~yx-%(1f?-F|vt z*Nn`6>)c{y&AnJpc3MQ*1KY8uXt{^7;SPUAG``dMHgj;XKdhS-#27xJr$^hna?2+^ zHzJlLz1(D67HPg?cJ7LUER@GT`1BiUjq$06f!?gauP3W)?e3`NLdOJjtc1HGGt=p? zK*?UQe~Pm7S2{=wu&GpNRjiFzP-NzkEI_qBQ{wHoYw|d_{ z1!Dx`x*f5>S!rCf!g*Bsdu^RN26|#d*LP@U*!6t1puwkQ#%=-k5&)5!LmV(({fgzZ zqYDTf4o7SD-|Cr>`;0J{uoRaT5%aR?%+0T4p?0=cnckAxd$I4IH{t!sE3)c$2tsv)++Q-|VmA<3PCrvs+@91vYjYR!_v=R%)>i06bQW-n z7EetY*2<|?3I2s)7Zw}!Jl|lIyvDpW9#`IjL2a42Luy5H0M}Fg>dqv1V%`uKn-|yo z>HOkJQ2Pit(|~wM>WPf5*7&`^0Kzvhr$dq_Z43|Ds`Y^}JRae38tFQ5Mfij<`+TT7 zt$d=vUGPg($hFrlf_3gzD|gd85!OMV03#I8ez8R0mQ9a4OziK}#Mvhc$Z^Rh2kCb< z(0{cd8?><3)w0i8o|@CkBp2nUvv_nGo=a}yPa~~LyFUWO+V)ufEy9@h&z@A05~q4* zb3*L>*F7FX9(afYUgPVu3&_2!L$O?aehKMhe#31o1FRCApq4X6S1gs5p3ybpN_nV` z_sth6UX2!cMK|45P6%CjNb^fUi8Y!woNqmcn$iE9aEO{EehN#%QghTnzveDdn7|ei z-6Gp`@r~`^ZlyYl?U=a_+A~K)#aGj4XkaX!CyUM{EgS<_+F{|v4-x0^*7WgA3f%q) zx2|>5?%qJKj`~LRX*(ffbr_4AHJKaQ!n}VYg3Ek*q=ZI4endk)j0iIZllCc2Y2>x?IOOYhX0Srbe6)%2QfbY zKavQ>9a+Bn{b}3zXI7*<2Vs3p!hA|R9OU;ugA8m?{qBC7iw#TYwr=(sntG$g7d43H z^~B3!Aid+IGi`?gUFw+)>&*EU>;JLn714gn#U$xh9d>u9$VtA4BDB>0q1Ry#4`Eg-Pr>GA`=v!!^g%_9WI>)k@s7ZH20n9X?(7i+8Tj#*r5nRf08j z8OL{?K^t?A`axdWKw^nYwx{P{Mu+=ucfOoPiX=X!r@6dp3BB^w%~Z_#Mv2- z&$q*}pWKvft;Q8c!NMcnPD0)7&Al8RvxA1-p2p|}Ivs7(lJaT#9M{`V?(Eq{PNpOk z-s|+Q@9LvkC-eEq-y@*DS?8I3v+EYK-EJeg*ls)~hs~BfEu|}LYqopsg~p|-Gu)YF zrf%p^ajqL$?Fzwno5GRtS_SFp3jQGL0a0)5lu+0Lt~pJ#F|TEV-ji@JkFj@b98&-74vVr zj*)JprgRn2RMFNi;YWA3X;k$k&fNNT-L!0F|S3qkjCYVMH zpR|)`;fg3%20@AS^n;L{PeK~gKN$k{qe~kbH_UyhsTU^ptl#=oD411L=Io>|)~|`p zBY_L4CqLXbjY$j*PeVo!kLr?f2!+zy>zhM(uu9mIkuJpr30k1y%fzw$k4Eobv#g2xHRaxot1z{?|Y^Zp$JjsxBM z%#vQ)T%@6}QL+b7$(n#i2Yi`;WWLwV=7X`dnV~=-r3OPNgUK(U?x?{?77o7}tfNPT zB2`6oj1<#1|L6U9u@2IHBXbOh)xf-Mtmmka;c86bbfNNo=T8toXZu4_d{Cp z(MNKM7V7w8fHOash#w1KJyv5_iicI7m}!Pxrdv&&{WtS_Y?y2Nl2?lA=-N|Ji@<8< zG*yCX)P>wP+MP--pdn8;l@=-4B-ofPn!`X!;J)F&>fh3Q-`PN>bYrRVcn=0cjNmWy z?X9ubtP$Ok0<@RACVQdCRJrz@Oq2Gc;zd{71!R&_vrN+%ai9oE&5-VpJH=M~7+r#G zJ&5@;XP&tiOb*JL#LZds#j`{1V)?%I$>#1b(%7L__(jM`tSW;XyMpw6Tnv(L9*)F{=Yl6u_TwVhcu>(nYYaJ*$S+$JRVMcj+_}RoSg5TX0brBGkg;|^e8cjir_^esS1T5Tb za}6QbG5fN5vE{wb-%jZ`lD&VPw_%wK2K}=eyXA(GsS@Idnb4F2Nk%>t>BD(7&%4md=u zzZhmvir0n8o<}+7!&XVse!tL^@FM%Y{RZRmc-?PLYj$$>k(l<-MIM(GAm9t z8_xJH*EH~^zA-YP+F>jX zyl#U#BESF}ac4)J&1!%`y#{NVG<`P* zB2g^2{^i!$NLA|6QxmFPz9E(9vx_jM%LvJZS~|sCd$A$>k|JPJ;QN%~cM(2TLnqUm zs%0W)OpYyGxe#NirhQKaV08;9PSmsDjMNkRiJG0x-zz9<++Ka2yUw+$vO zbfi?$V5v2$tw4zr&s{=5XP|FhRhhozIMT#p?03+6K)2Y*PP8z?1@ao7Zo1(T zff-uBlczhi>guW@_=el_T^LkN!W(P;xaBaU>yPRWV*{~;0q*31eYL~~VfT*fms=!m zi*e9co_cN5ci9SXwJfk&w-F~Y~t=E-1h+dwUq=P(#{AHqo z@-GX3Ugim_dX*p&pK;A^ z{VP{EHhML2MnuB)@nAmNr7F7JxYwiAQJ*ZG&@Qtj)DQZz)qzgWGglA6H3-M;Ev*Rb zF(%8`e+KId(mv8jsCq9wP`Bg~j~DCGi?Q$a!x zO73Tup~p>L;^P2ej8O7?_IHG;g*Fc?Rz{TK|8c`cl6GQu?#8^V2eJ%9~>-%`0Ey5;gfX7zwmv?0xvvEg%^~xCOHCq zM|?{-+&FWw=}Dp`#pVT>qEA9aYdi*VrTC5}rDgRao0t>Zr~O-9a5 z%{w*E;)HJJ)|Jh>ZY(1FF(A|EmB02O0aBn@y@E6D_&5z)&sM(7hg9?g_2n+@DdveF zr?fuvDua-$Xf#kt5oxaWPakU?#Epzj<@rkFJXxn-mr=pA- z>_aOlYIhDr`h_A6eVem!)nKJIVt#*vV1UB-edJy=E!< zgmsk#4eQfICttpiLJDJN5{#STdGd_G9t$Lc@;fvhbNiWJdDkf_+73^!y}QL++rklH z(*zk^^!u&ybMrRLBgbq7(Z2c$A)1}@~&Xv@WwtK$h{ zEs*eH=`m?GB_5M2@kAAV3M*V$c}ai7m23|R9}Pg(akEA{lx^4|7(fo_5WNx&m>u|I zvCpg*V6LuY!K)MZ?9Ix_(qGE$0jWh=i{gWDC?hRI8CN2^ni`pfO*s^JHr9{ruWXV) zStw@r`|x$JeWN*jRssjV-lK;rg?RGd_j&8S8eNmg4fc$3igjpLB6AP*dp%5tq}naCk-LXJ-*eHDMMIck|8>ZG@oyARo|M?NNe1mjqV(rh zL>aEM=H3|*<6LiVedjPv0iw4<%ok%ldP)m>)k6!P&lA=GNV}A@GmQRQ&q(d%Hy2^c{S{r?^lvAfaJ%Ruffac<`i^AZrCc#>ygf%@ z-?0YN(hn||U$Q`YdK?zv31FL)Zey&Y4-rrJ#kYP*reFTa?O|!73o8oZh`C*$)b90i zERFqQMLWPVV~djdf#jGy%1jiB5b^_uoI({z?Uxe4ZQ+B>;rLyBC!k8o6QNly-EJnA zmscA#f<~?6(wKOhRX4#ND{YTn{wlQq7S!ceNMfccTvq~JYd)H@->k$%Q{|Pj)NGz# zSp;iR3UezW1+BE^D$3Szu%C3rzsrB8)<#=MY0xz}^2IsZ6BHjc*-5GPb3ki1zw|M0 z1&Rx4@vMrrGNX&wYClY*x#$*t{}G<+*O8W4``6DRdwXOg?_AFYK{9zx%;)dX7Gb?n zvX=xlx21Mqiq#RpjLIs&ce}-LYQMf{u%lr2d&Xd$R6J2eQ1`~ax7arNTvlakNmUhu F{{yhqc@6*o diff --git a/design/package_structure/design.Rmd b/design/package_structure/design.Rmd index 948463b..788b467 100644 --- a/design/package_structure/design.Rmd +++ b/design/package_structure/design.Rmd @@ -99,7 +99,7 @@ pred <- predict_cf(fit, "trt") pred ``` -## Obtain the biasness +## Obtain the bias ```{r} bias <- function(fit, trt, strat, data) { @@ -123,9 +123,9 @@ bias.glm <- function(fit, trt, strat, data) { } ``` -## Ensure the unbiasness +## Ensure unbias -to ensure the unbiasness, the predict_cf need additional argument +to ensure unbias, the predict_cf need additional argument ```{r} predict_cf.glm <- function(fit, trt, data, ensure_unbias = TRUE, strata = NULL, ...) { @@ -207,10 +207,12 @@ vcov_robin.glm <- function(fit, trt, strat = NULL, preds = predict_cf(fit, trt, idx <- split(seq_len(nrow(raw_data)), raw_data[[trt]]) cov_ymu <- vapply(idx, function(is) stats::cov(fit$y[is], preds[is, ]), FUN.VALUE = rep(0, ncol(preds))) - vcov_sr <- vapply(idx, function(is) stats::var(fit$y[is]), FUN.VALUE = 0) / pit if (sr_decompose) { - vcov_sr <- vcov_sr + (diag(var_preds) - 2 * diag(cov_ymu)) / pit + vcov_sr <- (vapply(idx, function(is) stats::var(fit$y[is]), FUN.VALUE = 0) + diag(var_preds) - 2 * diag(cov_ymu)) / pit + } else { + vcov_sr <- vapply(idx, function(is) stats::var(resi[is]), FUN.VALUE = 0) / pit } + v <- diag(vcov_sr) + cov_ymu + t(cov_ymu) - var_preds if (!randomization %in% c("simple", "pocock-simon") && length(strat) > 0) { v <- v - get_erb(resi, raw_data[strat], raw_data[[trt]], pit, randomization) diff --git a/inst/WORDLIST b/inst/WORDLIST index b1dfb39..33d4d44 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,2 +1,4 @@ RobinCar initializer +pbo +trt diff --git a/man/RobinCar2-package.Rd b/man/RobinCar2-package.Rd new file mode 100644 index 0000000..e99d026 --- /dev/null +++ b/man/RobinCar2-package.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RobinCar2-package.R +\docType{package} +\name{RobinCar2-package} +\alias{RobinCar2} +\alias{RobinCar2-package} +\title{\code{RobinCar2} Package} +\description{ +\code{RobinCar2} implements unbiased prediction and robust inference of variance of a fit in R. +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://github.com/openpharma/RobinCar2/} + \item Report bugs at \url{https://github.com/openpharma/RobinCar2/issues} +} + +} +\author{ +\strong{Maintainer}: openpharma \email{openpharma@example.com} + +} diff --git a/man/bias.Rd b/man/bias.Rd new file mode 100644 index 0000000..4affe68 --- /dev/null +++ b/man/bias.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bias.R +\name{bias} +\alias{bias} +\title{Prediction Bias} +\usage{ +bias(residual, treatment, group_idx) +} +\arguments{ +\item{residual}{(\code{numeric}) residuals.} + +\item{treatment}{(\code{factor}) treatment.} + +\item{group_idx}{(\code{character}) stratum index.} +} +\value{ +Numeric matrix of bias in each stratum. +} +\description{ +Obtain prediction bias within each stratum. +} diff --git a/man/dummy_data.Rd b/man/dummy_data.Rd index f684a01..bac2655 100644 --- a/man/dummy_data.Rd +++ b/man/dummy_data.Rd @@ -5,13 +5,13 @@ \alias{dummy_data} \title{Dummy Trial Data with Permute-Block Randomization} \format{ -A data frame with 400 rows and 3 columns: +A data frame with 600 rows and 7 columns: \describe{ \item{id}{The ID of the patients.} -\item{treatment}{The treatment assignment.} -\item{s1}{The first stratification variable.} -\item{s2}{The second stratification variable.} -\item{covar}{The covariate.} +\item{treatment}{The treatment assignment, "pbo", "trt1" and "trt2"} +\item{s1}{The first stratification variable, "a" and "b".} +\item{s2}{The second stratification variable, "c" and "d".} +\item{covar}{The covariate following normal distribution.} \item{y}{The continuous response.} \item{y_b}{The binary response.} } diff --git a/man/h_get_vars.Rd b/man/h_get_vars.Rd new file mode 100644 index 0000000..7a63c73 --- /dev/null +++ b/man/h_get_vars.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{h_get_vars} +\alias{h_get_vars} +\title{Extract Variable Names} +\usage{ +h_get_vars(treatment) +} +\arguments{ +\item{treatment}{(\code{string} or \code{formula}) string name of the treatment, or a formula.} +} +\value{ +A list of two elements, \code{treatmetn} and \code{strata}. +} +\description{ +Extract Variable Names +} +\details{ +Extract the formula elements, including \code{treatment} and \code{strata}. +} diff --git a/man/hello.Rd b/man/hello.Rd deleted file mode 100644 index e7eb880..0000000 --- a/man/hello.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/hello.R -\name{hello} -\alias{hello} -\title{Personal greeting} -\usage{ -hello(name = "your name") -} -\arguments{ -\item{name}{Your name (character string; e.g. "john doe").} -} -\value{ -A character string, capitalized to title case. -} -\description{ -Greet a person and appropriately capitalize their name. -} -\examples{ -hello("james bond") -} diff --git a/man/predict_counterfactual.Rd b/man/predict_counterfactual.Rd new file mode 100644 index 0000000..581f5b8 --- /dev/null +++ b/man/predict_counterfactual.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/predict_couterfactual.R +\name{predict_counterfactual} +\alias{predict_counterfactual} +\title{Counterfactual Prediction} +\usage{ +predict_counterfactual(fit, treatment, data, unbiased) +} +\arguments{ +\item{fit}{fitted object.} + +\item{treatment}{(\code{string} or \code{formula}) treatment variable in string, or a formula of form +treatment ~ strata(s).} + +\item{data}{(\code{data.frame}) raw dataset.} + +\item{unbiased}{(\code{flag}) indicator of whether to remove potential bias of the prediction.} +} +\value{ +Numeric matrix of counter factual prediction. +} +\description{ +Obtain counterfactual prediction of a fit. +} diff --git a/man/prediction_cf_methods.Rd b/man/prediction_cf_methods.Rd new file mode 100644 index 0000000..97fda7b --- /dev/null +++ b/man/prediction_cf_methods.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prediction_cf.R +\name{prediction_cf_methods} +\alias{prediction_cf_methods} +\alias{print.prediction_cf} +\title{S3 Methods for \code{prediction_cf}} +\usage{ +\method{print}{prediction_cf}(x, ...) +} +\arguments{ +\item{x}{(\code{prediction_cf})\cr the obtained counter-factual prediction object.} +} +\description{ +S3 Methods for \code{prediction_cf} +} +\section{Functions}{ +\itemize{ +\item \code{print(prediction_cf)}: prints the prediction_cf object. + +}} +\keyword{internal} diff --git a/tests/testthat/_snaps/bias.md b/tests/testthat/_snaps/bias.md new file mode 100644 index 0000000..3fd3322 --- /dev/null +++ b/tests/testthat/_snaps/bias.md @@ -0,0 +1,608 @@ +# bias works for guassian + + Code + bias(residuals(fit_glm), treatment = dummy_data$treatment, group_idx = list( + seq_len(nrow(dummy_data)))) + Output + [,1] [,2] [,3] + [1,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [2,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [3,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [4,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [5,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [6,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [7,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [8,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [9,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [10,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [11,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [12,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [13,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [14,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [15,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [16,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [17,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [18,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [19,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [20,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [21,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [22,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [23,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [24,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [25,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [26,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [27,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [28,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [29,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [30,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [31,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [32,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [33,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [34,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [35,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [36,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [37,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [38,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [39,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [40,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [41,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [42,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [43,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [44,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [45,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [46,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [47,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [48,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [49,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [50,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [51,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [52,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [53,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [54,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [55,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [56,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [57,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [58,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [59,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [60,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [61,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [62,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [63,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [64,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [65,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [66,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [67,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [68,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [69,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [70,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [71,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [72,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [73,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [74,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [75,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [76,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [77,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [78,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [79,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [80,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [81,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [82,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [83,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [84,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [85,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [86,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [87,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [88,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [89,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [90,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [91,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [92,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [93,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [94,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [95,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [96,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [97,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [98,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [99,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [100,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [101,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [102,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [103,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [104,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [105,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [106,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [107,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [108,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [109,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [110,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [111,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [112,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [113,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [114,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [115,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [116,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [117,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [118,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [119,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [120,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [121,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [122,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [123,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [124,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [125,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [126,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [127,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [128,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [129,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [130,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [131,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [132,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [133,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [134,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [135,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [136,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [137,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [138,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [139,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [140,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [141,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [142,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [143,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [144,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [145,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [146,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [147,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [148,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [149,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [150,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [151,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [152,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [153,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [154,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [155,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [156,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [157,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [158,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [159,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [160,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [161,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [162,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [163,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [164,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [165,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [166,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [167,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [168,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [169,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [170,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [171,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [172,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [173,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [174,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [175,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [176,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [177,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [178,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [179,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [180,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [181,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [182,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [183,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [184,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [185,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [186,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [187,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [188,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [189,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [190,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [191,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [192,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [193,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [194,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [195,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [196,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [197,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [198,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [199,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [200,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [201,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [202,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [203,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [204,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [205,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [206,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [207,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [208,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [209,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [210,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [211,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [212,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [213,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [214,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [215,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [216,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [217,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [218,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [219,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [220,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [221,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [222,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [223,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [224,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [225,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [226,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [227,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [228,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [229,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [230,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [231,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [232,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [233,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [234,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [235,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [236,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [237,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [238,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [239,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [240,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [241,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [242,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [243,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [244,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [245,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [246,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [247,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [248,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [249,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [250,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [251,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [252,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [253,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [254,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [255,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [256,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [257,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [258,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [259,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [260,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [261,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [262,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [263,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [264,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [265,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [266,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [267,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [268,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [269,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [270,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [271,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [272,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [273,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [274,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [275,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [276,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [277,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [278,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [279,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [280,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [281,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [282,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [283,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [284,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [285,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [286,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [287,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [288,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [289,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [290,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [291,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [292,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [293,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [294,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [295,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [296,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [297,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [298,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [299,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [300,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [301,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [302,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [303,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [304,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [305,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [306,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [307,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [308,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [309,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [310,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [311,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [312,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [313,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [314,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [315,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [316,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [317,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [318,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [319,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [320,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [321,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [322,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [323,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [324,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [325,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [326,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [327,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [328,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [329,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [330,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [331,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [332,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [333,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [334,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [335,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [336,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [337,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [338,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [339,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [340,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [341,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [342,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [343,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [344,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [345,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [346,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [347,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [348,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [349,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [350,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [351,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [352,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [353,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [354,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [355,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [356,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [357,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [358,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [359,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [360,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [361,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [362,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [363,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [364,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [365,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [366,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [367,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [368,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [369,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [370,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [371,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [372,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [373,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [374,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [375,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [376,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [377,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [378,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [379,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [380,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [381,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [382,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [383,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [384,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [385,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [386,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [387,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [388,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [389,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [390,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [391,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [392,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [393,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [394,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [395,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [396,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [397,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [398,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [399,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [400,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [401,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [402,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [403,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [404,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [405,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [406,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [407,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [408,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [409,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [410,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [411,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [412,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [413,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [414,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [415,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [416,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [417,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [418,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [419,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [420,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [421,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [422,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [423,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [424,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [425,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [426,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [427,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [428,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [429,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [430,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [431,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [432,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [433,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [434,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [435,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [436,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [437,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [438,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [439,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [440,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [441,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [442,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [443,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [444,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [445,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [446,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [447,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [448,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [449,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [450,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [451,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [452,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [453,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [454,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [455,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [456,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [457,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [458,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [459,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [460,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [461,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [462,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [463,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [464,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [465,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [466,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [467,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [468,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [469,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [470,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [471,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [472,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [473,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [474,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [475,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [476,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [477,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [478,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [479,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [480,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [481,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [482,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [483,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [484,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [485,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [486,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [487,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [488,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [489,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [490,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [491,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [492,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [493,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [494,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [495,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [496,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [497,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [498,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [499,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [500,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [501,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [502,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [503,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [504,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [505,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [506,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [507,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [508,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [509,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [510,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [511,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [512,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [513,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [514,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [515,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [516,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [517,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [518,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [519,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [520,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [521,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [522,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [523,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [524,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [525,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [526,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [527,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [528,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [529,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [530,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [531,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [532,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [533,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [534,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [535,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [536,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [537,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [538,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [539,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [540,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [541,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [542,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [543,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [544,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [545,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [546,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [547,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [548,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [549,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [550,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [551,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [552,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [553,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [554,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [555,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [556,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [557,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [558,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [559,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [560,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [561,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [562,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [563,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [564,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [565,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [566,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [567,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [568,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [569,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [570,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [571,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [572,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [573,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [574,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [575,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [576,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [577,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [578,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [579,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [580,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [581,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [582,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [583,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [584,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [585,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [586,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [587,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [588,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [589,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [590,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [591,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [592,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [593,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [594,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [595,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [596,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [597,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [598,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [599,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + [600,] -2.074669e-16 -8.721366e-17 -1.63251e-16 + diff --git a/tests/testthat/_snaps/predict_counterfactual.md b/tests/testthat/_snaps/predict_counterfactual.md new file mode 100644 index 0000000..dbce9d3 --- /dev/null +++ b/tests/testthat/_snaps/predict_counterfactual.md @@ -0,0 +1,33 @@ +# predict_counterfactual works for guassian + + Code + predict_counterfactual(fit_glm, "treatment") + Output + counter-factual prediction + + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + +# predict_counterfactual works for guassian with lm + + Code + predict_counterfactual(fit_lm, "treatment", data = dummy_data) + Output + counter-factual prediction + + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + +# predict_counterfactual works for binomial + + Code + predict_counterfactual(fit_binom, "treatment") + Output + counter-factual prediction + + pbo trt1 trt2 + 0.3560965 0.5806957 0.6213865 + + diff --git a/tests/testthat/helper-fit.R b/tests/testthat/helper-fit.R new file mode 100644 index 0000000..b7fe21f --- /dev/null +++ b/tests/testthat/helper-fit.R @@ -0,0 +1,8 @@ +# glm fit on iris +fit_glm <- glm(y ~ treatment * s1 + covar, data = dummy_data) + +# lm fit on iris +fit_lm <- lm(y ~ treatment * s1 + covar, data = dummy_data) + +# glm fit on ToothGrowth +fit_binom <- glm(y_b ~ treatment * s1 + covar, data = dummy_data, family = binomial()) diff --git a/tests/testthat/test-bias.R b/tests/testthat/test-bias.R new file mode 100644 index 0000000..83fd119 --- /dev/null +++ b/tests/testthat/test-bias.R @@ -0,0 +1,5 @@ +test_that("bias works for guassian", { + expect_snapshot( + bias(residuals(fit_glm), treatment = dummy_data$treatment, group_idx = list(seq_len(nrow(dummy_data)))) + ) +}) diff --git a/tests/testthat/test-hello.R b/tests/testthat/test-hello.R deleted file mode 100644 index 4dccc45..0000000 --- a/tests/testthat/test-hello.R +++ /dev/null @@ -1,5 +0,0 @@ -test_that("hello greets the entity", { - result <- hello("foo") - expected <- "Hello, Foo" - expect_identical(result, expected) -}) diff --git a/tests/testthat/test-predict_counterfactual.R b/tests/testthat/test-predict_counterfactual.R new file mode 100644 index 0000000..03c3ed6 --- /dev/null +++ b/tests/testthat/test-predict_counterfactual.R @@ -0,0 +1,11 @@ +test_that("predict_counterfactual works for guassian", { + expect_snapshot(predict_counterfactual(fit_glm, "treatment")) +}) + +test_that("predict_counterfactual works for guassian with lm", { + expect_snapshot(predict_counterfactual(fit_lm, "treatment", data = dummy_data)) +}) + +test_that("predict_counterfactual works for binomial", { + expect_snapshot(predict_counterfactual(fit_binom, "treatment")) +}) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R new file mode 100644 index 0000000..ef50a98 --- /dev/null +++ b/tests/testthat/test-utils.R @@ -0,0 +1,34 @@ +test_that("h_get_vars works for single character", { + res <- expect_silent(h_get_vars("abc")) + expect_identical(res, list(treatment = "abc", strata = character(0))) + + expect_error( + h_get_vars(c("abc", "def")), + "length 1 character" + ) + expect_error( + h_get_vars(NULL), + "length 1 character" + ) +}) + +test_that("h_get_vars works for formula", { + res <- expect_silent(h_get_vars(a ~ b + c)) + expect_identical(res, list(treatment = "a", strata = c("b", "c"))) + + res <- expect_silent(h_get_vars(`~`(a, ))) + expect_identical(res, list(treatment = "a", strata = character(0))) + + expect_error( + h_get_vars(log(a) ~ strata(b)), + "left hand side of the treatment formula should be a single name!" + ) + + expect_error( + h_get_vars(a + b ~ strata(b)), + "left hand side of the treatment formula should be a single name!" + ) + + res <- expect_silent(h_get_vars(a ~ strata(b))) + expect_identical(res, list(treatment = "a", strata = "b")) +}) diff --git a/vignettes/hello.Rmd b/vignettes/hello.Rmd deleted file mode 100644 index 4b917c8..0000000 --- a/vignettes/hello.Rmd +++ /dev/null @@ -1,12 +0,0 @@ - ---- -title: "Hello" -date: "`r Sys.Date()`" ---- - -Hello World! - -```{r} -library(RobinCar2) -hello("RobinCar2!") -``` diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd new file mode 100644 index 0000000..ffa063f --- /dev/null +++ b/vignettes/intro.Rmd @@ -0,0 +1,20 @@ + +--- +title: "Introduction" +package: RobinCar2 +output: + rmarkdown::html_vignette: + toc: true +vignette: | + %\VignetteIndexEntry{Model Fitting Algorithm} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +# Brief Introduction of RobinCar2 + +```{r} +library(RobinCar2) +``` From 99a0c84e2c2c8862763882d8ba5d9494ccc0969c Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 19 May 2024 21:41:15 +0000 Subject: [PATCH 05/12] add robincar vcov calc functions and report function --- .Rbuildignore | 3 + R/report_vcovRobinCar.R | 83 ++++++++ R/variance.R | 190 +++++++++++++++++- .../package_structure/developer_note_xgc_v1.R | 1 - 4 files changed, 275 insertions(+), 2 deletions(-) create mode 100644 R/report_vcovRobinCar.R diff --git a/.Rbuildignore b/.Rbuildignore index 5df0095..5b712fe 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -19,3 +19,6 @@ workflows.md images design data-raw +^RobinCar2\.Rcheck$ +^RobinCar2.*\.tar\.gz$ +^RobinCar2.*\.tgz$ diff --git a/R/report_vcovRobinCar.R b/R/report_vcovRobinCar.R new file mode 100644 index 0000000..81a8692 --- /dev/null +++ b/R/report_vcovRobinCar.R @@ -0,0 +1,83 @@ +#' Reporting function for vcov_robincar +#' +#' Pretty reporting for the object returned by \code{\link{vcovRobinCar.prediction_cf}}. +#' Note, if eff_measure was not specified (i.e. NULL), no reporting table will be generated. +#' +#' @param x object returned by \code{\link{vcovRobinCar.prediction_cf}} +#' @param digits integer, number of decimal places of estimates in the reported table +#' +#' @return a data frame +#' @export + +report_vcov_robincar <- function(x, digits = 3) { + eff_measure <- x$eff_measure + + if(is.null(eff_measure)){ + cli::cli_alert_info("When eff_measure is NULL, nothing to report") + }else{ + report( + fit = x$fit, + trt_var_name = x$trt_var_name, + theta = x$theta, + f_vcov = x$f_vcov, + eff_measure = eff_measure, + digits = digits + ) + } + +} + +report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { + n <- nrow(fit$model) + treatments <- fit$xlevels[[trt_var_name]] + n.arm <- summary(fit$model[[trt_var_name]]) + n.arm.perc <- round(n.arm * 100 / n, 1) |> format(nsmall = 1) + + out_cases <- expand.grid(seq_along(theta), seq_along(theta)) + out_cases <- out_cases[out_cases[, 1] != out_cases[, 2], , drop = FALSE] + t1 <- out_cases[, 1] + t0 <- out_cases[, 2] + + outtmp <- data.frame( + Total = n, + Arm_1 = treatments[t1], + N_Arm_1 = paste0(n.arm[t1], "(", n.arm.perc[t1], "%)"), + Arm_0 = treatments[t0], + N_Arm_0 = paste0(n.arm[t0], "(", n.arm.perc[t0], "%)"), + Eff_measure = eff_measure, + Estimate = calculate_f(theta[t1], theta[t0], eff_measure) + ) + outtmp[["S.E."]] <- sqrt(f_vcov[cbind(t1, t0)]) + outtmp[["95%CI_lower"]] <- outtmp$Estimate + qnorm(0.025) * outtmp[["S.E."]] + outtmp[["95%CI_upper"]] <- outtmp$Estimate + qnorm(0.975) * outtmp[["S.E."]] + + if (eff_measure == "diff") { + outtmp["Estimate"] <- round(outtmp["Estimate"], digits) |> format(nsmall = digits) + outtmp["S.E."] <- round(outtmp["S.E."], digits) |> format(nsmall = digits) + outtmp["95%CI_lower"] <- round(outtmp["95%CI_lower"], digits) |> format(nsmall = digits) + outtmp["95%CI_upper"] <- round(outtmp["95%CI_upper"], digits) |> format(nsmall = digits) + } else { + outtmp["S.E."] <- ((exp(outtmp["S.E."]) - 1) * exp(2 * outtmp["Estimate"] + outtmp["S.E."]^2)) |> + round(digits) |> + format(nsmall = digits) + outtmp["Estimate"] <- exp(outtmp["Estimate"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_lower"] <- exp(outtmp["95%CI_lower"]) |> + round(digits) |> + format(nsmall = digits) + outtmp["95%CI_upper"] <- exp(outtmp["95%CI_upper"]) |> + round(digits) |> + format(nsmall = digits) + } + outtmp +} + +calculate_f <- function(theta_t, theta_s, eff_measure) { + switch(eff_measure, + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) + ) +} + diff --git a/R/variance.R b/R/variance.R index be5f633..3d15e66 100644 --- a/R/variance.R +++ b/R/variance.R @@ -1,5 +1,7 @@ +#### Main functions #### + #' @exportS3Method -vcovHC.prediction_cf <- function(x, type = "HC3",...) { +vcovHC.prediction_cf <- function(x, type = "HC3", ...) { fit <- attr(x, "fit") vc <- vcovHC(fit, type = type) mm <- attr(x, "model_matrix") @@ -10,3 +12,189 @@ vcovHC.prediction_cf <- function(x, type = "HC3",...) { dimnames(ret) <- list(names(x), names(x)) ret } + +#' Covariance matrix of Counterfactual prediction +#' +#' @param prediction_cf R object returned by \code{\line{predict_counterfactual}} +#' @param eff_measure default is NULL, the covariance matrix of the mean counterfactual +#' predictions of the response. If not NULL, a valid option is "diff", "odds ratio", or "risk ratio", +#' in which case covariance matrix of the estimated effect measure will also be returned. +#' +#' @return an R object of class 'vcov_robincar', essentially a list consisting +#' \itemize{ +#' \item theta: mean counterfactual prediction of the reponse, as the calculation result of +#' \code{\line{predict_counterfactual}} +#' \item V: covariance matrix of the mean counterfactual predictions of the response +#' \item f_vcov: covariance matrix of the estimated effect measures. It is NULL +#' if \code{eff_measure} is set to NULL (default) +#' \item n: sample size in \code{fit} of \code{\line{predict_counterfactual}} +#' \item eff_measure: \code{eff_measure} in the argument of this function +#' \item fit: \code{fit} of \code{\line{predict_counterfactual}} +#' \item trt_var_name: name of treatment variable in the data used in \code{fit} +#' of \code{\line{predict_counterfactual}} +#' } +#' @export + +vcovRobinCar.prediction_cf <- function(x, eff_measure = NULL, ...) { + target_model <- attr(x, "fit") + tm_use_data <- target_model$data + tm_use_trts <- attr(x, "treatment") + trt_idx_in_data <- sapply(tm_use_data, function(ik) identical(ik, tm_use_trts)) + tm_use_trtvar <- names(trt_idx_in_data)[trt_idx_in_data] + + interim_data <- prepare_interim_vcovrc(target_model, tm_use_trtvar) + + fit.fvcov <- calculate_f_vcov( + fit = target_model, + trt_var_name = tm_use_trtvar, + interim_data = interim_data, + eff_measure = eff_measure + ) + + return(fit.fvcov) +} + + +#### Helper functions #### + +#' helper function: wrap calculations from fit to interested covariance estimate +calculate_f_vcov <- function(fit, trt_var_name, interim_data, eff_measure) { + treatments <- fit$xlevels[[trt_var_name]] + + theta_t <- sapply(interim_data, function(dat_t) unique(dat_t$theta_t)) + + V <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_V", + counterfact_res = interim_data + ) + rownames(V) <- colnames(V) <- treatments + + out <- list( + theta = theta_t, + V = V, + f_vcov = NULL, + n = nrow(fit$model), + eff_measure = NULL, + trt_var_name = trt_var_name, + fit = fit + ) + + if (!is.null(eff_measure)) { + n_fVf <- outer( + X = treatments, + Y = treatments, + FUN = "vectorized_calculate_fVf", + theta = theta_t, + V = V, + n = nrow(fit$model), + eff_measure = eff_measure + ) + rownames(n_fVf) <- colnames(n_fVf) <- treatments + + out$eff_measure <- eff_measure + out$f_vcov <- n_fVf + } + + class(out) <- c(class(out), "vcov_robincar") + out +} + + +#' helper function: prepare interim_data for calculate_f_vcov +prepare_interim_vcovrc <- function(fit, trt_var_name) { + treatments <- fit$xlevels[[trt_var_name]] + use_dat <- fit$model + + dat_counterfact <- sapply(treatments, function(t) { + tmp <- use_dat + tmp$y_observed <- use_dat[, 1] + tmp$trt_observed <- use_dat[[trt_var_name]] + tmp$trt_counterfact <- t + tmp + }, simplify = FALSE) + + out_counterfact <- sapply(treatments, function(t) { + outtmp <- tmp <- dat_counterfact[[t]] + tmp[[trt_var_name]] <- tmp$trt_counterfact + outtmp$mu_t <- predict(fit, newdata = tmp, type = "response") + outtmp$theta_t <- mean(outtmp$mu_t) + outtmp$resid_mu_t <- outtmp[, 1] - outtmp$mu_t + outtmp$resid_theta_t <- outtmp[, 1] - outtmp$mu_t + outtmp + }, simplify = FALSE) + + out_counterfact +} + +#' helper function: calculate V (covariance of the mean counterfactual prediction of response) +calculate_V <- function(t, s, counterfact_res) { + if (t == s) { # : calculate v_tt + # v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + + # 2*cov(Y_j,mu_t_j | j \in I_t) - + # var(mu_t) + calres <- counterfact_res[[t]] + + n <- nrow(calres) + It_ind <- (calres$trt_observed == t) + pi_t <- sum(It_ind) / n + + var_y_mut_It <- var(calres$resid_mu_t[It_ind]) + cov_y_mut_It <- cov(calres$y_observed[It_ind], calres$mu_t[It_ind]) + var_mut <- var(calres$mu_t) + + v_out <- var_y_mut_It / pi_t + 2 * cov_y_mut_It - var_mut + } else { # : calculate v_ts + # v_ts = cov(Y_j,mu_s_j | j \in I_t) + + # cov(Y_j,mu_t_j | j \in I_s) - + # cov(mu_t,mu_s) + + calres_t <- counterfact_res[[t]] + calres_s <- counterfact_res[[s]] + + n <- nrow(calres_t) + It_ind <- (calres_s$trt_observed == t) + Is_ind <- (calres_t$trt_observed == s) + pi_t <- sum(It_ind) / n + pi_s <- sum(Is_ind) / n + + cov_y_mus_It <- cov(calres_s$y_observed[It_ind], calres_s$mu_t[It_ind]) + cov_y_mut_Is <- cov(calres_t$y_observed[Is_ind], calres_t$mu_t[Is_ind]) + cov_mut_mus <- cov(calres_s$mu_t, calres_t$mu_t) + + v_out <- cov_y_mus_It + cov_y_mut_Is - cov_mut_mus + } + v_out +} + +vectorized_calculate_V <- Vectorize(calculate_V, c("t", "s")) + +#' helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure) +op <- function(x, square = FALSE) (x * (1 - x))^ifelse(square, 2, 1) + +odds <- function(x) (x / (1 - x)) + +calculate_fVf <- function(t, s, theta, V, n, eff_measure) { + t_ind <- which(names(theta) == t) + s_ind <- which(names(theta) == s) + + # depneding on eff_measure (1.diff; 2.risk ratio; 3.odds ratio) + # 1. v_tt - 2*v_ts + v_ss + # 2. v_tt/(theta_t^2) - 2*v_ts/(theta_t*theta_s) + v_ss/(theta_s^2) + # 3. v_tt/(op(theta_t)^2) - 2*v_ts/(op(theta_t)*op(theta_s)) + v_ss/(op(theta_s)^2) + if (eff_measure == "diff") { + outtmp <- V[t_ind, t_ind] + V[s_ind, s_ind] - 2 * V[t_ind, s_ind] + } else if (eff_measure == "risk ratio") { + outtmp <- V[t_ind, t_ind] / (theta[t_ind]^2) + + V[s_ind, s_ind] / (theta[s_ind]^2) - + 2 * V[t_ind, s_ind] / (theta[t_ind] * theta[s_ind]) + } else if (eff_measure == "odds ratio") { + outtmp <- V[t_ind, t_ind] / op(theta[t_ind], square = T) + + V[s_ind, s_ind] / op(theta[s_ind], square = T) - + 2 * V[t_ind, s_ind] / (op(theta[t_ind]) * op(theta[s_ind])) + } + outtmp / sqrt(n) +} + +vectorized_calculate_fVf <- Vectorize(calculate_fVf, c("t", "s")) diff --git a/design/package_structure/developer_note_xgc_v1.R b/design/package_structure/developer_note_xgc_v1.R index 0bdbda5..53fd6e0 100644 --- a/design/package_structure/developer_note_xgc_v1.R +++ b/design/package_structure/developer_note_xgc_v1.R @@ -263,7 +263,6 @@ d <- sim_data(500) fit <- glm(y ~ trt:z1 + trt, family = binomial(link = "logit"), data = d) ## start of the pkg -source(file.path(getwd(), "design", "package_structure", "all_draft_functions.R")) pred_counterfact <- get_countfact_pred(fit, "trt") fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact) report_fvcov(fit.fvcov, 3) From 1fc1fd8b6a5bc43e7c008779b6be05961092fef4 Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 19 May 2024 22:09:02 +0000 Subject: [PATCH 06/12] adding all the tests --- R/variance.R | 3 +- tests/testthat/_snaps/report.md | 113 +++++++++ tests/testthat/_snaps/variance.md | 375 ++++++++++++++++++++++++++++++ tests/testthat/test-report.R | 41 ++++ tests/testthat/test-variance.R | 35 +++ 5 files changed, 566 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/_snaps/report.md create mode 100644 tests/testthat/_snaps/variance.md create mode 100644 tests/testthat/test-report.R create mode 100644 tests/testthat/test-variance.R diff --git a/R/variance.R b/R/variance.R index 3d15e66..82ff9bf 100644 --- a/R/variance.R +++ b/R/variance.R @@ -37,7 +37,8 @@ vcovHC.prediction_cf <- function(x, type = "HC3", ...) { vcovRobinCar.prediction_cf <- function(x, eff_measure = NULL, ...) { target_model <- attr(x, "fit") - tm_use_data <- target_model$data + data_idx <- ifelse(identical(class(target_model), "lm"), "model", "data") + tm_use_data <- target_model[[data_idx]] tm_use_trts <- attr(x, "treatment") trt_idx_in_data <- sapply(tm_use_data, function(ik) identical(ik, tm_use_trts)) tm_use_trtvar <- names(trt_idx_in_data)[trt_idx_in_data] diff --git a/tests/testthat/_snaps/report.md b/tests/testthat/_snaps/report.md new file mode 100644 index 0000000..e0cf4de --- /dev/null +++ b/tests/testthat/_snaps/report.md @@ -0,0 +1,113 @@ +# reporting function give info message when eff_measure is NULL + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, + "treatment", dummy_data)), digits = 3) + Message + i When eff_measure is NULL, nothing to report + +# reporting function works for gaussian fitted by lm + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, + "treatment", dummy_data), eff_measure = "diff"), digits = 3) + Output + Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. + 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.564 0.498 + 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.771 0.502 + 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.564 0.498 + 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.207 0.529 + 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.771 0.502 + 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.207 0.529 + 95%CI_lower 95%CI_upper + 1 -0.412 1.539 + 2 -0.213 1.755 + 3 -1.539 0.412 + 4 -0.829 1.244 + 5 -1.755 0.213 + 6 -1.244 0.829 + +# reporting function works for guassian fitted by glm + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, + "treatment", dummy_data), eff_measure = "diff"), digits = 3) + Output + Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. + 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.564 0.498 + 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.771 0.502 + 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.564 0.498 + 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.207 0.529 + 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.771 0.502 + 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.207 0.529 + 95%CI_lower 95%CI_upper + 1 -0.412 1.539 + 2 -0.213 1.755 + 3 -1.539 0.412 + 4 -0.829 1.244 + 5 -1.755 0.213 + 6 -1.244 0.829 + +# reporting function works for binomial fitted by glm + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( + fit_binom, "treatment", dummy_data), eff_measure = "diff"), digits = 3) + Output + Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. + 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.225 0.236 + 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.265 0.236 + 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.225 0.236 + 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.041 0.237 + 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.265 0.236 + 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.041 0.237 + 95%CI_lower 95%CI_upper + 1 -0.238 0.687 + 2 -0.196 0.727 + 3 -0.687 0.238 + 4 -0.425 0.506 + 5 -0.727 0.196 + 6 -0.506 0.425 + +--- + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( + fit_binom, "treatment", dummy_data), eff_measure = "risk ratio"), digits = 3) + Output + Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. + 1 600 trt1 198(33.0%) pbo 202(33.7%) risk ratio 1.631 2.610 + 2 600 trt2 200(33.3%) pbo 202(33.7%) risk ratio 1.745 2.890 + 3 600 pbo 202(33.7%) trt1 198(33.0%) risk ratio 0.613 0.369 + 4 600 trt2 200(33.3%) trt1 198(33.0%) risk ratio 1.070 0.650 + 5 600 pbo 202(33.7%) trt2 200(33.3%) risk ratio 0.573 0.312 + 6 600 trt1 198(33.0%) trt2 200(33.3%) risk ratio 0.935 0.496 + 95%CI_lower 95%CI_upper + 1 0.558 4.762 + 2 0.609 5.001 + 3 0.210 1.791 + 4 0.493 2.324 + 5 0.200 1.642 + 6 0.430 2.030 + +--- + + Code + report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( + fit_binom, "treatment", dummy_data), eff_measure = "odds ratio"), digits = 3) + Output + Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. + 1 600 trt1 198(33.0%) pbo 202(33.7%) odds ratio 2.504 29.088 + 2 600 trt2 200(33.3%) pbo 202(33.7%) odds ratio 2.968 43.252 + 3 600 pbo 202(33.7%) trt1 198(33.0%) odds ratio 0.399 0.740 + 4 600 trt2 200(33.3%) trt1 198(33.0%) odds ratio 1.185 6.372 + 5 600 pbo 202(33.7%) trt2 200(33.3%) odds ratio 0.337 0.558 + 6 600 trt1 198(33.0%) trt2 200(33.3%) odds ratio 0.844 3.231 + 95%CI_lower 95%CI_upper + 1 0.354 17.710 + 2 0.407 21.653 + 3 0.056 2.824 + 4 0.170 8.280 + 5 0.046 2.459 + 6 0.121 5.896 + diff --git a/tests/testthat/_snaps/variance.md b/tests/testthat/_snaps/variance.md new file mode 100644 index 0000000..8f16314 --- /dev/null +++ b/tests/testthat/_snaps/variance.md @@ -0,0 +1,375 @@ +# robincar convariance estimation works for guassian fitted by lm + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, "treatment", + dummy_data)) + Output + $theta + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + $V + pbo trt1 trt2 + pbo 2.74152393 0.05972904 0.05201403 + trt1 0.05972904 3.44436398 0.06333053 + trt2 0.05201403 0.06333053 3.53553897 + + $f_vcov + NULL + + $n + [1] 600 + + $eff_measure + NULL + + $trt_var_name + [1] "treatment" + + $fit + + Call: + lm(formula = y ~ treatment * s1 + covar, data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + 0.33108 0.54454 0.73872 -0.27656 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.20396 0.03796 0.06399 + + + attr(,"class") + [1] "list" "vcov_robincar" + +--- + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, "treatment", + dummy_data), eff_measure = "diff") + Output + $theta + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + $V + pbo trt1 trt2 + pbo 2.74152393 0.05972904 0.05201403 + trt1 0.05972904 3.44436398 0.06333053 + trt2 0.05201403 0.06333053 3.53553897 + + $f_vcov + pbo trt1 trt2 + pbo 0.0000000 0.2476610 0.2520131 + trt1 0.2476610 0.0000000 0.2797824 + trt2 0.2520131 0.2797824 0.0000000 + + $n + [1] 600 + + $eff_measure + [1] "diff" + + $trt_var_name + [1] "treatment" + + $fit + + Call: + lm(formula = y ~ treatment * s1 + covar, data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + 0.33108 0.54454 0.73872 -0.27656 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.20396 0.03796 0.06399 + + + attr(,"class") + [1] "list" "vcov_robincar" + +# robincar convariance estimation works for guassian fitted by glm + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, "treatment", + dummy_data)) + Output + $theta + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + $V + pbo trt1 trt2 + pbo 2.74152393 0.05972904 0.05201403 + trt1 0.05972904 3.44436398 0.06333053 + trt2 0.05201403 0.06333053 3.53553897 + + $f_vcov + NULL + + $n + [1] 600 + + $eff_measure + NULL + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y ~ treatment * s1 + covar, data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + 0.33108 0.54454 0.73872 -0.27656 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.20396 0.03796 0.06399 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 724.3 + Residual Deviance: 632.3 AIC: 1750 + + attr(,"class") + [1] "list" "vcov_robincar" + +--- + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, "treatment", + dummy_data), eff_measure = "diff") + Output + $theta + pbo trt1 trt2 + 0.2003208 0.7639709 0.9712499 + + $V + pbo trt1 trt2 + pbo 2.74152393 0.05972904 0.05201403 + trt1 0.05972904 3.44436398 0.06333053 + trt2 0.05201403 0.06333053 3.53553897 + + $f_vcov + pbo trt1 trt2 + pbo 0.0000000 0.2476610 0.2520131 + trt1 0.2476610 0.0000000 0.2797824 + trt2 0.2520131 0.2797824 0.0000000 + + $n + [1] 600 + + $eff_measure + [1] "diff" + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y ~ treatment * s1 + covar, data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + 0.33108 0.54454 0.73872 -0.27656 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.20396 0.03796 0.06399 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 724.3 + Residual Deviance: 632.3 AIC: 1750 + + attr(,"class") + [1] "list" "vcov_robincar" + +# robincar convariance estimation works for binomial fitted by glm + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", + dummy_data)) + Output + $theta + pbo trt1 trt2 + 0.3560965 0.5806957 0.6213865 + + $V + pbo trt1 trt2 + pbo 0.676245627 0.01113740 0.008003312 + trt1 0.011137402 0.70765788 0.013068673 + trt2 0.008003312 0.01306867 0.698427494 + + $f_vcov + NULL + + $n + [1] 600 + + $eff_measure + NULL + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), + data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + -0.42927 0.97415 1.10707 -0.42043 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.41049 -0.01731 0.06947 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 831 + Residual Deviance: 773.7 AIC: 787.7 + + attr(,"class") + [1] "list" "vcov_robincar" + +--- + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", + dummy_data), eff_measure = "diff") + Output + $theta + pbo trt1 trt2 + 0.3560965 0.5806957 0.6213865 + + $V + pbo trt1 trt2 + pbo 0.676245627 0.01113740 0.008003312 + trt1 0.011137402 0.70765788 0.013068673 + trt2 0.008003312 0.01306867 0.698427494 + + $f_vcov + pbo trt1 trt2 + pbo 0.00000000 0.05558826 0.05546733 + trt1 0.05558826 0.00000000 0.05633614 + trt2 0.05546733 0.05633614 0.00000000 + + $n + [1] 600 + + $eff_measure + [1] "diff" + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), + data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + -0.42927 0.97415 1.10707 -0.42043 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.41049 -0.01731 0.06947 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 831 + Residual Deviance: 773.7 AIC: 787.7 + + attr(,"class") + [1] "list" "vcov_robincar" + +--- + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", + dummy_data), eff_measure = "risk ratio") + Output + $theta + pbo trt1 trt2 + 0.3560965 0.5806957 0.6213865 + + $V + pbo trt1 trt2 + pbo 0.676245627 0.01113740 0.008003312 + trt1 0.011137402 0.70765788 0.013068673 + trt2 0.008003312 0.01306867 0.698427494 + + $f_vcov + pbo trt1 trt2 + pbo 0.0000000 0.2989941 0.2886095 + trt1 0.2989941 0.0000000 0.1565623 + trt2 0.2886095 0.1565623 0.0000000 + + $n + [1] 600 + + $eff_measure + [1] "risk ratio" + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), + data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + -0.42927 0.97415 1.10707 -0.42043 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.41049 -0.01731 0.06947 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 831 + Residual Deviance: 773.7 AIC: 787.7 + + attr(,"class") + [1] "list" "vcov_robincar" + +--- + + Code + vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", + dummy_data), eff_measure = "odds ratio") + Output + $theta + pbo trt1 trt2 + 0.3560965 0.5806957 0.6213865 + + $V + pbo trt1 trt2 + pbo 0.676245627 0.01113740 0.008003312 + trt1 0.011137402 0.70765788 0.013068673 + trt2 0.008003312 0.01306867 0.698427494 + + $f_vcov + pbo trt1 trt2 + pbo 0.0000000 0.9961186 1.0281436 + trt1 0.9961186 0.0000000 0.9838131 + trt2 1.0281436 0.9838131 0.0000000 + + $n + [1] 600 + + $eff_measure + [1] "odds ratio" + + $trt_var_name + [1] "treatment" + + $fit + + Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), + data = dummy_data) + + Coefficients: + (Intercept) treatmenttrt1 treatmenttrt2 s1b + -0.42927 0.97415 1.10707 -0.42043 + covar treatmenttrt1:s1b treatmenttrt2:s1b + 0.41049 -0.01731 0.06947 + + Degrees of Freedom: 599 Total (i.e. Null); 593 Residual + Null Deviance: 831 + Residual Deviance: 773.7 AIC: 787.7 + + attr(,"class") + [1] "list" "vcov_robincar" + diff --git a/tests/testthat/test-report.R b/tests/testthat/test-report.R new file mode 100644 index 0000000..2c7a2f8 --- /dev/null +++ b/tests/testthat/test-report.R @@ -0,0 +1,41 @@ +test_that("reporting function give info message when eff_measure is NULL", { + expect_snapshot( + predict_counterfactual(fit_lm, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf() |> + report_vcov_robincar(digits = 3) + ) +}) + +test_that("reporting function works for gaussian fitted by lm", { + expect_snapshot( + predict_counterfactual(fit_lm, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf(eff_measure = "diff") |> + report_vcov_robincar(digits = 3) + ) +}) + +test_that("reporting function works for guassian fitted by glm", { + expect_snapshot( + predict_counterfactual(fit_glm, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf(eff_measure = "diff") |> + report_vcov_robincar(digits = 3) + ) +}) + +test_that("reporting function works for binomial fitted by glm", { + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf(eff_measure = "diff") |> + report_vcov_robincar(digits = 3) + ) + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf(eff_measure = "risk ratio") |> + report_vcov_robincar(digits = 3) + ) + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> + vcovRobinCar.prediction_cf(eff_measure = "odds ratio") |> + report_vcov_robincar(digits = 3) + ) +}) diff --git a/tests/testthat/test-variance.R b/tests/testthat/test-variance.R new file mode 100644 index 0000000..41469d9 --- /dev/null +++ b/tests/testthat/test-variance.R @@ -0,0 +1,35 @@ +test_that("robincar convariance estimation works for guassian fitted by lm", { + expect_snapshot( + predict_counterfactual(fit_lm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() + ) + + expect_snapshot( + predict_counterfactual(fit_lm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") + ) +}) + +test_that("robincar convariance estimation works for guassian fitted by glm", { + expect_snapshot( + predict_counterfactual(fit_glm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() + ) + + expect_snapshot( + predict_counterfactual(fit_glm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") + ) +}) + +test_that("robincar convariance estimation works for binomial fitted by glm", { + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() + ) + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") + ) + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "risk ratio") + ) + expect_snapshot( + predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "odds ratio") + ) +}) + From ac97b374e9ab388cffe522b29a06c90eab6d32c5 Mon Sep 17 00:00:00 2001 From: "27856297+dependabot-preview[bot]@users.noreply.github.com" <27856297+dependabot-preview[bot]@users.noreply.github.com> Date: Sun, 19 May 2024 22:18:59 +0000 Subject: [PATCH 07/12] [skip roxygen] [skip vbump] Roxygen Man Pages Auto Update --- man/calculate_V.Rd | 11 +++++++++++ man/calculate_f_vcov.Rd | 11 +++++++++++ man/op.Rd | 11 +++++++++++ man/prepare_interim_vcovrc.Rd | 11 +++++++++++ man/report_vcov_robincar.Rd | 20 +++++++++++++++++++ man/vcovRobinCar.prediction_cf.Rd | 33 +++++++++++++++++++++++++++++++ 6 files changed, 97 insertions(+) create mode 100644 man/calculate_V.Rd create mode 100644 man/calculate_f_vcov.Rd create mode 100644 man/op.Rd create mode 100644 man/prepare_interim_vcovrc.Rd create mode 100644 man/report_vcov_robincar.Rd create mode 100644 man/vcovRobinCar.prediction_cf.Rd diff --git a/man/calculate_V.Rd b/man/calculate_V.Rd new file mode 100644 index 0000000..5f08c36 --- /dev/null +++ b/man/calculate_V.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{calculate_V} +\alias{calculate_V} +\title{helper function: calculate V (covariance of the mean counterfactual prediction of response)} +\usage{ +calculate_V(t, s, counterfact_res) +} +\description{ +helper function: calculate V (covariance of the mean counterfactual prediction of response) +} diff --git a/man/calculate_f_vcov.Rd b/man/calculate_f_vcov.Rd new file mode 100644 index 0000000..63f0a0d --- /dev/null +++ b/man/calculate_f_vcov.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{calculate_f_vcov} +\alias{calculate_f_vcov} +\title{helper function: wrap calculations from fit to interested covariance estimate} +\usage{ +calculate_f_vcov(fit, trt_var_name, interim_data, eff_measure) +} +\description{ +helper function: wrap calculations from fit to interested covariance estimate +} diff --git a/man/op.Rd b/man/op.Rd new file mode 100644 index 0000000..c2c6202 --- /dev/null +++ b/man/op.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{op} +\alias{op} +\title{helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure)} +\usage{ +op(x, square = FALSE) +} +\description{ +helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure) +} diff --git a/man/prepare_interim_vcovrc.Rd b/man/prepare_interim_vcovrc.Rd new file mode 100644 index 0000000..ba2b275 --- /dev/null +++ b/man/prepare_interim_vcovrc.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{prepare_interim_vcovrc} +\alias{prepare_interim_vcovrc} +\title{helper function: prepare interim_data for calculate_f_vcov} +\usage{ +prepare_interim_vcovrc(fit, trt_var_name) +} +\description{ +helper function: prepare interim_data for calculate_f_vcov +} diff --git a/man/report_vcov_robincar.Rd b/man/report_vcov_robincar.Rd new file mode 100644 index 0000000..1176ce5 --- /dev/null +++ b/man/report_vcov_robincar.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/report_vcovRobinCar.R +\name{report_vcov_robincar} +\alias{report_vcov_robincar} +\title{Reporting function for vcov_robincar} +\usage{ +report_vcov_robincar(x, digits = 3) +} +\arguments{ +\item{x}{object returned by \code{\link{vcovRobinCar.prediction_cf}}} + +\item{digits}{integer, number of decimal places of estimates in the reported table} +} +\value{ +a data frame +} +\description{ +Pretty reporting for the object returned by \code{\link{vcovRobinCar.prediction_cf}}. +Note, if eff_measure was not specified (i.e. NULL), no reporting table will be generated. +} diff --git a/man/vcovRobinCar.prediction_cf.Rd b/man/vcovRobinCar.prediction_cf.Rd new file mode 100644 index 0000000..5ba4d88 --- /dev/null +++ b/man/vcovRobinCar.prediction_cf.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{vcovRobinCar.prediction_cf} +\alias{vcovRobinCar.prediction_cf} +\title{Covariance matrix of Counterfactual prediction} +\usage{ +vcovRobinCar.prediction_cf(x, eff_measure = NULL, ...) +} +\arguments{ +\item{eff_measure}{default is NULL, the covariance matrix of the mean counterfactual +predictions of the response. If not NULL, a valid option is "diff", "odds ratio", or "risk ratio", +in which case covariance matrix of the estimated effect measure will also be returned.} + +\item{prediction_cf}{R object returned by \code{\line{predict_counterfactual}}} +} +\value{ +an R object of class 'vcov_robincar', essentially a list consisting +\itemize{ +\item theta: mean counterfactual prediction of the reponse, as the calculation result of +\code{\line{predict_counterfactual}} +\item V: covariance matrix of the mean counterfactual predictions of the response +\item f_vcov: covariance matrix of the estimated effect measures. It is NULL +if \code{eff_measure} is set to NULL (default) +\item n: sample size in \code{fit} of \code{\line{predict_counterfactual}} +\item eff_measure: \code{eff_measure} in the argument of this function +\item fit: \code{fit} of \code{\line{predict_counterfactual}} +\item trt_var_name: name of treatment variable in the data used in \code{fit} +of \code{\line{predict_counterfactual}} +} +} +\description{ +Covariance matrix of Counterfactual prediction +} From 4acaaad1f8cc68ef4acaf0ba2352f6973b7458a1 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Sun, 19 May 2024 22:19:11 +0000 Subject: [PATCH 08/12] [skip style] [skip vbump] Restyle files --- R/report_vcovRobinCar.R | 12 +++++------- tests/testthat/test-variance.R | 1 - 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/R/report_vcovRobinCar.R b/R/report_vcovRobinCar.R index 81a8692..f810e4c 100644 --- a/R/report_vcovRobinCar.R +++ b/R/report_vcovRobinCar.R @@ -12,9 +12,9 @@ report_vcov_robincar <- function(x, digits = 3) { eff_measure <- x$eff_measure - if(is.null(eff_measure)){ + if (is.null(eff_measure)) { cli::cli_alert_info("When eff_measure is NULL, nothing to report") - }else{ + } else { report( fit = x$fit, trt_var_name = x$trt_var_name, @@ -24,7 +24,6 @@ report_vcov_robincar <- function(x, digits = 3) { digits = digits ) } - } report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { @@ -75,9 +74,8 @@ report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { calculate_f <- function(theta_t, theta_s, eff_measure) { switch(eff_measure, - "diff" = theta_t - theta_s, - "risk ratio" = log(theta_t / theta_s), - "odds ratio" = log(odds(theta_t) / odds(theta_s)) + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) ) } - diff --git a/tests/testthat/test-variance.R b/tests/testthat/test-variance.R index 41469d9..b7a4014 100644 --- a/tests/testthat/test-variance.R +++ b/tests/testthat/test-variance.R @@ -32,4 +32,3 @@ test_that("robincar convariance estimation works for binomial fitted by glm", { predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "odds ratio") ) }) - From a8834eb5d8ebe14b306434b8d3e0f119e34935a7 Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 19 May 2024 22:27:40 +0000 Subject: [PATCH 09/12] add Rds --- NAMESPACE | 2 ++ R/report_vcovRobinCar.R | 13 ++++++------ R/variance.R | 4 ++++ man/report_vcov_robincar.Rd | 20 +++++++++++++++++++ man/vcovRobinCar.prediction_cf.Rd | 33 +++++++++++++++++++++++++++++++ 5 files changed, 65 insertions(+), 7 deletions(-) create mode 100644 man/report_vcov_robincar.Rd create mode 100644 man/vcovRobinCar.prediction_cf.Rd diff --git a/NAMESPACE b/NAMESPACE index 55f3734..e0729b2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,6 +6,8 @@ S3method(print,prediction_cf) S3method(vcovHC,prediction_cf) export(bias) export(predict_counterfactual) +export(report_vcov_robincar) +export(vcovRobinCar.prediction_cf) import(checkmate) importFrom(sandwich,vcovHC) importFrom(stats,coefficients) diff --git a/R/report_vcovRobinCar.R b/R/report_vcovRobinCar.R index 81a8692..617ead7 100644 --- a/R/report_vcovRobinCar.R +++ b/R/report_vcovRobinCar.R @@ -12,9 +12,9 @@ report_vcov_robincar <- function(x, digits = 3) { eff_measure <- x$eff_measure - if(is.null(eff_measure)){ + if (is.null(eff_measure)) { cli::cli_alert_info("When eff_measure is NULL, nothing to report") - }else{ + } else { report( fit = x$fit, trt_var_name = x$trt_var_name, @@ -24,7 +24,6 @@ report_vcov_robincar <- function(x, digits = 3) { digits = digits ) } - } report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { @@ -73,11 +72,11 @@ report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { outtmp } +#' @noRd calculate_f <- function(theta_t, theta_s, eff_measure) { switch(eff_measure, - "diff" = theta_t - theta_s, - "risk ratio" = log(theta_t / theta_s), - "odds ratio" = log(odds(theta_t) / odds(theta_s)) + "diff" = theta_t - theta_s, + "risk ratio" = log(theta_t / theta_s), + "odds ratio" = log(odds(theta_t) / odds(theta_s)) ) } - diff --git a/R/variance.R b/R/variance.R index 82ff9bf..50331a0 100644 --- a/R/variance.R +++ b/R/variance.R @@ -59,6 +59,7 @@ vcovRobinCar.prediction_cf <- function(x, eff_measure = NULL, ...) { #### Helper functions #### #' helper function: wrap calculations from fit to interested covariance estimate +#' @noRd calculate_f_vcov <- function(fit, trt_var_name, interim_data, eff_measure) { treatments <- fit$xlevels[[trt_var_name]] @@ -104,6 +105,7 @@ calculate_f_vcov <- function(fit, trt_var_name, interim_data, eff_measure) { #' helper function: prepare interim_data for calculate_f_vcov +#' @noRd prepare_interim_vcovrc <- function(fit, trt_var_name) { treatments <- fit$xlevels[[trt_var_name]] use_dat <- fit$model @@ -130,6 +132,7 @@ prepare_interim_vcovrc <- function(fit, trt_var_name) { } #' helper function: calculate V (covariance of the mean counterfactual prediction of response) +#' @noRd calculate_V <- function(t, s, counterfact_res) { if (t == s) { # : calculate v_tt # v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + @@ -172,6 +175,7 @@ calculate_V <- function(t, s, counterfact_res) { vectorized_calculate_V <- Vectorize(calculate_V, c("t", "s")) #' helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure) +#' @noRd op <- function(x, square = FALSE) (x * (1 - x))^ifelse(square, 2, 1) odds <- function(x) (x / (1 - x)) diff --git a/man/report_vcov_robincar.Rd b/man/report_vcov_robincar.Rd new file mode 100644 index 0000000..1176ce5 --- /dev/null +++ b/man/report_vcov_robincar.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/report_vcovRobinCar.R +\name{report_vcov_robincar} +\alias{report_vcov_robincar} +\title{Reporting function for vcov_robincar} +\usage{ +report_vcov_robincar(x, digits = 3) +} +\arguments{ +\item{x}{object returned by \code{\link{vcovRobinCar.prediction_cf}}} + +\item{digits}{integer, number of decimal places of estimates in the reported table} +} +\value{ +a data frame +} +\description{ +Pretty reporting for the object returned by \code{\link{vcovRobinCar.prediction_cf}}. +Note, if eff_measure was not specified (i.e. NULL), no reporting table will be generated. +} diff --git a/man/vcovRobinCar.prediction_cf.Rd b/man/vcovRobinCar.prediction_cf.Rd new file mode 100644 index 0000000..5ba4d88 --- /dev/null +++ b/man/vcovRobinCar.prediction_cf.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance.R +\name{vcovRobinCar.prediction_cf} +\alias{vcovRobinCar.prediction_cf} +\title{Covariance matrix of Counterfactual prediction} +\usage{ +vcovRobinCar.prediction_cf(x, eff_measure = NULL, ...) +} +\arguments{ +\item{eff_measure}{default is NULL, the covariance matrix of the mean counterfactual +predictions of the response. If not NULL, a valid option is "diff", "odds ratio", or "risk ratio", +in which case covariance matrix of the estimated effect measure will also be returned.} + +\item{prediction_cf}{R object returned by \code{\line{predict_counterfactual}}} +} +\value{ +an R object of class 'vcov_robincar', essentially a list consisting +\itemize{ +\item theta: mean counterfactual prediction of the reponse, as the calculation result of +\code{\line{predict_counterfactual}} +\item V: covariance matrix of the mean counterfactual predictions of the response +\item f_vcov: covariance matrix of the estimated effect measures. It is NULL +if \code{eff_measure} is set to NULL (default) +\item n: sample size in \code{fit} of \code{\line{predict_counterfactual}} +\item eff_measure: \code{eff_measure} in the argument of this function +\item fit: \code{fit} of \code{\line{predict_counterfactual}} +\item trt_var_name: name of treatment variable in the data used in \code{fit} +of \code{\line{predict_counterfactual}} +} +} +\description{ +Covariance matrix of Counterfactual prediction +} From cbb8f70e8717909863ec5b55108ca6985b2d29ff Mon Sep 17 00:00:00 2001 From: Gregory Chen Date: Sun, 19 May 2024 22:29:33 +0000 Subject: [PATCH 10/12] delete unneeded RdS --- man/calculate_V.Rd | 11 ----------- man/calculate_f_vcov.Rd | 11 ----------- man/op.Rd | 11 ----------- man/prepare_interim_vcovrc.Rd | 11 ----------- 4 files changed, 44 deletions(-) delete mode 100644 man/calculate_V.Rd delete mode 100644 man/calculate_f_vcov.Rd delete mode 100644 man/op.Rd delete mode 100644 man/prepare_interim_vcovrc.Rd diff --git a/man/calculate_V.Rd b/man/calculate_V.Rd deleted file mode 100644 index 5f08c36..0000000 --- a/man/calculate_V.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/variance.R -\name{calculate_V} -\alias{calculate_V} -\title{helper function: calculate V (covariance of the mean counterfactual prediction of response)} -\usage{ -calculate_V(t, s, counterfact_res) -} -\description{ -helper function: calculate V (covariance of the mean counterfactual prediction of response) -} diff --git a/man/calculate_f_vcov.Rd b/man/calculate_f_vcov.Rd deleted file mode 100644 index 63f0a0d..0000000 --- a/man/calculate_f_vcov.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/variance.R -\name{calculate_f_vcov} -\alias{calculate_f_vcov} -\title{helper function: wrap calculations from fit to interested covariance estimate} -\usage{ -calculate_f_vcov(fit, trt_var_name, interim_data, eff_measure) -} -\description{ -helper function: wrap calculations from fit to interested covariance estimate -} diff --git a/man/op.Rd b/man/op.Rd deleted file mode 100644 index c2c6202..0000000 --- a/man/op.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/variance.R -\name{op} -\alias{op} -\title{helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure)} -\usage{ -op(x, square = FALSE) -} -\description{ -helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure) -} diff --git a/man/prepare_interim_vcovrc.Rd b/man/prepare_interim_vcovrc.Rd deleted file mode 100644 index ba2b275..0000000 --- a/man/prepare_interim_vcovrc.Rd +++ /dev/null @@ -1,11 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/variance.R -\name{prepare_interim_vcovrc} -\alias{prepare_interim_vcovrc} -\title{helper function: prepare interim_data for calculate_f_vcov} -\usage{ -prepare_interim_vcovrc(fit, trt_var_name) -} -\description{ -helper function: prepare interim_data for calculate_f_vcov -} From a4f6519f2aebbbdc991c73c043e215f6f4010d41 Mon Sep 17 00:00:00 2001 From: Liming <36079400+clarkliming@users.noreply.github.com> Date: Thu, 30 May 2024 19:52:04 +0800 Subject: [PATCH 11/12] add anhecova covariance (#24) * add anhecova covariance * fix error --- NAMESPACE | 1 + R/variance.R | 2 +- R/variance_anhecova.R | 75 +++++++++++++++++++++++++++++++++++++++++++ man/vcovAHECOVA.Rd | 23 +++++++++++++ 4 files changed, 100 insertions(+), 1 deletion(-) create mode 100644 R/variance_anhecova.R create mode 100644 man/vcovAHECOVA.Rd diff --git a/NAMESPACE b/NAMESPACE index e0729b2..6491f4e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ S3method(vcovHC,prediction_cf) export(bias) export(predict_counterfactual) export(report_vcov_robincar) +export(vcovAHECOVA) export(vcovRobinCar.prediction_cf) import(checkmate) importFrom(sandwich,vcovHC) diff --git a/R/variance.R b/R/variance.R index 50331a0..49f3d37 100644 --- a/R/variance.R +++ b/R/variance.R @@ -5,7 +5,7 @@ vcovHC.prediction_cf <- function(x, type = "HC3", ...) { fit <- attr(x, "fit") vc <- vcovHC(fit, type = type) mm <- attr(x, "model_matrix") - n <- length(md) / length(names(x)) + n <- nrow(mm) / length(names(x)) md <- family(fit)$mu.eta(attr(x, "predictions_linear")) / n z <- block_sum(as.vector(md) * mm, n) ret <- z %*% vc %*% t(z) diff --git a/R/variance_anhecova.R b/R/variance_anhecova.R new file mode 100644 index 0000000..051584a --- /dev/null +++ b/R/variance_anhecova.R @@ -0,0 +1,75 @@ +#' ANHECOVA Covariance +#' +#' @param x (`prediction_cf`) Counter-factual prediction. +#' @param decompose (`flag`) whether to use decompose method to calculate the variance. +#' @param randomization (`string`) randomization method. +#' @param ... Not used. +#' +#' @return Named covariance matrix. +#' @export +vcovAHECOVA <- function(x, decompose = TRUE, randomization = "simple", ...) { + assert_class(x, "prediction_cf") + assert_flag(decompose) + assert_string(randomization) + resi <- attr(x, "residual") + est <- as.numeric(x) + preds <- attr(x, "predictions") + var_preds <- var(preds) + y <- attr(x, "response") + trt <- attr(x, "treatment") + pi_t <- as.numeric(table(trt) / length(trt)) + trt_lvls <- levels(trt) + group_idx <- attr(x, "group_idx") + + idx <- split(seq_len(length(trt)), trt) + cov_ymu <- vapply(idx, function(is) stats::cov(y[is], preds[is, ]), FUN.VALUE = rep(0, ncol(preds))) + + if (decompose) { + vcov_sr <- (vapply(idx, function(is) stats::var(y[is]), FUN.VALUE = 0) + diag(var_preds) - 2 * diag(cov_ymu)) / pi_t + } else { + vcov_sr <- vapply(idx, function(is) stats::var(resi[is]), FUN.VALUE = 0) / pi_t + } + + v <- diag(vcov_sr) + cov_ymu + t(cov_ymu) - var_preds + v <- v - h_get_erb(resi, group_idx, trt, pi_t, randomization) + ret <- v / length(resi) + dimnames(ret) <- list(trt_lvls, trt_lvls) + return(ret) +} + +#' @keywords internal +h_adjust_pi <- function(pi_t) { + diag(pi_t) - pi_t %*% t(pi_t) +} + +#' @keywords internal +h_get_erb <- function(resi, group_idx, trt, pi_t, randomization) { + if (length(group_idx) == 1L) { + return(0) + } + if (randomization %in% c("simple", "pocock-simon")) { + return(0) + } + omegaz_sr <- h_adjust_pi(pi_t) + resi_per_strata <- vapply( + group_idx, + function(ii) vapply(split(resi[ii], trt[ii]), mean, FUN.VALUE = 0), + FUN.VALUE = rep(0, length(levels(trt))) + ) + strata_props <- vapply( + group_idx, + length, + FUN.VALUE = 0L + ) + strata_props <- strata_props / sum(strata_props) + rb_z <- resi_per_strata / as.numeric(pi_t) + strata_levels <- length(resi_per_strata) + n_trt <- length(pi_t) + ind <- matrix(seq_len(strata_levels), byrow = TRUE, ncol = n_trt) + rb_z_sum <- lapply( + seq_len(ncol(rb_z)), + function(x) rb_z[, x] %*% t(rb_z[, x]) * strata_props[x] + ) + rb_z_sum <- Reduce(`+`, rb_z_sum) + rb_z_sum * omegaz_sr +} diff --git a/man/vcovAHECOVA.Rd b/man/vcovAHECOVA.Rd new file mode 100644 index 0000000..e4c8ec5 --- /dev/null +++ b/man/vcovAHECOVA.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{vcovAHECOVA} +\alias{vcovAHECOVA} +\title{ANHECOVA Covariance} +\usage{ +vcovAHECOVA(x, decompose = TRUE, randomization = "simple", ...) +} +\arguments{ +\item{x}{(\code{prediction_cf}) Counter-factual prediction.} + +\item{decompose}{(\code{flag}) whether to use decompose method to calculate the variance.} + +\item{randomization}{(\code{string}) randomization method.} + +\item{...}{Not used.} +} +\value{ +Named covariance matrix. +} +\description{ +ANHECOVA Covariance +} From 033de5a01d19b91aa170a019754225940d56d937 Mon Sep 17 00:00:00 2001 From: Liming <36079400+clarkliming@users.noreply.github.com> Date: Wed, 26 Jun 2024 08:50:33 +0800 Subject: [PATCH 12/12] use new implementation of anhecova variance (#25) * replace original function with faster version * update anhecova variance and remove report functions --- NAMESPACE | 4 +- R/report_vcovRobinCar.R | 82 ---- R/variance.R | 205 ---------- R/variance_anhecova.R | 50 ++- R/variance_hc.R | 14 + .../design_note.R} | 0 .../developer_note.R} | 0 man/h_adjust_pi.Rd | 18 + man/h_get_erb.Rd | 23 ++ man/report_vcov_robincar.Rd | 20 - man/{vcovAHECOVA.Rd => vcovANHECOVA.Rd} | 6 +- man/vcovRobinCar.prediction_cf.Rd | 33 -- tests/testthat/_snaps/report.md | 113 ------ tests/testthat/_snaps/variance.md | 379 +----------------- tests/testthat/test-report.R | 41 -- tests/testthat/test-variance.R | 32 +- 16 files changed, 111 insertions(+), 909 deletions(-) delete mode 100644 R/report_vcovRobinCar.R delete mode 100644 R/variance.R create mode 100644 R/variance_hc.R rename design/{package_structure/design_note_xgc_v1.R => code_replicate/design_note.R} (100%) rename design/{package_structure/developer_note_xgc_v1.R => code_replicate/developer_note.R} (100%) create mode 100644 man/h_adjust_pi.Rd create mode 100644 man/h_get_erb.Rd delete mode 100644 man/report_vcov_robincar.Rd rename man/{vcovAHECOVA.Rd => vcovANHECOVA.Rd} (81%) delete mode 100644 man/vcovRobinCar.prediction_cf.Rd delete mode 100644 tests/testthat/_snaps/report.md delete mode 100644 tests/testthat/test-report.R diff --git a/NAMESPACE b/NAMESPACE index 6491f4e..2fce319 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,9 +6,7 @@ S3method(print,prediction_cf) S3method(vcovHC,prediction_cf) export(bias) export(predict_counterfactual) -export(report_vcov_robincar) -export(vcovAHECOVA) -export(vcovRobinCar.prediction_cf) +export(vcovANHECOVA) import(checkmate) importFrom(sandwich,vcovHC) importFrom(stats,coefficients) diff --git a/R/report_vcovRobinCar.R b/R/report_vcovRobinCar.R deleted file mode 100644 index 617ead7..0000000 --- a/R/report_vcovRobinCar.R +++ /dev/null @@ -1,82 +0,0 @@ -#' Reporting function for vcov_robincar -#' -#' Pretty reporting for the object returned by \code{\link{vcovRobinCar.prediction_cf}}. -#' Note, if eff_measure was not specified (i.e. NULL), no reporting table will be generated. -#' -#' @param x object returned by \code{\link{vcovRobinCar.prediction_cf}} -#' @param digits integer, number of decimal places of estimates in the reported table -#' -#' @return a data frame -#' @export - -report_vcov_robincar <- function(x, digits = 3) { - eff_measure <- x$eff_measure - - if (is.null(eff_measure)) { - cli::cli_alert_info("When eff_measure is NULL, nothing to report") - } else { - report( - fit = x$fit, - trt_var_name = x$trt_var_name, - theta = x$theta, - f_vcov = x$f_vcov, - eff_measure = eff_measure, - digits = digits - ) - } -} - -report <- function(fit, trt_var_name, theta, f_vcov, eff_measure, digits = 3) { - n <- nrow(fit$model) - treatments <- fit$xlevels[[trt_var_name]] - n.arm <- summary(fit$model[[trt_var_name]]) - n.arm.perc <- round(n.arm * 100 / n, 1) |> format(nsmall = 1) - - out_cases <- expand.grid(seq_along(theta), seq_along(theta)) - out_cases <- out_cases[out_cases[, 1] != out_cases[, 2], , drop = FALSE] - t1 <- out_cases[, 1] - t0 <- out_cases[, 2] - - outtmp <- data.frame( - Total = n, - Arm_1 = treatments[t1], - N_Arm_1 = paste0(n.arm[t1], "(", n.arm.perc[t1], "%)"), - Arm_0 = treatments[t0], - N_Arm_0 = paste0(n.arm[t0], "(", n.arm.perc[t0], "%)"), - Eff_measure = eff_measure, - Estimate = calculate_f(theta[t1], theta[t0], eff_measure) - ) - outtmp[["S.E."]] <- sqrt(f_vcov[cbind(t1, t0)]) - outtmp[["95%CI_lower"]] <- outtmp$Estimate + qnorm(0.025) * outtmp[["S.E."]] - outtmp[["95%CI_upper"]] <- outtmp$Estimate + qnorm(0.975) * outtmp[["S.E."]] - - if (eff_measure == "diff") { - outtmp["Estimate"] <- round(outtmp["Estimate"], digits) |> format(nsmall = digits) - outtmp["S.E."] <- round(outtmp["S.E."], digits) |> format(nsmall = digits) - outtmp["95%CI_lower"] <- round(outtmp["95%CI_lower"], digits) |> format(nsmall = digits) - outtmp["95%CI_upper"] <- round(outtmp["95%CI_upper"], digits) |> format(nsmall = digits) - } else { - outtmp["S.E."] <- ((exp(outtmp["S.E."]) - 1) * exp(2 * outtmp["Estimate"] + outtmp["S.E."]^2)) |> - round(digits) |> - format(nsmall = digits) - outtmp["Estimate"] <- exp(outtmp["Estimate"]) |> - round(digits) |> - format(nsmall = digits) - outtmp["95%CI_lower"] <- exp(outtmp["95%CI_lower"]) |> - round(digits) |> - format(nsmall = digits) - outtmp["95%CI_upper"] <- exp(outtmp["95%CI_upper"]) |> - round(digits) |> - format(nsmall = digits) - } - outtmp -} - -#' @noRd -calculate_f <- function(theta_t, theta_s, eff_measure) { - switch(eff_measure, - "diff" = theta_t - theta_s, - "risk ratio" = log(theta_t / theta_s), - "odds ratio" = log(odds(theta_t) / odds(theta_s)) - ) -} diff --git a/R/variance.R b/R/variance.R deleted file mode 100644 index 49f3d37..0000000 --- a/R/variance.R +++ /dev/null @@ -1,205 +0,0 @@ -#### Main functions #### - -#' @exportS3Method -vcovHC.prediction_cf <- function(x, type = "HC3", ...) { - fit <- attr(x, "fit") - vc <- vcovHC(fit, type = type) - mm <- attr(x, "model_matrix") - n <- nrow(mm) / length(names(x)) - md <- family(fit)$mu.eta(attr(x, "predictions_linear")) / n - z <- block_sum(as.vector(md) * mm, n) - ret <- z %*% vc %*% t(z) - dimnames(ret) <- list(names(x), names(x)) - ret -} - -#' Covariance matrix of Counterfactual prediction -#' -#' @param prediction_cf R object returned by \code{\line{predict_counterfactual}} -#' @param eff_measure default is NULL, the covariance matrix of the mean counterfactual -#' predictions of the response. If not NULL, a valid option is "diff", "odds ratio", or "risk ratio", -#' in which case covariance matrix of the estimated effect measure will also be returned. -#' -#' @return an R object of class 'vcov_robincar', essentially a list consisting -#' \itemize{ -#' \item theta: mean counterfactual prediction of the reponse, as the calculation result of -#' \code{\line{predict_counterfactual}} -#' \item V: covariance matrix of the mean counterfactual predictions of the response -#' \item f_vcov: covariance matrix of the estimated effect measures. It is NULL -#' if \code{eff_measure} is set to NULL (default) -#' \item n: sample size in \code{fit} of \code{\line{predict_counterfactual}} -#' \item eff_measure: \code{eff_measure} in the argument of this function -#' \item fit: \code{fit} of \code{\line{predict_counterfactual}} -#' \item trt_var_name: name of treatment variable in the data used in \code{fit} -#' of \code{\line{predict_counterfactual}} -#' } -#' @export - -vcovRobinCar.prediction_cf <- function(x, eff_measure = NULL, ...) { - target_model <- attr(x, "fit") - data_idx <- ifelse(identical(class(target_model), "lm"), "model", "data") - tm_use_data <- target_model[[data_idx]] - tm_use_trts <- attr(x, "treatment") - trt_idx_in_data <- sapply(tm_use_data, function(ik) identical(ik, tm_use_trts)) - tm_use_trtvar <- names(trt_idx_in_data)[trt_idx_in_data] - - interim_data <- prepare_interim_vcovrc(target_model, tm_use_trtvar) - - fit.fvcov <- calculate_f_vcov( - fit = target_model, - trt_var_name = tm_use_trtvar, - interim_data = interim_data, - eff_measure = eff_measure - ) - - return(fit.fvcov) -} - - -#### Helper functions #### - -#' helper function: wrap calculations from fit to interested covariance estimate -#' @noRd -calculate_f_vcov <- function(fit, trt_var_name, interim_data, eff_measure) { - treatments <- fit$xlevels[[trt_var_name]] - - theta_t <- sapply(interim_data, function(dat_t) unique(dat_t$theta_t)) - - V <- outer( - X = treatments, - Y = treatments, - FUN = "vectorized_calculate_V", - counterfact_res = interim_data - ) - rownames(V) <- colnames(V) <- treatments - - out <- list( - theta = theta_t, - V = V, - f_vcov = NULL, - n = nrow(fit$model), - eff_measure = NULL, - trt_var_name = trt_var_name, - fit = fit - ) - - if (!is.null(eff_measure)) { - n_fVf <- outer( - X = treatments, - Y = treatments, - FUN = "vectorized_calculate_fVf", - theta = theta_t, - V = V, - n = nrow(fit$model), - eff_measure = eff_measure - ) - rownames(n_fVf) <- colnames(n_fVf) <- treatments - - out$eff_measure <- eff_measure - out$f_vcov <- n_fVf - } - - class(out) <- c(class(out), "vcov_robincar") - out -} - - -#' helper function: prepare interim_data for calculate_f_vcov -#' @noRd -prepare_interim_vcovrc <- function(fit, trt_var_name) { - treatments <- fit$xlevels[[trt_var_name]] - use_dat <- fit$model - - dat_counterfact <- sapply(treatments, function(t) { - tmp <- use_dat - tmp$y_observed <- use_dat[, 1] - tmp$trt_observed <- use_dat[[trt_var_name]] - tmp$trt_counterfact <- t - tmp - }, simplify = FALSE) - - out_counterfact <- sapply(treatments, function(t) { - outtmp <- tmp <- dat_counterfact[[t]] - tmp[[trt_var_name]] <- tmp$trt_counterfact - outtmp$mu_t <- predict(fit, newdata = tmp, type = "response") - outtmp$theta_t <- mean(outtmp$mu_t) - outtmp$resid_mu_t <- outtmp[, 1] - outtmp$mu_t - outtmp$resid_theta_t <- outtmp[, 1] - outtmp$mu_t - outtmp - }, simplify = FALSE) - - out_counterfact -} - -#' helper function: calculate V (covariance of the mean counterfactual prediction of response) -#' @noRd -calculate_V <- function(t, s, counterfact_res) { - if (t == s) { # : calculate v_tt - # v_tt = pi_t^(-1) * var(Y_j - mu_t_j | j \in I_t) + - # 2*cov(Y_j,mu_t_j | j \in I_t) - - # var(mu_t) - calres <- counterfact_res[[t]] - - n <- nrow(calres) - It_ind <- (calres$trt_observed == t) - pi_t <- sum(It_ind) / n - - var_y_mut_It <- var(calres$resid_mu_t[It_ind]) - cov_y_mut_It <- cov(calres$y_observed[It_ind], calres$mu_t[It_ind]) - var_mut <- var(calres$mu_t) - - v_out <- var_y_mut_It / pi_t + 2 * cov_y_mut_It - var_mut - } else { # : calculate v_ts - # v_ts = cov(Y_j,mu_s_j | j \in I_t) + - # cov(Y_j,mu_t_j | j \in I_s) - - # cov(mu_t,mu_s) - - calres_t <- counterfact_res[[t]] - calres_s <- counterfact_res[[s]] - - n <- nrow(calres_t) - It_ind <- (calres_s$trt_observed == t) - Is_ind <- (calres_t$trt_observed == s) - pi_t <- sum(It_ind) / n - pi_s <- sum(Is_ind) / n - - cov_y_mus_It <- cov(calres_s$y_observed[It_ind], calres_s$mu_t[It_ind]) - cov_y_mut_Is <- cov(calres_t$y_observed[Is_ind], calres_t$mu_t[Is_ind]) - cov_mut_mus <- cov(calres_s$mu_t, calres_t$mu_t) - - v_out <- cov_y_mus_It + cov_y_mut_Is - cov_mut_mus - } - v_out -} - -vectorized_calculate_V <- Vectorize(calculate_V, c("t", "s")) - -#' helper function: calculate n^(-1)f'Vf (covariance of estimated effect measure) -#' @noRd -op <- function(x, square = FALSE) (x * (1 - x))^ifelse(square, 2, 1) - -odds <- function(x) (x / (1 - x)) - -calculate_fVf <- function(t, s, theta, V, n, eff_measure) { - t_ind <- which(names(theta) == t) - s_ind <- which(names(theta) == s) - - # depneding on eff_measure (1.diff; 2.risk ratio; 3.odds ratio) - # 1. v_tt - 2*v_ts + v_ss - # 2. v_tt/(theta_t^2) - 2*v_ts/(theta_t*theta_s) + v_ss/(theta_s^2) - # 3. v_tt/(op(theta_t)^2) - 2*v_ts/(op(theta_t)*op(theta_s)) + v_ss/(op(theta_s)^2) - if (eff_measure == "diff") { - outtmp <- V[t_ind, t_ind] + V[s_ind, s_ind] - 2 * V[t_ind, s_ind] - } else if (eff_measure == "risk ratio") { - outtmp <- V[t_ind, t_ind] / (theta[t_ind]^2) + - V[s_ind, s_ind] / (theta[s_ind]^2) - - 2 * V[t_ind, s_ind] / (theta[t_ind] * theta[s_ind]) - } else if (eff_measure == "odds ratio") { - outtmp <- V[t_ind, t_ind] / op(theta[t_ind], square = T) + - V[s_ind, s_ind] / op(theta[s_ind], square = T) - - 2 * V[t_ind, s_ind] / (op(theta[t_ind]) * op(theta[s_ind])) - } - outtmp / sqrt(n) -} - -vectorized_calculate_fVf <- Vectorize(calculate_fVf, c("t", "s")) diff --git a/R/variance_anhecova.R b/R/variance_anhecova.R index 051584a..670f685 100644 --- a/R/variance_anhecova.R +++ b/R/variance_anhecova.R @@ -7,7 +7,7 @@ #' #' @return Named covariance matrix. #' @export -vcovAHECOVA <- function(x, decompose = TRUE, randomization = "simple", ...) { +vcovANHECOVA <- function(x, decompose = TRUE, randomization = "simple", ...) { assert_class(x, "prediction_cf") assert_flag(decompose) assert_string(randomization) @@ -37,39 +37,45 @@ vcovAHECOVA <- function(x, decompose = TRUE, randomization = "simple", ...) { return(ret) } +#' Obtain Adjustment for Proportion of Treatment Assignment #' @keywords internal -h_adjust_pi <- function(pi_t) { - diag(pi_t) - pi_t %*% t(pi_t) +#' @param pi (`numeric`) vector of proportions. +#' @return Numeric matrix. +h_adjust_pi <- function(pi) { + assert_numeric(pi) + diag(pi, nrow = length(pi)) - pi %*% t(pi) } +#' Obtain Adjustment for Covariance Matrix +#' @param resi (`numeric`) vector of residuals. +#' @param group_idx (`list` of `integer`) index for each groups. +#' @param trt (`factor`) of treatment assignment. +#' @param pi (`numeric`) proportion of treatment assignment. +#' @param randomization (`string`) name of the randomization schema. #' @keywords internal -h_get_erb <- function(resi, group_idx, trt, pi_t, randomization) { +h_get_erb <- function(resi, group_idx, trt, pi, randomization) { + assert_list(group_idx, types = "integer") if (length(group_idx) == 1L) { return(0) } + assert_string(randomization) if (randomization %in% c("simple", "pocock-simon")) { return(0) } - omegaz_sr <- h_adjust_pi(pi_t) + assert_numeric(resi) + assert_factor(trt) + assert_numeric(pi) + + omegaz_sr <- h_adjust_pi(pi) + n_tot <- length(resi) resi_per_strata <- vapply( group_idx, - function(ii) vapply(split(resi[ii], trt[ii]), mean, FUN.VALUE = 0), + function(ii) { + v <- vapply(split(resi[ii], trt[ii]), mean, FUN.VALUE = 0) + v * sqrt(length(ii) / n_tot) + }, FUN.VALUE = rep(0, length(levels(trt))) ) - strata_props <- vapply( - group_idx, - length, - FUN.VALUE = 0L - ) - strata_props <- strata_props / sum(strata_props) - rb_z <- resi_per_strata / as.numeric(pi_t) - strata_levels <- length(resi_per_strata) - n_trt <- length(pi_t) - ind <- matrix(seq_len(strata_levels), byrow = TRUE, ncol = n_trt) - rb_z_sum <- lapply( - seq_len(ncol(rb_z)), - function(x) rb_z[, x] %*% t(rb_z[, x]) * strata_props[x] - ) - rb_z_sum <- Reduce(`+`, rb_z_sum) - rb_z_sum * omegaz_sr + rb_z <- resi_per_strata / as.numeric(pi) + tcrossprod(rb_z) * omegaz_sr } diff --git a/R/variance_hc.R b/R/variance_hc.R new file mode 100644 index 0000000..7f0be20 --- /dev/null +++ b/R/variance_hc.R @@ -0,0 +1,14 @@ +#### Main functions #### + +#' @exportS3Method +vcovHC.prediction_cf <- function(x, type = "HC3", ...) { + fit <- attr(x, "fit") + vc <- vcovHC(fit, type = type) + mm <- attr(x, "model_matrix") + n <- nrow(mm) / length(names(x)) + md <- family(fit)$mu.eta(attr(x, "predictions_linear")) / n + z <- block_sum(as.vector(md) * mm, n) + ret <- z %*% vc %*% t(z) + dimnames(ret) <- list(names(x), names(x)) + ret +} diff --git a/design/package_structure/design_note_xgc_v1.R b/design/code_replicate/design_note.R similarity index 100% rename from design/package_structure/design_note_xgc_v1.R rename to design/code_replicate/design_note.R diff --git a/design/package_structure/developer_note_xgc_v1.R b/design/code_replicate/developer_note.R similarity index 100% rename from design/package_structure/developer_note_xgc_v1.R rename to design/code_replicate/developer_note.R diff --git a/man/h_adjust_pi.Rd b/man/h_adjust_pi.Rd new file mode 100644 index 0000000..c32e099 --- /dev/null +++ b/man/h_adjust_pi.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{h_adjust_pi} +\alias{h_adjust_pi} +\title{Obtain Adjustment for Proportion of Treatment Assignment} +\usage{ +h_adjust_pi(pi) +} +\arguments{ +\item{pi}{(\code{numeric}) vector of proportions.} +} +\value{ +Numeric matrix. +} +\description{ +Obtain Adjustment for Proportion of Treatment Assignment +} +\keyword{internal} diff --git a/man/h_get_erb.Rd b/man/h_get_erb.Rd new file mode 100644 index 0000000..7714c0b --- /dev/null +++ b/man/h_get_erb.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{h_get_erb} +\alias{h_get_erb} +\title{Obtain Adjustment for Covariance Matrix} +\usage{ +h_get_erb(resi, group_idx, trt, pi, randomization) +} +\arguments{ +\item{resi}{(\code{numeric}) vector of residuals.} + +\item{group_idx}{(\code{list} of \code{integer}) index for each groups.} + +\item{trt}{(\code{factor}) of treatment assignment.} + +\item{pi}{(\code{numeric}) proportion of treatment assignment.} + +\item{randomization}{(\code{string}) name of the randomization schema.} +} +\description{ +Obtain Adjustment for Covariance Matrix +} +\keyword{internal} diff --git a/man/report_vcov_robincar.Rd b/man/report_vcov_robincar.Rd deleted file mode 100644 index 1176ce5..0000000 --- a/man/report_vcov_robincar.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/report_vcovRobinCar.R -\name{report_vcov_robincar} -\alias{report_vcov_robincar} -\title{Reporting function for vcov_robincar} -\usage{ -report_vcov_robincar(x, digits = 3) -} -\arguments{ -\item{x}{object returned by \code{\link{vcovRobinCar.prediction_cf}}} - -\item{digits}{integer, number of decimal places of estimates in the reported table} -} -\value{ -a data frame -} -\description{ -Pretty reporting for the object returned by \code{\link{vcovRobinCar.prediction_cf}}. -Note, if eff_measure was not specified (i.e. NULL), no reporting table will be generated. -} diff --git a/man/vcovAHECOVA.Rd b/man/vcovANHECOVA.Rd similarity index 81% rename from man/vcovAHECOVA.Rd rename to man/vcovANHECOVA.Rd index e4c8ec5..37b859d 100644 --- a/man/vcovAHECOVA.Rd +++ b/man/vcovANHECOVA.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/variance_anhecova.R -\name{vcovAHECOVA} -\alias{vcovAHECOVA} +\name{vcovANHECOVA} +\alias{vcovANHECOVA} \title{ANHECOVA Covariance} \usage{ -vcovAHECOVA(x, decompose = TRUE, randomization = "simple", ...) +vcovANHECOVA(x, decompose = TRUE, randomization = "simple", ...) } \arguments{ \item{x}{(\code{prediction_cf}) Counter-factual prediction.} diff --git a/man/vcovRobinCar.prediction_cf.Rd b/man/vcovRobinCar.prediction_cf.Rd deleted file mode 100644 index 5ba4d88..0000000 --- a/man/vcovRobinCar.prediction_cf.Rd +++ /dev/null @@ -1,33 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/variance.R -\name{vcovRobinCar.prediction_cf} -\alias{vcovRobinCar.prediction_cf} -\title{Covariance matrix of Counterfactual prediction} -\usage{ -vcovRobinCar.prediction_cf(x, eff_measure = NULL, ...) -} -\arguments{ -\item{eff_measure}{default is NULL, the covariance matrix of the mean counterfactual -predictions of the response. If not NULL, a valid option is "diff", "odds ratio", or "risk ratio", -in which case covariance matrix of the estimated effect measure will also be returned.} - -\item{prediction_cf}{R object returned by \code{\line{predict_counterfactual}}} -} -\value{ -an R object of class 'vcov_robincar', essentially a list consisting -\itemize{ -\item theta: mean counterfactual prediction of the reponse, as the calculation result of -\code{\line{predict_counterfactual}} -\item V: covariance matrix of the mean counterfactual predictions of the response -\item f_vcov: covariance matrix of the estimated effect measures. It is NULL -if \code{eff_measure} is set to NULL (default) -\item n: sample size in \code{fit} of \code{\line{predict_counterfactual}} -\item eff_measure: \code{eff_measure} in the argument of this function -\item fit: \code{fit} of \code{\line{predict_counterfactual}} -\item trt_var_name: name of treatment variable in the data used in \code{fit} -of \code{\line{predict_counterfactual}} -} -} -\description{ -Covariance matrix of Counterfactual prediction -} diff --git a/tests/testthat/_snaps/report.md b/tests/testthat/_snaps/report.md deleted file mode 100644 index e0cf4de..0000000 --- a/tests/testthat/_snaps/report.md +++ /dev/null @@ -1,113 +0,0 @@ -# reporting function give info message when eff_measure is NULL - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, - "treatment", dummy_data)), digits = 3) - Message - i When eff_measure is NULL, nothing to report - -# reporting function works for gaussian fitted by lm - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, - "treatment", dummy_data), eff_measure = "diff"), digits = 3) - Output - Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. - 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.564 0.498 - 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.771 0.502 - 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.564 0.498 - 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.207 0.529 - 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.771 0.502 - 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.207 0.529 - 95%CI_lower 95%CI_upper - 1 -0.412 1.539 - 2 -0.213 1.755 - 3 -1.539 0.412 - 4 -0.829 1.244 - 5 -1.755 0.213 - 6 -1.244 0.829 - -# reporting function works for guassian fitted by glm - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, - "treatment", dummy_data), eff_measure = "diff"), digits = 3) - Output - Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. - 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.564 0.498 - 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.771 0.502 - 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.564 0.498 - 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.207 0.529 - 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.771 0.502 - 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.207 0.529 - 95%CI_lower 95%CI_upper - 1 -0.412 1.539 - 2 -0.213 1.755 - 3 -1.539 0.412 - 4 -0.829 1.244 - 5 -1.755 0.213 - 6 -1.244 0.829 - -# reporting function works for binomial fitted by glm - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( - fit_binom, "treatment", dummy_data), eff_measure = "diff"), digits = 3) - Output - Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. - 1 600 trt1 198(33.0%) pbo 202(33.7%) diff 0.225 0.236 - 2 600 trt2 200(33.3%) pbo 202(33.7%) diff 0.265 0.236 - 3 600 pbo 202(33.7%) trt1 198(33.0%) diff -0.225 0.236 - 4 600 trt2 200(33.3%) trt1 198(33.0%) diff 0.041 0.237 - 5 600 pbo 202(33.7%) trt2 200(33.3%) diff -0.265 0.236 - 6 600 trt1 198(33.0%) trt2 200(33.3%) diff -0.041 0.237 - 95%CI_lower 95%CI_upper - 1 -0.238 0.687 - 2 -0.196 0.727 - 3 -0.687 0.238 - 4 -0.425 0.506 - 5 -0.727 0.196 - 6 -0.506 0.425 - ---- - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( - fit_binom, "treatment", dummy_data), eff_measure = "risk ratio"), digits = 3) - Output - Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. - 1 600 trt1 198(33.0%) pbo 202(33.7%) risk ratio 1.631 2.610 - 2 600 trt2 200(33.3%) pbo 202(33.7%) risk ratio 1.745 2.890 - 3 600 pbo 202(33.7%) trt1 198(33.0%) risk ratio 0.613 0.369 - 4 600 trt2 200(33.3%) trt1 198(33.0%) risk ratio 1.070 0.650 - 5 600 pbo 202(33.7%) trt2 200(33.3%) risk ratio 0.573 0.312 - 6 600 trt1 198(33.0%) trt2 200(33.3%) risk ratio 0.935 0.496 - 95%CI_lower 95%CI_upper - 1 0.558 4.762 - 2 0.609 5.001 - 3 0.210 1.791 - 4 0.493 2.324 - 5 0.200 1.642 - 6 0.430 2.030 - ---- - - Code - report_vcov_robincar(vcovRobinCar.prediction_cf(predict_counterfactual( - fit_binom, "treatment", dummy_data), eff_measure = "odds ratio"), digits = 3) - Output - Total Arm_1 N_Arm_1 Arm_0 N_Arm_0 Eff_measure Estimate S.E. - 1 600 trt1 198(33.0%) pbo 202(33.7%) odds ratio 2.504 29.088 - 2 600 trt2 200(33.3%) pbo 202(33.7%) odds ratio 2.968 43.252 - 3 600 pbo 202(33.7%) trt1 198(33.0%) odds ratio 0.399 0.740 - 4 600 trt2 200(33.3%) trt1 198(33.0%) odds ratio 1.185 6.372 - 5 600 pbo 202(33.7%) trt2 200(33.3%) odds ratio 0.337 0.558 - 6 600 trt1 198(33.0%) trt2 200(33.3%) odds ratio 0.844 3.231 - 95%CI_lower 95%CI_upper - 1 0.354 17.710 - 2 0.407 21.653 - 3 0.056 2.824 - 4 0.170 8.280 - 5 0.046 2.459 - 6 0.121 5.896 - diff --git a/tests/testthat/_snaps/variance.md b/tests/testthat/_snaps/variance.md index 8f16314..5c5c8c3 100644 --- a/tests/testthat/_snaps/variance.md +++ b/tests/testthat/_snaps/variance.md @@ -1,375 +1,30 @@ -# robincar convariance estimation works for guassian fitted by lm +# vcovHC works Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, "treatment", - dummy_data)) + vcovHC(pc) Output - $theta - pbo trt1 trt2 - 0.2003208 0.7639709 0.9712499 - - $V - pbo trt1 trt2 - pbo 2.74152393 0.05972904 0.05201403 - trt1 0.05972904 3.44436398 0.06333053 - trt2 0.05201403 0.06333053 3.53553897 - - $f_vcov - NULL - - $n - [1] 600 - - $eff_measure - NULL - - $trt_var_name - [1] "treatment" - - $fit - - Call: - lm(formula = y ~ treatment * s1 + covar, data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - 0.33108 0.54454 0.73872 -0.27656 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.20396 0.03796 0.06399 - - - attr(,"class") - [1] "list" "vcov_robincar" + pbo trt1 trt2 + pbo 1.106024e-03 4.523445e-07 -9.709004e-06 + trt1 4.523445e-07 1.164889e-03 -7.709031e-07 + trt2 -9.709004e-06 -7.709031e-07 1.170214e-03 ---- - - Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_lm, "treatment", - dummy_data), eff_measure = "diff") - Output - $theta - pbo trt1 trt2 - 0.2003208 0.7639709 0.9712499 - - $V - pbo trt1 trt2 - pbo 2.74152393 0.05972904 0.05201403 - trt1 0.05972904 3.44436398 0.06333053 - trt2 0.05201403 0.06333053 3.53553897 - - $f_vcov - pbo trt1 trt2 - pbo 0.0000000 0.2476610 0.2520131 - trt1 0.2476610 0.0000000 0.2797824 - trt2 0.2520131 0.2797824 0.0000000 - - $n - [1] 600 - - $eff_measure - [1] "diff" - - $trt_var_name - [1] "treatment" - - $fit - - Call: - lm(formula = y ~ treatment * s1 + covar, data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - 0.33108 0.54454 0.73872 -0.27656 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.20396 0.03796 0.06399 - - - attr(,"class") - [1] "list" "vcov_robincar" - -# robincar convariance estimation works for guassian fitted by glm - - Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, "treatment", - dummy_data)) - Output - $theta - pbo trt1 trt2 - 0.2003208 0.7639709 0.9712499 - - $V - pbo trt1 trt2 - pbo 2.74152393 0.05972904 0.05201403 - trt1 0.05972904 3.44436398 0.06333053 - trt2 0.05201403 0.06333053 3.53553897 - - $f_vcov - NULL - - $n - [1] 600 - - $eff_measure - NULL - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y ~ treatment * s1 + covar, data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - 0.33108 0.54454 0.73872 -0.27656 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.20396 0.03796 0.06399 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 724.3 - Residual Deviance: 632.3 AIC: 1750 - - attr(,"class") - [1] "list" "vcov_robincar" - ---- - - Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_glm, "treatment", - dummy_data), eff_measure = "diff") - Output - $theta - pbo trt1 trt2 - 0.2003208 0.7639709 0.9712499 - - $V - pbo trt1 trt2 - pbo 2.74152393 0.05972904 0.05201403 - trt1 0.05972904 3.44436398 0.06333053 - trt2 0.05201403 0.06333053 3.53553897 - - $f_vcov - pbo trt1 trt2 - pbo 0.0000000 0.2476610 0.2520131 - trt1 0.2476610 0.0000000 0.2797824 - trt2 0.2520131 0.2797824 0.0000000 - - $n - [1] 600 - - $eff_measure - [1] "diff" - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y ~ treatment * s1 + covar, data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - 0.33108 0.54454 0.73872 -0.27656 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.20396 0.03796 0.06399 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 724.3 - Residual Deviance: 632.3 AIC: 1750 - - attr(,"class") - [1] "list" "vcov_robincar" - -# robincar convariance estimation works for binomial fitted by glm - - Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", - dummy_data)) - Output - $theta - pbo trt1 trt2 - 0.3560965 0.5806957 0.6213865 - - $V - pbo trt1 trt2 - pbo 0.676245627 0.01113740 0.008003312 - trt1 0.011137402 0.70765788 0.013068673 - trt2 0.008003312 0.01306867 0.698427494 - - $f_vcov - NULL - - $n - [1] 600 - - $eff_measure - NULL - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), - data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - -0.42927 0.97415 1.10707 -0.42043 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.41049 -0.01731 0.06947 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 831 - Residual Deviance: 773.7 AIC: 787.7 - - attr(,"class") - [1] "list" "vcov_robincar" - ---- - - Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", - dummy_data), eff_measure = "diff") - Output - $theta - pbo trt1 trt2 - 0.3560965 0.5806957 0.6213865 - - $V - pbo trt1 trt2 - pbo 0.676245627 0.01113740 0.008003312 - trt1 0.011137402 0.70765788 0.013068673 - trt2 0.008003312 0.01306867 0.698427494 - - $f_vcov - pbo trt1 trt2 - pbo 0.00000000 0.05558826 0.05546733 - trt1 0.05558826 0.00000000 0.05633614 - trt2 0.05546733 0.05633614 0.00000000 - - $n - [1] 600 - - $eff_measure - [1] "diff" - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), - data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - -0.42927 0.97415 1.10707 -0.42043 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.41049 -0.01731 0.06947 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 831 - Residual Deviance: 773.7 AIC: 787.7 - - attr(,"class") - [1] "list" "vcov_robincar" - ---- +# vcovANHECOVA works Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", - dummy_data), eff_measure = "risk ratio") + vcovANHECOVA(pc) Output - $theta - pbo trt1 trt2 - 0.3560965 0.5806957 0.6213865 - - $V - pbo trt1 trt2 - pbo 0.676245627 0.01113740 0.008003312 - trt1 0.011137402 0.70765788 0.013068673 - trt2 0.008003312 0.01306867 0.698427494 - - $f_vcov - pbo trt1 trt2 - pbo 0.0000000 0.2989941 0.2886095 - trt1 0.2989941 0.0000000 0.1565623 - trt2 0.2886095 0.1565623 0.0000000 - - $n - [1] 600 - - $eff_measure - [1] "risk ratio" - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), - data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - -0.42927 0.97415 1.10707 -0.42043 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.41049 -0.01731 0.06947 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 831 - Residual Deviance: 773.7 AIC: 787.7 - - attr(,"class") - [1] "list" "vcov_robincar" + pbo trt1 trt2 + pbo 1.128902e-03 1.856234e-05 1.333885e-05 + trt1 1.856234e-05 1.184599e-03 2.178112e-05 + trt2 1.333885e-05 2.178112e-05 1.157268e-03 --- Code - vcovRobinCar.prediction_cf(predict_counterfactual(fit_binom, "treatment", - dummy_data), eff_measure = "odds ratio") + vcovANHECOVA(pc, randomization = "permute_block") Output - $theta - pbo trt1 trt2 - 0.3560965 0.5806957 0.6213865 - - $V - pbo trt1 trt2 - pbo 0.676245627 0.01113740 0.008003312 - trt1 0.011137402 0.70765788 0.013068673 - trt2 0.008003312 0.01306867 0.698427494 - - $f_vcov - pbo trt1 trt2 - pbo 0.0000000 0.9961186 1.0281436 - trt1 0.9961186 0.0000000 0.9838131 - trt2 1.0281436 0.9838131 0.0000000 - - $n - [1] 600 - - $eff_measure - [1] "odds ratio" - - $trt_var_name - [1] "treatment" - - $fit - - Call: glm(formula = y_b ~ treatment * s1 + covar, family = binomial(), - data = dummy_data) - - Coefficients: - (Intercept) treatmenttrt1 treatmenttrt2 s1b - -0.42927 0.97415 1.10707 -0.42043 - covar treatmenttrt1:s1b treatmenttrt2:s1b - 0.41049 -0.01731 0.06947 - - Degrees of Freedom: 599 Total (i.e. Null); 593 Residual - Null Deviance: 831 - Residual Deviance: 773.7 AIC: 787.7 - - attr(,"class") - [1] "list" "vcov_robincar" + pbo trt1 trt2 + pbo 1.128902e-03 1.856234e-05 1.333885e-05 + trt1 1.856234e-05 1.184599e-03 2.178112e-05 + trt2 1.333885e-05 2.178112e-05 1.157268e-03 diff --git a/tests/testthat/test-report.R b/tests/testthat/test-report.R deleted file mode 100644 index 2c7a2f8..0000000 --- a/tests/testthat/test-report.R +++ /dev/null @@ -1,41 +0,0 @@ -test_that("reporting function give info message when eff_measure is NULL", { - expect_snapshot( - predict_counterfactual(fit_lm, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf() |> - report_vcov_robincar(digits = 3) - ) -}) - -test_that("reporting function works for gaussian fitted by lm", { - expect_snapshot( - predict_counterfactual(fit_lm, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf(eff_measure = "diff") |> - report_vcov_robincar(digits = 3) - ) -}) - -test_that("reporting function works for guassian fitted by glm", { - expect_snapshot( - predict_counterfactual(fit_glm, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf(eff_measure = "diff") |> - report_vcov_robincar(digits = 3) - ) -}) - -test_that("reporting function works for binomial fitted by glm", { - expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf(eff_measure = "diff") |> - report_vcov_robincar(digits = 3) - ) - expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf(eff_measure = "risk ratio") |> - report_vcov_robincar(digits = 3) - ) - expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> - vcovRobinCar.prediction_cf(eff_measure = "odds ratio") |> - report_vcov_robincar(digits = 3) - ) -}) diff --git a/tests/testthat/test-variance.R b/tests/testthat/test-variance.R index b7a4014..217f294 100644 --- a/tests/testthat/test-variance.R +++ b/tests/testthat/test-variance.R @@ -1,34 +1,16 @@ -test_that("robincar convariance estimation works for guassian fitted by lm", { +test_that("vcovHC works", { + pc <- predict_counterfactual(fit_binom, treatment ~ s1) expect_snapshot( - predict_counterfactual(fit_lm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() - ) - - expect_snapshot( - predict_counterfactual(fit_lm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") - ) -}) - -test_that("robincar convariance estimation works for guassian fitted by glm", { - expect_snapshot( - predict_counterfactual(fit_glm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() - ) - - expect_snapshot( - predict_counterfactual(fit_glm, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") + vcovHC(pc) ) }) -test_that("robincar convariance estimation works for binomial fitted by glm", { - expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf() - ) - expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "diff") - ) +test_that("vcovANHECOVA works", { + pc <- predict_counterfactual(fit_binom, treatment ~ s1) expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "risk ratio") + vcovANHECOVA(pc) ) expect_snapshot( - predict_counterfactual(fit_binom, "treatment", dummy_data) |> vcovRobinCar.prediction_cf(eff_measure = "odds ratio") + vcovANHECOVA(pc, randomization = "permute_block") ) })