updates for merge

James-Thorson-NOAA · Jan 10, 2025 · e143457 · e143457
1 parent 6c349a5
commit e143457
Show file tree

Hide file tree

Showing 10 changed files with 105 additions and 25 deletions.
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
@@ -15,7 +15,7 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - uses: r-lib/actions/setup-r@v2
         with:
@@ -56,7 +56,7 @@ jobs:
 
       - name: Upload test results
         if: failure()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: coverage-test-failures
           path: ${{ runner.temp }}/package
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: dsem
 Type: Package
 Title: Fit Dynamic Structural Equation Models
 Version: 1.4.0
-Date: 2024-10-07
+Date: 2025-01-10
 Authors@R: 
     c(person(given = "James",
              family = "Thorson",
@@ -14,8 +14,9 @@ Imports:
     Matrix,
     sem,
     igraph,
-    RTMB,
+    RTMB (>= 1.7.0),
     ggraph,
+    ggplot2,
     grid,
     methods
 Depends:
@@ -29,7 +30,6 @@ Suggests:
     gridExtra,
     dynlm,
     MARSS,
-    ggplot2,
     ggpubr,
     grid,
     vars,

diff --git a/NAMESPACE b/NAMESPACE
@@ -34,11 +34,17 @@ importFrom(Matrix,solve)
 importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,t)
 importFrom(Matrix,tcrossprod)
+importFrom(RTMB,AD)
+importFrom(RTMB,ADREPORT)
+importFrom(RTMB,ADoverload)
+importFrom(RTMB,REPORT)
+importFrom(RTMB,dgmrf)
 importFrom(TMB,MakeADFun)
 importFrom(TMB,compile)
 importFrom(TMB,dynlib)
 importFrom(TMB,sdreport)
 importFrom(TMB,summary.sdreport)
+importFrom(ggplot2,aes)
 importFrom(ggraph,create_layout)
 importFrom(ggraph,geom_edge_arc)
 importFrom(ggraph,geom_node_text)
@@ -60,12 +66,12 @@ importFrom(stats,na.omit)
 importFrom(stats,nlminb)
 importFrom(stats,optimHess)
 importFrom(stats,pnorm)
+importFrom(stats,rbinom)
+importFrom(stats,rgamma)
 importFrom(stats,rnorm)
+importFrom(stats,rpois)
 importFrom(stats,sd)
 importFrom(stats,simulate)
 importFrom(stats,time)
 importFrom(stats,vcov)
 useDynLib(dsem, .registration = TRUE)
-
-# Edited by hand
-import("RTMB")
diff --git a/NEWS.md b/NEWS.md
@@ -7,7 +7,9 @@
 * Adding `convert_equations` to extend `sem::specifyEquations` and simplify
   specification for large models
 * Adding argument `prior_negloglike` as interface to specify Bayesian priors
-  and/or likelihood penalties
+  and/or likelihood penalties in `dsem`
+* Adding `dsemRTMB` using RTMB as interchangeable alternative to `dsem` using
+  TMB
 
 # dsem 1.3.0
 

diff --git a/R/compute_nll.R b/R/compute_nll.R
@@ -185,14 +185,14 @@ function( parlist,
     }
     # familycode = 2 :  binomial
     if( family[j]=="binomial" ){
-      mu_tj[t,j] = invlogit(z_tj[t,j]);
+      mu_tj[t,j] = plogis(z_tj[t,j]);
       if(!is.na(y_tj[t,j])){
-        loglik_tj[t,j] = dbinom( y_tj[t,j], Type(1.0), mu_tj[t,j], TRUE );
+        loglik_tj[t,j] = dbinom( y_tj[t,j], 1.0, mu_tj[t,j], TRUE );
         if( isTRUE(simulate_data) ){
           y_tj[t,j] = rbinom( n=1, size=1, prob=mu_tj[t,j] );
         }
       }
-      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(-2*(((1-y_tj[t,j])*log(1-mu_tj[t,j]) + y_tj[t,j]*log(mu_tj[t,j]))), 0.5);
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * sqrt( -2*((1-y_tj[t,j])*log(1-mu_tj[t,j]) + y_tj[t,j]*log(mu_tj[t,j])) );
     }
     # familycode = 3 :  Poisson
     if( family[j]=="poisson" ){
@@ -203,7 +203,7 @@ function( parlist,
           y_tj[t,j] = rpois( n=1, lambda=mu_tj[t,j] );
         }
       }
-      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(2*(y_tj[t,j]*log((Type(1e-10) + y_tj[t,j])/mu_tj[t,j]) - (y_tj[t,j]-mu_tj[t,j])), 0.5);
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * sqrt( 2*(y_tj[t,j]*log((1e-10 + y_tj[t,j])/mu_tj[t,j]) - (y_tj[t,j]-mu_tj[t,j])) );
     }
     # familycode = 4 :  Gamma:   shape = 1/CV^2; scale = mean*CV^2
     if( family[j]=="gamma" ){
@@ -214,7 +214,7 @@ function( parlist,
           y_tj[t,j] = rgamma( n=1, shape=pow(sigma_j[j],-2), scale=mu_tj[t,j]*pow(sigma_j[j],2) );
         }
       }
