R/hill.R

In oro.pet: Rigorous - Positron Emission Tomography

##
## Copyright (c) 2012, Brandon Whitcher
## All rights reserved.
## 
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:
## 
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer. 
##     * Redistributions in binary form must reproduce the above
##       copyright notice, this list of conditions and the following
##       disclaimer in the documentation and/or other materials provided
##       with the distribution.
##     * Neither the name of Rigorous Analytics Ltd. nor the names of
##       its contributors may be used to endorse or promote products 
##       derived from this software without specific prior written 
##       permission.
## 
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##
## $Id: $
##

hillEquation <- function(conc, occ, guess=c(1,100), control=minpack.lm::nls.lm.control()) {
    func <- function(x, conc, occ) {
    IC50 <- x[1]
    rmax <- x[2]
    occ - rmax * conc / (conc + IC50)
  }
  out <- minpack.lm::nls.lm(par=guess, fn=func, control=control, conc=conc, occ=occ)
  rdf <- length(out$fvec) - length(coef(out))
  varcovmat <- (out$deviance / rdf) * chol2inv(chol(out$hessian))
  list(IC50=out$par[1],
       rmax=out$par[2],
       IC50SE=sqrt(varcovmat[1,1]),
       rmaxSE=sqrt(varcovmat[2,2]),
       hessian=out$hessian, info=out$info, deviance=out$deviance,
       message=out$message)
}