add Gaussian spatial covariance

.gitignore

@@ -47,3 +47,8 @@ renv.lock
 *.pyc
 nohup.out
 inst/doc
+*xpt
+.python-version
+pyproject.toml
+sascfg_personal.py
+uv.lock

NEWS.md

@@ -8,6 +8,7 @@
 ### New Features
 
 - `mmrm` now returns score per subject in empirical covariance. It can be accessed by `component(obj, name = "score_per_subject")`. 
+- `mmrm` now support spatial Gaussian covariance matrix.
 
 # mmrm 0.3.14
 

R/cov_struct.R

@@ -40,7 +40,8 @@ COV_TYPES <- local({ # nolint
       type("auto-regressive order one", "ar1", "ar1h", TRUE, FALSE),
       type("ante-dependence", "ad", "adh", TRUE, FALSE),
       type("compound symmetry", "cs", "csh", TRUE, FALSE),
-      type("spatial exponential", "sp_exp", NA, FALSE, TRUE)
+      type("spatial exponential", "sp_exp", NA, FALSE, TRUE),
+      type("spatial Gaussian", "sp_gau", NA, FALSE, TRUE)
     )
   )
 })
@@ -150,7 +151,10 @@ COV_TYPES <- local({ # nolint
 #' \tab spatial exponential
 #' \tab \eqn{2}
 #' \tab \eqn{\sigma^{2}\rho^{-d_{ij}}}
-#'
+#' sp_gau
+#' \tab spatial Guassian
+#' \tab \eqn{2}
+#' \tab \eqn{\sigma^{2}\rho^{-d_{ij}^2}}
 #' }
 #'
 #' where \eqn{d_{ij}} denotes the Euclidean distance between time points
@@ -378,7 +382,7 @@ print.cov_struct <- function(x, ...) {
 #' sp_exp(coord1, coord2 | group / subject)
 #' ```
 #'
-#' Note that only `sp_exp` (spatial) covariance structures may provide multiple
+#' Note that only spatial covariance structures may provide multiple
 #' coordinates, which identify the Euclidean distance between the time points.
 #'
 #' @param x an object from which to derive a covariance structure. See object

R/mmrm-methods.R

@@ -162,7 +162,8 @@ h_print_cov <- function(cov_type, n_theta, n_groups) {
     adh = "heterogeneous ante-dependence",
     cs = "compound symmetry",
     csh = "heterogeneous compound symmetry",
-    sp_exp = "spatial exponential"
+    sp_exp = "spatial exponential",
+    sp_gau = "spatial Gaussian"
   )
 
   catstr <- sprintf(

R/tmb.R

@@ -35,7 +35,7 @@ h_mmrm_tmb_formula_parts <- function(
       model_formula = model_formula,
       full_formula = h_add_covariance_terms(model_formula, covariance),
       cov_type = tmb_cov_type(covariance),
-      is_spatial = covariance$type == "sp_exp",
+      is_spatial = grepl("sp_", covariance$type),
       visit_var = covariance$visits,
       subject_var = covariance$subject,
       group_var = if (length(covariance$group) < 1) NULL else covariance$group,

R/utils.R

@@ -377,7 +377,8 @@ std_start <- function(cov_type, n_visits, n_groups, ...) {
     adh = rep(0, 2 * n_visits - 1),
     cs = rep(0, 2),
     csh = rep(0, n_visits + 1),
-    sp_exp = rep(0, 2)
+    sp_exp = rep(0, 2),
+    sp_gau = rep(0, 2)
   )
   rep(start_value, n_groups)
 }

design/SAS/sas_sp_gau.Rmd

@@ -0,0 +1,31 @@
+---
+title: Spatial Gaussian
+---
+
+
+```{sas}
+libname xptfile xport "fev.xpt";
+proc copy inlib = xptfile outlib = work;
+run;
+
+ods output CovParms = covparm Rcorr = rcorr SolutionF = fixed FitStatistics= fit;
+PROC MIXED DATA = fev cl method=ml;
+  CLASS USUBJID;
+  MODEL FEV1 = / ddfm=satterthwaite solution chisq;
+  REPEATED / subject=USUBJID type=sp(gau)(visitn) rcorr;
+RUN;
+
+libname xptout xport "out.xpt";
+
+proc copy inlib = work outlib = xptout memtype = data;
+  select covparm rcorr fixed fit;
+run;
+```
+
+
+```{r}
+fit <- mmrm(FEV1 ~ sp_gau(VISITN | USUBJID), data = fev_data, reml = FALSE, method = "Satterthwaite")
+summary(fit)
+
+# degree of freedom seems not working right now.
+```

src/covariance.h

@@ -129,14 +129,26 @@ matrix<T> get_compound_symmetry_heterogeneous(const vector<T>& theta, int n_visi
 // Spatial Exponential Cholesky factor.
 template <class T>
 matrix<T> get_spatial_exponential(const vector<T>& theta, const matrix<T>& distance) {
+  T const_sd = exp(theta(0));
+  T rho = invlogit(theta(1)); // theta (-Inf, Inf), rho (0, 1), log(rho) (-Inf, 0)
+  matrix<T> expdist = exp(distance.array() * log(rho)); 
+  matrix<T> result = expdist * const_sd;
+  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > cov_i_chol(result);
+  return cov_i_chol.matrixL();
+}
+
+// Spatial Gaussian Cholesky factor.
+template <class T>
+matrix<T> get_spatial_gaussian(const vector<T>& theta, const matrix<T>& distance) {
   T const_sd = exp(theta(0));
   T rho = invlogit(theta(1));
-  matrix<T> expdist = exp(distance.array() * log(rho));
+  matrix<T> expdist = exp(distance.array().square() * log(rho)); // only take square of the dist
   matrix<T> result = expdist * const_sd;
   Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > cov_i_chol(result);
   return cov_i_chol.matrixL();
 }
 
+
 // Creates a new correlation object dynamically.
 template <class T>
 matrix<T> get_covariance_lower_chol(const vector<T>& theta, int n_visits, std::string cov_type) {
@@ -173,7 +185,9 @@ matrix<T> get_spatial_covariance_lower_chol(const vector<T>& theta, const matrix
   matrix<T> result;
   if (cov_type == "sp_exp") {
     result = get_spatial_exponential<T>(theta, distance);
-  } else {
+  } else if (cov_type == "sp_gau") {
+    result = get_spatial_gaussian<T>(theta, distance);
+  }else {
     Rf_error("%s", ("Unknown spatial covariance type '" + cov_type + "'.").c_str());
   }
   return result;

src/derivatives.h

@@ -221,4 +221,45 @@ struct derivatives_sp_exp: public lower_chol_spatial<Type>, virtual derivatives_
   }
 };
 
+
+// derivatives_sp_gau struct is created to obtain the exact derivatives of spatial Guassian
+// covariance structure, and its inverse.
+// Inherit from sp_exp because of similar parameterization
+template <class Type>
+struct derivatives_sp_gau: public derivatives_sp_exp<Type> {
+  Type const_sd;
+  Type rho;
+  Type logrho;
+  derivatives_sp_gau() {
+    // This default constructor is needed because the use of `[]` in maps.
+  }
+  // Initialize the theta values; the reason to have theta is that for a fit, the theta
+  // is the same for all subjects, while the distance between each visits for each subject
+  // can be different.
+  derivatives_sp_gau(vector<Type> theta, std::string cov_type): lower_chol_spatial<Type>(theta, cov_type) ,const_sd(exp(theta(0))), rho(invlogit(theta(1))) {
+    this->logrho = log(this->rho);
+  }
+  // Obtain first order derivatives
+  matrix<Type> get_sigma_derivative1(std::vector<int> visits, matrix<Type> dist) override {
+    matrix<Type> ret(2 * dist.rows(), dist.cols());
+    // partial sigma / partial theta_1 = sigma.
+    auto sigma = this->get_sigma(visits, dist);
+    ret.block(0, 0, dist.rows(), dist.cols()) = sigma;
+    ret.block(dist.rows(), 0, dist.rows(), dist.cols()) = sigma.array() * dist.array().suqare() * (1 - this->rho);
+    return ret;
+  }
+  // Obtain second order derivatives.
+  matrix<Type> get_sigma_derivative2(std::vector<int> visits, matrix<Type> dist) override {
+    matrix<Type> ret(4 * dist.rows(), dist.cols());
+    auto sigma = this->get_sigma(visits, dist);
+    ret.block(0, 0, dist.rows(), dist.cols()) = sigma;
+    Type rho_r = 1 - this->rho;
+    auto dtheta1dtheta2 = sigma.array() * dist.array().suqare() * rho_r;
+    ret.block(dist.rows(), 0, dist.rows(), dist.cols()) =  dtheta1dtheta2;
+    ret.block(dist.rows() * 2, 0, dist.rows(), dist.cols()) = dtheta1dtheta2;
+    matrix<Type> dtheta2s = dtheta1dtheta2 * (dist.array().suqare() * rho_r - this->rho);
+    ret.block(dist.rows() * 3, 0, dist.rows(), dist.cols()) = dtheta2s;
+    return ret;
+  }
+};
 #endif

tests/testthat/test-fit.R

@@ -71,7 +71,7 @@ test_that("fit_single_optimizer gives error messages", {
     paste(
       "Covariance structure must be specified in formula.",
       "Possible covariance structures include:",
-      "us, toep, toeph, ar1, ar1h, ad, adh, cs, csh, sp_exp"
+      "us, toep, toeph, ar1, ar1h, ad, adh, cs, csh, sp_exp, sp_gau"
     ),
     fixed = TRUE
   )

tests/testthat/test-tmb.R

@@ -1847,6 +1847,20 @@ test_that("fit_mmrm works with sp_exp covariance structure and REML(2-dimension)
   expect_equal(exp(result$theta_est[1]), 87.8870, tolerance = 1e-3)
 })
 
+
+## spatial gaussian ----
+
+test_that("fit_mmrm works with sp_exp covariance structure and ML", {
+  formula <- FEV1 ~ sp_gau(VISITN | USUBJID)
+  result <- fit_mmrm(formula, fev_data, weights = rep(1, nrow(fev_data)), reml = FALSE)
+  expect_class(result, "mmrm_tmb")
+  # See design/SAS/sas_sp_exp_ml.txt for the source of numbers.
+  expect_equal(deviance(result), 3871.058334)
+  expect_equal(as.numeric(result$beta_est[1]), 42.35391569, tolerance = 1e-4)
+  expect_equal(sqrt(-1 / plogis(result$theta_est[2], log.p = TRUE)), 1.08533042962185, tolerance = 1e-3)
+  expect_equal(exp(result$theta_est[1]), 88.48735528483263, tolerance = 1e-3)
+})
+
 ## misc ----
 
 test_that("fit_mmrm also works with character ID variable", {

tests/testthat/test-utils.R

@@ -231,6 +231,10 @@ test_that("std_start works", {
     std_start("sp_exp", 5, 5),
     rep(0, 10)
   )
+  expect_identical(
+    std_start("sp_gau", 5, 5),
+    rep(0, 10)
+  )
 })
 
 # h_get_theta_from_cov ----

vignettes/covariance.Rmd

@@ -239,4 +239,16 @@ Hence we have the following parametrization form.
   \theta = (\log(\sigma), \logit(\rho))
 \]
 
+## Spatial Gaussian (`sp_gau`)
+
+For spatial Gaussian, the covariance structure is defined as follows:
+
+\[
+  \rho_{ij} = \rho^{d_{ij}^2}
+\]
+
+where $d_{ij}$ is the distance between time point $i$ and time point $j$,
+
+The parameterization is the same as spatial exponential.
+
 # References