Merge pull request #8 from ATpoint/master

Support sample weights in circa_single() and circacompare()

DESCRIPTION

@@ -1,27 +1,31 @@
-Package: circacompare
-Title: Analyses of Circadian Data
-Version: 0.1.1.9000
-Authors@R: 
-    person(given = "Rex",
-           family = "Parsons",
-           role = c("aut", "cre"),
-           email = "[email protected]",
-           comment = c(ORCID = "0000-0002-6053-8174"))
-Description: Uses non-linear regression to statistically compare two circadian rhythms.
-    Groups are only compared if both are rhythmic (amplitude is non-zero).
-    Performs analyses regarding mesor, phase, and amplitude, reporting on estimates and statistical differences, for each, between groups.
-    Details can be found in Parsons et al (2020) <doi:10.1093/bioinformatics/btz730>.
-License: MIT + file LICENSE
-Encoding: UTF-8
-Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.1
-Imports: 
-    ggplot2 (>= 2.2.1),
-    stats
-Suggests: 
-    testthat (>= 3.0.0),
-    nlme,
-    knitr,
-    rmarkdown
-Config/testthat/edition: 3
-VignetteBuilder: knitr
+Package: circacompare
+Title: Analyses of Circadian Data
+Version: 0.1.1.9000
+Authors@R: 
+    c(person(given = "Rex",
+             family = "Parsons",
+             role = c("aut", "cre"),
+             email = "[email protected]",
+             comment = c(ORCID = "0000-0002-6053-8174")),
+       person(given = "Alexander",
+              family = "Bender",
+              role = c("ctb"),
+              email = "[email protected]"))
+Description: Uses non-linear regression to statistically compare two circadian rhythms.
+    Groups are only compared if both are rhythmic (amplitude is non-zero).
+    Performs analyses regarding mesor, phase, and amplitude, reporting on estimates and statistical differences, for each, between groups.
+    Details can be found in Parsons et al (2020) <doi:10.1093/bioinformatics/btz730>.
+License: MIT + file LICENSE
+Encoding: UTF-8
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.2.3
+Imports: 
+    ggplot2 (>= 2.2.1),
+    stats
+Suggests: 
+    testthat (>= 3.0.0),
+    nlme,
+    knitr,
+    rmarkdown
+Config/testthat/edition: 3
+VignetteBuilder: knitr

NEWS.md

@@ -2,6 +2,8 @@
 
 ## Improvements
 
+* support per-sample weights in `circacompare()`, `circa_single()`, `circa_single_mixed()` and `circacompare_mixed()`
+
 * better error messages - also makes the shiny app more verbose for users
 
 * allow the user to use `suppress_all` argument to show/hide messages during model fitting across all circacompare functions.

R/circa_single.R

@@ -11,6 +11,7 @@
 #' @param timeout_n The upper limit for the model fitting attempts. Default is 10,000.
 #' @param return_figure Whether or not to return a ggplot graph of the rhythm and cosine model.
 #' @param control \code{list}. Used to control the parameterization of the model.
+#' @param weights An optional numeric vector of (fixed) weights. When present, the objective function is weighted least squares.
 #' @param suppress_all Logical. Set to \code{TRUE} to avoid seeing errors or messages during model fitting procedure. Default is \code{FALSE}.
 #'
 #' @return list
@@ -19,7 +20,20 @@
 #' @examples
 #' df <- make_data()
 #' df <- df[df$group == "g1", ]
-#' circa_single(x = df, col_time = "time", col_outcome = "measure")
+#' out <- circa_single(x = df, col_time = "time", col_outcome = "measure")
+#' out
+#'
+#' # with sample weights (arbitrary weights for demonstration)
+#' sw <- runif(n = nrow(df))
+#' out2 <- circa_single(
+#'   x = df,
+#'   col_time = "time",
+#'   col_outcome = "measure",
+#'   weights = sw,
+#'   suppress_all = TRUE
+#' )
+#' out2
+#'
 circa_single <- function(x,
                          col_time,
                          col_outcome,
@@ -28,6 +42,7 @@ circa_single <- function(x,
                          timeout_n = 10000,
                          return_figure = TRUE,
                          control = list(),
+                         weights = NULL,
                          suppress_all = FALSE) {
   controlVals <- circa_single_control()
   controlVals[names(control)] <- control
@@ -81,6 +96,14 @@ circa_single <- function(x,
       ))
     }
   }
+
+  if (!is.null(weights)) {
+    check_weights(x, weights)
+    x$weights <- weights
+  } else {
+    x$weights <- rep(1, nrow(x))
+  }
+
   success <- FALSE
   n <- 0
   form <- create_formula(main_params = controlVals$main_params, decay_params = controlVals$decay_params)$formula
@@ -91,7 +114,8 @@ circa_single <- function(x,
         stats::nls(
           formula = form,
           data = x,
-          start = start_list(outcome = x$measure, controlVals = controlVals)
+          start = start_list(outcome = x$measure, controlVals = controlVals),
+          weights = weights
         )
       },
       silent = suppress_all

