Dev (#62)

* add varying paths

* Update dsem_control.Rd

* fix option for TMBad

* Using TMB@dev for now

* simplify CPP

* fix last commit

* updates

* add logistic-regression demo

* Update test_deterministic_with_measurements.R

* add control$stabilize_Q

* Update NEWS.md

* augment estimate_mu

* Update DESCRIPTION

* Update model-description.Rmd

* update vignettes

* Update model-description.Rmd

* Update model-description.Rmd

* add data

* add vignettes

* Update statistical_interactions.Rmd

* Create utils.hpp

* fix defaults for estimate_mu = NULL...

... i.e., estimating mu for all variables with >= 1 observation

* update to use `gmrf_project` ...

... by default

* renaming options for `gmrf_parameterization`

* try to make CI happy

* try again to make CI happy

* Update dynamic_factor_analysis.Rmd

* still trying to make CI happy

* Create dsemRTMB.Rd

* another attempt to make CI happy

* updates

* Update dynamic_factor_analysis.Rmd

* loosen test tol to 0.001

* updates to vignettes

* add informative error ...

... when using gmrf_parameterization = "full" and a zero-exogenous variance variable

* fix bug in previous error-check

* add integrated tests for varying-slopes

* make check_win_devel happy

* Create explore_varying_slopes.R

* Update model-description.Rmd

* add rk4sys for vignettes

* add GLM example

* improved `family` interface

* fix docs for family-list

* add branch to CI

* fix vignettes and bugs

* make CI happy

* adding `gaussian_fixed_sd`

* Update test-distributions.R

* make CI happy

* add lognormal family

* fix typo

* make CI happy

* Update beverton-holt.R

* Create explore_GARCH.R

* add error check

* Update dsem.R

* add informative error for projecting moderators

* Update dsem.R

* add GARCH vignette

* add informative error ...

... for when `sem` includes variables not in `tsdata` (i.e., due to typos)

* adding tweeide and make_msv

* change version number for dev

* fix bug from in merge

* Update model-description.Rmd

* Update _pkgdown.yml

* Update _pkgdown.yml

* Update make_dsem_ram.R

* make R_CMD_check happy

* update docs

* fix typo in nonlinear vignette ...

... h/t Gabriel Leiva-Gomez

* Create github_notes.R

* Update loo_residuals.Rd

* update Version to prep for CRAN

* Update NEWS.md

* add checkmate assertions

---------

Co-authored-by: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com>
Co-authored-by: James Thorson <James.Thorson@noaa.gov>

Author: James-Thorson

.github/workflows/R-CMD-check.yaml

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, dev]
+    branches: [main, dev, dev_add-family]
   pull_request:
-    branches: [main, dev]
+    branches: [main, dev, dev_add-family]
 
 name: R-CMD-check
 

.github/workflows/test-coverage.yaml

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, dev]
+    branches: [main, dev, dev_add-family]
   pull_request:
-    branches: [main, dev]
+    branches: [main, dev, dev_add-family]
 
 name: test-coverage
 

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: dsem
 Type: Package
 Title: Dynamic Structural Equation Models
