aboutsummaryrefslogtreecommitdiff
path: root/R/mkinplot.R
diff options
context:
space:
mode:
authorjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-18 12:58:38 +0000
committerjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2010-05-18 12:58:38 +0000
commit16f5b1d3c0136413e92b2be0f20d365e92e9cd1c (patch)
treeedab7b27ba4253b9dd45520e504c54851f65f26f /R/mkinplot.R
parent30cbb4092f6d2d3beff5800603374a0d009ad770 (diff)
Much more complete version that was just submitted to CRAN.
git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@9 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
Diffstat (limited to 'R/mkinplot.R')
-rw-r--r--R/mkinplot.R49
1 files changed, 49 insertions, 0 deletions
diff --git a/R/mkinplot.R b/R/mkinplot.R
new file mode 100644
index 0000000..d2880c0
--- /dev/null
+++ b/R/mkinplot.R
@@ -0,0 +1,49 @@
+mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), ...)
+{
+ fixed <- fit$fixed$value
+ names(fixed) <- rownames(fit$fixed)
+ parms.all <- c(fit$par, fixed)
+ ininames <- c(
+ rownames(subset(fit$start, type == "state")),
+ rownames(subset(fit$fixed, type == "state")))
+ odeini <- parms.all[ininames]
+ names(odeini) <- names(fit$diffs)
+
+ outtimes <- seq(xlim[1], xlim[2], length.out=100)
+
+ odenames <- c(
+ rownames(subset(fit$start, type == "deparm")),
+ rownames(subset(fit$fixed, type == "deparm")))
+ odeparms <- parms.all[odenames]
+
+ # Solve the ode
+ out <- ode(
+ y = odeini,
+ times = outtimes,
+ func = fit$mkindiff,
+ parms = odeparms)
+
+ # Output transformation for models with unobserved compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(fit$map)) {
+ if(length(fit$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, fit$map[[var]]])
+ }
+ }
+
+ # Plot the data and model output
+ plot(0, type="n",
+ xlim = xlim, ylim = ylim,
+ xlab = xlab, ylab = ylab, ...)
+ col_obs <- pch_obs <- 1:length(fit$map)
+ names(col_obs) <- names(pch_obs) <- names(fit$map)
+ for (obs_var in names(fit$map)) {
+ points(subset(fit$data, variable == obs_var, c(time, observed)),
+ pch = pch_obs[obs_var], col = col_obs[obs_var])
+ }
+ matlines(out_transformed$time, out_transformed[-1])
+ legend("topright", inset=c(0.05, 0.05), legend=names(fit$map),
+ col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
+}

Contact - Imprint