Merge pull request #14 from nuvanrath/master

use unweighted variance for denominator for weighted smd
bsaul · May 1, 2024 · e6a993b · e6a993b
2 parents 29cc482 + a99df33
commit e6a993b
Show file tree

Hide file tree

Showing 8 changed files with 139 additions and 79 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: smd
 Type: Package
 Title: Compute Standardized Mean Differences
-Version: 0.6.8
+Version: 0.7.0
 Authors@R: c(person("Bradley", "Saul", role = c("aut", "cre"), 
                     email = "[email protected]"),
              person("Alex", "Breskin", role = c("ctb"), 
@@ -13,11 +13,11 @@ Authors@R: c(person("Bradley", "Saul", role = c("aut", "cre"),
              person("Daniel", "Sjoberg", role = c("ctb"),
                     email = "[email protected]"),
              person("Nuvan", "Rathnayaka", role = c("ctb"),
-                    email = "nuvanrath@mail.proton.me")
+                    email = "[email protected]")
              )
 Description: Computes standardized mean differences and confidence intervals for 
     multiple data types based on Yang, D., & Dalton, J. E. (2012) 
-    <http://www.lerner.ccf.org/qhs/software/lib/stddiff.pdf>. 
+    <https://support.sas.com/resources/papers/proceedings12/335-2012.pdf>. 
 Imports:
   MASS (>= 7.3-50),
   methods (>= 3.5.1)
@@ -35,6 +35,6 @@ URL: https://bsaul.github.io/smd/
 BugReports: https://github.com/bsaul/smd/issues
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 VignetteBuilder: knitr
 Repository: cran.novisci.com
diff --git a/R/mean_var.R b/R/mean_var.R
@@ -21,27 +21,28 @@ multinom_var <- function(p){
 #' @param x a vector of values
 #' @param w an optional vector of \code{numeric} weights
 #' @param na.rm passed to \code{sum}
+#' @param unwgt.var Use unweighted or weighted covariance matrix
 #' @importFrom stats var
 #' @importFrom methods setGeneric setMethod
 #' @return a list containing \code{mean} and \code{var}
 #' @keywords internal
-setGeneric("n_mean_var", def = function(x, w = NULL, na.rm = FALSE){
+setGeneric("n_mean_var", def = function(x, w = NULL, na.rm = FALSE, unwgt.var = TRUE){
   standardGeneric("n_mean_var")
 })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("numeric", "missing"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       x <- stats::na.omit(x)
     }
 
     n    <- length(x)
     mean <- sum(x)/n
-
+    
     list(
       n    = n,
       mean = mean,
@@ -53,7 +54,7 @@ setMethod(
 setMethod(
   f          = "n_mean_var",
   signature  = c("numeric", "numeric"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       kp <- !is.na(x)
@@ -67,13 +68,25 @@ setMethod(
 
     xw   <- x * w
     n    <- sum(w)
+
     # Handle case were sum of weights is 0
-    mean <- if(n == 0) 0 else sum(xw)/n
+    if(n == 0){
+      mean = 0
+      var = 0
+    } else if(unwgt.var == TRUE){
+      mean = sum(xw)/n
+      unwgt_n = length(x)
+      unwgt_mean = sum(x)/unwgt_n
+      var = sum((x - unwgt_mean)^2)/unwgt_n
+    } else {
+      mean = sum(xw)/n
+      var = sum(w*(x - mean)^2)/n
+    }
 
     list(
       n    = n,
       mean = mean,
-      var  = if(n == 0) 0 else sum(w*(x - mean)^2)/n
+      var  = var
     )
 })
 
@@ -82,43 +95,43 @@ setMethod(
 setMethod(
   f          = "n_mean_var",
   signature  = c("integer", "missing"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     check_for_two_levels(x)
-    n_mean_var(x = as.numeric(x), na.rm = na.rm)
+    n_mean_var(x = as.numeric(x), na.rm = na.rm, unwgt.var = unwgt.var)
 })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("integer", "numeric"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     check_for_two_levels(x)
-    n_mean_var(x = as.numeric(x), w = w, na.rm = na.rm)
+    n_mean_var(x = as.numeric(x), w = w, na.rm = na.rm, unwgt.var = unwgt.var)
 })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("logical", "missing"),
-  definition = function(x, na.rm = FALSE){
-    n_mean_var(x = as.numeric(x), na.rm = na.rm)
+  definition = function(x, na.rm = FALSE, unwgt.var = TRUE){
+    n_mean_var(x = as.numeric(x), na.rm = na.rm, unwgt.var = unwgt.var)
 })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("logical", "numeric"),
-  definition = function(x, w, na.rm = FALSE){
-    n_mean_var(x = as.numeric(x), w = w, na.rm = na.rm)
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
+    n_mean_var(x = as.numeric(x), w = w, na.rm = na.rm, unwgt.var = unwgt.var)
 })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("factor", "missing"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       x <- stats::na.omit(x)
@@ -132,24 +145,31 @@ setMethod(
 setMethod(
   f          = "n_mean_var",
   signature  = c("factor", "numeric"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       kp <- !is.na(x)
       w <- w[kp]
       x <- x[kp]
     }
-
+    
     n <- sum(w)
     p <- tapply(w, x, function(r) if(n == 0) 0 else sum(r)/n, default = 0)
-    list(n = n, mean = p, var = multinom_var(p))
+    if(unwgt.var == TRUE){
+      unwt_p <- prop.table(table(x)) #for unweighted variance
+      var <- multinom_var(unwt_p)
+    } else {
+      var = multinom_var(p)
+    }
+
+    list(n = n, mean = p, var = var)
   })
 
 #' @rdname n_mean_var
 setMethod(
   f          = "n_mean_var",
   signature  = c("character", "missing"),
-  definition = function(x, w, na.rm = FALSE){
+  definition = function(x, w, na.rm = FALSE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       x <- stats::na.omit(x)
@@ -168,7 +188,7 @@ setMethod(
 setMethod(
   f          = "n_mean_var",
   signature  = c("character", "numeric"),
-  definition = function(x, w, na.rm = TRUE){
+  definition = function(x, w, na.rm = TRUE, unwgt.var = TRUE){
 
     if(na.rm == TRUE){
       kp <- !is.na(x)
@@ -182,7 +202,7 @@ setMethod(
       warning("x has more than 50 levels. Are you sure you meant for this?")
     }
 
-    n_mean_var(x, w)
+    n_mean_var(x, w, unwgt.var = unwgt.var)
 })