diff --git a/R/matching.R b/R/matching.R index 85294ff3..81282b72 100644 --- a/R/matching.R +++ b/R/matching.R @@ -421,8 +421,11 @@ plot.maicplus_estimate_weights <- function(x, ggplot = FALSE, #' Check to see if weights are optimized correctly #' -#' This function checks to see if the optimization is done properly by checking the covariate averages -#' before and after adjustment. +#' This function checks to see if the optimization is done properly +#' by checking the covariate averages before and after adjustment. +#' In case of ties when calculating median, +#' we return the mean of the two numbers. For more details, +#' see `ties` parameter in [matrixStats::weightedMedian]. #' #' @param weighted_data object returned after calculating weights using \code{\link{estimate_weights}} #' @param processed_agd a data frame, object returned after using \code{\link{process_agd}} or @@ -487,6 +490,8 @@ check_weights <- function(weighted_data, processed_agd) { outdata$internal_trial_after_weighted[ii] <- weightedMedian( x = ipd_with_weights[[covname]], w = ipd_with_weights$weights, + interpolate = FALSE, + ties = "mean", na.rm = TRUE ) # no IPD equals to reported AgD median diff --git a/man/check_weights.Rd b/man/check_weights.Rd index 6da46394..630af0e3 100644 --- a/man/check_weights.Rd +++ b/man/check_weights.Rd @@ -39,8 +39,11 @@ data.frame of weighted and unweighted covariate averages of the IPD, average of aggregate data, and sum of inner products of covariate \eqn{x_i} and the weights (\eqn{exp(x_i\beta)}) } \description{ -This function checks to see if the optimization is done properly by checking the covariate averages -before and after adjustment. +This function checks to see if the optimization is done properly +by checking the covariate averages before and after adjustment. +In case of ties when calculating median, +we return the mean of the two numbers. For more details, +see \code{ties} parameter in \link[matrixStats:weightedMedian]{matrixStats::weightedMedian}. } \section{Methods (by generic)}{ \itemize{ diff --git a/tests/testthat/_snaps/maic_unanchored.md b/tests/testthat/_snaps/maic_unanchored.md index 490a4819..38f5ecad 100644 --- a/tests/testthat/_snaps/maic_unanchored.md +++ b/tests/testthat/_snaps/maic_unanchored.md @@ -287,7 +287,7 @@ --- Code - print(testout_boot_RR$inferential$fit, digits = 5) + testout_boot_RR$inferential$fit Output $model_before @@ -295,11 +295,11 @@ Coefficients: (Intercept) ARMA - -0.35667 0.10821 + -0.3567 0.1082 Degrees of Freedom: 899 Total (i.e. Null); 898 Residual Null Deviance: 1023 - Residual Deviance: 1015.6 AIC: 1019.6 + Residual Deviance: 1016 AIC: 1020 $model_after @@ -308,44 +308,44 @@ Coefficients: (Intercept) ARMA - -0.356675 0.023352 + -0.35667 0.02335 Degrees of Freedom: 899 Total (i.e. Null); 898 Residual - Null Deviance: 726.66 - Residual Deviance: 726.48 AIC: 712.47 + Null Deviance: 726.7 + Residual Deviance: 726.5 AIC: 712.5 $res_AB $res_AB$est - [1] 1.0236 + [1] 1.023627 $res_AB$se - [1] 0.060252 + [1] 0.06025155 $res_AB$ci_l - [1] 0.91236 + [1] 0.9123647 $res_AB$ci_u - [1] 1.1485 + [1] 1.148457 $res_AB$pval - [1] 0.69081 + [1] 0.6908096 $res_AB_unadj $res_AB_unadj$est - [1] 1.1143 + [1] 1.114286 $res_AB_unadj$se - [1] 0.045119 + [1] 0.04511891 $res_AB_unadj$ci_l - [1] 1.0294 + [1] 1.029372 $res_AB_unadj$ci_u - [1] 1.2062 + [1] 1.206204 $res_AB_unadj$pval - [1] 0.0074553 + [1] 0.007455267 $boot_res @@ -359,22 +359,22 @@ Bootstrap Statistics : - original bias std. error - t1* 0.0233518 0.01366708 0.05380528 - t2* 0.0030551 -0.00004315 0.00038235 + original bias std. error + t1* 0.023351847 1.366708e-02 0.0538052765 + t2* 0.003055105 -4.315049e-05 0.0003823512 $boot_res_AB $boot_res_AB$est - [1] 1.0236 + [1] 1.023627 $boot_res_AB$se [1] NA $boot_res_AB$ci_l - [1] 0.90867 + [1] 0.9086715 $boot_res_AB$ci_u - [1] 1.122 + [1] 1.122032 $boot_res_AB$pval [1] NA @@ -384,11 +384,11 @@ --- Code - print(testout_boot_RR$inferential$summary, digits = 5) + testout_boot_RR$inferential$summary Output - case RR LCL UCL pval - 1 AB 1.1143 1.02937 1.2062 0.0074553 - 2 adjusted_AB 1.0236 0.91236 1.1485 0.6908096 + case RR LCL UCL pval + 1 AB 1.114286 1.0293724 1.206204 0.007455267 + 2 adjusted_AB 1.023627 0.9123647 1.148457 0.690809607 # test time to event case diff --git a/tests/testthat/_snaps/matching.md b/tests/testthat/_snaps/matching.md new file mode 100644 index 00000000..38e51764 --- /dev/null +++ b/tests/testthat/_snaps/matching.md @@ -0,0 +1,20 @@ +# weighted median is correct + + Code + check_weights_result + Output + covariate match_stat internal_trial internal_trial_after_weighted + 1 AGE Mean 59.850 51.00 + 2 AGE SD 9.011 3.25 + 3 SEX_MALE Prop 0.380 0.49 + 4 ECOG0 Prop 0.410 0.40 + 5 SMOKE Prop 0.320 0.20 + 6 N_PR_THER Median 3.000 2.00 + external_trial sum_centered_IPD_with_weights + 1 51.00 0.0001 + 2 3.25 0.0125 + 3 0.49 0.0000 + 4 0.40 0.0000 + 5 0.20 0.0000 + 6 2.00 0.0000 + diff --git a/tests/testthat/test-matching.R b/tests/testthat/test-matching.R index 1cbdacc1..a644e4d4 100644 --- a/tests/testthat/test-matching.R +++ b/tests/testthat/test-matching.R @@ -253,3 +253,30 @@ test_that("check_weights works as expected", { expect_equal(result[, "match_stat"], c("Mean", "Median", "SD", "Prop", "Prop", "Prop")) expect_equal(result[, "sum_centered_IPD_with_weights"], c(0, 0, -0.0045, 0, 0, 0)) }) + +test_that("weighted median is correct", { + # based on calculating_weights vignette + data <- adsl_sat + data$SEX_MALE <- ifelse(data$SEX == "Male", 1, 0) + data$SAGE_SQUARED <- data$AGE^2 + agd <- data.frame( + AGE_MEAN = 51, + AGE_SD = 3.25, + SEX_MALE_PROP = 147 / 300, + ECOG0_PROP = 0.40, + SMOKE_PROP = 58 / (300 - 5), + N_PR_THER_MEDIAN = 2 + ) + + ipd <- adsl_sat + ipd_centered <- center_ipd(ipd = ipd, agd = agd) + centered_colnames <- c("AGE", "AGE_SQUARED", "SEX_MALE", "ECOG0", "SMOKE", "N_PR_THER_MEDIAN") + centered_colnames <- paste0(centered_colnames, "_CENTERED") + + weighted_sat <- estimate_weights( + data = ipd_centered, + centered_colnames = centered_colnames + ) + check_weights_result <- check_weights(weighted_sat, agd) + expect_snapshot(check_weights_result) +})