mhunter1 / dynr

Dynamic Modeling in R
5 stars 6 forks source link

NaNs produced during model fitting in dynr leading to errors in dynr.taste2. #279

Open chriscmr opened 2 months ago

chriscmr commented 2 months ago

Good day. I am working on a model using the dynr package and encountering several issues related to NaNs being produced during the fitting process. Below is the code I am using, along with the warnings and errors I keep encountering.

// Load necessary libraries
library(dynr)
dat <- temp

mean_x <- mean(dat)
sd_x <- sd(dat)
var_x <- var(dat)
cov_x <- cov(dat,dat)

// Prepare the data for dynr
data_shk <- dynr.data(dat)

// Measurement preparation
meas_shk <- prep.measurement(
  values.load = matrix(c(1), ncol = 1, byrow = TRUE),
  params.load = matrix(c('l_11'), ncol = 1, byrow = TRUE),
  state.names = c('eta_1'),
  obs.names = c('x')
)

// Noise preparation using the variance of the observed variable
nois_shk <- prep.noise(
  values.latent = matrix(c(var_x), ncol = 1, byrow = TRUE), # Variance of latent variable
  params.latent = matrix(c('psi_11'), ncol = 1, byrow = TRUE),
  values.observed = diag(var_x*0.5, ncol = 1, nrow = 1),        # Variance of observed variable
  params.observed = diag(paste0('e_', 1:1), 1)
)

// Initial conditions for state and latent using the mean and variance of 'x'
init_shk <- prep.initial(
  values.inistate = c(mean_x),  # Mean as starting value
  params.inistate = c('mu_1'),
  values.inicov = matrix(c(var_x*), ncol = 1, byrow = TRUE),  # Variance as initial covariance
  params.inicov = matrix(c('c_11'), ncol = 1, byrow = TRUE)
)

// Dynamics preparation 
dynm_shk <- prep.matrixDynamics(
  values.dyn = matrix(c(1), ncol = 1, byrow = TRUE),  # Transition matrix
  params.dyn = matrix(c('b_11'), ncol = 1, byrow = TRUE),
  isContinuousTime = FALSE
)

// Prepare the model with the updated initial conditions and noise parameters
model_shk <- dynr.model(
  dynamics = dynm_shk,
  measurement = meas_shk,
  noise = nois_shk,
  initial = init_shk,
  data = data_shk,
  outfile = paste0("tempk.c")
)

// Update options to increase the number of evaluations
opt <- model_shk$options
opt$maxeval <- 10000
model_shk@options <- opt

// Fit the model (run KL and FIS)
cook_shk <- dynr.cook(model_shk, debug_flag = TRUE)

//Print the model results or debug information
print(cook_shk)

When running the above code, I consistently receive the following warning messages:

Warning message in sqrt(diag(iHess)):
"NaNs produced"
Warning message in sqrt(diag(x$inv.hessian)):
"NaNs produced"
Warning message:
"These parameters may have untrustworthy standard errors: e_1, c_11."

I proceed with the second part of the analysis, which involves computing shocks and chi-square using dynr.taste, followed by re-fitting the state-space model using dynr.taste2.

taste_shk <- dynr.taste(model_shk, cook_shk, conf.level = .99) #compute shocks and chi-square

taste_plot <- autoplot(taste_shk)

taste2_shk <- dynr.taste2(#Re-fit state-space model using the estimated outliers.
  model_shk, cook_shk, taste_shk,
  newOutfile = "taste2_shk.c"
)

When I run this second part, I get the following warning and error messages:

Warning message in qt(conf_level, nrow(t_df_inn_plot_i) - obs_n):
"NaNs produced"
Warning message in qt(conf_level, t_add_df_i - obs_n):
"NaNs produced"
Error in deltaLat[, begT:endT] <- cbind(rep(0, lat_n), deltaLat[, begT:(endT - : number of items to replace is not a multiple of replacement length
Traceback:
1. dynr.taste2(model_shk, cook_shk, taste_shk, newOutfile = "taste2_shk.c")

What I Have Tried:

  1. Adjusting starting values: I've changed the starting values several times, especially for e_1 and c_11, but the warnings and errors persist.
  2. Testing with simple data: I even tried using a simple dataset generated from a normal distribution to check if it was an issue with the data, but I still encountered the same warnings and errors.
  3. Modifying noise and initial condition values: I've adjusted the variance scaling, but no success so far.
mhunter1 commented 1 week ago

The following could be issues with your model.

// Initial conditions for state and latent using the mean and variance of 'x'
init_shk <- prep.initial(
  values.inistate = c(mean_x),  # Mean as starting value
  params.inistate = c('mu_1'),
  values.inicov = matrix(c(var_x*), ncol = 1, byrow = TRUE),  # Variance as initial covariance
  params.inicov = matrix(c('c_11'), ncol = 1, byrow = TRUE)
)

The above code freely estimates the initial mean and variance. For single-subject data, these parameters are not identified. You cannot estimate them for single-subject data. For multi-subject data you can estimate them. A typical solution for single-subject data is to have a fixed initial mean of zero and a fixed initial variance of some reasonably large value (i.e., diffuse initial conditions).

// Dynamics preparation 
dynm_shk <- prep.matrixDynamics(
  values.dyn = matrix(c(1), ncol = 1, byrow = TRUE),  # Transition matrix
  params.dyn = matrix(c('b_11'), ncol = 1, byrow = TRUE),
  isContinuousTime = FALSE
)

The above code starts the autoregressive dynamics at the boundary of stability. The AR coefficient is 1.0 (i.e. a random walk). This is usually a bad spot to start your parameter optimization. I suggest an AR starting value in the range of -.3 to +.3.

I'd also recommend setting lower and upper bounds on most of your free parameters.

Finally for dynr.taste, if the initial model does not yield reasonable parameter estimates and standard errors, then I would not proceed with dynr.taste. I'm not sure what exactly is causing your errors, but the warnings are because the standard errors are all messed up.