update residual functions (#267)

* update residual functions

* [skip actions] Restyle files

* update tests

* [skip actions] Roxygen Man Pages Auto Update

* simplify

* allow normal residual for spatial mmrm

* update pearson residuals

* style update

---------

Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com>

NEWS.md

@@ -7,8 +7,7 @@
   customization through the `tidymodels` interface.
 - Add Kenward-Roger support for spatial covariance structures.
 - Add support for `residuals` method with a `type` argument allowing for 
-  raw (the default, and only option for models with a spatial covariance structure),
-  Pearson and normalized residuals to be calculated from an `mmrm` fit.
+  raw, Pearson and normalized residuals to be calculated from an `mmrm` fit.
 - Add empirical, empirical Jackknife and empirical bias-reduced adjusted coefficients covariance matrix.
   In addition, the argument `method` now only specifies the method used
   for the degrees of freedom, another argument `vcov` is added to specify the

R/tmb-methods.R

@@ -244,98 +244,49 @@ print.mmrm_tmb <- function(x,
 #' - \insertRef{galecki2013linear}{mmrm}
 residuals.mmrm_tmb <- function(object, type = c("response", "pearson", "normalized"), ...) {
   type <- match.arg(type)
-  resids_unscaled <- component(object, "y_vector") - unname(fitted(object))
-  if (type == "response") {
-    resids_unscaled
-  } else {
-    if (object$formula_parts$is_spatial) {
-      stop("Only 'response' residuals are available for models with spatial covariance structures.")
-    }
-    if (type == "pearson") {
-      h_residuals_pearson(object, resids_unscaled)
-    } else if (type == "normalized") {
-      h_residuals_normalized(object, resids_unscaled)
-    }
-  }
+  switch(type,
+    "response" = h_residuals_response(object),
+    "pearson" = h_residuals_pearson(object),
+    "normalized" = h_residuals_normalized(object)
+  )
 }
-
 #' Calculate Pearson Residuals
 #'
 #' This is used by [residuals.mmrm_tmb()] to calculate Pearson residuals.
 #'
 #' @param object (`mmrm_tmb`)\cr the fitted MMRM.
-#' @param resids_unscaled (`numeric`)\cr the response residuals.
 #'
 #' @return Vector of residuals.
 #'
 #' @keywords internal
-h_residuals_pearson <- function(object, resids_unscaled) {
+h_residuals_pearson <- function(object) {
   assert_class(object, "mmrm_tmb")
-  assert_numeric(resids_unscaled)
-  visits <- as.numeric(object$tmb_data$full_frame[[object$formula_parts$visit_var]])
-  cov_list <- if (component(object, "n_groups") == 1) {
-    list(object$cov)
-  } else {
-    object$cov
-  }
-  visit_sigmas <- lapply(cov_list, function(x) sqrt(diag(x, names = FALSE)))
-  nobs <- nrow(object$tmb_data$full_frame)
-  subject_grps <- if (component(object, "n_groups") == 1) {
-    rep(1, times = nobs)
-  } else {
-    object$tmb_data$full_frame[[object$formula_parts$group_var]]
-  }
-  sapply(1:nobs, function(x) {
-    resids_unscaled[x] / visit_sigmas[[subject_grps[x]]][visits[x]] * sqrt(object$tmb_data$weights_vector[x])
-  })
+  h_residuals_response(object) * object$tmb_object$report()$diag_cov_inv_sqrt
 }
 
 #' Calculate normalized residuals
 #'
 #' This is used by [residuals.mmrm_tmb()] to calculate normalized / scaled residuals.
 #'
 #' @param object (`mmrm_tmb`)\cr the fitted MMRM.
-#' @param resids_unscaled (`numeric`)\cr the raw/response residuals.
 #'
 #' @return Vector of residuals
 #'
 #' @keywords internal
-h_residuals_normalized <- function(object, resids_unscaled) {
+h_residuals_normalized <- function(object) {
   assert_class(object, "mmrm_tmb")
-  assert_numeric(resids_unscaled)
-
-  resid_df <- data.frame(
-    subject = object$tmb_data$full_frame[[object$formula_parts$subject_var]],
-    time = as.numeric(object$tmb_data$full_frame[[object$formula_parts$visit_var]]),
-    residual = resids_unscaled,
-    weights = object$tmb_data$weights_vector
-  )
-
-  subject_list <- split(resid_df, resid_df$subject)
-
-  lower_chol_list <- if (component(object, "n_groups") == 1) {
-    lapply(seq_along(subject_list), function(x) {
-      weighted_cov <- object$cov[subject_list[[x]]$time, subject_list[[x]]$time] /
-        sqrt(tcrossprod(matrix(subject_list[[x]]$weights, ncol = 1)))
-
-      solve(t(chol(weighted_cov)))
-    })
-  } else {
-    groups <- data.frame(
-      subject = object$tmb_data$full_frame[[object$formula_parts$subject_var]],
-      group = object$tmb_data$full_frame[[object$formula_parts$group_var]]
-    )
-    groups <- groups[!duplicated(groups), ]
-    lapply(seq_along(subject_list), function(x) {
-      this_cov <- object$cov[[groups$group[x]]]
-
-      weighted_cov <- this_cov[subject_list[[x]]$time, subject_list[[x]]$time] /
-        sqrt(tcrossprod(matrix(subject_list[[x]]$weights, ncol = 1)))
-
-      solve(t(chol(weighted_cov)))
-    })
-  }
-  unlist(lapply(seq_along(subject_list), function(x) {
-    lower_chol_list[[x]] %*% matrix(subject_list[[x]]$residual, ncol = 1)
-  }))
+  object$tmb_object$report()$epsilonTilde
+}
+#' Calculate response residuals.
+#'
+#' This is used by [residuals.mmrm_tmb()] to calculate response residuals.
+#'
+#' @param object (`mmrm_tmb`)\cr the fitted MMRM.
+#'
+#' @return Vector of residuals
+#'
+#' @keywords internal
+h_residuals_response <- function(object) {
+  assert_class(object, "mmrm_tmb")
+  component(object, "y_vector") - unname(fitted(object))
 }

src/mmrm.cpp

@@ -59,6 +59,8 @@ Type objective_function<Type>::operator() ()
   int theta_one_group_size = theta.size() / n_groups;
   // Convert is_spatial_int to bool.
   bool is_spatial = (is_spatial_int == 1);
+  // Diagonal of weighted covariance
+  vector<Type> diag_cov_inv_sqrt(x_matrix.rows());
   // Create the lower triangular Cholesky factor of the visit x visit covariance matrix.
   std::map<int, lower_chol_base<Type>*> chols_by_group;
   for (int r = 0; r < n_groups; r++) {
@@ -84,14 +86,11 @@ Type objective_function<Type>::operator() ()
     } else {
       dist_i = euclidean(matrix<Type>(coordinates.block(start_i, 0, n_visits_i, coordinates.cols())));
     }
-
     // Obtain Cholesky factor Li.
-    matrix<Type> Li = chols_by_group[subject_groups[i]]->get_chol(visit_i, dist_i);
-
+    matrix<Type> Li = chols_by_group[subject_groups[i]]->get_chol(visit_i, dist_i);  
     // Calculate weighted Cholesky factor for this subject.
     Eigen::DiagonalMatrix<Type,Eigen::Dynamic,Eigen::Dynamic> Gi_inv_sqrt = weights_vector.segment(start_i, n_visits_i).cwiseInverse().sqrt().matrix().asDiagonal();
     Li = Gi_inv_sqrt * Li;
-
     // Calculate scaled design matrix and response vector for this subject.
     matrix<Type> Xi = x_matrix.block(start_i, 0, n_visits_i, x_matrix.cols());
     matrix<Type> XiTilde = Li.template triangularView<Eigen::Lower>().solve(Xi);
@@ -104,7 +103,8 @@ Type objective_function<Type>::operator() ()
     XtWY += XiTilde.transpose() * YiTilde;
     vector<Type> LiDiag = Li.diagonal();
     sum_log_det += sum(log(LiDiag));
-
+    // Cache the reciprocal of square root of diagonal of covariance
+    diag_cov_inv_sqrt.segment(start_i, n_visits_i) = vector<Type>(tcrossprod(Li).diagonal()).rsqrt();
     // Save stuff.
     x_mat_tilde.block(start_i, 0, n_visits_i, x_matrix.cols()) = XiTilde;
     y_vec_tilde.segment(start_i, n_visits_i) = YiTilde.col(0);
@@ -146,7 +146,10 @@ Type objective_function<Type>::operator() ()
   Identity.setIdentity();
   matrix<Type> beta_vcov = XtWX_decomposition.solve(Identity);
   REPORT(beta_vcov);
-
+  // normalized residual
+  REPORT(epsilonTilde);
+  // inverse square root of diagonal of covariance
+  REPORT(diag_cov_inv_sqrt);
   matrix<Type> covariance_lower_chol = get_chol_and_clean(chols_by_group, is_spatial, n_visits);
   REPORT(covariance_lower_chol);
 

tests/testthat/test-tmb-methods.R

@@ -261,25 +261,31 @@ test_that("residuals works as expected with a model using spatial covariance str
   result_resp <- expect_silent(residuals(object, type = "response"))
   expect_double(result_resp, len = length(object$tmb_data$y_vector))
   expect_equal(head(result_resp, 5), c(-4.5428, -24.0301, -8.8329, -3.4092, 8.5200), tolerance = 1e-4)
-
-  expect_error(residuals(object, type = "pearson"))
-  expect_error(residuals(object, type = "normalized"))
+  result_normal <- residuals(object, type = "normalized")
+  expect_equal(head(result_normal, 5), c(-0.4943, -2.5698, -0.9613, 0.0046, 1.1645), tolerance = 1e-4)
+  result_pearson <- residuals(object, type = "normalized")
+  expect_equal(head(result_pearson, 5), c(-0.4944, -2.5699, -0.9613, 0.0046, 1.1645), tolerance = 1e-4)
 })
 
 test_that("pearson residuals helper function works as expected", {
   object <- get_mmrm()
   resid_response <- residuals(object, type = "response")
-
-  result_pearson <- h_residuals_pearson(object, resid_response)
+  result_pearson <- h_residuals_pearson(object)
   expect_double(result_pearson, len = length(object$tmb_data$y_vector))
   expect_equal(tail(result_pearson, 5), c(2.22057, 1.79050, 0.53322, 0.87243, 0.70477), tolerance = 1e-4)
 })
 
 test_that("normalized residuals helper function works as expected", {
   object <- get_mmrm()
-  resid_response <- residuals(object, type = "response")
-
-  result_norm <- h_residuals_normalized(object, resid_response)
+  result_norm <- h_residuals_normalized(object)
   expect_double(result_norm, len = length(object$tmb_data$y_vector))
   expect_equal(tail(result_norm, 5), c(1.99092, 1.49689, 0.53322, 0.71055, 0.56152), tolerance = 1e-4)
 })
+
+test_that("response residuals helper function works as expected", {
+  object <- get_mmrm()
+  result_rsp <- h_residuals_response(object)
+  expect_double(result_rsp, len = length(object$tmb_data$y_vector))
+  expected <- component(object, "y_vector") - as.vector(fitted(object))
+  expect_equal(result_rsp, expected)
+})