return meat per subject when computing empirical covariance matrix of coefficients

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: mmrm
 Title: Mixed Models for Repeated Measures
-Version: 0.3.14.9001
+Version: 0.3.15
 Authors@R: c(
     person("Daniel", "Sabanes Bove", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-0176-9239")),

NEWS.md

@@ -1,3 +1,9 @@
+# mmrm 0.3.15
+
+### New Features
+
+- `mmrm` now returns score per subject in empirical covariance. It can be accessed by `component(obj, name = "score_per_subject")`. 
+
 # mmrm 0.3.14.9001
 
 ### Bug Fixes

R/component.R

@@ -31,6 +31,8 @@
 #' - `varcor`: estimated covariance matrix for residuals. If there are multiple
 #'      groups, a named list of estimated covariance matrices for residuals will be
 #'      returned. The names are the group levels.
+#' - `score_per_subject`: score per subject in empirical covariance.
+#'      See the vignette \code{vignette("coef_vcov", package = "mmrm")}.
 #' - `theta_est`: estimated variance parameters.
 #' - `beta_est`: estimated coefficients (excluding aliased coefficients).
 #' - `beta_est_complete`: estimated coefficients including aliased coefficients
@@ -64,7 +66,7 @@ component <- function(object,
                       name = c(
                         "cov_type", "subject_var", "n_theta", "n_subjects", "n_timepoints",
                         "n_obs", "beta_vcov", "beta_vcov_complete",
-                        "varcor", "formula", "dataset", "n_groups",
+                        "varcor", "score_per_subject", "formula", "dataset", "n_groups",
                         "reml", "convergence", "evaluations", "method", "optimizer",
                         "conv_message", "call", "theta_est",
                         "beta_est", "beta_est_complete", "beta_aliased",
@@ -132,6 +134,7 @@ component <- function(object,
         object$beta_vcov
       },
     "varcor" = object$cov,
+    "meat_per_subject" = object$meat_per_subject,
     "x_matrix" = object$tmb_data$x_matrix,
     "xlev" = stats::.getXlevels(terms(object), object$tmb_data$full_frame),
     "contrasts" = attr(object$tmb_data$x_matrix, "contrasts"),

R/fit.R

@@ -501,7 +501,7 @@ mmrm <- function(formula,
     )
     fit$beta_vcov_adj <- empirical_comp$cov
     fit$empirical_df_mat <- empirical_comp$df_mat
-    fit$meat_per_subject <- empirical_comp$meat_per_subject
+    fit$score_per_subject <- empirical_comp$score_per_subject
     dimnames(fit$beta_vcov_adj) <- dimnames(fit$beta_vcov)
   } else if (identical(control$vcov, "Asymptotic")) {
     # Note that we only need the Jacobian list under Asymptotic covariance method,

src/empirical.cpp

@@ -29,7 +29,7 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
   vector<double> G_sqrt = as_num_vector_tmb(sqrt(weights_vector));
   int p = x.cols();
 
-  NumericMatrix meat_per_subject(n_subjects, p);
+  matrix<double> score_per_subject = matrix<double>::Zero(n_subjects, p);
 
   // Use map to hold these base class pointers (can also work for child class objects).
   auto derivatives_by_group = cache_obj<double, derivatives_base<double>, derivatives_sp_exp<double>, derivatives_nonspatial<double>>(theta_v, n_groups, is_spatial, cov_type, n_visits);
@@ -70,9 +70,7 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
     xt_g_simga_inv_chol.block(0, start_i, p, n_visits_i) = xt_gi_simga_inv_chol;
     ax.block(start_i, 0, n_visits_i, p) = xta.transpose();
 
-    for(int j = 0; j < p; ++j){
-      meat_per_subject(i, j) = z(j, 0);
-    }
+    score_per_subject.row(i) = z.transpose();
   }
   matrix<double> h = xt_g_simga_inv_chol.transpose() * beta_vcov_matrix * xt_g_simga_inv_chol;
   matrix<double> imh = matrix<double>::Identity(n_observations, n_observations) - h;
@@ -92,7 +90,7 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
   //  ret = ret * (n_subjects - 1) / n_subjects;
   //}
   return List::create(
-    Named("meat_per_subject") = meat_per_subject,
+    Named("score_per_subject") = as_num_matrix_rcpp(score_per_subject),
     Named("cov") = as_num_matrix_rcpp(ret),
     Named("df_mat") = as_num_matrix_rcpp(gtvg)
   );

tests/testthat/test-component.R

@@ -9,7 +9,8 @@ test_that("component works as expected for mmrm_tmb objects", {
     names(component(object_mmrm_tmb)),
     c(
       "cov_type", "subject_var", "n_theta", "n_subjects", "n_timepoints",
-      "n_obs", "beta_vcov", "beta_vcov_complete", "varcor", "formula", "dataset",
+      "n_obs", "beta_vcov", "beta_vcov_complete", "varcor", "score_per_subject",
+      "formula", "dataset",
       "n_groups", "reml", "convergence", "evaluations", "method", "optimizer",
       "conv_message", "call", "theta_est",
       "beta_est", "beta_est_complete", "beta_aliased",

tests/testthat/test-empirical.R

@@ -37,7 +37,7 @@ test_that("empirical covariance are the same with SAS result for ar1", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -55,7 +55,7 @@ test_that("empirical covariance are the same with SAS result for ar1h", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -73,7 +73,7 @@ test_that("empirical covariance are the same with SAS result for cs", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -91,7 +91,7 @@ test_that("empirical covariance are the same with SAS result for csh", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -109,7 +109,7 @@ test_that("empirical covariance are the same with SAS result for toep", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -127,7 +127,7 @@ test_that("empirical covariance are the same with SAS result for toeph", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -145,7 +145,7 @@ test_that("empirical covariance are the same with SAS result for adh", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -163,7 +163,7 @@ test_that("empirical covariance are the same with SAS result for us", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -184,7 +184,7 @@ test_that("empirical covariance are the same with SAS result for sp_exp", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -209,7 +209,7 @@ test_that("empirical covariance are the same with SAS result for ar1", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -230,7 +230,7 @@ test_that("empirical covariance are the same with SAS result for ar1h", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -251,7 +251,7 @@ test_that("empirical covariance are the same with SAS result for cs", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -272,7 +272,7 @@ test_that("empirical covariance are the same with SAS result for csh", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -293,7 +293,7 @@ test_that("empirical covariance are the same with SAS result for toep", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -314,7 +314,7 @@ test_that("empirical covariance are the same with SAS result for toeph", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -335,7 +335,7 @@ test_that("empirical covariance are the same with SAS result for adh", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -356,7 +356,7 @@ test_that("empirical covariance are the same with SAS result for us", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4
@@ -377,7 +377,7 @@ test_that("empirical covariance are the same with SAS result for sp_exp", {
 
   expect_equal(
     fit$beta_vcov %*%
-      (t(fit$meat_per_subject) %*% fit$meat_per_subject) %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
       t(fit$beta_vcov),
     fit$beta_vcov_adj,
     tolerance = 1e-4

vignettes/coef_vcov.Rmd

@@ -58,7 +58,9 @@ covariance matrix.
 \]
 
 Where $W_i$ is the block diagonal part for subject $i$ of $W$ matrix, $\hat\epsilon_i$ is the observed residuals for subject i,
-$L_i$ is the Cholesky factor of $\Sigma_i^{-1}$ ($W_i = L_i L_i^\top$).
+$L_i$ is the Cholesky factor of $\Sigma_i^{-1}$ ($W_i = L_i L_i^\top$). 
+In the sandwich, the score $X_i^\top W_i \hat\epsilon_i$ computed for 
+subject $i$ can be accessed by `component(mmrm_obj, name = "score_per_subject")`. 
 
 See the detailed explanation of these formulas in the [Weighted Least Square Empirical Covariance](empirical_wls.html) vignette.