Bugfix, model shorthand, state.ini[[1]] from observed data

- The bug occurred when using transform_rates=FALSE for FOMC, DFOP or HS
- Make it possible to use mkinfit("SFO", ...)
- Take initial mean value at time zero for the variable with the highest
  value in the observed data
- Update of vignette/FOCUS_L
- Improve the Makefile to build single vignettes

2 files changed, 39 insertions, 10 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index a8fbfc78..39d084cb 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -23,7 +23,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value"))
 

 mkinfit <- function(mkinmod, observed,

   parms.ini = "auto",

-  state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), 

+  state.ini = "auto", 

   fixed_parms = NULL,

   fixed_initials = names(mkinmod$diffs)[-1],

   solution_type = "auto",

@@ -41,6 +41,21 @@ mkinfit <- function(mkinmod, observed,
   trace_parms = FALSE,

   ...)

 {

+  # Check mkinmod and generate a model for the variable whithe the highest value

+  # if a suitable string is given

+  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB") 

+  presumed_parent_name = observed[which.max(observed$value), "name"]

+  if (class(mkinmod) != "mkinmod") {

+    if (mkinmod[[1]] %in% parent_models_available) {

+      speclist <- list(list(type = mkinmod, sink = TRUE))

+      names(speclist) <- presumed_parent_name

+      mkinmod <- mkinmod(speclist = speclist)

+    } else {

+      stop("Argument mkinmod must be of class mkinmod or a string containing one of\n  ",

+           paste(parent_models_available, collapse = ", "))

+    } 

+  }

+

   # Check optimisation method and set maximum number of iterations if specified

   method.modFit = match.arg(method.modFit)

   if (maxit.modFit != "auto") {

@@ -55,7 +70,7 @@ mkinfit <- function(mkinmod, observed,
   mod_vars <- names(mkinmod$diffs)

 

   # Get the names of observed variables

-  obs_vars = names(mkinmod$spec)

+  obs_vars <- names(mkinmod$spec)

 

   # Subset observed data with names of observed data in the model

   observed <- subset(observed, name %in% obs_vars)

@@ -137,6 +152,12 @@ mkinfit <- function(mkinmod, observed,
     }

   }

 

+  # Set default for state.ini if appropriate

+  if (state.ini[1] == "auto") {

+    state.ini = c(mean(subset(observed, time == 0 & name == presumed_parent_name)$value), 

+                  rep(0, length(mkinmod$diffs) - 1))

+  }

+

   # Name the inital state variable values if they are not named yet

   if(is.null(names(state.ini))) names(state.ini) <- mod_vars

 

@@ -279,7 +300,7 @@ mkinfit <- function(mkinmod, observed,
   if (!transform_rates) {

     index_k <- grep("^k_", names(lower))

     lower[index_k] <- 0

-    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower))

+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))

     lower[other_rate_parms] <- 0

   }

 

diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 912a5c0a..f518ae32 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -69,7 +69,11 @@ transform_odeparms <- function(parms, mkinmod,
   # and HS parameter tb if transformation of rates is requested

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

     if (!is.na(parms[pname])) {

-      transparms[paste0("log_", pname)] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

+      if (transform_rates) {

+        transparms[paste0("log_", pname)] <- log(parms[pname])

+      } else {

+        transparms[pname] <- parms[pname]

+      }

     } 

   }

   if (!is.na(parms["g"])) {

@@ -130,12 +134,16 @@ backtransform_odeparms <- function(transparms, mkinmod,
 

   # Transform parameters also for FOMC, DFOP and HS models

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

-    pname_trans = paste0("log_", pname)

-    if (!is.na(transparms[pname_trans])) {

-      parms[pname] <- ifelse(transform_rates, 

-                             exp(transparms[pname_trans]), 

-                             transparms[pname])

-    } 

+    if (transform_rates) {

+      pname_trans = paste0("log_", pname)

+      if (!is.na(transparms[pname_trans])) {

+        parms[pname] <- exp(transparms[pname_trans])

+      }

+    } else {

+      if (!is.na(transparms[pname])) {

+        parms[pname] <- transparms[pname]

+      }

+    }

   }

   if (!is.na(transparms["g_ilr"])) {

     g_ilr <- transparms["g_ilr"]