aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2014-07-21 08:20:44 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2014-07-21 09:20:58 +0200
commita1567638a3ba9f4d62fa199525097a94ddfd7912 (patch)
treea223ff07ae669b61c44b97dcbfd5a48f5243ac24 /R
parent8def5006fc81c032c3fc99751e062cdb32a81cc1 (diff)
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
Diffstat (limited to 'R')
-rw-r--r--R/mkinfit.R27
-rw-r--r--R/transform_odeparms.R22
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"]

Contact - Imprint