Wrapper around trust allowing for multiple fits
from randomly chosen initial values.
mstrust(objfun, center, studyname, rinit = 0.1, rmax = 10, fits = 20,
cores = 1, samplefun = "rnorm", resultPath = ".", stats = FALSE, ...)Objective function, see trust.
The names of the study or fit. This name is used to determine filenames for interim and final results. See Details.
Starting trust region radius, see trust.
Maximum allowed trust region radius, see trust.
Number of fits (jobs).
Number of cores for job parallelization.
Function to sample random initial values. It is assumed,
that samplefun has a named parameter "n" which defines how many
random numbers are to be returned, such as for rnorm or
runif. By default rnorm is used. Parameteres
for samplefun are simply appended as named parameters to the mstrust call
and automatically handed to samplefun by matching parameter names.
character indicating the folder where the results should be stored. Defaults to ".".
If true, the same summary statistic as written to the logfile is printed to command line on mstrust completion.
Additional parameters handed to trust(), samplefun(), or the objective function by matching parameter names. All unmatched parameters are handed to the objective function objfun(). The log file starts with a table telling which parameter was assigend to which function.
A parlist holding errored and converged fits.
By running multiple fits starting at randomly chosen inital
parameters, the chisquare landscape can be explored using a deterministic
optimizer. Here, trust is used for optimization. The standard
procedure to obtain random initial values is to sample random variables
from a uniform distribution (rnorm) and adding these to
center. It is, however, possible, to employ any other sampling
strategy by handing the respective function to mstrust(),
samplefun.
In case a special sampling is required, a customized sampling function can be used. If, e.g., inital values leading to a non-physical systems are to be discarded upfront, the objective function can be addapted accordingly.
All started fits either lead to an error or complete converged or unconverged. A statistics about the return status of fits can be shown by setting stats to TRUE.
Fit final and intermediat results are stored under studyname. For
each run of mstrust for the same study name, a folder under
studyname of the form "trial-x-date" is created. "x" is the number
of the trial, date is the current time stamp. In this folder, the
intermediate results are stored. These intermediate results can be loaded
by load.parlist. These are removed on successfull completion
of mstrust. In this case, the final list of fit parameters
(parameterList.Rda) and the fit log (mstrust.log) are found instead.
# NOT RUN {
################################################################
## Walk through the most frequently used functions in
## connection with ODE models
################################################################
library(deSolve)
## Model definition (text-based, scripting part)
r <- NULL
r <- addReaction(r, "A", "B", "k1*A", "translation")
r <- addReaction(r, "B", "", "k2*B", "degradation")
f <- as.eqnvec(r)
observables <- eqnvec(Bobs = "s1*B")
innerpars <- getSymbols(c(names(f), f, observables))
trafo0 <- structure(innerpars, names = innerpars)
trafo0 <- replaceSymbols(innerpars,
paste0("exp(", innerpars, ")"),
trafo0)
conditions <- c("a", "b")
trafo <- list()
trafo$a <- replaceSymbols("s1", "s1_a", trafo0)
trafo$b <- replaceSymbols("s1", "s1_b", trafo0)
events <- list()
events$a <- data.frame(var = "A", time = 5, value = 1, method = "add")
events$b <- data.frame(var = "A", time = 5, value = 2, method = "add")
## Model definition (generate compiled objects
## and functions from above information)
# ODE model
model <- odemodel(f)
# Observation function
g <- Y(observables, f, compile = TRUE, modelname = "obsfn")
# Parameter transformation for steady states
pSS <- P(f, method = "implicit", condition = NULL)
# Condition-specific transformation and prediction functions
p0 <- x <- NULL
for (C in conditions) {
p0 <- p0 + P(trafo[[C]], condition = C)
x <- x + Xs(model, events = events[[C]], condition = C)
}
## Process data
# Data
data <- datalist(
a = data.frame(time = c(0, 2, 7),
name = "Bobs",
value = c(.3, .3, .3),
sigma = c(.03, .03, .03)),
b = data.frame(time = c(0, 2, 7),
name = "Bobs",
value = c(.1, .1, .2),
sigma = c(.01, .01, .02))
)
## Evaluate model/data
# Initialize times and parameters
times <- seq(0, 10, .1)
outerpars <- attr(p0, "parameters")
pouter <- structure(rnorm(length(outerpars)),
names = outerpars)
# Plot prediction
plot((g*x*p0)(times, pouter))
plot((g*x*pSS*p0)(times, pouter))
plotFluxes(pouter, g*x*p0, times, getFluxes(r)$B)
# Fit data with L2 constraint
constr <- constraintL2(mu = 0*pouter, sigma = 5)
myfit <- trust(normL2(data, g*x*p0) + constr,
pouter, rinit = 1, rmax = 10)
plot((g*x*p0)(times, myfit$argument), data)
# List of fits
obj <- normL2(data, g*x*p0) + constr
mylist <- mstrust(normL2(data, g*x*p0) + constr,
pouter, fits = 10, cores = 4, sd = 4)
plot(mylist)
pars <- as.parframe(mylist)
plotValues(pars)
bestfit <- as.parvec(pars)
# Compute profiles of the fit
profile <- profile(normL2(data, g*x*p0) + constr,
bestfit, "k1", limits = c(-5, 5))
plotProfile(profile)
# Add a validation objective function
vali <- datapointL2("A", 2, "mypoint", .1, condition = "a")
# Get validation point parameters
mypoint <- trust(normL2(data, g*x*p0) + constr + vali,
c(mypoint = 0), rinit = 1, rmax = 10,
fixed = bestfit)$argument
# Compoute profile and plot
profile <- profile(normL2(data, g*x*p0) + constr + vali,
c(bestfit, mypoint), "mypoint")
plotProfile(profile)
## Using the controls() function to manipulate objects ------------------------
# List the available controls of an object
controls(x)
# List a specific control
controls(x, condition = "a", name = "optionsSens")
# Change specific control
controls(x, condition = "a", name = "optionsSens") <-
list(method = "lsoda", rtol = 1e-8, atol = 1e-8)
# Almost every function (objfn, parfn, prdfn, obsfn) saves the arguments
# by which it was invoked in a list named "controls", which can be manipulated
# by the controls() function.
## New example: condition-specific observation parameters ---------------------
# Condition-specific observation parameters
f <- eqnvec(A = "-k1*A", B = "k1*A - k2*B")
observables <- eqnvec(Bobs = "s1*B")
conditions <- c("scale1", "scale2")
dynpars <- getSymbols(c(names(f), f))
obspars <- setdiff(getSymbols(observables), dynpars)
trafo0 <- structure(c(dynpars, obspars), names = c(dynpars, obspars))
obspars_all <- outer(obspars, conditions, paste, sep = "_")
trafo_all <- structure(c(dynpars, obspars_all),
names = c(dynpars, obspars_all))
trafo_all <- replaceSymbols(dynpars,
paste0("exp(", dynpars, ")"),
trafo_all)
trafo_all <- replaceSymbols(obspars_all,
paste0("exp(", obspars_all, ")"),
trafo_all)
pObs <- NULL
for (C in conditions) {
pObs <- pObs + P(replaceSymbols(obspars,
paste(obspars, C, sep = "_"),
trafo0),
condition = C)
}
x <- Xs(model)
p <- P(trafo_all)
y <- g*pObs*x*p
times <- seq(0, 10, .1)
outerpars <- attr(p, "parameters")
pouter <- structure(rnorm(length(outerpars)), names = outerpars)
# Plot prediction
plot((g*pObs*x*p)(times, pouter))
# }
Run the code above in your browser using DataLab