aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/confint.mkinfit.R3
-rw-r--r--R/logLik.mkinfit.R9
-rw-r--r--R/lrtest.mkinfit.R57
-rw-r--r--R/mkinfit.R52
-rw-r--r--R/residuals.mkinfit.R31
-rw-r--r--R/update.mkinfit.R57
6 files changed, 191 insertions, 18 deletions
diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R
index fadd14ae..5e1703d6 100644
--- a/R/confint.mkinfit.R
+++ b/R/confint.mkinfit.R
@@ -57,7 +57,8 @@
#' # c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 1))
#' # If we exclude parent_0 (the confidence of which is often of minor interest), we get a nice
#' # performance improvement from about 30 seconds to about 12 seconds
-#' # system.time(ci_profile_no_parent_0 <- confint(f_d_1, c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 4))
+#' # system.time(ci_profile_no_parent_0 <- confint(f_d_1,
+#' # c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 4))
#' ci_profile
#' ci_quadratic_transformed <- confint(f_d_1, method = "quadratic")
#' ci_quadratic_transformed
diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
index 4ec3d7d4..cadc0d0a 100644
--- a/R/logLik.mkinfit.R
+++ b/R/logLik.mkinfit.R
@@ -1,9 +1,10 @@
#' Calculated the log-likelihood of a fitted mkinfit object
#'
-#' This function simply calculates the product of the likelihood densities
-#' calculated using \code{\link{dnorm}}, i.e. assuming normal distribution,
-#' with of the mean predicted by the degradation model, and the standard
-#' deviation predicted by the error model.
+#' This function returns the product of the likelihood densities of each
+#' observed value, as calculated as part of the fitting procedure using
+#' \code{\link{dnorm}}, i.e. assuming normal distribution, and with the means
+#' predicted by the degradation model, and the standard deviations predicted by
+#' the error model.
#'
#' The total number of estimated parameters returned with the value of the
#' likelihood is calculated as the sum of fitted degradation model parameters
diff --git a/R/lrtest.mkinfit.R b/R/lrtest.mkinfit.R
new file mode 100644
index 00000000..9c0a9039
--- /dev/null
+++ b/R/lrtest.mkinfit.R
@@ -0,0 +1,57 @@
+#' @importFrom lmtest lrtest
+#' @export
+lmtest::lrtest
+
+#' Likelihood ratio test for mkinfit models
+#'
+#' Compare two mkinfit models based on their likelihood. If two fitted
+#' mkinfit objects are given as arguments, it is checked if they have been
+#' fitted to the same data. It is the responsibility of the user to make sure
+#' that the models are nested, i.e. one of them has less degrees of freedom
+#' and can be expressed by fixing the parameters of the other.
+#'
+#' Alternatively, an argument to mkinfit can be given which is then passed
+#' to \code{\link{update.mkinfit}} to obtain the alternative model.
+#'
+#' The comparison is then made by the \code{\link[lmtest]{lrtest.default}}
+#' method from the lmtest package. The model with the higher number of fitted
+#' parameters (alternative hypothesis) is listed first, then the model with the
+#' lower number of fitted parameters (null hypothesis).
+#'
+#' @importFrom stats logLik update
+#' @param object An \code{\link{mkinfit}} object
+#' @param object_2 Optionally, another mkinfit object fitted to the same data.
+#' @param \dots Argument to \code{\link{mkinfit}}, passed to
+#' \code{\link{update.mkinfit}} for creating the alternative fitted object.
+#' @examples
+#' \dontrun{
+#' test_data <- subset(synthetic_data_for_UBA_2014[[12]]$data, name == "parent")
+#' sfo_fit <- mkinfit("SFO", test_data, quiet = TRUE)
+#' dfop_fit <- mkinfit("DFOP", test_data, quiet = TRUE)
+#' lrtest(dfop_fit, sfo_fit)
+#' lrtest(sfo_fit, dfop_fit)
+#' lrtest(dfop_fit, error_model = "tc")
+#' lrtest(dfop_fit, fixed_parms = c(k2 = 0))
+#' }
+#' @export
+lrtest.mkinfit <- function(object, object_2 = NULL, ...) {
+
+ name_function <- function(x) {
+ paste(x$mkinmod$name, "with error model", x$err_mod)
+ }
+
+ if (is.null(object_2)) {
+ object_2 <- update(object, ...)
+ } else {
+ data_object <- object$data[c("time", "variable", "observed")]
+ data_object_2 <- object_2$data[c("time", "variable", "observed")]
+ if (!identical(data_object, data_object_2)) {
+ stop("It seems that the mkinfit objects have not been fitted to the same data")
+ }
+ }
+ if (attr(logLik(object), "df") > attr(logLik(object_2), "df")) {
+ lmtest::lrtest.default(object, object_2, name = name_function)
+ } else {
+ lmtest::lrtest.default(object_2, object, name = name_function)
+ }
+}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index a3e39053..27c769db 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -52,7 +52,9 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
#' default values. Otherwise, inital values for all error model parameters
#' must be given.
#' @param fixed_parms The names of parameters that should not be optimised but
-#' rather kept at the values specified in \code{parms.ini}.
+#' rather kept at the values specified in \code{parms.ini}. Alternatively,
+#' a named numeric vector of parameters to be fixed, regardless of the values
+#' in parms.ini.
#' @param fixed_initials The names of model variables for which the initial
#' state at time 0 should be excluded from the optimisation. Defaults to all
#' state variables except for the first one.
@@ -253,6 +255,12 @@ mkinfit <- function(mkinmod, observed,
trace_parms = FALSE,
...)
{
+ call <- match.call()
+
+ # Derive the name used for the model
+ if (is.character(mkinmod)) mkinmod_name <- mkinmod
+ else mkinmod_name <- deparse(substitute(mkinmod))
+
# Check mkinmod and generate a model for the variable whith the highest value
# if a suitable string is given
parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
@@ -302,6 +310,14 @@ mkinfit <- function(mkinmod, observed,
# Define starting values for parameters where not specified by the user
if (parms.ini[[1]] == "auto") parms.ini = vector()
+ # Override parms.ini for parameters given as a numeric vector in
+ # fixed_parms
+ if (is.numeric(fixed_parms)) {
+ fixed_parm_names <- names(fixed_parms)
+ parms.ini[fixed_parm_names] <- fixed_parms
+ fixed_parms <- fixed_parm_names
+ }
+
# Warn for inital parameter specifications that are not in the model
wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
if (length(wrongpar.names) > 0) {
@@ -384,15 +400,22 @@ mkinfit <- function(mkinmod, observed,
# Set default for state.ini if appropriate
parent_name = names(mkinmod$spec)[[1]]
+ parent_time_0 = subset(observed, time == 0 & name == parent_name)$value
+ parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)
+ if (is.na(parent_time_0_mean)) {
+ state.ini_auto = c(100, rep(0, length(mkinmod$diffs) - 1))
+ } else {
+ state.ini_auto = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1))
+ }
+ names(state.ini_auto) <- mod_vars
+
if (state.ini[1] == "auto") {
- parent_time_0 = subset(observed, time == 0 & name == parent_name)$value
- parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)
- if (is.na(parent_time_0_mean)) {
- state.ini = c(100, rep(0, length(mkinmod$diffs) - 1))
- } else {
- state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1))
- }
+ state.ini_used <- state.ini_auto
+ } else {
+ state.ini_used <- state.ini_auto
+ state.ini_used[names(state.ini)] <- state.ini
}
+ state.ini <- state.ini_used
# Name the inital state variable values if they are not named yet
if(is.null(names(state.ini))) names(state.ini) <- mod_vars
@@ -799,19 +822,21 @@ mkinfit <- function(mkinmod, observed,
names(c(degparms, errparms)))
# Backtransform parameters
- bparms.optim = backtransform_odeparms(fit$par, mkinmod,
+ bparms.optim = backtransform_odeparms(degparms, mkinmod,
transform_rates = transform_rates,
transform_fractions = transform_fractions)
bparms.fixed = c(state.ini.fixed, parms.fixed)
bparms.all = c(bparms.optim, parms.fixed)
- fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms),
+ fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.optim, errparms),
OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE)
- dimnames(fit$hessian_notrans) <- list(names(c(bparms.all, errparms)),
- names(c(bparms.all, errparms)))
+ dimnames(fit$hessian_notrans) <- list(names(c(bparms.optim, errparms)),
+ names(c(bparms.optim, errparms)))
})
+ fit$call <- call
+
fit$error_model_algorithm <- error_model_algorithm
if (fit$convergence != 0) {
@@ -831,8 +856,9 @@ mkinfit <- function(mkinmod, observed,
fit$calls <- calls
fit$time <- fit_time
- # We also need the model for summary and plotting
+ # We also need the model and a model name for summary and plotting
fit$mkinmod <- mkinmod
+ fit$mkinmod$name <- mkinmod_name
# We need data and predictions for summary and plotting
fit$observed <- observed
diff --git a/R/residuals.mkinfit.R b/R/residuals.mkinfit.R
new file mode 100644
index 00000000..96bcf01c
--- /dev/null
+++ b/R/residuals.mkinfit.R
@@ -0,0 +1,31 @@
+#' Extract residuals from an mkinfit model
+#'
+#' @param object An \code{\link{mkinfit}} object
+#' @param standardized Should the residuals be standardized by dividing by the
+#' standard deviation obtained from the fitted error model?
+#' @param \dots Not used
+#' @export
+#' @examples
+#' f <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
+#' residuals(f)
+#' residuals(f, standardized = TRUE)
+residuals.mkinfit <- function(object, standardized = FALSE, ...) {
+ res <- object$data[["residual"]]
+ if (standardized) {
+ if (object$err_mod == "const") {
+ sigma_fitted <- object$errparms["sigma"]
+ }
+ if (object$err_mod == "obs") {
+ sigma_names = paste0("sigma_", object$data[["variable"]])
+ sigma_fitted <- object$errparms[sigma_names]
+ }
+ if (object$err_mod == "tc") {
+ sigma_fitted <- sigma_twocomp(object$data[["predicted"]],
+ sigma_low = object$errparms[1],
+ rsd_high = object$errparms[2])
+ }
+ return(res / sigma_fitted)
+ }
+ return(res)
+}
+
diff --git a/R/update.mkinfit.R b/R/update.mkinfit.R
new file mode 100644
index 00000000..2f0814e0
--- /dev/null
+++ b/R/update.mkinfit.R
@@ -0,0 +1,57 @@
+#' Update an mkinfit model with different arguments
+#'
+#' This function will return an updated mkinfit object. The fitted degradation
+#' model parameters from the old fit are used as starting values for the
+#' updated fit. Values specified as 'parms.ini' and/or 'state.ini' will
+#' override these starting values.
+#'
+#' @param object An mkinfit object to be updated
+#' @param \dots Arguments to \code{\link{mkinfit}} that should replace
+#' the arguments from the original call. Arguments set to NULL will
+#' remove arguments given in the original call
+#' @param evaluate Should the call be evaluated or returned as a call
+#' @examples
+#' \dontrun{
+#' fit <- mkinfit("DFOP", subset(FOCUS_2006_D, value != 0), quiet = TRUE)
+#' update(fit, error_model = "tc")
+#' }
+#' @export
+update.mkinfit <- function(object, ..., evaluate = TRUE)
+{
+ call <- object$call
+
+ update_arguments <- match.call(expand.dots = FALSE)$...
+
+ # Get optimised ODE parameters and let parms.ini override them
+ ode_optim_names <- intersect(names(object$bparms.optim), names(object$bparms.ode))
+ ode_start <- object$bparms.optim[ode_optim_names]
+ if ("parms.ini" %in% names(update_arguments)) {
+ ode_start[names(update_arguments["parms.ini"])] <- update_arguments["parms.ini"]
+ }
+ if (length(ode_start)) update_arguments[["parms.ini"]] <- ode_start
+
+ # Get optimised values for initial states and let state.ini override them
+ state_optim_names <- intersect(names(object$bparms.optim), paste0(names(object$bparms.state), "_0"))
+ state_start <- object$bparms.optim[state_optim_names]
+ names(state_start) <- gsub("_0$", "", names(state_start))
+ if ("state.ini" %in% names(update_arguments)) {
+ state_start[names(update_arguments["state.ini"])] <- update_arguments["state.ini"]
+ }
+ if (length(state_start)) update_arguments[["state.ini"]] <- state_start
+
+ if (length(update_arguments) > 0) {
+ update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
+
+ for (a in names(update_arguments)[update_arguments_in_call]) {
+ call[[a]] <- update_arguments[[a]]
+ }
+
+ update_arguments_not_in_call <- !update_arguments_in_call
+ if(any(update_arguments_not_in_call)) {
+ call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
+ call <- as.call(call)
+ }
+ }
+ if(evaluate) eval(call, parent.frame())
+ else call
+}

Contact - Imprint