R/kinresults.R

In kinfit: Routines for Fitting Kinetic Models to Chemical Degradation Data

# $Id: kinresults.R 124 2011-11-09 10:58:27Z jranke $

# Copyright (C) 2008-2011 Johannes Ranke
# Contact: mkin-devel@lists.berlios.de

# This file is part of the R package kinfit

# kinfit is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.

# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.

# You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>

kinresults <- function(kinfits, alpha = 0.05, DTmax = 1000, SFORB=TRUE)
{
	kindata <- data.frame(
    		t = kinfits[[1]]$model$t, 
    		parent = kinfits[[1]]$model$parent)
 	kindata.means <- aggregate(kindata, list(kindata$t), mean)
	kindata.means.mean <- mean(kindata.means$parent, na.rm=TRUE)
	n.times <- length(kindata.means$parent)
	parms <- list()
	df <- err.min <- R2 <- RSS <- TSS <- RSS.means <- TSS.means <- vector()
	DT50 <- DT90 <- vector()
	f <- list()
	
	for (kinmodel in names(kinfits))
	{
		m = kinfits[[kinmodel]]
		if(class(m) == "nls") {
			kindata.means$est <- predict(m, kindata.means)
			
      		if (!"parent.0" %in% names(coef(m))) {
				
       			parms[[kinmodel]] <- switch(kinmodel,
          			SFO = list(k = coef(m)[["k"]]),
          			FOMC = list(alpha = coef(m)[["alpha"]],
              			beta = coef(m)[["beta"]]),
	    			HS = list(k1 = coef(m)[["k1"]],
              			k2 = coef(m)[["k2"]], 
              			tb = coef(m)[["tb"]]),
          			DFOP= list(k1 = coef(m)[["k1"]],
              			k2 = coef(m)[["k2"]], 
              			g = coef(m)[["g"]]))
      		} else {
        			parms[[kinmodel]] <- switch(kinmodel,
          				SFO = list(parent.0 = coef(m)[["parent.0"]], 
              				k = coef(m)[["k"]]),
          				FOMC = list(parent.0 = coef(m)[["parent.0"]],
              				alpha = coef(m)[["alpha"]],
             				 beta = coef(m)[["beta"]]),
          				HS = list(parent.0 = coef(m)[["parent.0"]], 
              				k1 = coef(m)[["k1"]],
             			 	k2 = coef(m)[["k2"]], 
              				tb = coef(m)[["tb"]]),
          				DFOP = list(parent.0 = coef(m)[["parent.0"]],
              				k1 = coef(m)[["k1"]],
              				k2 = coef(m)[["k2"]], 
              				g = coef(m)[["g"]]))
        			if(kinmodel == "DFOP" & SFORB) {
          				k1 = coef(m)[["k1"]]
          				k2 = coef(m)[["k2"]]
          				g = coef(m)[["g"]]
          				parms[["SFORB"]] = 
            				list(parent.0 = coef(m)[["parent.0"]],
              				k1out = g * k1 + (1 - g) * k2,
             			 	k21 = k1 * k2 / (g * k1 + (1 - g) * k2),
              				k12 = (g * (1 - g) * (k1 - k2)^2) / (g * k1 + (1 - g) * k2))
        			}
      		}

  			n.parms = length(coef(m))
			if (!"parent.0" %in% names(coef(m))) {

				f[[kinmodel]] = switch(kinmodel,
					HS = function(t, x) {
						(HS(t, kinfits[[kinmodel]]$model[["parent.0.user"]], 
							coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["tb"]]) - 
						(1 - x/100) * kinfits[[kinmodel]]$model[["parent.0.user"]])^2
					},
					DFOP = function(t, x) {
						(DFOP(t, kinfits[[kinmodel]]$model[["parent.0.user"]], 
							coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["g"]]) - 
						(1 - x/100) * kinfits[[kinmodel]]$model[["parent.0.user"]])^2
					}
				)
			}else{
				f[[kinmodel]] = switch(kinmodel,
					HS = function(t, x) {
						(HS(t, coef(m)[["parent.0"]], 
							coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["tb"]]) - 
						(1 - x/100) * coef(m)[["parent.0"]])^2
					},
					DFOP = function(t, x) {
						(DFOP(t, coef(m)[["parent.0"]], 
							coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["g"]]) - 
						(1 - x/100) * coef(m)[["parent.0"]])^2
					}
				)
			}
				
			coef(m)

			df[[kinmodel]] = n.times - n.parms
			RSS[[kinmodel]] = sum(summary(m)$residuals^2)
      TSS[[kinmodel]] = sum((m$model$parent - mean(m$model$parent))^2)
	
			DT50.o = switch(kinmodel,
        SFO = log(2)/coef(m)[["k"]],
				FOMC = coef(m)[["beta"]] * (2^(1/coef(m)[["alpha"]]) - 1),
				HS = optimize(f[[kinmodel]], c(0, DTmax), x=50)$minimum,
				DFOP = optimize(f[[kinmodel]], c(0, DTmax), x=50)$minimum)
			
			
      DT50[[kinmodel]] = ifelse(abs(DT50.o - DTmax) < 0.1, NA, DT50.o)
      DT90.o = switch(kinmodel,
          SFO = log(10)/coef(m)[["k"]],
          FOMC = coef(m)[["beta"]] * (10^(1/coef(m)[["alpha"]]) - 1),
          HS = optimize(f[[kinmodel]], c(0, DTmax), x=90)$minimum,
          DFOP = optimize(f[[kinmodel]], c(0, DTmax), x=90)$minimum)
      DT90[[kinmodel]] = ifelse(abs(DT90.o - DTmax) < 0.1, NA, DT90.o)

      # Chi2 error level as defined in FOCUS kinetics (see ?kinerrmin)
      err.min[[kinmodel]] <- kinerrmin(kinfits, kinmodel)

      # Coefficient of determination calculated from residual sum of squares and totals sum of squares
      # so this r2 is what is called model efficiency in FOCUS kinetics (2006), p. 99
      R2[[kinmodel]] = 1 - RSS[[kinmodel]]/TSS[[kinmodel]]

		}
	}

	stats <- data.frame(n.times = n.times, df = df, 
   	mean.means = kindata.means.mean, 
		RSS = RSS, err.min = err.min, R2 = R2)
	results <- data.frame(DT50 = DT50, DT90 = DT90)
	list(parms = parms, stats = stats, results = results)
}