-Version: 2.0.1
-Date: 2026-05-12
+Version: 3.0.0
+Date: 2026-06-10
 Authors@R: 
     c(person(given = "James",
              family = "Thorson",
@@ -34,7 +34,8 @@ Imports:
     ggplot2,
     grid,
     methods,
-    stats
+    stats,
+    checkmate
 Depends:
     R (>= 4.0.0),
 Suggests:

NAMESPACE

@@ -15,17 +15,22 @@ export(classify_variables)
 export(convert_equations)
 export(dsem)
 export(dsem_control)
+export(fixed)
+export(gaussian_fixed_sd)
 export(list_parameters)
+export(lognormal)
 export(loo_residuals)
 export(make_dfa)
 export(make_dsem_ram)
 export(make_matrices)
+export(make_msv)
 export(parse_path)
 export(partition_variance)
 export(rk4sys)
 export(stepwise_selection)
 export(test_dsep)
 export(total_effect)
+export(tweedie)
 importFrom(Matrix,Cholesky)
 importFrom(Matrix,Diagonal)
 importFrom(Matrix,bandSparse)
@@ -53,6 +58,10 @@ importFrom(TMB,config)
 importFrom(TMB,dynlib)
 importFrom(TMB,sdreport)
 importFrom(TMB,summary.sdreport)
+importFrom(checkmate,assertClass)
+importFrom(checkmate,assertList)
+importFrom(checkmate,assertLogical)
+importFrom(checkmate,assertNames)
 importFrom(ggplot2,aes)
 importFrom(ggraph,create_layout)
 importFrom(ggraph,geom_edge_arc)
@@ -71,13 +80,17 @@ importFrom(methods,is)
 importFrom(stats,"tsp<-")
 importFrom(stats,.preformat.ts)
 importFrom(stats,AIC)
+importFrom(stats,Gamma)
+importFrom(stats,binomial)
+importFrom(stats,gaussian)
 importFrom(stats,logLik)
 importFrom(stats,na.omit)
 importFrom(stats,nlminb)
 importFrom(stats,optimHess)
 importFrom(stats,pchisq)
 importFrom(stats,plogis)
 importFrom(stats,pnorm)
+importFrom(stats,poisson)
 importFrom(stats,rbinom)
 importFrom(stats,rgamma)
 importFrom(stats,rnorm)

NEWS.md

@@ -1,3 +1,14 @@
+# dsem 3.0.0
+
+* Change to requiring `family` options that agree with `glm`
+* Adding `family = gaussian_fixed_sd()` option to pre-specify a known standard deviation
+* Adding `family = lognormal()` option (e.g., for linearized Beverton-Holt stock recruit analysis)
+* Add `family = tweedie()` option (e.g., for paleo Lake cores)
+* Adding option for moderator variables affecting exogenous variance in Gamma (not just paths in Rho)
+* adding `make_msv` for helping multivariate stochastic volatility models
+* Adding MGARCH vignette
+* Adding `checkmate` assertions for user inputs
+
 # dsem 2.0.1
 
 * Fix an "out-of-bounds read" error identified by SAN

R/dsem.R

@@ -7,14 +7,20 @@
 #'        \code{\link[dsem]{make_dsem_ram}} for more description
 #' @param tsdata time-series data, as outputted using \code{\link[stats]{ts}}, with \code{NA}
 #'        for missing values.
-#' @param family Character-vector listing the distribution used for each column of \code{tsdata}, where
-#'        each element must be \code{fixed} (for no measurement error), 
-#'        \code{normal} for normal measurement error using an identity link,
-#'        \code{gamma} for a gamma measurement error using a fixed CV and log-link, 
-#'        \code{bernoulli} for a Bernoulli measurement error using a logit-link, or
-#'        \code{poisson} for a Poisson measurement error using a log-link.
-#'        \code{family="fixed"} is default behavior and assumes that a given variable is measured exactly.
-#'        Other options correspond to different specifications of measurement error.
+#' @param family A named list of families, each returning a class \code{family},
+#'        including [fixed()], [gaussian()], [binomial()],  [Gamma()], [poisson()],
+#'        [lognormal()], [tweedie()], or [gaussian_fixed_sd()] with names that match levels of
+#'        \code{colnames(tsdata)} to allow different
+#'        families by variable.  Family [fixed()] specifies that states
+#'        are known (i.e., measurements for that variable have no error).
+#'        Other families allow users to supply a link function including `identity`,
+#'        `log`, `logit`, or `cloglog`.  For example
+#'        \code{family = list(y = binomial("logit"), x = fixed())}
+#'        would specify logit-linked Bernoulli distribution for variable `tsdata$y`
+#'        and a fixed (no measurement error) distribution for `tsdata$x`.
+#'        For many variables, it is convenient to do e.g.,
+#'        \code{family = Map(function(.) gaussian(), colnames(tsdata))} rather than writing
+#'        them all manually.
 #' @param estimate_delta0 Boolean indicating whether to estimate deviations from equilibrium in initial year
 #'        as fixed effects, or alternatively to assume that dynamics start at some stochastic draw away from
 #'        the stationary distribution
@@ -42,14 +48,17 @@
 #'        settings, and see documentation for that function for details.
 #'
 #' @importFrom TMB compile dynlib MakeADFun sdreport summary.sdreport config
-#' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rbinom rgamma rpois rnorm simulate time tsp<- plogis pchisq
+#' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rbinom
+#' @importFrom stats rgamma rpois rnorm simulate time tsp<- plogis pchisq
+#' @importFrom stats gaussian poisson Gamma binomial
 #' @importFrom Matrix solve Cholesky sparseMatrix mat2triplet drop0 t
 #' @importFrom igraph plot.igraph graph_from_data_frame with_sugiyama layout_ graph.adjacency clusters
 #' @importFrom ggraph ggraph geom_edge_arc create_layout rectangle geom_node_text theme_graph
 #' @importFrom ggplot2 aes
 #' @importFrom grid arrow
 #' @importFrom methods is
 #' @importFrom utils combn
+#' @importFrom checkmate assertList assertClass assertNames assertLogical
 #'
 #' @details
 #' A DSEM involves (at a minimum):
@@ -186,10 +195,12 @@
 #'                "gwage","invest","capital")]
 #'
 #' # Fit model
-#' fit = dsem( sem=sem,
-#'             tsdata = tsdata,
-#'             estimate_delta0 = TRUE,
-#'             control = dsem_control(quiet=TRUE) )
+#' fit = dsem(
+#'   sem = sem,
+#'   tsdata = tsdata,
+#'   estimate_delta0 = TRUE,
+#'   control = dsem_control(quiet=TRUE)
+#' )
 #' summary( fit )
 #' plot( fit )
 #' plot( fit, edge_label="value" )
@@ -199,22 +210,32 @@
 dsem <-
 function( sem,
           tsdata,
-          family = rep("fixed",ncol(tsdata)),
+          family = Map(function(.) fixed(), colnames(tsdata)),
           estimate_delta0 = FALSE,
           estimate_mu = NULL,
           prior_negloglike = NULL,
           control = dsem_control(),
           covs = colnames(tsdata) ){
   start_time = Sys.time()
 
+  # Temporary warning
+  if( isTRUE(is(family, "character")) ){
+    stop("starting with release 3.0.0, `family` must be a named list of ")
+  }
+
   # General error checks
-  if( isFALSE(is(control, "dsem_control")) ) stop("`control` must be made by `dsem_control()`")
+  #if( isFALSE(is(control, "dsem_control")) ) stop("`control` must be made by `dsem_control()`")
   if( isTRUE(control$gmrf_parameterization=="project") ){
     if( isTRUE(any(family=="fixed" & colSums(!is.na(tsdata))>0)) ){
       stop("`family` cannot be `fixed` using `gmrf_parameterization=projection` for any variable with data")
     }
   }
-  if( isFALSE(is(tsdata,"ts")) ) stop("`tsdata` must be a `ts` object")
+  #if( isFALSE(is(tsdata,"ts")) ) stop("`tsdata` must be a `ts` object")
+  assertClass(control, classes = "dsem_control")
+  assertClass(tsdata, classes = "ts")
+  assertNames( names(family), must.include = colnames(tsdata))
+  assertList(family, types = "family")
+  assertLogical(estimate_delta0, len = 1)
 
   # General warnings
   if( isFALSE(control$quiet) ){
@@ -249,7 +270,7 @@ function( sem,
       stop("If specifying `lower` or `upper`, please set `dsem_control('newton_loops'=0)`")
     }
   }
-  if( any(ram$heads==0) ){
+  if( any(ram$heads %in% c(3,4)) ){
     # Using moderating variables, their raw values are used to construct Rho_kk and raw values follow standard-normal distribution
     # so Gamma_kk and Rho_kk are not properly applied to moderating variables
     if( control$gmrf_parameterization %in% c("project") ){
@@ -261,7 +282,8 @@ function( sem,
   options = c(
     switch(control$gmrf_parameterization, "full" = 0, "project" = 1, "mvn_project" = 2, "gmrf_project" = 3, NA),
     switch(control$constant_variance, "conditional"=0, "marginal"=1, "diagonal"=2),
-    ifelse( isTRUE(control$stabilize_Q), 1, 0 )
+    ifelse( isTRUE(control$stabilize_Q), 1, 0 ),
+    ifelse( isTRUE(control$logscale_moderating_variance), 1, 0)
   )
 
   # define variables to project to `project_k` unless provided by user
@@ -304,14 +326,89 @@ function( sem,
     }
   }
 
+  # Check for errors
+  if( !all(colnames(tsdata) %in% names(family)) ){
+    stop("some variable in `colnames(tsdata)` not found in `names(family)`")
+  }
+
+  # distribution/link
+  build_distributions <-
+  function( variables ){
+
+    # Construct log_sigma based on family
+    remove_last = function(x) x[-length(x)]
+
+    # Fixed values
+    sigma_j = lapply( family, FUN=function(x){
+                       switch( x$family[length(x$family)],
+                         "fixed" = c(),
+                         "gaussian" = NA,
+                         "poisson" = c(),
+                         "binomial" = c(),
+                         "bernoulli" = c(),
+                         "Gamma" = NA,
+                         "gaussian_fixed_sd" = c(),
+                         "lognormal" = NA,
+                         "tweedie" = c(NA, NA)
+                       )} )
+    Nsigma_j = sapply(sigma_j, length)
+    sigmastart_j = remove_last(cumsum(c(0,Nsigma_j)))
+    names(sigmastart_j) = variables
+
+    # Extract known SDs
+    sd_tj = sapply( family,
+      FUN = function(x){
+        if(length(x$fixed_sd)==0){
+          return(rep(0, nrow(tsdata)))
+        }else{
+          sd = rep(0, nrow(tsdata))
+          sd[] = x$fixed_sd
+          return(sd)
+        }
+      }
+    )
+
+    #
+    family_code = sapply( family, FUN=function(x){
+                       c("fixed" = 0,
+                         "gaussian" = 1,
+                         "binomial" = 2,
+                         "bernoulli" = 2,
+                         "poisson" = 3,
+                         "Gamma" = 4,
+                         "gaussian_fixed_sd" = 5,
+                         "lognormal" = 6,
+                         "tweedie" = 7
+                       )[x$family]} )
+    link_code = sapply( family, FUN=function(x){
+                       c("identity" = 0,
+                         "log" = 1,
+                         "logit" = 2,
+                         "cloglog" = 3
+                       )[x$link]} )
+    out = list(
+      family_code = family_code,
+      link_code = link_code,
+      Nsigma_j = Nsigma_j,
+      sigmastart_j = sigmastart_j,
+      sigma_j = sigma_j,
+      sd_tj = sd_tj
+    )
+    return(out)
+  }
+  family = family[match(colnames(tsdata), names(family))]
+  distributions = build_distributions( colnames(tsdata) )
+
   #
   Data = list( 
-    "options" = options,
-    #"RAM" = as.matrix(na.omit(ram[,1:4])),
-    "RAM" = as.matrix(ram[,-5]),
-    "RAMstart" = as.numeric(ram[,5]),
-    "familycode_j" = sapply(family, FUN=switch, "fixed"=0, "normal"=1, "bernoulli"=2, "poisson"=3, "gamma"=4 ),
-    "y_tj" = tsdata 
+    options = options,
+    RAM = as.matrix(ram[,-5]),
+    RAMstart = as.numeric(ram[,5]),
+    familycode_j = distributions$family_code,
+    linkcode_j = distributions$link_code,
+    sigmastart_j = distributions$sigmastart_j,
+    eps_tj = distributions$sd_tj,
+    y_tj = tsdata
   )
 
   # Default ... unobs_idx for fixed variables with zero variance
@@ -323,10 +420,10 @@ function( sem,
   # Construct parameters
   if( is.null(control$parameters) ){
     Params = list(
-      beta_z = rep(0,max(ram[,4],na.rm=TRUE)),  # NA for spatially-varying paths in ram
-      lnsigma_j = rep(0,ncol(tsdata)),
-      mu_j = rep(0,ncol(tsdata)),
-      delta0_j = rep(0,ncol(tsdata)),
+      beta_z = rep(0, max(ram[,4],na.rm=TRUE)),  # NA for spatially-varying paths in ram
+      lnsigma_z = rep(0, sum(distributions$Nsigma_j)),
+      mu_j = rep(0, ncol(tsdata)),
+      delta0_j = rep(0, ncol(tsdata)),
       x_tj = ifelse( is.na(tsdata), 0, tsdata )
     )
     #if( control$gmrf_parameterization=="full" ){
@@ -368,9 +465,9 @@ function( sem,
   if( is.null(control$map) ){
     Map = list()
     # Map off x_tj for fixed when data is available
-    Map$x_tj = ifelse( is.na(as.vector(tsdata)) | (Data$familycode_j[col(tsdata)] %in% c(1,2,3,4)), seq_len(prod(dim(tsdata))), NA )
+    Map$x_tj = ifelse( is.na(as.vector(tsdata)) | (Data$familycode_j[col(tsdata)] %in% c(1,2,3,4,5,6,7)), seq_len(prod(dim(tsdata))), NA )
     # Map off sigma_j for fixed / bernoulli / Poisson
-    Map$lnsigma_j = factor( ifelse(Data$familycode_j %in% c(0,2,3), NA, seq_along(Params$lnsigma_j)) )
+    # Map$lnsigma_j = factor( ifelse(Data$familycode_j %in% c(0,2,3), NA, seq_along(Params$lnsigma_j)) )
 
     # Map off mean for latent variables
     Map$mu_j = factor( ifelse(colnames(tsdata) %in% estimate_mu, seq_len(ncol(tsdata)), NA) )
@@ -645,6 +742,9 @@ function( sem,
 #'        are then replaced with Inf and avoided during estimation
 #' @param stabilize_Q add \code{stability_eps = 1e-10} to stabilize precision
 #'        (experimental)
+#' @param logscale_moderating_variance When users supply a moderating variable
+#'        for an exogenous variance (i.e., double-headed arrow), whether to
+#'        exponentiate the variable prior to
 #'
 #' @return
 #' An S3 object of class "dsem_control" that specifies detailed model settings,
@@ -673,7 +773,8 @@ function( nlminb_loops = 1,
           upper = Inf,
           project_k = NULL,
           suppress_nlminb_warnings = TRUE,
-          stabilize_Q = FALSE ){
+          stabilize_Q = FALSE,
+          logscale_moderating_variance = FALSE ){
 
   gmrf_parameterization = match.arg(gmrf_parameterization)
   constant_variance = match.arg(constant_variance)
@@ -701,7 +802,8 @@ function( nlminb_loops = 1,
     upper = upper,
     project_k = project_k,
     suppress_nlminb_warnings = suppress_nlminb_warnings,
-    stabilize_Q = stabilize_Q
+    stabilize_Q = stabilize_Q,
+    logscale_moderating_variance = logscale_moderating_variance
   ), class = "dsem_control" )
 }