Nothing
simode_multi_worker <- function(...) {
arg_list <- list(...)
cl <- arg_list$cl
parallel::clusterExport(cl, "arg_list", envir=environment())
parallel::clusterExport(cl, ls("package:simode"), envir=environment())
obs_sets <- length(arg_list$obs)
simode_fits <- parallel::parLapply(cl, 1:obs_sets, function(i) {
with(arg_list, {
simode_ctrl$job_id <- i
args <-
c(list(equations=equations, pars=pars, time=time, obs=obs[[i]],
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl), extra_args)
do.call("simode", args)
})
})
return (simode_fits)
}
simode_multi <- function(equations, pars, time, obs,
nlin_pars, likelihood_pars,
fixed, start, lower, upper,
im_method, decouple_equations,
gen_obs, calc_nll, scale_pars, simode_ctrl, ...) {
if(simode_ctrl$obs_sets_fit == 'together')
return (simode_multi_together(equations=equations, pars=pars, time=time, obs=obs,
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...))
else if(simode_ctrl$obs_sets_fit == 'separate_x0')
return (simode_multi_separate_x0(equations=equations, pars=pars, time=time, obs=obs,
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...))
return (simode_multi_separate(equations=equations, pars=pars, time=time, obs=obs,
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...))
}
simode_multi_together <- function(equations, pars, time, obs,
nlin_pars, likelihood_pars,
fixed, start, lower, upper,
im_method, decouple_equations,
gen_obs, calc_nll, scale_pars,
simode_ctrl, ...) {
obs_mean <- list()
obs_sets <- length(obs)
for(i in 1:length(obs[[1]])) {
obs1 <- lapply(1:obs_sets, function(j) obs[[j]][i])
obs1 <- matrix(unlist(obs1), ncol = obs_sets)
obs_mean[[i]] <- unlist(lapply(1:nrow(obs1), function(k) mean(obs1[k,])))
}
names(obs_mean) <- names(obs[[1]])
obs <- obs_mean
return (simode(equations=equations, pars=pars, time=time, obs=obs,
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...))
}
simode_multi_separate_x0 <- function(equations, pars, time, obs,
nlin_pars, likelihood_pars,
fixed, start, lower, upper,
im_method, decouple_equations,
gen_obs, calc_nll, scale_pars, simode_ctrl, ...) {
sep <- simode_ctrl$decouple_sep
eq_names <- names(equations)
vars <- names(obs[[1]])
equations_m <- c()
obs_m <- c()
time_m <- c()
if(!is.list(time))
time_m <- time
obs_sets <- length(obs)
for(i in 1:obs_sets) {
vars1 <- paste0(vars,sep,i)
names(vars1) <- vars
obs1 <- obs[[i]]
names(obs1) <- vars1
obs_m <- c(obs_m,obs1)
if(is.list(time)) {
stopifnot(length(time)==length(obs1))
time1 <- time
names(time1) <- vars1
time_m <- c(time_m,time1)
}
equations1 <- fix_pars(equations,vars1)
names(equations1) <- paste0(eq_names,sep,i)
equations_m <- c(equations_m,equations1)
}
pars_m <- pars
pars_x0 <- intersect(pars,eq_names)
if(!pracma::isempty(pars_x0)) {
pars_x0_ex <- paste0(pars_x0,sep,unlist(lapply(1:obs_sets,function(i) rep(i,length(pars_x0)))))
pars_m <- c(setdiff(pars,pars_x0),pars_x0_ex)
}
nlin_pars_m <- nlin_pars
nlin_pars_x0 <- intersect(nlin_pars,eq_names)
if(!pracma::isempty(nlin_pars_x0)) {
nlin_pars_x0_ex <- paste0(nlin_pars_x0,sep,unlist(lapply(1:obs_sets,function(i) rep(i,length(nlin_pars_x0)))))
nlin_pars_m <- c(setdiff(nlin_pars,nlin_pars_x0),nlin_pars_x0_ex)
}
if(!is.list(fixed)) {
fixed_m <- fixed
fixed_x0_names <- intersect(names(fixed),eq_names)
if(!pracma::isempty(fixed_x0_names)) {
fixed_x0 <- rep(fixed[fixed_x0_names],obs_sets)
names(fixed_x0) <-
paste0(fixed_x0_names,sep,
unlist(lapply(1:obs_sets,function(i) rep(i,length(equations)))))
fixed_m <- c(fixed[setdiff(names(fixed),fixed_x0_names)],fixed_x0)
}
}
else {
stopifnot(length(fixed)==obs_sets)
fixed_m <- c()
for(i in 1:obs_sets) {
fixed_x0_names <- intersect(names(fixed[[i]]),eq_names)
if(!pracma::isempty(fixed_x0_names)) {
fixed_x0 <- fixed[[i]][fixed_x0_names]
names(fixed_x0) <- paste0(fixed_x0_names,sep,i)
fixed_m <- c(fixed_m,fixed[[i]][setdiff(names(fixed[[i]]),fixed_x0_names)],fixed_x0)
}
else
fixed_m <- c(fixed_m,fixed[[i]])
}
}
start_m <- start
start_x0_names <- intersect(names(start),eq_names)
if(!pracma::isempty(start_x0_names)) {
start_x0 <- rep(start[start_x0_names],obs_sets)
names(start_x0) <-
paste0(start_x0_names,sep,
unlist(lapply(1:obs_sets,function(i) rep(i,length(equations)))))
start_m <- c(start[setdiff(names(start),start_x0_names)],start_x0)
}
lower_m <- lower
lower_x0_names <- intersect(names(lower),eq_names)
if(!pracma::isempty(lower_x0_names)) {
lower_x0 <- rep(lower[lower_x0_names],obs_sets)
names(lower_x0) <-
paste0(lower_x0_names,sep,
unlist(lapply(1:obs_sets,function(i) rep(i,length(equations)))))
lower_m <- c(lower[setdiff(names(lower),lower_x0_names)],lower_x0)
}
upper_m <- upper
upper_x0_names <- intersect(names(upper),eq_names)
if(!pracma::isempty(upper_x0_names)) {
upper_x0 <- rep(upper[upper_x0_names],obs_sets)
names(upper_x0) <-
paste0(upper_x0_names,sep,
unlist(lapply(1:obs_sets,function(i) rep(i,length(equations)))))
upper_m <- c(upper[setdiff(names(upper),upper_x0_names)],upper_x0)
}
x <- simode(equations=equations_m, pars=pars_m, time=time_m, obs=obs_m,
nlin_pars=nlin_pars_m, likelihood_pars=likelihood_pars,
fixed=fixed_m, start=start_m, lower=lower_m, upper=upper_m,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...)
if(is.null(x))
return (NULL)
x$equations <- equations
x$pars <- pars
x$nlin_pars <- nlin_pars
x$fixed <- fixed
x$start <- start
x$lower <- lower
x$upper <- upper
x_list <- list()
for(i in 1:obs_sets) {
x_list[[i]] <- x
x_list[[i]]$obs <- obs[[i]]
x_list[[i]]$x0 <- NULL
x_list[[i]]$x0[pars_x0] <- x$x0[paste0(pars_x0,sep,i)]
fixed_x0_names <- intersect(names(fixed[[i]]),eq_names)
if(!pracma::isempty(fixed_x0_names)) {
x_list[[i]]$x0[fixed_x0_names] <- fixed[[i]][fixed_x0_names]
}
if(!is.null(x$im_pars_est)) {
x_list[[i]]$im_pars_est <- x$im_pars_est[setdiff(pars,pars_x0)]
x_list[[i]]$im_pars_est[pars_x0] <- x$im_pars_est[paste0(pars_x0,sep,i)]
x_list[[i]]$im_pars_est <- x_list[[i]]$im_pars_est[pars]
x_list[[i]]$im_smooth$val <- x$im_smooth$val[,paste0(eq_names,sep,i)]
colnames(x_list[[i]]$im_smooth$val) <- eq_names
}
if(!is.null(x$nls_pars_est)) {
x_list[[i]]$nls_pars_est <- x$nls_pars_est[setdiff(pars,pars_x0)]
x_list[[i]]$nls_pars_est[pars_x0] <- x$nls_pars_est[paste0(pars_x0,sep,i)]
x_list[[i]]$nls_pars_est <- x_list[[i]]$nls_pars_est[pars]
}
}
class(x_list) <- 'list.simode'
return (x_list)
}
simode_multi_separate <- function(equations, pars, time, obs,
nlin_pars, likelihood_pars,
fixed, start, lower, upper,
im_method, decouple_equations,
gen_obs, calc_nll, scale_pars, simode_ctrl, ...) {
if(simode_ctrl$parallel==TRUE && !requireNamespace("parallel", quietly=TRUE)) {
warning("parallel package not installed - running sequentially",
immediate.=TRUE)
simode_ctrl$parallel <- F
}
log_file <- NULL
if(simode_ctrl$save_to_log) {
log_file <- file.path(tempdir(),"simode.log")
simode_ctrl$save_to_log <- F
}
simode_fits <- NULL
obs_sets <- length(obs)
tryCatch(
{
if(simode_ctrl$parallel==FALSE) {
if(!is.null(log_file)) {
sink(file=log_file, append=FALSE)
cat(noquote((paste0("Call to simode on [", Sys.time(), "]:\n"))))
}
simode_fits <- list()
for(i in 1:obs_sets) {
simode_ctrl$job_id <- i
simode_fits[[i]] <-
simode(equations=equations, pars=pars, time=time, obs=obs[[i]],
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, ...)
}
}
else {
cl <- parallel::makeCluster(obs_sets, outfile=log_file)
extra_args <- list(...)
simode_fits <- simode_multi_worker(
cl=cl, equations=equations, pars=pars, time=time, obs=obs,
nlin_pars=nlin_pars, likelihood_pars=likelihood_pars,
fixed=fixed, start=start, lower=lower, upper=upper,
im_method=im_method, decouple_equations=decouple_equations,
gen_obs=gen_obs, calc_nll=calc_nll, scale_pars=scale_pars,
simode_ctrl=simode_ctrl, extra_args=extra_args)
}
},
error = function(e) { print(e) },
finally={
if(simode_ctrl$parallel==TRUE && !is.null(cl))
parallel::stopCluster(cl)
}
)
if(!is.null(simode_fits))
class(simode_fits) <- "list.simode"
return (simode_fits)
}
#' Summary method for \code{list.simode} objects
#'
#' @param object \code{list.simode} object returned by a call to \code{\link{simode}}
#' with \code{obs_sets>1}
#' @param sum_mean_sd Whether to calculate mean and standard deviation
#' for the parameter estimates in the fits included in the given \code{object}.
#' To be used when \code{object} is the result of fitting monte-carlo simulations.
#' @param pars_true The true parameter values (relevant only for
#' when \code{sum_mean_sd=TRUE}). When given, the summary will also include the
#' bias and RMSE for each parameter estimate.
#' @param digits The number of significant digits to use.
#' @param ... Additional argument(s) for methods.
#' @return The mean and standard deviation for the
#' loss values and parameter estimates obtained from the integral-matching
#' and nonlinear least squares optimizations. If \code{pars_true} is given, then
#' will also calculate bias and RMSE for the parameter estimates.
#' @export
#'
summary.list.simode <- function(object, sum_mean_sd=FALSE, pars_true=NULL,
digits=max(3, getOption("digits")-3), ...) {
pars <- object[[1]]$pars
obs_sets_fit <- object[[1]]$ctrl$obs_sets_fit
if(!is.null(pars_true)) {
if(!sum_mean_sd)
warning('parameter pars_true is only relevant when sum_mean_sd==TRUE')
else
stopifnot(!is.null(names(pars_true)),all(names(pars_true) %in% pars))
}
simode_num <- length(object)
im_est <- matrix(NA,nrow=simode_num,ncol=length(pars))
nls_est <- matrix(NA,nrow=simode_num,ncol=length(pars))
im_loss_vals <- rep(NA,simode_num)
nls_loss_vals <- rep(NA,simode_num)
for(i in 1:simode_num) {
if(!is.null(object[[i]]$im_pars_est))
im_est[i,] <- object[[i]]$im_pars_est
if(!is.null(object[[i]]$nls_pars_est))
nls_est[i,] <- object[[i]]$nls_pars_est
if(!is.null(object[[i]]$im_loss))
im_loss_vals[i] <- object[[i]]$im_loss
if(!is.null(object[[i]]$nls_loss))
nls_loss_vals[i] <- object[[i]]$nls_loss
}
summary <- list()
if(!sum_mean_sd) {
if(!all(is.na(im_est))) {
summary$im_est <- signif(im_est,digits)
colnames(summary$im_est) <- pars
}
if(!all(is.na(nls_est))) {
summary$nls_est <- signif(nls_est,digits)
colnames(summary$nls_est) <- pars
}
if(obs_sets_fit == "separate") {
if(any(!is.na(im_loss_vals)))
summary$im_loss <- signif(sum(im_loss_vals),digits)
if(any(!is.na(nls_loss_vals)))
summary$nls_loss <- signif(sum(nls_loss_vals),digits)
}
else
{
if(any(!is.na(im_loss_vals)))
summary$im_loss <- signif(im_loss_vals[1],digits)
if(any(!is.na(nls_loss_vals)))
summary$nls_loss <- signif(nls_loss_vals[1],digits)
}
}
else {
im_loss_mean <- signif(mean(im_loss_vals),digits)
im_loss_sd <- signif(stats::sd(im_loss_vals),digits)
nls_loss_mean <- signif(mean(nls_loss_vals),digits)
nls_loss_sd <- signif(stats::sd(nls_loss_vals),digits)
if(!is.na(im_loss_mean)) {
summary$loss$im_mean <- im_loss_mean
summary$loss$im_sd <- im_loss_sd
}
if(!is.na(nls_loss_mean)) {
summary$loss$nls_mean <- nls_loss_mean
summary$loss$nls_sd <- nls_loss_sd
}
summary$loss <- as.data.frame(summary$loss)
im_est_mean <- signif(apply(im_est,2,mean),digits)
im_est_sd <- signif(apply(im_est,2,stats::sd),digits)
names(im_est_mean) <- pars
nls_est_mean <- signif(apply(nls_est,2,mean),digits)
nls_est_sd <- signif(apply(nls_est,2,stats::sd),digits)
names(nls_est_mean) <- pars
summary$est <- data.frame(par=pars)
if(!is.null(pars_true)) {
summary$est$true <- signif(pars_true[pars], digits)
im_est_bias <- signif(im_est_mean - pars_true[pars],digits)
nls_est_bias <- signif(nls_est_mean - pars_true[pars],digits)
im_est_err <- im_est-matrix(rep(pars_true[pars],simode_num),nrow=simode_num,byrow=TRUE)
im_est_rmse <- signif(sqrt(apply(im_est_err^2,2,mean)),digits)
nls_est_err <- nls_est-matrix(rep(pars_true[pars],simode_num),nrow=simode_num,byrow=TRUE)
nls_est_rmse <- signif(sqrt(apply(nls_est_err^2,2,mean)),digits)
}
if(any(!is.na(im_est_mean))) {
summary$est$im_mean <- im_est_mean
summary$est$im_sd <- im_est_sd
if(!is.null(pars_true)) {
summary$est$im_bias <- im_est_bias
summary$est$im_rmse <- im_est_rmse
}
}
if(any(!is.na(nls_est_mean))) {
summary$est$nls_mean <- nls_est_mean
summary$est$nls_sd <- nls_est_sd
if(!is.null(pars_true)) {
summary$est$nls_bias <- nls_est_bias
summary$est$nls_rmse <- nls_est_rmse
}
}
}
class(summary) <- "summary.list.simode"
return(summary)
}
#' Plot the fit/estimates of a \code{list.simode} object
#'
#' Plot the fit or parameter estimates obtained from a call to \code{simode}
#' with \code{obs_sets>1}.
#'
#' @param x \code{list.simode} object.
#' @param type Type of plot - 'fit' to plot the fitted equations
#' and 'est' to plot the parameter estimates.
#' @param show Whether to plot the fit/estimates obtained
#' using nonlinear least squares ('nls'), integral-matching ('im') or both ('both').
#' @param which Which variables to plot in case \code{type='fit'}, or which
#' parameters to plot in case \code{type='est'}. If empty, the plot will include
#' all of the variables/parameters in x.
#' @param pars_true The true parameter values (if are known). Should be
#' named using the parameter names. If given, the true values for
#' the variables/parameters will be added to the plot.
#' @param time The time points to use for the fitted curves
#' (relevant only for \code{type='fit'}).
#' If not given then the time points of the observations in x will be used.
#' @param plot_mean_sd Plot the mean and standard deviation for the fit/estimates
#' of the \code{simode} objects in \code{x}. To be used when \code{x} is the
#' result of fitting monte-carlo simulations.
#' @param plot_im_smooth Whether or not to plot the smoothed curves created and
#' used by the integral-matching procedure (relevant only for \code{type='fit'}).
#' @param legend Whether or not to add a figure legend.
#' @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 plot points lines axis par
#' @importFrom pracma errorbar
#' @export
#'
plot.list.simode <-
function(x, type=c('fit','est'), show=c('nls','im','both'),
which=NULL, pars_true=NULL, time=NULL,
plot_mean_sd=FALSE, plot_im_smooth=FALSE,
legend=FALSE, mfrow=par('mfrow'),
cols=list(nls_fit="blue",im_fit="green", true="black",
obs="red", im_smooth="magenta"), ...) {
type <- match.arg(type)
show <- match.arg(show)
if((show=='nls' || show=='both') && is.null(x[[1]]$nls_pars_est)) {
warning("no nls estimates in this simode object - will show only \'im\' fit\\estimates")
show <- 'im'
}
if((show=='im' || show=='both') && is.null(x[[1]]$im_pars_est)) {
warning("no im estimates in this simode object - will show only \'nls\' fit\\estimates")
show <- 'nls'
}
if(length(x)>1 && !plot_mean_sd) {
for(i in 1:length(x)) {
xi <- x[i]
class(xi) <- 'list.simode'
if(type == 'fit')
plot_list_simode_fit(x=xi,show=show,which=which,pars_true=pars_true,
time=time,plot_im_smooth=plot_im_smooth,plot_mean_sd=plot_mean_sd,
legend=legend, mfrow=mfrow, cols=cols, fit_num=i, ...)
else #if(type == 'est')
plot_list_simode_est(x=xi,show=show,which=which,pars_true=pars_true,
plot_mean_sd=plot_mean_sd, legend=legend, cols=cols, fit_num=i, ...)
}
}
else {
if(type == 'fit')
plot_list_simode_fit(x=x,show=show,which=which,pars_true=pars_true,
time=time,plot_im_smooth=plot_im_smooth,plot_mean_sd=plot_mean_sd,
legend=legend, mfrow=mfrow, cols=cols, ...)
else #if(type == 'est')
plot_list_simode_est(x=x,show=show,which=which,pars_true=pars_true,
plot_mean_sd=plot_mean_sd,legend=legend, cols=cols, ...)
}
}
plot_list_simode_est <-
function(x, show=c('nls','im','both'),
which=NULL, pars_true=NULL, plot_mean_sd=FALSE, legend=FALSE,
cols=list(nls_fit="blue",im_fit="green", true="black"),
fit_num=NULL, ...) {
show <- match.arg(show)
plot_nls_est <- (show=='nls' || show=='both')
plot_im_est <- (show=='im' || show=='both')
pars <- which
if(is.null(pars))
pars <- x[[1]]$pars
if(length(x)==1)
plot_mean_sd <- F
if(!is.null(pars_true)) {
stopifnot(is.numeric(pars_true))
stopifnot(all(pars %in% names(pars_true)))
}
lp <- length(pars)
pind <- which(as.vector(x[[1]]$pars) %in% pars)
if(plot_mean_sd) {
sum <- summary(x, sum_mean_sd=TRUE)
im_est_mean <- apply(as.matrix(1:lp),1,function(i) sum$est[pind[i],]$im_mean)
im_est_sd <- apply(as.matrix(1:lp),1,function(i) sum$est[pind[i],]$im_sd)
nls_est_mean <- apply(as.matrix(1:lp),1,function(i) sum$est[pind[i],]$nls_mean)
nls_est_sd <- apply(as.matrix(1:lp),1,function(i) sum$est[pind[i],]$nls_sd)
names(im_est_mean) <- pars
names(im_est_sd) <- pars
names(nls_est_mean) <- pars
names(nls_est_sd) <- pars
im_est_sd[which(is.na(im_est_sd))] <- 0
nls_est_sd[which(is.na(nls_est_sd))] <- 0
}
y_min <- Inf
y_max <- -Inf
if(!is.null(pars_true)) {
y_min <- min(pars_true[pars],na.rm=TRUE)
y_max <- max(pars_true[pars],na.rm=TRUE)
}
if(plot_mean_sd) {
if(plot_nls_est) {
min_nls <- nls_est_mean[pars] - nls_est_sd[pars]
max_nls <- nls_est_mean[pars] + nls_est_sd[pars]
y_min <- min(c(y_min,min_nls),na.rm=TRUE)
y_max <- max(c(y_max,max_nls),na.rm=TRUE)
}
if(plot_im_est) {
min_im <- im_est_mean[pars] - im_est_sd[pars]
max_im <- im_est_mean[pars] + im_est_sd[pars]
y_min <- min(c(y_min,min_im),na.rm=TRUE)
y_max <- max(c(y_max,max_im),na.rm=TRUE)
}
}
else {
if(plot_nls_est) {
min_nls <- min(unlist(lapply(1:length(x),function(i) min(x[[i]]$nls_pars_est[pars]))))
max_nls <- max(unlist(lapply(1:length(x),function(i) max(x[[i]]$nls_pars_est[pars]))))
y_min <- min(c(y_min,min_nls),na.rm=TRUE)
y_max <- max(c(y_max,max_nls),na.rm=TRUE)
}
if(plot_im_est) {
min_im <- min(unlist(lapply(1:length(x),function(i) min(x[[i]]$im_pars_est[pars]))))
max_im <- max(unlist(lapply(1:length(x),function(i) max(x[[i]]$im_pars_est[pars]))))
y_min <- min(c(y_min,min_im),na.rm=TRUE)
y_max <- max(c(y_max,max_im),na.rm=TRUE)
}
}
legend_text <- c()
legend_cols <- c()
legend_lty <- c()
legend_pch <- c()
if(legend) {
old.par <- par(xpd=TRUE, mar=c(4,4,2,6))
on.exit(par(old.par))
}
ylab <- ifelse(is.null(fit_num),'value',paste0('value (fit-',fit_num,')'))
plot(-100, -100, xlim=c(0.5,lp+0.5), ylim=c(y_min,y_max), xaxt="n",
xlab='parameter', ylab=ylab, ...)
if(!is.null(pars_true)) { # && !is.list(pars_true)) {
points(x=1:lp, pars_true[pars], 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)
}
if(plot_mean_sd) {
if(plot_im_est && !is.null(im_est_mean)) {
x_vals <- 1:lp
if(plot_nls_est) # || is.list(pars_true))
x_vals <- x_vals-0.1
points(x=x_vals, im_est_mean[pars], col=cols[['im_fit']], xaxt="n",...)
errorbar(x=x_vals, im_est_mean[pars], yerr=im_est_sd[pars],
bar.col=cols[['im_fit']], grid=FALSE, add=TRUE, xaxt="n")
legend_text <- c(legend_text,'im_mean','im_sd')
legend_cols <- c(legend_cols,cols[['im_fit']],cols[['im_fit']])
legend_lty <- c(legend_lty,0,1)
legend_pch <- c(legend_pch,1,NA)
}
if(plot_nls_est && !is.null(nls_est_mean)) {
x_vals <- 1:lp
if(plot_im_est)
x_vals <- x_vals+0.1
points(x=x_vals, nls_est_mean[pars], col=cols[['nls_fit']], xaxt="n",...)
errorbar(x=x_vals, nls_est_mean[pars], yerr=nls_est_sd[pars],
bar.col=cols[['nls_fit']], grid=FALSE, add=TRUE, xaxt="n")
legend_text <- c(legend_text,'nls_mean','nls_sd')
legend_cols <- c(legend_cols,cols[['nls_fit']],cols[['nls_fit']])
legend_lty <- c(legend_lty,0,1)
legend_pch <- c(legend_pch,1,NA)
}
}
else {
for(i in 1:length(x)) {
if(plot_im_est) {
x_vals <- 1:lp
if(plot_nls_est)
x_vals <- x_vals-0.1
points(x=x_vals, x[[i]]$im_pars_est[pars], col=cols[['im_fit']], xaxt="n",...)
if(i==1) {
legend_text <- c(legend_text,'im_est')
legend_cols <- c(legend_cols,cols[['im_fit']])
legend_lty <- c(legend_lty,0)
legend_pch <- c(legend_pch,1)
}
}
if(plot_nls_est) {
x_vals <- 1:lp
if(plot_im_est)
x_vals <- x_vals+0.1
points(x=x_vals, x[[i]]$nls_pars_est[pars], col=cols[['nls_fit']], xaxt="n",...)
if(i==1) {
legend_text <- c(legend_text,'nls_est')
legend_cols <- c(legend_cols,cols[['nls_fit']])
legend_lty <- c(legend_lty,0)
legend_pch <- c(legend_pch,1)
}
}
}
}
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)
}
}
plot_list_simode_fit <-
function(x, show=c('nls','im','both'),
which=NULL, pars_true=NULL, time=NULL,
plot_mean_sd=FALSE, plot_im_smooth=FALSE, legend=FALSE, mfrow=par('mfrow'),
cols=list(nls_fit="blue",im_fit="green", true="black",
obs="red", im_smooth="magenta"),
fit_num=NULL, main=NULL,...) {
show <- match.arg(show)
plot_nls_fit <- (show=='nls' || show=='both')
plot_im_fit <- (show=='im' || show=='both')
x0.na <- names(x[[1]]$x0[is.na(x[[1]]$x0)])
xvars <- setdiff(names(x[[1]]$obs),names(x[[1]]$equations))
vars <- which
if(is.null(vars))
vars <- setdiff(names(x[[1]]$obs),xvars)
if(!is.null(main)) {
if(length(main)==1) {
main <- rep(main,length(vars))
}
else {
stopifnot(length(main)==length(vars))
}
names(main) <- vars
}
if(length(x)==1)
plot_mean_sd <- F
if(!is.null(pars_true)) {
stopifnot(is.numeric(pars_true))
stopifnot(all(setdiff(x[[1]]$pars,x[[1]]$likelihood_pars) %in% names(pars_true)))
}
if(!is.list(x[[1]]$time)) {
if(is.null(time))
time <- x[[1]]$time
time_obs <- rep(list(x[[1]]$time),length(x[[1]]$obs[vars]))
}
else {
if(is.null(time))
time <- sort(unique(unlist(x[[1]]$time)))
time_obs <- x[[1]]$time[which(names(x[[1]]$obs) %in% vars)]
}
names(time_obs) <- vars
y_min <- unlist(lapply(vars, function(i) min(x[[1]]$obs[[i]])))
y_max <- unlist(lapply(vars, function(i) max(x[[1]]$obs[[i]])))
names(y_min) <- names(x[[1]]$obs[vars])
names(y_max) <- names(x[[1]]$obs[vars])
legend_text <- c('obs')
legend_cols <- cols[['obs']]
if(!is.null(pars_true)) {
x0_true <- x[[1]]$x0
x0_true[x0.na] <- pars_true[x0.na]
pars_true <- pars_true[setdiff(names(pars_true),x0.na)]
model_true <- solve_ode(x[[1]]$equations, pars_true, x0_true, time, x[[1]]$obs[xvars])
if(is.null(model_true)) {
stop('Error running model with pars_true values')
}
legend_text <- c(legend_text, 'true')
legend_cols <- c(legend_cols, cols[['true']])
y_min1 <- apply(as.matrix(model_true[,vars]), 2, function(col) min(col))
y_max1 <- apply(as.matrix(model_true[,vars]), 2, function(col) max(col))
y_min <- pmin(y_min,y_min1)
y_max <- pmax(y_max,y_max1)
}
model_nls_est <- list()
if(plot_nls_fit) {
for(i in 1:length(x)) {
simode_obj <- x[[i]]
simode_obj$time <- time
model_nls_est[[i]] <- solve_ode2(simode_obj,'nls')$nls
if(is.null(model_nls_est[[i]])) {
fn <- ifelse(!is.null(fit_num),fit_num,i)
stop(paste0('Error running model with nls estimates of fit #',fn))
}
}
model_nls_est_mean <- matrix(0,nrow=length(time), ncol=length(vars))
model_nls_est_sd <- matrix(0,nrow=length(time), ncol=length(vars))
colnames(model_nls_est_mean) <- vars
colnames(model_nls_est_sd) <- vars
for(var in vars) {
nls_est_mat <- matrix(0,nrow=length(time), ncol=length(x))
for(i in 1:length(x)) {
nls_est_mat[,i] <- model_nls_est[[i]][,var]
}
model_nls_est_mean[,var] <- apply(nls_est_mat,1,mean)
if(length(x)>1)
model_nls_est_sd[,var] <- apply(nls_est_mat,1,stats::sd)
}
text <- ifelse(plot_mean_sd,'nls_mean_sd','nls_fit')
legend_text <- c(legend_text, text)
legend_cols <- c(legend_cols, cols[['nls_fit']])
y_min1 <- apply(as.matrix(model_nls_est_mean[,vars]-model_nls_est_sd[,var]), 2, function(col) min(col))
y_max1 <- apply(as.matrix(model_nls_est_mean[,vars]+model_nls_est_sd[,var]), 2, function(col) max(col))
y_min <- pmin(y_min,y_min1)
y_max <- pmax(y_max,y_max1)
}
model_im_est <- list()
if(plot_im_fit) {
for(i in 1:length(x)) {
simode_obj <- x[[i]]
simode_obj$time <- time
model_im_est[[i]] <- solve_ode2(simode_obj,'im')$im
if(is.null(model_im_est[[i]])) {
fn <- ifelse(!is.null(fit_num),fit_num,i)
stop(paste0('Error running model with im estimates of fit #',fn))
}
}
model_im_est_mean <- matrix(0,nrow=length(time), ncol=length(vars))
model_im_est_sd <- matrix(0,nrow=length(time), ncol=length(vars))
colnames(model_im_est_mean) <- vars
colnames(model_im_est_sd) <- vars
for(var in vars) {
im_est_mat <- matrix(0,nrow=length(time), ncol=length(x))
for(i in 1:length(x)) {
im_est_mat[,i] <- model_im_est[[i]][,var]
}
model_im_est_mean[,var] <- apply(im_est_mat,1,mean)
if(length(x)>1)
model_im_est_sd[,var] <- apply(im_est_mat,1,stats::sd)
}
text <- ifelse(plot_mean_sd,'im_mean_sd','im_fit')
legend_text <- c(legend_text, text)
legend_cols <- c(legend_cols, cols[['im_fit']])
y_min1 <- apply(model_im_est_mean[,vars]-model_im_est_sd[,var], 2, function(col) min(col))
y_max1 <- apply(model_im_est_mean[,vars]+model_im_est_sd[,var], 2, function(col) max(col))
y_min <- pmin(y_min,y_min1)
y_max <- pmax(y_max,y_max1)
}
if(plot_im_smooth) {
im_time <- x[[1]]$im_smooth$time
model_im_smooth_mean <- matrix(0,nrow=length(im_time), ncol=length(vars))
model_im_smooth_sd <- matrix(0,nrow=length(im_time), ncol=length(vars))
colnames(model_im_smooth_mean) <- vars
colnames(model_im_smooth_sd) <- vars
im_smooth_mat <- matrix(0,nrow=length(im_time), ncol=length(x))
for(var in vars) {
for(i in 1:length(x)) {
im_smooth_mat[,i] <- x[[i]]$im_smooth$val[,var]
}
model_im_smooth_mean[,var] <- apply(im_smooth_mat,1,mean)
if(length(x)>1)
model_im_smooth_sd[,var] <- apply(im_smooth_mat,1,stats::sd)
}
text <- ifelse(plot_mean_sd,'im_smooth_mean_sd','im_smooth')
legend_text <- c(legend_text, text)
legend_cols <- c(legend_cols, cols[['im_smooth']])
}
if(legend)
old.par <- par(mfrow=mfrow, mar=c(4,4,2,6), xpd=TRUE)
else
old.par <- par(mfrow=mfrow, mar=c(4,4,2,6))
on.exit(par(old.par))
for(var in vars) {
ylab <- ifelse(is.null(fit_num),var,paste0(var,' (fit ',fit_num,')'))
time_var <- time_obs[[var]]
plot(-100,-100,xlab='time',ylab=ylab,
xlim=c(time_var[1],time_var[length(time_var)]),
ylim=c(y_min[var],y_max[var]), main=main[var], ...)
for(i in 1:length(x)) {
points(time_var,x[[i]]$obs[[var]],col=cols[['obs']],...)
}
if(!is.null(pars_true)) {
lines(time,model_true[,var],col=cols[['true']],...)
}
if(plot_mean_sd) {
if(plot_im_smooth) {
lines(im_time,model_im_smooth_mean[,var],col=cols[['im_smooth']],...)
errorbar(x=im_time, model_im_smooth_mean[,var], yerr=model_im_smooth_sd[,var],
bar.col=cols[['im_smooth']], grid=FALSE, add=TRUE, xaxt="n",...)
}
if(plot_im_fit) {
lines(time,model_im_est_mean[,var],col=cols[['im_fit']],...)
errorbar(x=time, model_im_est_mean[,var], yerr=model_im_est_sd[,var],
bar.col=cols[['im_fit']], grid=FALSE, add=TRUE, xaxt="n")
}
if(plot_nls_fit) {
lines(time,model_nls_est_mean[,var],col=cols[['nls_fit']],...)
errorbar(x=time, model_nls_est_mean[,var], yerr=model_nls_est_sd[,var],
bar.col=cols[['nls_fit']], grid=FALSE, add=TRUE, xaxt="n")
}
}
else {
for(i in 1:length(x)) {
if(plot_im_smooth)
lines(im_time,x[[i]]$im_smooth$val[,var],col=cols[['im_smooth']],...)
if(plot_im_fit)
lines(time,model_im_est[[i]][,var],col=cols[['im_fit']],...)
if(plot_nls_fit)
lines(time,model_nls_est[[i]][,var],col=cols[['nls_fit']],...)
}
}
if(legend) {
legend(1.05*time[length(time)],
y_max[[var]],
cex=0.8,
bty="n",
legend_text ,col=legend_cols,
lty=c(0,rep(1,length(legend_text)-1)),
pch=c(1,rep(NA,length(legend_text)-1)))
}
}
}
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.