Nothing
###
### $Id: sc3-2-mainFit.R 1000 2017-11-14 14:11:30Z krcoombes $
###
##
## Module Variables
##
## :TODO: Migrate this to .onLoad since it's singleton-like
.ModelEnv <- new.env(hash=TRUE)
attr(.ModelEnv, "name") <- "SuperCurveModels"
##
## Private Methods
##
##-----------------------------------------------------------------------------
## Returns private environment for storing registered models
modelenv <- function() {
return(.ModelEnv)
}
##-----------------------------------------------------------------------------
registerPkgFitModels <- function() {
registerModel("logistic", "LogisticFitClass", "Logistic")
registerModel("cobs", "CobsFitClass", "Monotone Increasing B-spline")
registerModel("loess", "LoessFitClass", "Loess")
}
##
## Public Methods
##
##-----------------------------------------------------------------------------
## REQUIRED INPUTS:
## rppa - an RPPA object
## design - a RPPADesign object
## measure - name of column in rppa@data to view as intensity
##
##-----------------------------------------------------------------------------
## All optional inputs are handled through RPPAFitParams
##
## OPTIONAL INPUTS CONTROLLING THE ALGORITHM:
## model - which statistical model to fit
## method - how to fit the alpha and beta values.
## xform - function used to transform the values before fitting
## trim - if 0, concentrations will not be trimmed / bounded; otherwise,
## trim the concentrations using this value as the trim level.
## ci - if true, calculate a 90% confidence interval for the LC50
## supercurve values
## ignoreNegative - if true, then values below zero are converted to NA before
## estimating alpha and beta. Only matters for the 'quantiles' method.
##
## OPTIONAL INPUTS CONTROLLING VERBOSITY:
## trace - passed to nls in the 'bayesian' portion
## verbose - if true, have the function tell you what it is doing
## veryVerbose - if true, have the function overwhelm you telling you what it
## is doing
## warnLevel - used to set the 'warn' option before calling 'rlm'.
## since this is wrapped in a 'try', it won't cause failure but will
## give us a chance to figure out which dilution series failed.
## Setting warnLevel to two or greater may change computed results.
##-----------------------------------------------------------------------------
## Refactored this to improve maintainability; everything now happens once.
## Retain this generator for backwards compatibility.
RPPAFit <- function(rppa,
design,
measure,
model="logistic",
xform=NULL,
method=c("nls", "nlrob", "nlrq"),
trim=2,
ci=FALSE,
ignoreNegative=TRUE,
trace=FALSE,
verbose=FALSE,
veryVerbose=FALSE,
warnLevel=0) {
## Check arguments
## [1] 'rppa' and 'design' are checked in the call to RPPAFitFromParams
## [2] Everything else is checked in RPPAFitParams
params <- RPPAFitParams(measure,
model,
xform,
method,
trim,
ci,
ignoreNegative,
trace,
verbose,
veryVerbose,
warnLevel)
RPPAFitFromParams(rppa, design, params)
}
##-----------------------------------------------------------------------------
## Collect parameters in one place and make sure they are reasonably sensible
RPPAFitParams <- function(measure,
model="logistic",
xform=NULL,
method=c("nls", "nlrob", "nlrq"),
trim=2,
ci=FALSE,
ignoreNegative=TRUE,
trace=FALSE,
verbose=FALSE,
veryVerbose=FALSE,
warnLevel=0) {
## Check arguments
## Start with the critical parameter, 'model', since it tells us which
## statistical model we are supposed to fit. At the time the curve is
## fit, 'model' has to be the name of a registered FitClass. You could
## conceivably define your own class and register it after setting up the
## fit parameters. Because of that fact, we defer some checking until
## later (in RPPAFitFromParams).
if (!is.character(model)) {
stop(sprintf("argument %s must be a character string",
sQuote("model")))
}
model <- model[1]
## This parameter is also critical, and we cannot complete checking it
## until later, when we know whether it refers to an actual column in
## the RPPA object being fit.
if (!is.character(measure)) {
stop(sprintf("argument %s must be character",
sQuote("measure")))
} else if (!(length(measure) == 1)) {
stop(sprintf("argument %s must be of length 1",
sQuote("measure")))
}
if (!is.null(xform)) {
if (!is.function(xform)) {
stop(sprintf("argument %s must be function, if specified",
sQuote("xform")))
}
}
## Remaining parameters are less important, so we mostly just transform
## them into the correct type and length.
method <- match.arg(method)
ci <- as.logical(ci)[1]
ignoreNegative <- as.logical(ignoreNegative)[1]
trace <- as.logical(trace)[1]
verbose <- as.logical(verbose)[1]
veryVerbose <- as.logical(veryVerbose)[1]
warnLevel <- as.integer(warnLevel)[1]
if (!(is.numeric(trim) || is.logical(trim))) {
stop(sprintf("argument %s must be numeric (or logical)",
sQuote("trim")))
}
trim <- if (is.logical(trim) && isTRUE(trim[1])) {
## Substitute numeric equivalent for original behavior
formals(RPPAFitParams)$trim
} else {
as.numeric(trim)[1]
}
if (is.na(trim) || trim < 0) {
stop(sprintf("argument %s must be a non-negative quantity",
sQuote("trim")))
}
## Create new class
new("RPPAFitParams",
measure=measure,
xform=xform,
method=method,
ci=ci,
ignoreNegative=ignoreNegative,
trace=trace,
verbose=verbose,
veryVerbose=veryVerbose,
warnLevel=warnLevel,
trim=trim,
model=model)
}
##-----------------------------------------------------------------------------
## Returns a string representation of this instance. The content and format of
## the returned string may vary between versions. Returned string may be
## empty, but never null.
setMethod("paramString", signature(object="RPPAFitParams"),
function(object,
slots=slotNames(object),
...) {
## :TODO: Implementation currently ignores the 'slots' argument and
## returns string containing parameters from various slots as though:
##
## slotsToDisplay <- c("measure", "model", "method", "trim",
## "ci", "ignoreNegative", "warnLevel")
## paramString(fp, slotsToDisplay)
##
paste(paste("measure:", shQuote(object@measure)), "\n",
paste("model:", shQuote(object@model)), "\n",
paste("method:", shQuote(object@method)), "\n",
paste("trim:", object@trim), "\n",
paste("ci:", object@ci), "\n",
paste("ignoreNegative:", object@ignoreNegative), "\n",
paste("warnLevel:", object@warnLevel), "\n",
sep="")
})
##-----------------------------------------------------------------------------
## Compute a multiple of the logit transform.
## in practice, epsilon never changes.
.calcLogitz <- function(data,
alpha,
beta,
gamma,
epsilon=0.001) {
z <- (data - alpha) / beta
## Truncate points at min and max for dilution curve
z[z < epsilon] <- epsilon
z[z > (1-epsilon)] <- 1-epsilon
log(z / (1-z)) / gamma
}
##-----------------------------------------------------------------------------
## Computes crude estimates of the parameters (i.e., alpha, beta, gamma, and
## the EC50s) so routine can be initialized.
.firstPass <- function(yval,
design,
ignoreNegative,
epsilon=1e-4) {
## 'yval' is a vector of intensity values selected from the MicroVigene
## file; its length of 'yval' should match the number of rows in
## design@layout, and be in the same order.
## in practice, 'epsilon' never changes
## Check arguments
stopifnot(is.numeric(yval))
stopifnot(is.RPPADesign(design))
stopifnot(is.logical(ignoreNegative) && length(ignoreNegative) == 1)
stopifnot(is.numeric(epsilon))
## Begin processing
if (ignoreNegative) {
yval[yval < 0] <- NA
}
temp <- yval[!.controlVector(design)]
## Compute a robust estimate for alpha and beta
## Tried to use robust 'rnls' from 'sfsmisc' to fit alpha and beta below
## but found it was still too sensitive to outliers.
## This is a crude way to get alpha and beta but it seems to be robust.
## HERE IS THE PROBLEM!
## Quantiles forces a double trim in conjunction with calcLogitz
## and epsilon above, artificially shrinking alpha and beta
## Try to match version 0.12
# lBot <- quantile(temp, probs=0.05, na.rm=TRUE)
# lTop <- quantile(temp, probs=0.95, na.rm=TRUE)
lBot <- min(temp, na.rm=TRUE)
lTop <- max(temp, na.rm=TRUE)
lindata <- .calcLogitz(yval, lBot, lTop-lBot, 1)
series <- seriesNames(design) # omits the control spots automagically
##-------------------------------------------------------------------------
.estimateLogisticSlope <- function(ser,
dsn,
ld) {
## Get the items in this dilution series
items <- dsn@layout$Series == ser
## Get the log2 steps for this series
steps <- dsn@layout$Steps[items]
## Get the transformed intensity values
x <- ld[items]
## Only use valid values
y <- x[!is.na(x) & !is.infinite(x)]
if (all(is.na(y))) {
NA
} else {
## ychange / xchange
(max(y) - min(y)) / (max(steps) - min(steps))
}
}
## Initial estimate of the logistic slope across the dilution steps.
## There should be a more comprehensible way to write this:
dat <- unlist(lapply(series,
.estimateLogisticSlope,
design,
lindata))
## Force a nonzero slope for starting estimate
minslope <- 0.001
## Filter out zero slopes from blank sample
dat <- dat[dat > minslope]
## Use mean, not median. We may have many blank samples on a slide,
## in which case median slope will be 0
gamma <- mean(dat, trim=0.01)
gamma <- max(gamma, minslope) # force nonzero slope for starting estimate
## For each sample, estimate the offset
passer <- rep(NA, length(series))
names(passer) <- series
for (this in series) {
layout <- design@layout
items <- layout$Series == this
passer[this] <- median(lindata[items] / gamma - layout$Steps[items],
na.rm=TRUE)
}
list(lBot=lBot,
lTop=lTop,
gamma=gamma,
passer=passer)
}
##-----------------------------------------------------------------------------
## Returns model associated with key for invocation.
getRegisteredModel <- function(key) {
return(classname <- getRegisteredObject(key,
envir=modelenv(),
"classname")$classname)
}
##-----------------------------------------------------------------------------
## Returns label associated with key for display by user interface.
getRegisteredModelLabel <- function(key) {
return(ui.label <- getRegisteredObject(key,
envir=modelenv(),
"classname")$ui.label)
}
##-----------------------------------------------------------------------------
## Returns vector containing "keys" for all registered models.
getRegisteredModelKeys <- function() {
keys <- getRegisteredObjectKeys(envir=modelenv())
if (length(keys) == 0) {
stop("no registered models exist")
}
return(keys)
}
##-----------------------------------------------------------------------------
## Registers specific model for use by SuperCurve package.
## Register names of what code will consider "valid" models.
registerModel <- function(key,
classname,
ui.label=names(key)) {
if (is.null(ui.label)) {
ui.label <- key
}
ui.label <- as.character(ui.label)[1]
## Verify class
tryCatch(new(classname),
error=function(e) {
stop(sprintf("cannot create instance of classname %s",
sQuote(classname)))
})
if (!extends(classname, "FitClass")) {
stop(sprintf("argument %s must be name of subclass of class %s",
sQuote("classname"),
sQuote("FitClass")))
}
registerClassname(key, classname, ui.label=ui.label, envir=modelenv())
}
##-----------------------------------------------------------------------------
RPPAFitFromParams <- function(rppa,
design,
fitparams,
progmethod=NULL) {
## Check arguments
if (!is.RPPA(rppa)) {
stop(sprintf("argument %s must be object of class %s",
sQuote("rppa"), "RPPA"))
}
if (!is.RPPADesign(design)) {
stop(sprintf("argument %s must be object of class %s",
sQuote("design"), "RPPADesign"))
}
if (!is.RPPAFitParams(fitparams)) {
stop(sprintf("argument %s must be object of class %s",
sQuote("fitparams"), "RPPAFitParams"))
}
if (!is.null(progmethod)) {
if (!is.function(progmethod)) {
stop(sprintf("argument %s must be function, if specified",
sQuote("progmethod")))
}
} else {
## Create a placeholder function
progmethod <- function(phase) {}
}
## :WORKAROUND: codetools (via R CMD check) complains unnecessarily about
## "no visible binding" as code below uses assign() rather than "<-".
measure <- model <- xform <- method <- trim <- warnLevel <-
ignoreNegative <- trace <- verbose <- veryVerbose <- ci <- NULL
## Create variables from 'fitparams' slots
for (slotname in slotNames(fitparams)) {
assign(slotname, slot(fitparams, slotname))
}
## Need to make certain that the 'model' is a registered FitClass
modelClass <- tryCatch(getRegisteredModel(model),
error=function(e) {
stop(sprintf("argument %s must be name of registered fit class",
sQuote("model")))
})
## Need to make sure that 'measure' refers to an actual data column
if (missing("measure")) {
stop("missing name of the measurement column to fit")
}
dn <- dimnames(rppa@data)[[2]]
temp <- pmatch(measure, dn)
if (is.na(temp)) {
stop("supply the name of a valid measurement column to fit")
} else if (length(temp) > 1) {
stop(sprintf("argument %s must identify unique column of argument %s",
sQuote("measure"),
sQuote("rppa")))
}
measure <- dn[temp]
## Begin processing
l1 <- levels(design@layout$Sample)
l2 <- levels(rppa@data$Sample)
if (!(length(l1) == length(l2)) || sum(l1 != l2) > 0) {
warning(sprintf("number of sample labels in design (%d) does not match number of sample labels given in raw RPPA quantification file (%d).",
length(l1),
length(l2)))
}
call <- match.call()
intensity <- if (!is.null(xform)) {
xform(rppa@data[, measure])
} else {
rppa@data[, measure]
}
silent <- warnLevel < 0
## Perform the first pass to initialize the estimates
progmethod("firstpass")
first <- .firstPass(intensity, design, ignoreNegative)
passer <- first$passer
gamma.est <- first$gamma
if (verbose) {
cat(paste("Completed first pass. Parameters:",
paste("\t", "lBot =", first$lBot),
paste("\t", "lTop =", first$lTop),
paste("\t", "G =", first$gamma),
"\n",
sep="\n"))
flush.console()
}
## Put our current guess at the x and y values into vectors
## Here we must finally explicitly poke into the internals of RPPADesign
## class in order to omit the control spots.
yval <- intensity[!.controlVector(design)]
## Create new class
fc <- switch(modelClass,
LogisticFitClass = new(modelClass,
coefficients=c(alpha=first[[1]],
beta=first[[2]],
gamma=first[[3]])),
new(modelClass)
)
## Do a two pass estimation, first using rough conc. estimates,
## then using better ones
steps <- getSteps(design)
series <- seriesNames(design)
for (pass in seq_len(2)) {
pass.name <- if (pass == 1) "coarse" else "fine"
xval <- if (pass == 1) {
steps + passer[names(design)]
} else {
steps + pass2[names(design)]
}
## Fit a response curve for the slide of the form yval = f(xval)
progmethod(sprintf("%s fit slide", pass.name))
fc <- fitSlide(fc,
conc=xval,
intensity=yval,
method=method)
## Conditional on the response curve fit for the slide
## Perform a separate fit of EC50 values for each dilution series.
pass2 <- rep(NA, length(series))
names(pass2) <- series
ss.ratio <- pass2
warn2 <- rep("", length(series))
names(warn2) <- series
series.len <- length(series)
report.incr <- as.integer(5)
i.this <- as.integer(1)
progmethod(sprintf("%s fit series", pass.name))
for (this in series) {
items <- names(design) == this
## Report in 5 percent increments rather than iterations (speed)
percent.this <- as.integer((i.this / series.len) * 100)
if (!(percent.this %% report.incr)) {
## Skip 0 and 100 when reporting
if (percent.this %% as.integer(100)) {
## Notify progress update
progmethod(sprintf("%s fit series (%d%%)",
pass.name,
percent.this))
}
}
fs <- fitSeries(fc,
diln=steps[items],
intensity=yval[items],
est.conc=passer[this],
method=method,
silent=silent,
trace=trace)
pass2[this] <- fs$est.conc
warn2[this] <- fs$warn
## Compute R^2 as sum(r[i]^2) / sum((y[i]-mean(y))^2),
## the fraction of variance explained for this series
resids <- fs$resids
sse <- sum(resids*resids)
sst <- var(yval[items]) * (length(yval[items])-1)
rsquared <- 1 - sse/sst
ss.ratio[this] <- rsquared
i.this <- i.this + as.integer(1)
} # end: this in series
progmethod(sprintf("%s fit series complete", pass.name))
if (verbose) {
cat("Finished estimating EC50 values. Coefficients:", "\n")
print(summary(pass2))
cat("SS Ratio:", "\n")
print(summary(ss.ratio))
flush.console()
}
}
## Create new class
result <- new("RPPAFit",
call=call,
rppa=rppa,
design=design,
measure=measure,
method=method,
trimset=c(lo.intensity=-100000,
hi.intensity=100000,
lo.conc=-1000,
hi.conc=1000),
model=fc,
concentrations=pass2,
lower=pass2,
upper=pass2,
intensities=fitted(fc, pass2),
ss.ratio=ss.ratio,
conf.width=0,
warn=warn2,
version=packageDescription("SuperCurve", fields="Version"))
if (trim > 0) {
if (verbose) {
cat("Trimming concentrations...", "\n")
flush.console()
}
progmethod("trim")
tc <- trimConc(fc,
conc=fitted(result, "X"),
intensity=intensity,
design=design,
trimLevel=trim)
if (any(is.nan(c(tc$lo.conc, tc$hi.conc)))) {
warning("concentration values zero'd out due to NaNs",
immediate.=TRUE)
series.conc <- rep(0, length(result@concentrations))
names(series.conc) <- names(result@concentrations)
} else {
series.conc <- result@concentrations
series.conc[series.conc < tc$lo.conc] <- tc$lo.conc
series.conc[series.conc > tc$hi.conc] <- tc$hi.conc
}
minunique <- 5 ## :TBD: Magic# for undetermined limit
nunique <- length(unique(sort(series.conc)))
if (!(nunique > minunique)) {
warning(sprintf("trim level %s is too high: #unique conc. values=%d",
as.character(trim),
nunique),
immediate.=TRUE)
}
result@concentrations <- series.conc
result@trimset <- unlist(tc)
}
## Should confidence intervals for the estimates be computed?
if (ci) {
if (verbose) {
cat("Computing confidence intervals...", "\n")
flush.console()
}
result <- getConfidenceInterval(result,
progmethod=progmethod)
}
result
}
##-----------------------------------------------------------------------------
## Extracts model coefficients from objects returned by modeling functions.
setMethod("coef", signature(object="RPPAFit"),
function(object,
...) {
callGeneric(object@model)
})
setMethod("coefficients", signature(object="RPPAFit"),
function(object,
...) {
coef(object, ...)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.