R/circa_single_mixed.R

@@ -12,6 +12,8 @@
 #' @param alpha_threshold The level of alpha for which the presence of rhythmicity is considered. Default is to \code{0.05}.
 #' @param nlme_control A list of control values for the estimation algorithm to replace the default values returned by the function nlme::nlmeControl. Defaults to an empty list.
 #' @param nlme_method A character string. If "REML" the model is fit by maximizing the restricted log-likelihood. If "ML" the log-likelihood is maximized. Defaults to "ML".
+#' @param weights An optional numeric vector of (fixed) weights internally passed to \code{nlme::nlme()} via \code{nlme::varPower()}.
+#' When present, the objective function is weighted least squares.
 #' @param suppress_all Logical. Set to \code{TRUE} to avoid seeing errors or messages during model fitting procedure. Default is \code{FALSE}. If \code{FALSE}, also runs \code{nlme()} with \code{verbose = TRUE}.
 #' @param timeout_n The upper limit for the model fitting attempts. Default is \code{10000}.
 #' @param return_figure Whether or not to return a ggplot graph of the rhythm and cosine model.
@@ -41,6 +43,14 @@
 #'   x = df, col_time = "time", col_outcome = "measure",
 #'   col_id = "id", randomeffects = c("k")
 #' )
+#'
+#' # with sample weights (arbitrary weights for demonstration)
+#' sw <- runif(n = nrow(df))
+#' out2 <- circa_single_mixed(
+#'   x = df, col_time = "time", col_outcome = "measure",
+#'   col_id = "id", randomeffects = c("k"), weights = sw
+#' )
+#'
 circa_single_mixed <- function(x,
                                col_time,
                                col_outcome,
@@ -50,6 +60,7 @@ circa_single_mixed <- function(x,
                                alpha_threshold = 0.05,
                                nlme_control = list(),
                                nlme_method = "ML",
+                               weights = NULL,
                                suppress_all = FALSE,
                                timeout_n = 10000,
                                return_figure = TRUE,
@@ -124,6 +135,13 @@ circa_single_mixed <- function(x,
     }
   }
 
+  if (!is.null(weights)) {
+    check_weights(x, weights)
+    x$weights <- weights
+  } else {
+    x$weights <- rep(1, nrow(x))
+  }
+
   success <- FALSE
   n <- 0
   form <- create_formula(main_params = controlVals$main_params, decay_params = controlVals$decay_params)$formula
@@ -141,7 +159,8 @@ circa_single_mixed <- function(x,
           start = unlist(start_list(outcome = x$measure, controlVals = controlVals)),
           control = nlme_control,
           method = nlme_method,
-          verbose = !suppress_all
+          verbose = !suppress_all,
+          weights = nlme::varPower(form = ~weights)
         )
       },
       silent = ifelse(suppress_all, TRUE, FALSE)

R/circacompare.R

@@ -11,6 +11,7 @@
 #' @param alpha_threshold The level of alpha for which the presence of rhythmicity is considered. Default is 0.05.
 #' @param timeout_n The upper limit for the model fitting attempts. Default is 10,000.
 #' @param control \code{list}. Used to control the parameterization of the model.
+#' @param weights An optional numeric vector of (fixed) weights. When present, the objective function is weighted least squares.
 #' @param suppress_all Logical. Set to \code{TRUE} to avoid seeing errors or messages during model fitting procedure. Default is \code{FALSE}.
 #'
 #' @return list
@@ -23,6 +24,15 @@
 #'   col_outcome = "measure"
 #' )
 #' out
+#'
+#' # with sample weights (arbitrary weights for demonstration)
+#' sw <- runif(n = nrow(df))
+#' out2 <- circacompare(
+#'   x = df, col_time = "time", col_group = "group",
+#'   col_outcome = "measure", weights = sw
+#' )
+#' out2
+#'
 circacompare <- function(x,
                          col_time,
                          col_group,
@@ -31,6 +41,7 @@ circacompare <- function(x,
                          alpha_threshold = 0.05,
                          timeout_n = 10000,
                          control = list(),
+                         weights = NULL,
                          suppress_all = FALSE) {
   controlVals <- circacompare_control()
   controlVals[names(control)] <- control
@@ -91,6 +102,13 @@ circacompare <- function(x,
     }
   }
 
+  if (!is.null(weights)) {
+    check_weights(x, weights)
+    x$weights <- weights
+  } else {
+    x$weights <- rep(1, nrow(x))
+  }
+
   group_1_text <- levels(as.factor(x$group))[1]
   group_2_text <- levels(as.factor(x$group))[2]
   x$x_group <- ifelse(x$group == group_1_text, 0, 1)
@@ -149,7 +167,8 @@ circacompare <- function(x,
           formula = form_group,
           data = x,
           start = start_list_grouped(g1 = g1_model$model, g2 = g2_model$model, grouped_params = controlVals$grouped_params),
-          control = stats::nls.control(maxiter = 100, minFactor = 1 / 10000)
+          control = stats::nls.control(maxiter = 100, minFactor = 1 / 10000),
+          weights = weights
         )
       },
       silent = suppress_all

