Nothing
#' Calculate likelihood profiles for the model parameters
#'
#' @param fitted \code{simode} object returned by a call to \code{\link{simode}}.
#' @param which Which parameters to estimate the profile for.
#' @param optim_type Whether to calculate the profiles based on
#' maximum-likelihood optimization only ('nls') or based on integral-matching
#' followed by maximum-likelihood optimization ('both').
#' @param step_size Step size for profiling (one value for all parameters or
#' a value for each parameter in \code{which}).
#' @param max_steps Maximum number of steps to take in each direction.
#' @param alpha Maximum (two-sided) likelihood ratio test confidence level to find.
#' @param skip_err Whether on not to stop the calculation if encountering a problem
#' with one point in the profile.
#' @param trace Report level (0-4), with higher values producing more tracing information.
#' @param save_to_log Whether to redirect output to log file. The log file
#' will be saved to tempdir().
#' @param ... Additional argument(s) for methods.
#' @details If the call to \code{\link{simode}}, which returned the fitted object given to
#' this method, included a user-defined likelihood function (with the \code{calc_nll} argument),
#' then the likelihood profiles will be calculated using this function. Otherwise,
#' the profiles will be calculated using a likelihood based on a Gaussian distribuion
#' with fixed sigma, where sigma will be estimated in the background together with
#' the rest of the model parameters.
#' @return The likelihood profiles.
#' @importFrom stats qchisq
#' @export
#'
profile.simode <- function(
fitted, which=NULL, optim_type=c("nls","both"),
step_size, max_steps=100, alpha=0.05, skip_err=TRUE, trace=0, save_to_log=FALSE, ...) {
optim_type <- match.arg(optim_type)
if(is.null(fitted$nls_pars_est)) {
stop(paste0('missing nls parameter estimates in given simode object - ',
'run simode again with \'optim_type\' set to \'nls\' or \'both\''))
}
stopifnot(is.null(which) || is.character(which))
if(is.null(which))
which <- fitted$pars
if(!pracma::isempty(setdiff(which,fitted$pars))) {
extra <- paste('[',setdiff(which,fitted$pars),']',collapse='')
stop(paste0("\'which\' contain parameters that do not appear in fitted$pars: ",
extra))
}
stopifnot(is.numeric(alpha),alpha>0,alpha<1)
stopifnot(is.numeric(step_size))
if(length(step_size)==1)
step_size <- rep(step_size,length(which))
stopifnot(length(step_size)==length(which))
stopifnot(is.numeric(max_steps),max_steps>=1)
max_delta <- qchisq(1-alpha,1)/2
if(is.null(fitted$calc_nll)) {
x0.na <- names(fitted$x0[is.na(fitted$x0)])
x0_nls <- fitted$x0
x0_nls[x0.na] <- fitted$nls_pars_est[x0.na]
pars_est_nls <- fitted$nls_pars_est[setdiff(names(fitted$nls_pars_est),x0.na)]
model_out <- solve_ode(fitted$equations, pars_est_nls, x0_nls, fitted$time)
vars <- names(fitted$equations)
model_all_vars <- unlist(lapply(vars, function(var) model_out[,var]))
obs_all_vars <- unlist(lapply(vars, function(var) fitted$obs[[var]]))
sigma <- stats::sd(model_all_vars-obs_all_vars)
sigma_name <- 'lik_normal_sigma'
# fitted$fixed[sigma_name] <- sigma
fitted$start[sigma_name] <- sigma
if(is.null(fitted$lower)) {
fitted$lower[fitted$pars] <- -Inf
}
fitted$lower[sigma_name] <- 0
if(is.null(fitted$upper)) {
fitted$upper[fitted$pars] <- Inf
}
fitted$upper[sigma_name] <- sigma*10
fitted$nls_pars_est[sigma_name] <- sigma
fitted$im_pars_est[sigma_name] <- NA
fitted$pars <- c(fitted$pars,sigma_name)
fitted$likelihood_pars <- c(fitted$likelihood_pars,sigma_name)
fitted$calc_nll <- calc_nll_normal
fitted$extra_args <- c(fitted$extra_args,sigma_name=sigma_name)
}
if(save_to_log) {
log_file <- file.path(tempdir(),"simode.log")
sink(file=log_file, append=TRUE)
cat(noquote((paste0("Call to profile on [", Sys.time(), "]:\n"))))
}
on.exit(expr=simode_cleanup(save_to_log))
res <- list(call=match.call(), alpha=alpha, pars=which,
nls_pars_est=fitted$nls_pars_est[which],
lower=fitted$lower[which],
upper=fitted$upper[which])
class(res) <- 'profile.simode'
res$profiles <- list()
tryCatch({
for(i in 1:length(which)) {
par <- which[i]
res$profiles[[par]] <- calc_profile(fitted, par, optim_type, step_size[i],
max_steps, max_delta, skip_err, trace)
}
},
warning=function(w) {
print(w)
},
error=function(e) {
print(e)
res <- NULL
}
)
return (res)
}
calc_profile <- function(x, par, optim_type, step_size, max_steps,
max_delta, skip_err, trace) {
if(trace > 0)
cat(noquote((paste0("Calculating profile for parameter [", par, "]\n"))))
is_lik_par <- F
if(!pracma::isempty(intersect(x$likelihood_pars,par)))
is_lik_par <- T
x$ctrl$optim_type <- optim_type
x$pars <- setdiff(x$pars,par)
x$nlin_pars <- setdiff(x$nlin_pars,par)
x$likelihood_pars <- setdiff(x$likelihood_pars,par)
par_min <- unname(x$lower[par])
par_max <- unname(x$upper[par])
x$lower <- x$lower[x$pars]
x$upper <- x$upper[x$pars]
x$start <- x$start[setdiff(names(x$start),par)]
x$ctrl$job_id <- NULL
x$ctrl$trace <- trace-1
nls_par_est <- x$nls_pars_est[par]
par_vals <- c(unname(nls_par_est))
init_pars_est_l <- x$nls_pars_est[x$pars]
init_pars_est_r <- x$nls_pars_est[x$pars]
step_res <- calc_profile_step(x, par, nls_par_est, is_lik_par,
init_pars_est_l, skip_err, trace)
min_nll <- step_res$nll
nll_vals <- min_nll
profile_l_done <- F
profile_r_done <- F
i <- 0
while(T)
{
i <- i+1
if(trace > 1)
cat(noquote((paste0("Performing step [", i, "] in calculating profile for parameter [", par ,"]\n"))))
if(!profile_l_done) {
par_val <- nls_par_est-i*step_size
if(!is.null(par_min) && par_val <= par_min) {
par_val <- par_min
profile_l_done <- T
}
step_res <-
calc_profile_step(x, par, par_val, is_lik_par,
init_pars_est_l, skip_err, trace)
if(!is.null(step_res)) {
par_vals <- c(unname(par_val),par_vals)
nll <- step_res$nll
nll_vals <- c(nll,nll_vals)
if(min_nll > nll)
min_nll <- nll
if(nll - min_nll > max_delta)
profile_l_done <- T
init_pars_est_l <- step_res$nls_pars_est
if(trace > 1)
cat(noquote((paste0("Parameter [", par, "], Step [", i, "] L, Delta NLL [",
round(nll-min_nll,2) ,"]\n"))))
}
}
if(!profile_r_done) {
par_val <- nls_par_est+i*step_size
if(!is.null(par_max) && par_val >= par_max) {
par_val <- par_max
profile_r_done <- T
}
step_res <-
calc_profile_step(x, par, par_val, is_lik_par,
init_pars_est_r, skip_err, trace)
if(!is.null(step_res)) {
par_vals <- c(par_vals,unname(par_val))
nll <- step_res$nll
nll_vals <- c(nll_vals,nll)
if(min_nll > nll)
min_nll <- nll
if(nll - min_nll > max_delta)
profile_r_done <- T
init_pars_est_r <- step_res$nls_pars_est
if(trace > 1)
cat(noquote((paste0("Parameter [", par, "], Step [", i, "] R, Delta NLL [",
round(nll-min_nll,2) ,"]\n"))))
}
}
if(profile_l_done && profile_r_done)
break
if(i >= max_steps) {
if(trace > 0) {
cat(noquote((paste0("Calculating profile for parameter [", par,
"] - max step has been reached\n"))))
}
break
}
}
return (list(par_vals=par_vals, nll_vals=nll_vals))
}
calc_profile_step <- function(x, par, par_val, is_lik_par,
init_pars_est, skip_err, trace) {
if(is_lik_par)
x$fixed<- c(x$fixed,par_val)
else if(par %in% names(x$x0))
x$x0[par] <- par_val
else
x$equations <- fix_pars(x$equations, par_val)
x$start <- init_pars_est
if(x$ctrl$optim_type=='both' && x$decouple_equations)
x <- simode_decouple(x)
else
x <- simode_impl(x)
if(is.null(x)) {
if(!skip_err)
stop(paste0('Error occurred during profile calculation for parameter [', par, ']\n'))
return (NULL)
}
res <- list(nll=x$nls_loss, nls_pars_est=x$nls_pars_est)
return (res)
}
#' Plot the likelihood profiles for the model parameters
#'
#' @param x \code{profile.simode} object returned by a call to \code{\link{profile}}.
#' @param which Which parameters to plot the likelihood profiles for.
#' If empty, the plot will include all of the parameters in x.
#' @param mfrow A vector of the form c(nr,nc) setting the layout of
#' subplots in one plot (see also \code{\link{par}}).
#' @param cols List of colors for each element of the plot.
#' @param ... Additional argument(s) for methods.
#' @importFrom graphics abline
#' @export
#'
plot.profile.simode <- function(x, which=NULL, mfrow=par('mfrow'),
cols=list(fit="blue",threshold="red"), ...) {
max_delta <- qchisq(1-x$alpha,1)/2
if(is.null(which))
which <- x$pars
if(!pracma::isempty(setdiff(which,x$pars))) {
extra <- paste('[',setdiff(which,x$pars),']',collapse='')
stop(paste0("\'which\' contain parameters that do not appear in x: ",extra))
}
old.par <- par(mar=c(4,4,2,2),mfrow=mfrow)
on.exit(par(old.par))
for(i in 1:length(which)) {
par <- which[i]
plot(x$profiles[[par]]$par_vals,
x$profiles[[par]]$nll_vals-min(x$profiles[[par]]$nll_vals,na.rm=TRUE),
xlab=par,ylab=expression(paste(Delta,' NLL')),
col=cols[['fit']], type='o', ...)
abline(h=max_delta, col=cols[['threshold']])
}
}
#' Print method for \code{profile.simode} objects
#'
#' @param x The \code{profile.simode} object.
#' @param ... Additional argument(s)for methods.
#' @export
print.profile.simode <- function(x, ...) {
cat("call:\n")
print(x$call)
cat("alpha:\n")
cat(x$alpha,"\n")
cat("profiles:\n")
print(x$profiles)
}
#' Calculates confidence intervals for the model parameters
#'
#' Calculates confidence intervals for the model parameters,
#' based on the given likelihood profiles
#' @param object A fitted model object: \code{profile.simode} object returned by a call to \code{\link{profile}}.
#' @param parm A specification of which parameters are to be given confidence intervals (named vector).
#' If missing, all parameters are considered.
#' @param level The confidence level required.
#' @param ... Additional argument(s) for methods.
#' @return The confidence intervals.
#' @importFrom stats approx
#' @export
#'
confint.profile.simode<- function(object, parm=NULL, level=0.95, ...) {
stopifnot(is.null(parm) || (is.character(parm) && is.vector(parm)))
pars <- parm
if(is.null(pars))
pars <- c(object$pars)
if(!pracma::isempty(setdiff(pars,object$pars))) {
extra <- paste('[',setdiff(pars,object$pars),']',collapse='')
stop(paste0("\'parm\' contain parameters that do not appear in object$pars: ",
extra))
}
stopifnot(is.numeric(level),level>0,level<1)
if(1-level<object$alpha)
stop('1-level must be greater than value of alpha used in profile call')
max_delta <- qchisq(level,1)/2
res <- list(call=match.call(), level=level)
intervals <- list()
for(par in pars) {
par_vals <- object$profiles[[par]]$par_vals
nll_vals <- object$profiles[[par]]$nll_vals
par_min <- object$lower[[par]]
par_max <- object$upper[[par]]
intervals[[par]] <- calc_confint(par_vals,nll_vals,max_delta,par_min,par_max)
}
ci_df <- data.frame(par=pars)
if(!is.null(object$nls_pars_est))
ci_df$nls_est <- object$nls_pars_est
lp <- length(pars)
ci_df$lower <- apply(as.matrix(1:lp),1,function(i) intervals[[i]]$lower)
ci_df$upper <- apply(as.matrix(1:lp),1,function(i) intervals[[i]]$upper)
res$intervals <- ci_df
class(res) <- 'confint.simode'
return(res)
}
calc_confint <- function(par_vals, nll_vals, max_delta, par_min, par_max) {
ci_lower <- NA
ci_upper <- NA
delta_nll <- nll_vals-min(nll_vals)
interp1 <- approx(x=par_vals, y=delta_nll, n=100)
x_vals <- interp1$x
y_vals <- interp1$y
t_vals <- as.numeric(y_vals < max_delta)
b_vals <- which(t_vals==1)
if(length(b_vals) >= 1) {
ind.min <- min(b_vals)
ind.max <- max(b_vals)
if(ind.min > 1)
ci_lower <- mean(x_vals[(ind.min-1):ind.min])
if(ind.max < length(x_vals))
ci_upper <- mean(x_vals[ind.max:(ind.max+1)])
}
if(is.na(ci_lower) && !is.null(par_min) && par_vals[1]==par_min)
ci_lower <- par_min
if(is.na(ci_upper) && !is.null(par_max) && par_vals[length(par_vals)]==par_max)
ci_upper <- par_max
res <- list(lower=ci_lower,upper=ci_upper)
return (res)
}
#' Plot confidence intervals for the model parameters
#'
#' Plot confidence intervals for the model parameters of a \code{simode} object,
#' calculated based on likelihood profiles
#' @param x \code{confint.simode} object returned by a call to \code{\link{confint}}
#' @param which Which parameters to plot the confidence intervals for.
#' If empty, the plot will include all of the parameters in x.
#' @param pars_true The true parameter values (if are known).
#' @param legend Whether or not to add a figure legend.
#' @param cols List of colors for each element of the plot.
#' @param ... Additional argument(s) for methods.
#' @export
#'
plot.confint.simode<- function(x, which=NULL, pars_true=NULL, legend=FALSE,
cols=list(fit="blue",true="black"), ...) {
pars <- which
if(is.null(pars))
pars <- as.vector(x$intervals$par)
if(!is.null(pars_true))
stopifnot(is.numeric(pars_true), length(pars_true)==length(pars))
lp <- length(pars)
pind <- which(as.vector(x$intervals$par) %in% pars)
y_min <- Inf
y_max <- -Inf
nls_pars_est <- apply(as.matrix(1:lp),1,function(i) x$intervals[pind[i],]$nls_est)
lower_vals <- apply(as.matrix(1:lp),1,function(i) x$intervals[pind[i],]$lower)
upper_vals <- apply(as.matrix(1:lp),1,function(i) x$intervals[pind[i],]$upper)
names(lower_vals) <- pars
names(upper_vals) <- pars
y_min <- min(c(y_min,lower_vals),na.rm=TRUE)
y_max <- max(c(y_max,upper_vals),na.rm=TRUE)
if(!is.null(pars_true)) {
y_min <- min(c(y_min,pars_true),na.rm=TRUE)
y_max <- max(c(y_max,pars_true),na.rm=TRUE)
}
legend_text <- c()
legend_cols <- c()
legend_lty <- c()
legend_pch <- c()
if(legend) {
old.par <- par(xpd = T, mar = par()$mar + c(0,0,0,4))
on.exit(par(old.par))
}
plot(-100, -100, xlim=c(0.5,lp+0.5), ylim=c(y_min,y_max), xaxt="n",
xlab='parameter', ylab='value', ...)
if(!is.null(pars_true)) {
points(x=1:lp, pars_true, col=cols[['true']], pch=4, xaxt="n")
legend_text <- c(legend_text,'true')
legend_cols <- c(legend_cols,cols[['true']])
legend_lty <- c(legend_lty,0)
legend_pch <- c(legend_pch,4)
}
x_vals <- 1:lp
points(x=x_vals, nls_pars_est, col=cols[['fit']], xaxt="n")
for(i in 1:lp) {
if(!is.na(lower_vals[[i]]) && !is.na(upper_vals[[i]])) {
lines(x=c(x_vals[i],x_vals[i]),
y=c(lower_vals[[i]], upper_vals[[i]]),
col=cols[['fit']])
}
}
legend_text <- c(legend_text,'nls_est','ci')
legend_cols <- c(legend_cols,cols[['fit']],cols[['fit']])
legend_lty <- c(legend_lty,0,1)
legend_pch <- c(legend_pch,1,NA)
axis(side=1,at=1:lp,labels=pars)
if(legend) {
legend(lp+1, y_max, cex=0.8, bty="n",
legend_text ,col=legend_cols, lty=legend_lty, pch=legend_pch)
}
}
#' Print method for \code{confint.simode} objects
#'
#' @param x The \code{confint.simode} object.
#' @param ... Additional argument(s) for methods.
#' @export
print.confint.simode <- function(x, ...) {
cat("call:\n")
print(x$call)
cat("level:\n")
cat(x$level,"\n")
cat("intervals:\n")
print(x$intervals)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.