R/fitCellGrowth.R

In cellGrowth: Fitting cell population growth models

#' Fit a cell growth curve
#' 
#' For the non-parametric "locfit" model, local regression is done by a call to
#' \code{\link[locfit:locfit]{locfit}}. The returned maximum growth rate values
#' the maximum value of the fitted derivative over the data points.
#' For the parametric models "logistic", "gompertz", "rosso" and "baranyi", the function
#' does a non-least square fit by calling \code{\link{nls}}. Initial parameters values are
#' generated by \code{\link{guessCellGrowthParams}}. The returned maximum growth
#' rate values the \code{mu} parameter of these models.
#' 
#' @title Fit growth curves
#' @param x \code{numeric} vector: time points  
#' @param z \code{numeric} vector: log(OD)
#' @param model which model to fit.
#' @param locfit.h \code{numeric}: \code{h} parameter (window size) in call to
#' \code{\link[locfit:locfit]{locfit}}. The default value is set to three hours assuming
#' \code{x} given in seconds. You can detect a better bandwidth by calling \code{\link{bandwidthCV}}
#' @param locfit.deg \code{numeric}: \code{deg} parameter (polynomials degree) in call to
#' \code{\link[locfit:locfit]{locfit}}
#' @param relative.height.at.lag Parameter used by \code{\link{guessCellGrowthParams}}
#' @return Fit as returned by \code{\link[locfit:locfit]{locfit}} for the
#' "locfit" model and as returned by \code{\link{nls}} for the "logistic", "gompertz", "rosso"
#' and "baranyi" models. The returned value also has an attribute \code{maxGrowthRate} valueing
#' the inferred maximum growth rate as well as an attribute \code{pointOfMaxGrowthRate}
#' valuing the datapoint at which the growth rate is maximal. Also, it has an attribute
#' \code{max} valuing the inferred maximum among the time points as well as \code{pointOfMax}
#' valuing the datapoint of max fitted value. It gets the additional class \code{cellCurveFit}
#' assigned.
#' @author Julien Gagneur and Moritz Matthey
#' @seealso \code{\link{nls}}, \code{\link[locfit:locfit]{locfit}}, \code{\link{baranyi}}, \code{\link{gompertz}}, \code{\link{logistic}}, \code{\link{rosso}}, \code{\link{guessCellGrowthParams}}, \code{\link{fitCellGrowths}}
#' @export
#' @examples x = 1:1000
#'           z = gompertz(x, mu=0.01, l=200, z0=1, zmax=5) + rnorm(length(x),sd=0.1)
#'           f = fitCellGrowth(x, z, model = "gompertz")
#'           floc = fitCellGrowth(x, z, model = "locfit", locfit.h=500)
#'           	plot(x,z, main="simulated data\nGompertz model")
#'           	lines(x, predict(f,x), lwd=2, col="red")
#'           	lines(x, predict(floc,x), lwd=2, col="blue")
#'           	legend( "right", legend=c("gompertz fit", "locfit"), lwd=1, col=c("red","blue") )
fitCellGrowth = function(x,
		z,						
        model = "locfit",						
        locfit.h=3*60*60,
		locfit.deg=2,
		relative.height.at.lag=0.1){													
    															
    if(model!="locfit"){																	
        # first guess values
        guess =  guessCellGrowthParams(x,z,relative.height.at.lag=relative.height.at.lag)															
        
        # make formula
        param.str = do.call(paste, c(as.list(names(guess)), sep="," )) 						
        expr = paste("z ~", model, "(x,", param.str, ")" )									
        
        # fit with least square
		error = 0
        rv = tryCatch(nls(as.formula(expr), start = guess),error=function(e){
					print("Fitting did not succeed. Try with other parameters for relative.height.at.lag ( See guessCellGrowthParams )")
					NA
				});
		if ( is.null(attr(rv,"class")) )
			return(NA)
		
		param = summary( rv )$parameters
		mu = param[ rownames(param) == "mu", "Estimate" ]
		l =  param[ rownames(param) == "l", "Estimate" ]
		z0 =  param[ rownames(param) == "z0", "Estimate" ]
		zmax  =  param[ rownames(param) == "zmax", "Estimate" ]
		
		attr(rv, "maxGrowthRate") = mu
		
		if ( model == 'gompertz')
			attr(rv,"pointOfMaxGrowthRate") = (exp(1)*l*mu-z0+zmax)/(exp(1)*mu);
		if ( model == 'logistic')
			attr(rv,"pointOfMaxGrowthRate") = (2*l*mu-z0+zmax)/(2*mu);
		if (model == 'rosso' )
			attr(rv,"pointOfMaxGrowthRate") = l;
		
		## numeric derivative
		# xNumeric <- as.numeric(x)	# in case x has integer values
		# numDrv <- numericDeriv(quote(predict(rv,data.frame(x=xNumeric))),c("xNumeric")) # Numeric derivative
		# numDrv <- diag(attr(numDrv,"gradient"))

		## max fitted value
		attr(rv, "max") = max( predict(rv, x) )
		attr(rv, "pointOfMax") = which.max( predict(rv, x) )
 			
    }else{																					
        require(locfit)																		
        
        ## main fit
        rv = locfit( z ~ lp(x, h=locfit.h, deg=locfit.deg) )
        
		## numeric derivative
		# xNumeric <- as.numeric(x)	# in case x has integer values
		# numDrv <- numericDeriv(quote(predict(rv,xNumeric)),c("xNumeric")) # Numeric derivative
		# numDrv <- diag(attr(numDrv,"gradient"))
		#
		## max growth rate
		# attr(rv, "maxGrowthRate") = max( numDrv )
		# attr(rv, "pointOfMaxGrowthRate") = which.max( numDrv )
		#####
		
		## get first derivative for the *same* fit
		drv = locfit( z ~ lp(x, h=locfit.h, deg=locfit.deg), deriv=1 )
	
		## max growth rate as the max value of fitted derivative on the input points
		attr(rv, "maxGrowthRate") = max( predict(drv, x) )
		attr(rv, "pointOfMaxGrowthRate") = which.max( predict(drv, x) )
	
		## max fitted value
		attr(rv, "max") = max( predict(rv, x) )
		attr(rv, "pointOfMax") = which.max( predict(rv, x) )
		
		## max growth reached too soon
		if ( x[attr(rv,"pointOfMaxGrowthRate")]-min(x) < locfit.h && max(x) - x[attr(rv,"pointOfMaxGrowthRate")] < locfit.h)
			warning("Maximal growth reached close to the border. You can identify the well by testing for pointOfMaxGrowthRate.")
		if ( max(x)-x[attr(rv,"pointOfMaxGrowthRate")] < locfit.h )
			warning("Maximal population reached close to the border. You can identify the well by testing for pointOfMax")
    }

	## time points / values
	attr(rv, "x") = x;
	attr(rv, "z") = z;
	
	attr(rv, "model") = model;
	
    class(rv) = c("cellGrowthFit", class(rv) )
    return(rv)																				 
}