hta-pharma · gravesti · Feb 13, 2025 · Feb 10, 2025 · Feb 10, 2025 · Feb 10, 2025
diff --git 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
diff --git a/tests/testthat/_snaps/maic_unanchored.md b/tests/testthat/_snaps/maic_unanchored.md
@@ -287,19 +287,19 @@
 ---
 
     Code
-      print(testout_boot_RR$inferential$fit, digits = 5)
+      testout_boot_RR$inferential$fit
     Output
       $model_before
 
       Call:  glm(formula = RESPONSE ~ ARM, family = glm_link, data = dat)
 
       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
@@ -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
@@ -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)
+})