-      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(2 * ( (y_tj[t,j]-mu_tj[t,j])/mu_tj[t,j] - log(y_tj[t,j]/mu_tj[t,j]) ), 0.5);
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * sqrt( 2*((y_tj[t,j]-mu_tj[t,j])/mu_tj[t,j] - log(y_tj[t,j]/mu_tj[t,j])) );
     }
   }}
 
@@ -253,9 +253,3 @@ function( parlist,
   }
 }
 
-if( FALSE ){
-  family = rep("fixed",ncol(tsdata))
-  y_tj = tsdata
-  estimate_delta0 = FALSE
-}
-
diff --git a/R/dsem.R b/R/dsem.R
@@ -24,16 +24,21 @@
 #'        returns the negative log-prior probability. For example
 #'        \code{prior_negloglike = function(obj) -1 * dnorm( obj$par[1], mean=0, sd=0.1, log=TRUE)}
 #'        specifies a normal prior probability
-#'        for the for the first fixed effect with mean of zero and logsd of 0.1
+#'        for the for the first fixed effect with mean of zero and logsd of 0.1.
+#'        NOTE:  this implementation does not work well with \code{tmbstan} and
+#'        is highly experimental.  If using priors, considering using \code{\link{dsemRTMB}}
+#'        instead.  The option in \code{dsem} is mainly intended to validate its
+#'        use in \code{dsemRTMB}
 #' @param control Output from \code{\link{dsem_control}}, used to define user
 #'        settings, and see documentation for that function for details.
 #'
 #' @importFrom TMB compile dynlib MakeADFun sdreport summary.sdreport
-#' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rnorm simulate time tsp<-
+#' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rnorm simulate time tsp<- rbinom rgamma rpois
 #' @importFrom Matrix solve Cholesky sparseMatrix mat2triplet drop0
 #' @importFrom sem sem
 #' @importFrom igraph plot.igraph graph_from_data_frame with_sugiyama layout_
 #' @importFrom ggraph ggraph geom_edge_arc create_layout rectangle geom_node_text theme_graph
+#' @importFrom ggplot2 aes
 #' @importFrom grid arrow
 #' @importFrom methods is
 #'

diff --git a/R/dsemRTMB.R b/R/dsemRTMB.R
@@ -12,13 +12,44 @@
 #'        with mean of zero and sd of 0.1
 #'
 #' @importFrom Matrix solve Diagonal sparseMatrix drop0 kronecker crossprod tcrossprod t diag
+#' @importFrom RTMB ADoverload AD dgmrf REPORT ADREPORT
 #'
 #' @details
-#' See \code{\link{dsem}} for details
+#' \code{dsemRTMB} is interchangeable with \code{\link{dsem}}, but uses RTMB
+#' instead of TMB for estimation.  Both are provided for comparison and
+#' real-world comparison.
 #'
 #' @return
 #' An object (list) of class `dsem`, fitted using RTMB
 #'
+#' @examples
+#' # Define model
+#' sem = "
+#'   # Link, lag, param_name
+#'   cprofits -> consumption, 0, a1
+#'   cprofits -> consumption, 1, a2
+#'   pwage -> consumption, 0, a3
+#'   gwage -> consumption, 0, a3
+#'   cprofits -> invest, 0, b1
+#'   cprofits -> invest, 1, b2
+#'   capital -> invest, 0, b3
+#'   gnp -> pwage, 0, c2
+#'   gnp -> pwage, 1, c3
+#'   time -> pwage, 0, c1
+#' "
+#'
+#' # Load data
+#' data(KleinI, package="AER")
+#' TS = ts(data.frame(KleinI, "time"=time(KleinI) - 1931))
+#' tsdata = TS[,c("time","gnp","pwage","cprofits",'consumption',
+#'                "gwage","invest","capital")]
+#'
+#' # Fit model
+#' fit = dsemRTMB( sem=sem,
+#'             tsdata = tsdata,
+#'             estimate_delta0 = TRUE,
+#'             control = dsem_control(quiet=TRUE) )
+#'
 #' @export
 dsemRTMB <-
 function( sem,

diff --git a/man/dsem.Rd b/man/dsem.Rd
diff --git a/man/dsemRTMB.Rd b/man/dsemRTMB.Rd
diff --git a/tests/testthat/test-platform.R b/tests/testthat/test-platform.R
@@ -49,6 +49,13 @@ test_that("dsem example is working ", {
   # Check objective function
   expect_equal( as.numeric(fit$opt$obj), 587.4755, tolerance=1e-2 )
 
+  #
+  fitRTMB = dsemRTMB( sem=sem,
+              tsdata=Z,
+              control = dsem_control(getsd=FALSE) )
+  # Check objective function
+  expect_equal( as.numeric(fit$opt$obj), as.numeric(fitRTMB$opt$obj), tolerance=1e-2 )
+
   # Convert and plot using phylopath
   as_fitted_DAG(fit)