From a2e4fc306dd4e477ffdfad5cb96d09bf0b5619fc Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Fri, 10 Jan 2025 17:26:06 -0800 Subject: [PATCH] add CI for fixed and mapped parameters and fix issue in dsemRTMB with these --- R/dsem.R | 4 ++-- R/dsemRTMB.R | 13 +++++++------ R/make_matrices.R | 22 +++++++++++++++------- R/read_model.R | 6 ++++-- tests/testthat/test-platform.R | 28 ++++++++++++++++++++++++++++ 5 files changed, 56 insertions(+), 17 deletions(-) diff --git a/R/dsem.R b/R/dsem.R index dabbe09..929e419 100644 --- 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 index 2c033dc..198030a 100644 --- 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 index 121889e..d8f87f1 100644 --- 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 index 4c2c084..cb851af 100644 --- 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 index f8c79cc..923193e 100644 --- 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 ) + +}) \ No newline at end of file