From e1434571bc986922a7cd8b292526ad1f570bb10f Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Fri, 10 Jan 2025 07:53:51 -0800 Subject: [PATCH] updates for merge --- .github/workflows/test-coverage.yaml | 4 ++-- DESCRIPTION | 6 ++--- NAMESPACE | 12 +++++++--- NEWS.md | 4 +++- R/compute_nll.R | 16 +++++--------- R/dsem.R | 9 ++++++-- R/dsemRTMB.R | 33 +++++++++++++++++++++++++++- man/dsem.Rd | 6 ++++- man/dsemRTMB.Rd | 33 +++++++++++++++++++++++++++- tests/testthat/test-platform.R | 7 ++++++ 10 files changed, 105 insertions(+), 25 deletions(-) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 35a8cf7..01ee697 100644 --- 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 index df6ebfb..e7fd4c4 100644 --- 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 index 6e56623..1eea5f2 100644 --- 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") \ No newline at end of file diff --git a/NEWS.md b/NEWS.md index d5a4724..77a6820 100644 --- 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 index bce50f9..ddd83e9 100644 --- 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 index 7394bd3..5780b6d 100644 --- 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 index a5a7d41..baf9797 100644 --- 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 index 74a990c..58ae123 100644 --- a/man/dsem.Rd +++ b/man/dsem.Rd @@ -34,7 +34,11 @@ the stationary distribution} 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}} \item{control}{Output from \code{\link{dsem_control}}, used to define user settings, and see documentation for that function for details.} diff --git a/man/dsemRTMB.Rd b/man/dsemRTMB.Rd index 929b193..909f084 100644 --- a/man/dsemRTMB.Rd +++ b/man/dsemRTMB.Rd @@ -55,5 +55,36 @@ An object (list) of class `dsem`, fitted using RTMB Fits a dynamic structural equation model } \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. +} +\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) ) + } diff --git a/tests/testthat/test-platform.R b/tests/testthat/test-platform.R index 5460773..2c9e6ac 100644 --- 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)