adding option for constant marginal variance ...

... which is still experimental

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: dsem
 Type: Package
 Title: Fit Dynamic Structural Equation Models
-Version: 1.2.1
-Date: 2024-04-02
+Version: 1.3.0
+Date: 2024-07-18
 Authors@R: 
     c(person(given = "James",
              family = "Thorson",
@@ -45,7 +45,7 @@ Description: Applies dynamic structural equation models to time-series data
 	constrained by ecological mechanisms."  
 License: GPL-3
 Encoding: UTF-8
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 LazyData: true
 URL: https://james-thorson-noaa.github.io/dsem/

NEWS.md

@@ -1,3 +1,9 @@
+# dsem 1.3.0
+
+* Adding option to specify constant marginal variance, as alternative to existing
+  calculation which results in a constant conditional variance but a changing marginal 
+  variance along the causal graph
+
 # dsem 1.2.1
 
 * removing `checkDepPackageVersion(dep_pkg="Matrix", this_pkg="TMB")` from `.onLoad()`

R/dsem.R

@@ -160,7 +160,13 @@ function( sem,
   }
 
   #
-  Data = list( "options" = ifelse(control$gmrf_parameterization=="separable",0,1),
+  options = c(
+    ifelse(control$gmrf_parameterization=="separable", 0, 1),
+    ifelse(control$constant_variance=="conditional", 0, 1)
+  )
+
+  #
+  Data = list( "options" = options,
                "RAM" = as.matrix(na.omit(ram[,1:4])),
                "RAMstart" = as.numeric(ram[,5]),
                "familycode_j" = sapply(family, FUN=switch, "fixed"=0, "normal"=1, "gamma"=4 ),
@@ -337,6 +343,10 @@ function( sem,
 #'        a full-rank and IID precision for variables over time, and then projects
 #'        this using the inverse-cholesky of the precision, where this projection
 #'        can be rank-deficient.
+#' @param constant_variance Whether to specify a constant conditional variance 
+#'        \eqn{ \mathbf{Gamma Gamma}^t}, which results in a changing marginal variance      
+#'        along the specified causal graph, or instead specify a constant marginal variance
+#'        such that \eqn{ \mathbf{Gamma}} is rescaled to achieve this constraint.
 #' @param quiet Boolean indicating whether to run model printing messages to terminal or not;
 #' @param use_REML Boolean indicating whether to treat non-variance fixed effects as random,
 #'        either to motigate bias in estimated variance parameters or improve efficiency for
@@ -364,7 +374,8 @@ function( nlminb_loops = 1,
           getsd = TRUE,
           quiet = FALSE,
           run_model = TRUE,
-          gmrf_parameterization = c("separable","projection"),
+          gmrf_parameterization = c("separable", "projection"),
+          constant_variance = c("conditional", "marginal"),
           use_REML = TRUE,
           profile = NULL,
           parameters = NULL,
@@ -373,6 +384,7 @@ function( nlminb_loops = 1,
           extra_convergence_checks = TRUE ){
 
   gmrf_parameterization = match.arg(gmrf_parameterization)
+  constant_variance = match.arg(constant_variance)
 
   # Return
   structure( list(
@@ -385,6 +397,7 @@ function( nlminb_loops = 1,
     quiet = quiet,
     run_model = run_model,
     gmrf_parameterization = gmrf_parameterization,
+    constant_variance = constant_variance,
     use_REML = use_REML,
     profile = profile,
     parameters = parameters,

src/dsem.cpp

@@ -16,7 +16,9 @@ Type objective_function<Type>::operator() ()
   using namespace density;
 
   // Data
-  DATA_IVECTOR( options ); // options(0) -> full rank or rank-reduced GMRF
+  DATA_IVECTOR( options ); 
+  // options(0) -> 0: full rank;  1: rank-reduced GMRF
+  // options(1) -> 0: constant conditional variance;  1: constant marginal variance
   //DATA_INTEGER( resimulate_gmrf );
   DATA_IMATRIX( RAM );
   DATA_VECTOR( RAMstart );
@@ -68,11 +70,49 @@ Type objective_function<Type>::operator() ()
     }
     //if(RAM(r,0)==2) Gammainv_kk.coeffRef( RAM(r,1)-1, RAM(r,2)-1 ) = 1 / tmp;
   }
-
-  // solve(I - Rho) %*% x
   Eigen::SparseMatrix<Type> IminusRho_kk = I_kk - Rho_kk;
+
+  // Compute inverse
   Eigen::SparseLU< Eigen::SparseMatrix<Type>, Eigen::COLAMDOrdering<int> > inverseIminusRho_kk;
   inverseIminusRho_kk.compute(IminusRho_kk);
+
+  // Rescale I-Rho if using constant marginal variance
+  if( options(1)==1 ){
+    //matrix<Type> invIminusRho_kk = atomic::matinv(); // Requires dense Gamma and Rho from above
+    Eigen::SparseMatrix<Type> invIminusRho_kk;
+
+    // Doesn't work
+    //invIminusRho_kk = IminusRho_kk.inverse();
+    // WORKS:  Based on: https://github.com/kaskr/adcomp/issues/74
+    invIminusRho_kk = inverseIminusRho_kk.solve(I_kk);
+    REPORT( invIminusRho_kk );
+
+    // Hadamard squared
+    // See: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html
+    Eigen::SparseMatrix<Type> squared_invIminusRho_kk(n_k, n_k);
+    squared_invIminusRho_kk = invIminusRho_kk.cwiseProduct(invIminusRho_kk);
+    REPORT( squared_invIminusRho_kk );
+
+    // Rowsums
+    Eigen::SparseLU< Eigen::SparseMatrix<Type>, Eigen::COLAMDOrdering<int> > invsquared_invIminusRho_kk;
+    invsquared_invIminusRho_kk.compute(squared_invIminusRho_kk);
+    matrix<Type> ones_k1( n_k, 1 );
+    ones_k1.setOnes();
+    matrix<Type> margvar_k1 = invsquared_invIminusRho_kk.solve(ones_k1);
+    REPORT( ones_k1 );
+    REPORT( margvar_k1 );
+
+    // Operations
+    Eigen::SparseMatrix<Type> invmargsd_kk(n_k, n_k);
+    for( int k=0; k<n_k; k++ ){
+      invmargsd_kk.coeffRef(k,k) = pow( margvar_k1(k,0), -0.5 );
+    }
+    IminusRho_kk = invmargsd_kk * IminusRho_kk;
+    REPORT( invmargsd_kk );
+
+    // Recompute inverse
+    inverseIminusRho_kk.compute(IminusRho_kk);
+  }
 
   // Calculate effect of initial condition -- SPARSE version
   vector<Type> delta_k( n_k );