diff --git a/NEWS.md b/NEWS.md
index a676cf671..5127ab3be 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,10 @@
 
 - Previously, when compiling `mmrm` from source using a `TMB` version below 1.9.15, and installing a newer `TMB` of version 1.9.15 or above, would render the `mmrm` package unusable. This is fixed now, by checking in the dynamic library of `mmrm` whether the version of `TMB` has been sufficient.
 
+### 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
 
 ### Bug Fixes
diff --git a/R/component.R b/R/component.R
index b8cebfc89..5112b0146 100644
--- a/R/component.R
+++ b/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,
+    "score_per_subject" = object$score_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"),
diff --git a/R/fit.R b/R/fit.R
index a738a813d..836cb6b8e 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -501,6 +501,7 @@ mmrm <- function(formula,
     )
     fit$beta_vcov_adj <- empirical_comp$cov
     fit$empirical_df_mat <- empirical_comp$df_mat
+    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,
diff --git a/man/component.Rd b/man/component.Rd
index 7a061712b..436aabd5e 100644
--- a/man/component.Rd
+++ b/man/component.Rd
@@ -7,10 +7,11 @@
 component(
   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",
-    "reml", "convergence", "evaluations", "method", "optimizer", "conv_message", "call",
-    "theta_est", "beta_est", "beta_est_complete", "beta_aliased", "x_matrix", "y_vector",
-    "neg_log_lik", "jac_list", "theta_vcov", "full_frame", "xlev", "contrasts")
+    "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",
+    "x_matrix", "y_vector", "neg_log_lik", "jac_list", "theta_vcov", "full_frame",
+    "xlev", "contrasts")
 )
 }
 \arguments{
@@ -51,6 +52,8 @@ aliased coefficients with entries set to \code{NA}.
 \item \code{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.
+\item \code{score_per_subject}: score per subject in empirical covariance.
+See the vignette \code{vignette("coef_vcov", package = "mmrm")}.
 \item \code{theta_est}: estimated variance parameters.
 \item \code{beta_est}: estimated coefficients (excluding aliased coefficients).
 \item \code{beta_est_complete}: estimated coefficients including aliased coefficients
diff --git a/src/empirical.cpp b/src/empirical.cpp
index a3bb83902..69c4f68d8 100644
--- a/src/empirical.cpp
+++ b/src/empirical.cpp
@@ -28,6 +28,9 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
   matrix<double> residual = y_matrix - fitted;
   vector<double> G_sqrt = as_num_vector_tmb(sqrt(weights_vector));
   int p = x.cols();
+
+  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);
   matrix<double> meat = matrix<double>::Zero(p, p);
@@ -66,6 +69,8 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
     meat = meat + z * z.transpose();
     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();
+
+    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;
@@ -85,6 +90,7 @@ List get_empirical(List mmrm_data, NumericVector theta, NumericVector beta, Nume
   //  ret = ret * (n_subjects - 1) / n_subjects;
   //}
   return List::create(
+    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)
   );
diff --git a/tests/testthat/_snaps/tmb-methods.md b/tests/testthat/_snaps/tmb-methods.md
index 6f1f3a19c..a05bf8c3c 100644
--- a/tests/testthat/_snaps/tmb-methods.md
+++ b/tests/testthat/_snaps/tmb-methods.md
@@ -1,3 +1,8 @@
+# predict works for data different factor levels
+
+    c("1" = 36.4082672191324, "2" = 41.2059294235457, "3" = 46.0566947323669
+    )
+
 # predict works for unconditional prediction if response does not exist
 
     c("1" = 36.4082672191322, "2" = 41.2059294235456, "3" = 46.0566947323668, 
diff --git a/tests/testthat/test-component.R b/tests/testthat/test-component.R
index afa6e2c27..3472da023 100644
--- a/tests/testthat/test-component.R
+++ b/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",
diff --git a/tests/testthat/test-empirical.R b/tests/testthat/test-empirical.R
index d98e71b0d..44974f21c 100644
--- a/tests/testthat/test-empirical.R
+++ b/tests/testthat/test-empirical.R
@@ -34,6 +34,14 @@ test_that("empirical covariance are the same with SAS result for ar1", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for ar1h", {
@@ -44,6 +52,14 @@ test_that("empirical covariance are the same with SAS result for ar1h", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for cs", {
@@ -54,6 +70,14 @@ test_that("empirical covariance are the same with SAS result for cs", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for csh", {
@@ -64,6 +88,14 @@ test_that("empirical covariance are the same with SAS result for csh", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for toep", {
@@ -74,6 +106,14 @@ test_that("empirical covariance are the same with SAS result for toep", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for toeph", {
@@ -84,6 +124,14 @@ test_that("empirical covariance are the same with SAS result for toeph", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for adh", {
@@ -94,6 +142,14 @@ test_that("empirical covariance are the same with SAS result for adh", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for us", {
@@ -104,6 +160,14 @@ test_that("empirical covariance are the same with SAS result for us", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for sp_exp", {
@@ -117,6 +181,14 @@ test_that("empirical covariance are the same with SAS result for sp_exp", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 # weighted mmrm ----
@@ -134,6 +206,14 @@ test_that("empirical covariance are the same with SAS result for ar1", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for ar1h", {
@@ -147,6 +227,14 @@ test_that("empirical covariance are the same with SAS result for ar1h", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for cs", {
@@ -160,6 +248,14 @@ test_that("empirical covariance are the same with SAS result for cs", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for csh", {
@@ -173,6 +269,14 @@ test_that("empirical covariance are the same with SAS result for csh", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for toep", {
@@ -186,6 +290,14 @@ test_that("empirical covariance are the same with SAS result for toep", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for toeph", {
@@ -199,6 +311,14 @@ test_that("empirical covariance are the same with SAS result for toeph", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for adh", {
@@ -212,6 +332,14 @@ test_that("empirical covariance are the same with SAS result for adh", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for us", {
@@ -225,6 +353,14 @@ test_that("empirical covariance are the same with SAS result for us", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 test_that("empirical covariance are the same with SAS result for sp_exp", {
@@ -238,6 +374,14 @@ test_that("empirical covariance are the same with SAS result for sp_exp", {
     dimnames = rep(list(c("(Intercept)", "ARMCDTRT")), 2)
   )
   expect_equal(component(fit, "beta_vcov"), expected, tolerance = 1e-4)
+
+  expect_equal(
+    fit$beta_vcov %*%
+      (t(fit$score_per_subject) %*% fit$score_per_subject) %*%
+      t(fit$beta_vcov),
+    fit$beta_vcov_adj,
+    tolerance = 1e-4
+  )
 })
 
 ## Empirical Satterthwaite vs gls/clubSandwich ----
diff --git a/vignettes/coef_vcov.Rmd b/vignettes/coef_vcov.Rmd
index 5709137d6..544aafa52 100644
--- a/vignettes/coef_vcov.Rmd
+++ b/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.