From a1567638a3ba9f4d62fa199525097a94ddfd7912 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 21 Jul 2014 08:20:44 +0200
Subject: 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
---
 R/mkinfit.R            | 27 ++++++++++++++++++++++++---
 R/transform_odeparms.R | 22 +++++++++++++++-------
 2 files changed, 39 insertions(+), 10 deletions(-)

(limited to 'R')

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"]
-- 
cgit v1.2.1