R/circacompare_mixed.R

@@ -13,6 +13,8 @@
 #' @param alpha_threshold The level of alpha for which the presence of rhythmicity is considered. Default is to \code{0.05}.
 #' @param nlme_control A list of control values for the estimation algorithm to replace the default values returned by the function nlme::nlmeControl. Defaults to an empty list.
 #' @param nlme_method A character string. If "REML" the model is fit by maximizing the restricted log-likelihood. If "ML" the log-likelihood is maximized. Defaults to "REML".
+#' @param weights An optional numeric vector of (fixed) weights internally passed to \code{nlme::nlme()} via \code{nlme::varPower()}.
+#' When present, the objective function is weighted least squares.
 #' @param suppress_all Logical. Set to \code{TRUE} to avoid seeing errors or messages during model fitting procedure. Default is \code{FALSE}. If \code{FALSE}, also runs \code{nlme()} with \code{verbose = TRUE}.
 #' @param timeout_n The upper limit for the model fitting attempts. Default is \code{10000}.
 #' @param control \code{list}. Used to control the parameterization of the model.
@@ -48,7 +50,18 @@
 #'   col_id = "id",
 #'   control = list(grouped_params = c("phi"), random_params = c("phi1"))
 #' )
-#' out
+#'
+#' # with sample weights (arbitrary weights for demonstration)
+#' sw <- runif(n = nrow(df))
+#' out2 <- circacompare_mixed(
+#'   x = df,
+#'   col_time = "time",
+#'   col_group = "group",
+#'   col_outcome = "measure",
+#'   col_id = "id",
+#'   control = list(grouped_params = c("phi"), random_params = c("phi1")),
+#'   weights = sw
+#' )
 #'
 circacompare_mixed <- function(x,
                                col_time,
@@ -60,6 +73,7 @@ circacompare_mixed <- function(x,
                                alpha_threshold = 0.05,
                                nlme_control = list(),
                                nlme_method = "REML",
+                               weights = NULL,
                                suppress_all = FALSE,
                                timeout_n = 10000,
                                control = list()) {
@@ -137,6 +151,13 @@ circacompare_mixed <- function(x,
     }
   }
 
+  if (!is.null(weights)) {
+    check_weights(x, weights)
+    x$weights <- weights
+  } else {
+    x$weights <- rep(1, nrow(x))
+  }
+
   group_1_text <- levels(as.factor(x$group))[1]
   group_2_text <- levels(as.factor(x$group))[2]
 
@@ -210,7 +231,8 @@ circacompare_mixed <- function(x,
           start = unlist(start_list_grouped(g1 = g1_model$model, g2 = g2_model$model, grouped_params = controlVals$grouped_params)),
           method = nlme_method,
           control = nlme_control,
-          verbose = !suppress_all
+          verbose = !suppress_all,
+          weights = nlme::varPower(form = ~weights)
         )
       },
       silent = ifelse(suppress_all, TRUE, FALSE)

R/utils.R

@@ -1,4 +1,6 @@
-utils::globalVariables(c("time", "measure", "group", "eq", "eq_1", "eq_2"))
+utils::globalVariables(
+  c("time", "measure", "group", "eq", "eq_1", "eq_2", "weights")
+)
 
 extract_model_coefs <- function(model) {
   if ("nls" %in% class(model)) {
@@ -36,7 +38,8 @@ model_each_group <- function(data, type, form = stats::as.formula("measure~k+alp
             start = unlist(starting_params),
             method = args$nlme_method,
             control = args$nlme_control,
-            verbose = args$verbose
+            verbose = args$verbose,
+            weights = nlme::varPower(form = ~weights)
           )
         },
         silent = ifelse(args$verbose, FALSE, TRUE)
@@ -47,7 +50,8 @@ model_each_group <- function(data, type, form = stats::as.formula("measure~k+alp
           stats::nls(
             formula = form,
             data = data,
-            start = starting_params
+            start = starting_params,
+            weights = weights
           )
         },
         silent = args$suppress_all
@@ -450,3 +454,19 @@ circa_summary <- function(model, period, control,
 
   return(res)
 }
+
+check_weights <- function(x, weights) {
+  len_weights <- length(weights)
+  len_x <- nrow(x)
+
+  if (len_weights != len_x) {
+    stop("weights must have the same length as the number of rows in x")
+  }
+
+  contains_na_or_negative <- sum(is.na(weights) | weights < 0) > 0
+  non_numeric <- !is.numeric(weights)
+
+  if (contains_na_or_negative | non_numeric) {
+    stop("weights be not be negative or contain NAs/missing values")
+  }
+}