add CI for fixed and mapped parameters

and fix issue in dsemRTMB with these
James-Thorson-NOAA · Jan 11, 2025 · a2e4fc3 · a2e4fc3
1 parent bf47302
commit a2e4fc3
Show file tree

Hide file tree

Showing 5 changed files with 56 additions and 17 deletions.
diff --git a/R/dsem.R b/R/dsem.R
@@ -152,8 +152,8 @@ function( sem,
   # General warnings
   if( isFALSE(control$quiet) ){
     tsdata_SD = apply( tsdata, MARGIN=2, FUN=sd, na.rm=TRUE )
-    if( any((max(tsdata_SD)/min(tsdata_SD)) > 100, rm.na=TRUE) ){
-      warning("Some variables in `tsdata` have much higher variance than others. Please consider rescaling variables to prevent issues with numerical convergence.")
+    if( any((max(tsdata_SD,rm.na=TRUE)/min(tsdata_SD,rm.na=TRUE)) > 100) ){
+       warning("Some variables in `tsdata` have much higher variance than others. Please consider rescaling variables to prevent issues with numerical convergence.")
     }
   }
 

diff --git a/R/dsemRTMB.R b/R/dsemRTMB.R
@@ -74,7 +74,7 @@ function( sem,
   # General warnings
   if( isFALSE(control$quiet) ){
     tsdata_SD = apply( tsdata, MARGIN=2, FUN=sd, na.rm=TRUE )
-    if( any((max(tsdata_SD)/min(tsdata_SD)) > 100, rm.na=TRUE) ){
+    if( any((max(tsdata_SD,rm.na=TRUE)/min(tsdata_SD,rm.na=TRUE)) > 100) ){
       warning("Some variables in `tsdata` have much higher variance than others. Please consider rescaling variables to prevent issues with numerical convergence.")
     }
   }
@@ -85,10 +85,11 @@ function( sem,
             variables = colnames(tsdata),
             covs = covs,
             quiet = FALSE )
-  ram = make_matrices(
-            model = model,
-            times = as.numeric(time(tsdata)),
-            variables = colnames(tsdata) )
+  #ram = make_matrices(
+  #          beta_p = rnorm(nrow(model)),
+  #          model = model,
+  #          times = as.numeric(time(tsdata)),
+  #          variables = colnames(tsdata) )
 
   #
   options = c(
@@ -189,7 +190,7 @@ function( sem,
 
   # Further bundle
   out = list( "obj" = obj,
-              "ram" = ram,
+              #"ram" = ram, # Not useful using RTMB structure
               "sem_full" = model,
               "tmb_inputs"=list("parameters"=Params, "random"=Random, "map"=Map),
               #"call" = match.call(),

diff --git a/R/make_matrices.R b/R/make_matrices.R
@@ -4,10 +4,19 @@ function( beta_p,
           times,
           variables ){
 
+  "c" <- ADoverload("c")
+  "[<-" <- ADoverload("[<-")
   model <- as.data.frame(model)
-  if(missing(beta_p)){
-    model_unique = model[match(unique(model$parameter),model$parameter),]
-    beta_p =  as.numeric(model_unique$start)
+
+  # Combine fixed, estimated, and mapped parameters into vector
+  beta_i = rep(0, nrow(model))
+  off = which(model[,'parameter'] == 0)
+  if( length(off) > 0 ){
+    beta_i[off] = as.numeric(model[off,'start'])
+  }
+  not_off = which(model[,'parameter'] > 0)
+  if( length(not_off) > 0 ){
+    beta_i[not_off] = beta_p[model[not_off,'parameter']]
   }
 
   # Loop through paths
@@ -27,12 +36,11 @@ function( beta_p,
                          dims = rep(length(variables),2) )
 
     # Assemble
+    tmp_kk = kronecker(P_jj, L_tt)
     if(abs(as.numeric(model[i,'direction']))==1){
-      tmp_kk = (kronecker(P_jj, L_tt))
-      P_kk = P_kk + beta_p[model$parameter[i]] * tmp_kk # AD(tmp_kk)
+      P_kk = P_kk + beta_i[i] * tmp_kk # AD(tmp_kk)
     }else{
-      tmp_kk = (kronecker(P_jj, L_tt))
-      G_kk = G_kk + beta_p[model$parameter[i]] * tmp_kk # AD(tmp_kk)
+      G_kk = G_kk + beta_i[i] * tmp_kk # AD(tmp_kk)
     }
   }
 

diff --git a/R/read_model.R b/R/read_model.R
@@ -85,8 +85,10 @@ function( sem,
                 quiet = quiet)
   model$path <- gsub("\\t", " ", model$path)
   model$par[model$par == ""] <- NA
-  model <- data.frame( "path"=model$path, "lag"=model$lag,
-                       "name"=model$par, "start"=model$start)
+  model <- data.frame( "path" = model$path,
+                       "lag" = model$lag,
+                       "name" = model$par,
+                       "start" = model$start )
 
   # Adding a SD automatically
   if( !is.null(covs) ){

diff --git a/tests/testthat/test-platform.R b/tests/testthat/test-platform.R
@@ -137,3 +137,31 @@ test_that("bering sea example is stable ", {
   expect_equal( as.numeric(fit$opt$obj), 189.3005, tolerance=1e-2 )
 })
 
+test_that("Fixing parameters works ", {
+  sem = "
+    A -> B, 0, NA, 1
+    B -> C, 0, NA, 0.5
+    # Extra links
+    A -> D, 1, beta
+    B -> D, 1, beta
+  "
+  set.seed(101)
+  nobs = 40
+  A = rnorm(nobs)
+  B = A + rnorm(nobs)
+  C = 0.5 * B + rnorm(nobs)
+  D = rnorm(nobs)
+  tsdata = ts(cbind(A=A, B=B, C=C, D=D))
+
+  # Run models
+  fit = dsem( sem=sem,
+               tsdata=tsdata,
+               control = dsem_control(getsd=FALSE) )
+  fitRTMB = dsemRTMB( sem=sem,
+               tsdata=tsdata,
+               control = dsem_control(getsd=FALSE) )
+  # Check objective function
+  expect_equal( as.numeric(fit$opt$obj), 224.2993, tolerance=1e-2 )
+  expect_equal( as.numeric(fit$opt$obj), as.numeric(fitRTMB$opt$obj), tolerance=1e-2 )
+
+})