Nothing
#######################################################################
## Copyright 2008, Novartis Pharma AG
##
## This program is Open Source 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/.
## functions to analyze dose finding data according to MCP-Mod
## t-stats for MCP step
getTstat <-
function(data, n, contMat, resp = "resp", dose = "dose")
{
if (any(is.na(match(c(resp, dose), names(data))))) {
stop(resp," and/or ", dose, " not found in data")
}
dose <- data[, dose]
resp <- data[, resp]
## means per dose
mn <- tapply(resp, dose, mean)
## pooled standard deviation
sdv <- tapply(resp, dose, sd)
if (length(n) == 1) n <- rep(n, length(sdv))
# remove NAs for dose groups with only 1 obs.
if(any(n==1)){
sdv[n==1] <- 0
}
n1 <- n - 1
sdv <- sqrt(sum(n1 * sdv^2)/sum(n1))
## contrasts
ct <- as.vector(mn %*% contMat)
den <- sdv * sqrt(colSums((contMat^2)/n))
## max t-stat
val <- ct/den
names(val) <- dimnames(contMat)[[2]]
val
}
# function to calculate p-values
pValues <-
function(cMat, n, alpha = 0.025, Tstats, control = mvtnorm.control(),
twoSide = FALSE, corMat = NULL, ...)
{
## function to calculate p-values
##
## cMat - model contrast matrix nDose x nMod
## n - scalar or vector with sample size(s)
## alpha - significance level
## control - see mvtnorm.control function
## twoSide - logical variable indicating if one or two-sided test
## corMat - correlation matrix
nD <- nrow(cMat)
nMod <- ncol(cMat)
if(length(Tstats) != nMod){
stop("Tstats needs to have length equal to the number of models")
}
if (length(n) == 1) {
nDF <- nD * (n - 1)
} else {
if (length(n) != nD) {
stop("'n' must have length as number of doses")
}
nDF <- sum(n) - nD
}
if(is.null(corMat)){
corMat <- t(cMat)%*%(cMat/n)
den <- sqrt(crossprod(t(colSums(cMat^2/n))))
corMat <- corMat / den
}
ctrl <- mvtnorm.control()
if (!missing(control)) {
control <- as.list(control)
ctrl[names(control)] <- control
}
if(twoSide){
lower <- matrix(rep(-Tstats, each = nMod), nrow = nMod)
}
else lower <- matrix(rep(-Inf, nMod^2), nrow = nMod)
upper <- matrix(rep(Tstats, each = nMod), nrow = nMod)
pVals <- numeric(nMod)
for(i in 1:nMod){
pVals[i] <- 1 - pmvt(lower[,i], upper[,i], df = nDF, corr = corMat,
algorithm = ctrl, ...)
}
pVals
}
### functions to fit models, and functions to get starting values
### for the nls function
### Initial estimates for built-in models
### two auxiliary functions to get starting values for the nls algorithm
getInitLRasymp <-
## left and right asymptotes
function(meanVal, dlt = 0.1)
{
rg <- range(meanVal)
dlt <- dlt * diff(rg)
c(e0 = rg[1] - dlt, eMax = rg[2] + dlt)
}
getInitP <-
## dose that gives certain percentage of max effect
function(meanVal, eVec = getInitLRasymp(meanVal), p = 0.5, dose, ve=FALSE)
{
e0 <- eVec[1]
eMax <- eVec[2]
targ <- e0 * (1 - p) + eMax * p
ind <- meanVal <= targ
ind1 <- !ind
if (!any(ind) || all(ind)) stop("cannot find initial values")
ind <- ind & (dose <= max(dose[ind1]))
d1 <- max(dose[ind]); p1 <- meanVal[dose == d1]
ind1 <- ind1 & dose > d1
d2 <- min(dose[ind1]); p2 <- meanVal[dose == d2]
res <- d1 + (d2 - d1) * (targ - p1)/(p2 - p1)
names(res) <- NULL
res
}
getInitBeta <-
function(m, scal, dose)
{
## calculates starting values for the beta model
## for the nls algorithm
## dose and m assumed to be ordered
## in the same order (according to dose)
dmax <- dose[which(m==max(m))]
e0 <- m[dose==0]
emax <- max(m) - e0
ds <- dose[order(m)[length(m)-1]] # select as 2nd dose the one
# with 2nd highest resp.
z <- dmax/(scal-dmax)
f <- (m[dose==ds] - e0)/emax
beta <- try(log(f)/(log(ds^z*(scal-ds)) - log(dmax^z*(scal-dmax))))
alpha <- z*beta
if(is.na(alpha)) alpha <- beta <- 1 # may occur if f < 0; (1,1) as
# initial values seems always to
# work quite well
res <- c(alpha, beta)
names(res) <- NULL
res
}
getInit <-
## Initial estimates for built-in models, if needed
## non-linear fitting is done with p-linear, but Gauss-Newton used
## later due to some problems with se.fit of the predictions of the
## p-linear fit
function(data, model = c("emax", "exponential", "logistic", "betaMod",
"sigEmax"), scal)
{
model <- match.arg(model) # only built-in allowed
meanVal <- data$respM
dose <- data$doseM
switch(model,
emax = {
c(ed50 = getInitP(meanVal, dose = dose))
},
exponential = {
e0 <- getInitLRasymp(meanVal)[1]
meanVal <- meanVal - e0
aux <- coef(lm(log(meanVal) ~ dose, na.action = na.omit))
names(aux) <- NULL
c(delta = 1/aux[2])
},
logistic = {
ed50 <- getInitP(meanVal, dose = dose)
delta <- getInitP(meanVal, p = 0.75, dose = dose) - ed50
c(ed50 = ed50, delta = delta)
},
betaMod = {
init <- getInitBeta(meanVal, scal, dose)
c(delta1 = init[1], delta2 = init[2])
},
sigEmax = {
ed50 <- getInitP(meanVal, dose = dose)
ed10 <- getInitP(meanVal, p = 0.1, dose = dose)
ed90 <- getInitP(meanVal, p = 0.9, dose = dose)
h <- 1.91/log10(ed90/ed10) # From Ch. 14, Ting (2006),
# Dose finding in drug development
c(ed50 = ed50, h = h)
})
}
fitModel <-
## fits dose-response model from list of built-in models,
## or according to model provided by user
function(data, model, resp = "resp", dose = "dose", start,
control = NULL, off = 1, scal = 1, uModPars = NULL,
addArgs = NULL, na.action = na.fail)
{
## data - data.frame with information on dose and response
## and possibly other covariates used in model
## model - string with name of model function to be used
## resp, dose - characters with names of columns in data
## corresponding to response and dose
## control - optional list with control parameters for fit
## off - offset for dealing with dose = 0 in lin-log model
## uModPars - optional character vector with names/expressions
## of user-defined model parameters (names(start) used by
## default)
## addArgs - optional character vector with names of additional
## arguments (variables) to user-defined model
## na.action - function specifying how to handle NAs in data
## ensuring resp and dose are defined in data
data$dose <- data[, dose]
data$resp <- data[, resp]
## reacting to possible NAs
## it would be better to only take into account NAs in
## varibles actually used in the model.
data <- na.action(data)
## checking if built-in model and assigning it number, if so
builtIn <- c("linlog", "linear", "quadratic", "emax",
"exponential","logistic", "betaMod", "sigEmax")
modelNum <- match(model, builtIn)
## under the built-in models with no covariates, the mean
## modeling can be used in all cases
respM <- tapply(data$resp, data$dose, mean) # fit models based on means
doseM <- sort(unique(data$dose))
n <- as.vector(table(data$dose))
dataM <- data.frame(doseM = doseM, respM = respM, n = n)
vars <- tapply(data$resp, data$dose, var)
# allow for dose groups with only one patient
vars[n==1] <- 0 # replace NAs with 0
S2 <- sum((n - 1) * vars)
if (!is.na(modelNum)) { # built-in model
if (is.element(modelNum, 1:3)) { # linear models
if (modelNum == 1) {
dataM$off <- rep(off, nrow(dataM))
## linear-log model
fm <- lm(respM ~ I(log(doseM + off)), dataM, qr=TRUE, weights = n)
base <- coef(fm)[1]+coef(fm)[2]*log(off)
}
if (modelNum == 2) {
## linear model
fm <- lm(respM ~ doseM, dataM, qr=TRUE, weights = n)
base <- coef(fm)[1]
}
if (modelNum == 3) {
## quadratic model
fm <- lm(respM ~ doseM + I(doseM^2), dataM, qr=TRUE, weights = n)
base <- coef(fm)[1] # baseline
}
} else { # nonlinear built-in model
if (is.null(start)) { # need to derive initial values
start <- getInit(dataM, model, scal)
}
if (modelNum == 4) { # Emax model
names(start) <- NULL
start <- c(led50 = log(start))
fm <- try(nls(respM ~ cbind(1, emax(doseM, 0, 1, exp(led50))),
start = start, data = dataM, model = TRUE,
algorithm = "plinear", control = control,
weights = n), silent = TRUE)
if (!inherits(fm, "nls")) {
fm <- NA
} else {
base <- coef(fm)[2]
}
}
if (modelNum == 5) { # exponential model
names(start) <- NULL
start <- c(ldelta = log(start))
fm <-
try(nls(respM ~ cbind(1, exponential(doseM, 0, 1, exp(ldelta))),
start = start, weights = n, data = dataM, model = TRUE,
control = control, algorithm = "plinear"), silent = TRUE)
if (!inherits(fm, "nls")) {
fm <- NA
} else {
base <- coef(fm)[2]
}
}
if (modelNum == 6) { # logistic model
start <- c(log(start["ed50"]), start["delta"])
names(start) <- c("led50", "delta")
fm <- try(nls(respM ~ cbind(1, logistic(doseM, 0, 1, exp(led50),
delta)),
start = start, algorithm = "plinear", data = dataM,
model = TRUE, control = control, weights = n),
silent = TRUE)
if (!inherits(fm, "nls")) {
fm <- NA
} else {
base <- predict(fm, newdata = list(doseM = 0))
}
}
if (modelNum == 7) { # beta model
dataM$scal <- rep(scal, nrow(dataM))
fm <-
try(nls(respM ~ cbind(1, betaMod(doseM, 0, 1, delta1, delta2, scal)),
start = start, weights = n, data = dataM, model = TRUE,
algorithm = "plinear", control = control), silent = TRUE)
if (!inherits(fm, "nls")) {
fm <- NA
} else {
base <- coef(fm)[3]
}
}
if (modelNum == 8) { # sigEmax model
start <- c(log(start["ed50"]), start["h"])
names(start) <- c("led50", "h")
fm <-
try(nls(respM ~ cbind(1, sigEmax(doseM, 0, 1, exp(led50), h)),
start = start, model = TRUE,
algorithm = "plinear", control = control, weights = n),
silent = TRUE)
if (!inherits(fm, "nls")) {
fm <- NA
} else {
base <- coef(fm)[3]
}
}
}
if (length(fm) > 1) ResidSS <- S2 + deviance(fm)
} else { # user defined model
## first check for starting estimates, parameter names, etc
if (is.null(start)) {
stop("must provide stating estimates for user-defined model")
}
namStart <- names(start)
if (is.null(namStart)) {
stop("'start' must have names for user-defined models")
}
if (is.null(uModPars)) uModPars <- namStart # default parameters
## create the model call
modForm <- paste("resp ~ ", model, "(dose,",
paste(uModPars, collapse = ","), sep = "")
if (!is.null(addArgs)) {
modForm <- paste(modForm, ",", paste(addArgs, collapse = ","))
}
modForm <- paste(modForm, ")")
modForm <- eval(parse(text = modForm))
## will assume non-linear model
fm <- try(do.call("nls", list(modForm, data, start, control)))
if (!inherits(fm, "nls")) {
fm <- NA
} else {
ResidSS <- deviance(fm)
dataM <- NULL
## following will not work when covariates are used in nls model
base <- predict(fm, newdata = list(dose = 0))
}
}
if (length(fm) > 1) {
res <- list(fm = fm, base = base)
## saving information for other methods
attr(res, "ResidSS") <- ResidSS
attr(res, "model") <- model
attr(res, "scal") <- scal
attr(res, "off") <- off
attr(res, "dataM") <- dataM
} else {
res <- NA
}
oldClass(res) <- "fitMCPMod" # allow use of methods later
res
}
# function to return analyt. gradient for built-in or
# user defined models
getGrad <-
function(obj, dose, uGrad = NULL)
## obj - an object inheriting from class fitMCPMod
## dose - vector with doses at which to calculate the gradient
## uGrad - optionally, a character string with the name of a
## user-defined gradient function
{
if(!inherits(obj, "fitMCPMod")) {
stop("obj must inherit from class fitMCPMod")
}
model <- attr(obj, "model")
fm <- obj$fm
if (!inherits(fm, "nls")) {
## linear models are left out
stop("fitted object must inherit from class nls")
}
## only used internally in getGrad
lg2 <- function(x) ifelse(x == 0, 0, log(x))
cf <- coef(obj) # estimated parameters
res <-
switch(model,
"emax" = {
eMax <- cf[2]; ed50 <- cf[3]
cbind(1, dose/(ed50 + dose), -eMax*dose/(dose + ed50)^2)
},
"logistic" = {
eMax <- cf[2]; ed50 <- cf[3]; delta <- cf[4]
den <- 1 + exp((ed50 - dose)/delta)
g1 <- -eMax*(den-1)/(delta*den^2)
g2 <- eMax*(den-1)*(ed50-dose)/(delta^2*den^2)
cbind(1, 1/den, g1, g2)
},
"sigEmax" = {
eMax <- cf[2]; ed50 <- cf[3]; h <- cf[4]
den <- (ed50^h + dose^h)
g1 <- dose^h/den
g2 <- -ed50^(h-1)*dose^h*h*eMax/den^2
g3 <- eMax*dose^h*ed50^h*lg2(dose/ed50)/den^2
cbind(1, g1, g2, g3)
},
"betaMod" = {
scal <- attr(obj, "scal")
dose <- dose/scal
if (any(dose > 1)) {
stop("doses cannot be larger than scal in betaModel")
}
delta1 <- cf[3]; delta2 <- cf[4]; eMax <- cf[2]
maxDens <- (delta1^delta1)*(delta2^delta2)/
((delta1 + delta2)^(delta1+delta2))
g1 <- ((dose^delta1) * (1 - dose)^delta2)/maxDens
g2 <- g1*eMax*(lg2(dose)+lg2(delta1+delta2)-lg2(delta1))
g3 <- g1*eMax*(lg2(1-dose)+lg2(delta1+delta2)-lg2(delta2))
cbind(1, g1, g2, g3)
},
"exponential" = {
delta <- cf[3]
e1 <- cf[2]
cbind(1, exp(dose/delta) - 1, -exp(dose/delta)*dose*e1/delta^2 )
},
{
## user defined gradient function
if(is.null(uGrad)) {
stop("user-defined gradient needs to be specified")
}
do.call(uGrad, c(list(dose), cf))
})
res
}
# function to return coefficients from fit
coef.fitMCPMod <-
function(object, ...)
{
## object - object inheriting from class MCPMod
fm <- object$fm
cf <- coef(fm)
if(inherits(fm, "lm")) return(cf)
temp <- names(cf) # save names for user-models
names(cf) <- NULL
model <- attr(object, "model")
cf <- switch(model,
"emax" = c(e0 = cf[2], eMax = cf[3], ed50 = exp(cf[1])),
"logistic" = c(e0=cf[3], eMax=cf[4], ed50 = exp(cf[1]),
delta = cf[2]),
"exponential" = c(e0 = cf[2], e1 = cf[3], delta = exp(cf[1])),
"betaMod" = c(e0=cf[3], eMax=cf[4], delta1=cf[1], delta2=cf[2]),
"sigEmax" = c(e0=cf[3], eMax=cf[4], ed50=exp(cf[1]), h=cf[2]),
{
names(cf) <- temp
cf
}
)
cf
}
# predict method
predict.fitMCPMod <-
function(object, doseSeq, se.fit = TRUE, uGrad = NULL, ...)
{
# object - either a lm object or a nls object in the latter case
# the code distinguishes between user-defined models
# and built-in models
# doseSeq - sequence for predictions
# se.fit - logical indicating whether sd of the mean should be
# calculated
model <- attr(object, "model")
scal <- attr(object, "scal")
off <- attr(object, "off")
fm <- object$fm
if(inherits(fm, "lm")){ # in this case use the implemented predict method
newdata <- data.frame(doseM = doseSeq, off = rep(off, length(doseSeq)))
res <- predict(fm, newdata, se.fit = se.fit)
## Need to correct se.fit, if requested
if(se.fit) {
df <- sum(fm$weights) - length(coef(fm))
sigCorr <- sqrt(attr(object, "ResidSS")/df)
sigOrig <- summary(fm)$sigma
res$se.fit <- (sigCorr/sigOrig) * res$se.fit
res$df <- df
res$residual.scale <- sigCorr
}
return(res)
}
builtIn <- is.element(model, c("emax", "exponential", "logistic",
"betaMod", "sigEmax"))
if(inherits(fm, "nls") & !builtIn){ # user defined model
mn <- predict(fm, data.frame(dose = doseSeq))
if(se.fit) {
## first, check if model includes a gradient attribute
grd <- attr(mn, "gradient")
if(is.null(grd)) {
if(is.null(uGrad)) {
stop("neither gradient attribute, nor uGrad defined")
}
grd <- getGrad(object, doseSeq, uGrad)
}
Rinv <- solve(fm$m$Rmat())
sig <- summary(fm)$sigma
seFit <- sig * sqrt(rowSums((grd %*% Rinv)^2))
df <- df.residual(fm)
res <- list(fit = mn, se.fit = as.vector(seFit),
residual.scale = sig, df = df)
return(res)
} else {
return(mn)
}
} else { # built in model
dd <- data.frame(doseM = doseSeq) # dose levels to be predicted
cf <- coef(object) # extract coefficients
if(model == "betaMod"){
cf <- c(cf, scal) # add scale parameter for betaModel
dd$scal <- scal
}
mn <- predict(fm, dd) # predict mean response
if(!se.fit){
return(mn)
}
RSS <- attr(object, "ResidSS") # get residual sum of squares
doseV <- attr(object, "dataM")$doseM # recover dose-vector
n <- fm$weights
df <- sum(n) - length(cf)
sig <- sqrt(RSS/df) # residual variance estimate
## Gradient matrix: see Bates/Watts p. 58/59
V <- getGrad(object, doseV)
if(all(!is.infinite(V)) & all(!is.nan(V))){
## checking for infinities (can happen because the analytic
## gradient is not used for fitting)
V <- t(crossprod(V, sqrt(diag(n)))) # weights for unequal sample sizes
R <- qr.R(qr(V))
Rinv <- try(solve(R))
if (!inherits(Rinv, "matrix")){
## sometimes R is singular (typically for very unequally
## spaced dose designs)
warning("Cannot cannot calculate standard deviation for ",
model, " model.\n")
seFit <- rep(NA, length(doseSeq))
} else {
v <- getGrad(object, doseSeq) # calc. gradient
seFit <- sig * sqrt(rowSums((v%*%Rinv)^2))
}
} else {
warning("Cannot cannot calculate standard deviation for ",
model, " model.\n")
seFit <- rep(NA, length(doseSeq))
}
res <- list(fit = mn, se.fit = as.vector(seFit),
residual.scale = sig, df = df)
}
res
}
# Calculate information criteria for nls and lm objects
AIC.fitMCPMod <-
function(object, ..., k = 2)
{
# object - an object of class fitMCPMod
# k - penalty term
fm <- object$fm
if (is.null(fm$weights)) {
## user defined model or non-averaged modeling
return(AIC(fm, k = k))
}
RSS <- attr(object, "ResidSS")
n <- sum(fm$weights)
sig2 <- RSS/n
logL <- -n/2*(log(2*pi) + 1 + log(sig2))
-2*logL + k*(length(coef(object)) + 1) # "+ 1" because of sigma parameter
}
## model selection
modelSelect <-
function(data, namSigMod, selMethod, pW, resp, dose, start, nlsControl,
off, scal, uModPars, addArgs)
{
fm <- NA
warn <- NULL
nSigMod <- length(namSigMod)
if (selMethod == "maxT") { # first maxT option
i <- 1
while((length(fm) == 1) && (i <= nSigMod)) {
nam <- namSigMod[i]
fm <- fitModel(data, nam, resp, dose, start[[nam]], nlsControl,
off, scal, uModPars[[nam]], addArgs[[nam]])
if(length(fm) == 1) # model didn't converge
warning(nam, " model did not converge\n")
i <- i + 1
}
if (length(fm) > 1) { # model converged
fm <- list(fit = list(fm), base = fm$base)
attr(fm$fit, "model2") <- names(fm$fit) <- nam
} else {
fm <- list(fit = NA, base = NA) # no model converged
}
} else { # AIC, BIC, aveAIC or aveBIC
fm <- vector("list", nSigMod)
crit <- fmBase <- rep(NA, nSigMod)
if (selMethod == "AIC"| selMethod == "aveAIC") {
pen <- 2
} else {
pen <- log(dim(data)[[1]])
}
names(fm) <- names(crit) <- names(fmBase) <- namSigMod
for(i in 1:nSigMod) {
nam <- namSigMod[i]
fitmod <- fitModel(data, nam, resp, dose, start[[nam]], nlsControl,
off, scal, uModPars[[nam]], addArgs[[nam]])
if(!is.list(fitmod)) { # model didn't converge
fm[[i]] <- NA
warning(nam, " model did not converge\n")
} else { # model converged
fm[[i]] <- fitmod
fmBase[i] <- fitmod$base
crit[i] <- AIC(fitmod, k = pen)
}
}
if (all(is.na(crit))) {
fm <- NA
} else {
if (selMethod == "AIC" | selMethod == "BIC") {
model2 <- namSigMod[which(crit == min(crit, na.rm = TRUE))]
fm <- list(fm[[model2]])
fmBase <- fmBase[model2]
attr(fm, "model2") <- names(fm) <- model2
attr(fm, "IC") <- crit
}
else { # model averaging
attr(fm, "model2") <- namSigMod[!is.na(fmBase)]
attr(fm, "IC") <- crit
crit <- crit[!is.na(crit)]
## subtract const from crit values to avoid numerically 0
## values (exp(-0.5*1500)=0!)
const <- mean(crit)
if(is.null(pW)){
pW <- rep(1, length(crit)) # standard 'noninformative' prior
names(pW) <- names(crit)
} else {
pW <- pW[names(crit)]
if(any(is.na(pW))) stop("pW needs to be a named vector with names equal to the models in the candidate set")
}
attr(fm, "weights") <-
pW*exp(-0.5*(crit-const))/sum(pW*exp(-0.5*(crit-const)))
attr(fm, "pweights") <- pW
}
}
fm <- list(fit = fm, base = fmBase)
}
fm
}
## dose estimation
getDose <-
function(dose, ind)
{
aa <- !is.na(ind)
if (!all(aa)) {
ind <- ind[aa]; dose <- dose[aa]
}
if (any(ind)) min(dose[ind])
else NA
}
getDoseEst <-
function(fmb, clinRel, rangeDose, doseEst, selMethod,
dePar = 0.1, lenDose = 101, uGrad = NULL)
{
## equally spaced points within dose range for prediction
doseSeq <- seq(rangeDose[1], rangeDose[2], len = lenDose)
off <- attr(fmb, "off"); scal <- attr(fmb, "scal")
newdata <- list(dose = doseSeq, off = rep(off, lenDose), scal = rep(scal, lenDose))
ldePar <- length(dePar)
val <- double(ldePar)
gma <- 1-dePar
base <- fmb$base
tDose <- matrix(ncol = length(base)-sum(is.na(base)), nrow = length(dePar))
z <- 1
for(m in which(!is.na(base))){ # only predict models that converged
if(doseEst[1] != "ED"){ # MED estimate
pred <- predict(fmb$fit[[m]], doseSeq, uGrad = uGrad)
fo <- data.frame(pred = pred$fit - base[m], se.fit = pred$se.fit)
df <- pred$df
## selects the desired dose according to MED method
for(i in 1:ldePar) {
crt <- qt(gma[i], df)
LL <- fo$pred - crt * fo$se.fit
UL <- fo$pred + crt * fo$se.fit
switch(doseEst,
"MED1" = ind <- LL >= 0 & UL >= clinRel,
"MED2" = ind <- LL >= 0 & fo$pred >= clinRel,
"MED3" = ind <- LL >= clinRel
)
val[i] <- getDose(doseSeq, ind)
}
} else { # ED estimate
pred <- predict(fmb$fit[[m]], doseSeq, se.fit = FALSE)
pEff <- dePar*(max(pred) - base[m])
for(i in 1:ldePar) {
ind <- pred >= pEff[i] + base[m]
val[i] <- getDose(doseSeq, ind)
}
}
tDose[,z] <- val
z <- z + 1
}
if(is.element(selMethod, c("aveAIC", "aveBIC"))){
doseAve <- as.vector(tDose%*%attr(fmb$fit, "weights"))
if(doseEst[1] == "ED")
names(doseAve) <- paste(doseEst, round(100 * dePar), "%", sep = "")
else
names(doseAve) <- paste(doseEst,"," , round(100 * (1-2*dePar)), "%",
sep = "")
dimnames(tDose) <- list(names(doseAve), attr(fmb$fit, "model2"))
attr(doseAve, "tdModels") <- tDose
tDose <- doseAve
} else {
tDose <- as.vector(tDose)
if(doseEst == "ED")
names(tDose) <- paste(doseEst, round(100 * dePar), "%", sep = "")
else
names(tDose) <- paste(doseEst,"," , round(100 * (1-2*dePar)), "%",
sep = "")
}
tDose
}
recovNames <-
function(names)
{
## function to recover model names (in case of multiple models from one class)
## just for use in MCPMod function
## example: recovNames(c("emax1", "betaMod", "emax2", "logistic", "usermodel"))
## returns: c("emax","betaMod","logistic","usermodel")
builtIn <- c("linlog", "linear", "quadratic", "emax",
"exponential","logistic", "betaMod", "sigEmax")
newnames <- character()
i <- 1
for(nam in names){
pm <- pmatch(builtIn, nam)
if(any(!is.na(pm))) {
newnames[i] <- builtIn[which(!is.na(pm))]
i <- i+1
}
if(all(is.na(pm))) {
newnames[i] <- nam
i <- i+1
}
}
unique(newnames)
}
## function to get reasonable defaults for, off, scal and dePar
getDef <-
function(off = NULL, scal = NULL, dePar = NULL, doses, doseEst)
{
mD <- max(doses)
if(is.null(dePar)){ # default for dePar depending on the dose estimate
dePar <- ifelse(doseEst=="ED", 0.5, 0.1)
}
if(is.null(scal)){ # default for scal parameter
scal <- 1.2*mD
} else { # check if valid scal is provided
if(scal < mD){
stop("'scal' should be >= maximum dose")
}
}
if(is.null(off)){ # default for off parameter
off <- 0.1*mD
}
list(scal = scal, off = off, dePar = dePar)
}
### main function implementing methodology, combining several others
MCPMod <-
function(data, models = NULL, contMat = NULL, critV = NULL,
resp = "resp", dose = "dose", off = NULL, scal = NULL,
alpha = 0.025, twoSide = FALSE,
selModel = c("maxT", "AIC", "BIC", "aveAIC", "aveBIC"),
doseEst = c("MED2", "MED1", "MED3", "ED"),
std = TRUE, start = NULL,
uModPars = NULL, addArgs = NULL,
dePar = NULL, clinRel = NULL, lenDose = 101, pW = NULL,
control = list(maxiter = 100, tol = 1e-6, minFactor = 1/1024),
signTtest = 1, pVal = FALSE, testOnly = FALSE,
mvtcontrol = mvtnorm.control(), na.action = na.fail,
uGrad = NULL)
{
## data - data.frame with, at least, columns "resp" and "dose"
## models - list with component names specifying models
## and optionally elements being used to calculate
## model contrasts; need to have named elements with
## corresponding models and be consistent with contMat,
## if that is non-null
## contMat - contrast matrix for candidate set and doses
## critV - critical value for MCP test
## alpha - significance level for calculating critVal (if = NULL)
## selModel - model selection criterion
## doseEst - type of dose estimator to be used
## dePar - gives "gamma" parameter for MED-type estimators or
## the "p" in the EDp estimate
## pW - named vector of prior probs. for different models
if (any(is.na(match(c(resp, dose), names(data))))) {
stop(resp," and/or ", dose, " not found in data")
}
data <- na.action(data)
ind <- match(dose, names(data))
data <- data[order(data[, ind]), ]
## target dose estimate type
doseEst <- match.arg(doseEst)
n <- as.vector(table(data[, dose])) # sample sizes per group
doses <- sort(unique(data[, dose]))
## getting defaults which depend on the DRdata
def <- getDef(off, scal, dePar, doses, doseEst)
scal <- def[[1]]; off <- def[[2]]; dePar <- def[[3]]
if (is.null(contMat)) {
## need to calculate it
mu <- modelMeans(models, doses, std, off, scal)
contMat <- modContr(mu, n)
}
## MCP test first
tStat <- signTtest * getTstat(data, n, contMat, resp, dose)
if(twoSide){
tStat <- abs(tStat)
}
if (is.null(critV)) {
if(pVal) {
pVals <- pValues(contMat, n, alpha, tStat, mvtcontrol, twoSide)
}
critV <- critVal(contMat, n, alpha, mvtcontrol, twoSide = twoSide)
attr(critV, "Calc") <- TRUE # determines whether cVal was calculated
} else {
pVal <- FALSE # pvals are not calculated if critV is supplied
attr(critV, "Calc") <- FALSE
}
indStat <- tStat > critV
if (!any(indStat) | testOnly) {
## only mcp-test or no significant t-stats
result <- list(signf = any(indStat), model1 = NA, model2 = NA)
} else {
## model selection method
selMethod <- match.arg(selModel)
control <- do.call("nls.control", control)
## at least one significant, select a model if possible
namMod <- names(tStat)
maxTstat <- max(tStat)
model1 <- namMod[which(tStat == maxTstat)] # model with most sig contrast
## significant models, in descending order of tstat
indSigMod <- 1:length(tStat)
indSigMod <- indSigMod[rev(order(tStat))][1:sum(indStat)]
namSigMod <- namMod[indSigMod] # significant models
namSigMod <- recovNames(namSigMod) # remove model nrs.
# (in case of multiple models)
fmb <-
modelSelect(data, namSigMod, selMethod, pW, resp, dose, start,
control, off, scal, uModPars, addArgs)
if (all(is.na(fmb$base))) { # none of sign. model converged
result <- list(signf = TRUE, model1 = model1, model2 = NA)
} else {
## dose estimation model(s) obtained
## move to final step, dose estimation
## dose range
rgDose <- range(doses)
tDose <- getDoseEst(fmb, clinRel, rgDose, doseEst, selMethod,
dePar, lenDose, uGrad)
result <- list(signf = TRUE, model1 = model1,
model2 = attr(fmb$fit, "model2"))
}
}
## add information to the object
result$input <- list(models=models, resp=resp, dose=dose, off=off,
scal=scal, alpha=alpha, twoSide=twoSide,
selModel=selModel[1], doseEst=doseEst,
std=std, dePar=dePar, uModPars=uModPars,
addArgs=addArgs, start = start, uGrad=uGrad,
clinRel=clinRel, lenDose=lenDose, signTtest=signTtest,
pVal=pVal, testOnly=testOnly)
result$data <- data
result$contMat <- contMat
result$corMat <- t(contMat)%*%(contMat/n) /
sqrt(crossprod(t(colSums(contMat^2/n))))
result$cVal <- critV
result$tStat <- tStat
if(pVal) attr(result$tStat, "pVal") <- pVals
if (!all(is.na(result$model2))){
result$fm <- fmb$fit
result$tdose <- tDose
}
oldClass(result) <- "MCPMod"
result
}
# print method for MCPMod objects
print.MCPMod <-
function(x, digits = 3,...)
{
cat("MCPMod\n\n")
if(x$input$twoSide) side <- "two-sided"
else side <- "one-sided"
if(!attr(x$cVal, "Calc")){
cat(paste("PoC:",sep=""),
paste(c("no", "yes")[x$signf+1]), "\n")
} else {
cat(paste("PoC (alpha = ", x$input$alpha,", ", side, "):",sep=""),
paste(c("no", "yes")[x$signf+1]), "\n")
}
if(x$signf & !x$input$testOnly){
cat("Model with highest t-statistic:", paste(x$model1),"\n")
if(!all(is.na(x$model2))){
modSel <- is.element(x$input$selModel, c("AIC", "BIC", "maxT"))
mo <- ifelse(modSel, "Model", "Models")
cat(paste(mo, "used for dose estimation:"), paste(x$model2),"\n")
cat("Dose estimate:","\n")
attr(x$tdose, "tdModels") <- NULL
print(round(x$tdose, digits))
} else if(!x$input$testOnly)
cat("\nNone of significant models converged")
}
cat("\n")
}
summary.MCPMod <-
function(object, digits = 3,...)
{
oldClass(object) <- "summary.MCPMod"
print(object, digits = digits)
}
print.summary.MCPMod <-
function(x, digits = 3,...)
{
cat("MCPMod\n\n")
if(x$input$twoSide) side <- "two-sided"
else side <- "one-sided"
cat("Input parameters:","\n")
if(attr(x$cVal, "Calc")){
cat(" alpha =", paste(x$input$alpha," (", side,")", sep=""), "\n")
}
if(!x$input$testOnly){
cat(" model selection:", paste(x$input$selModel[1]),"\n")
if(is.element(x$input$selModel[1], c("aveAIC", "aveBIC"))){
pW <- attr(x$fm, "pweights")
cat(" prior model weights:\n ")
print(round(pW/sum(pW), digits))
}
nr <- match(substr(x$input$doseEst,1,2), c("ME", "ED"))
if(nr == 1){
cat(" clinical relevance =",paste(x$input$clinRel),"\n")
}
dePar <- c("gamma", "p")[nr]
cat(paste(" dose estimator: ", x$input$doseEst, " (",dePar, " = ", x$input$dePar, ")", sep=""), "\n")
} # multiple contrast test information
cat("\n","Optimal Contrasts:","\n", sep="")
print(round(x$contMat, digits))
cat("\n","Contrast Correlation:","\n", sep="")
print(round(x$corMat, digits))
cat("\n","Multiple Contrast Test:","\n",sep="")
ord <- rev(order(x$tStat))
if(x$input$pVal){
print(data.frame(Tvalue = round(x$tStat, digits)[ord],
pValue = round(attr(x$tStat, "pVal"), digits)[ord]))
}
else {
print(data.frame(Tvalue = round(x$tStat, digits)[ord]))
}
cat("\n","Critical value: ", round(x$cVal, digits),"\n",sep="")
if(!x$signf) return(cat("")) # No significant model
if(all(is.na(x$model2))) {
if(!x$input$testOnly) cat("\nNone of significant models converged","\n")
return(cat(""))
}
else { # add IC values
if(x$input$selModel[1] != "maxT"){
if(x$input$selModel[1] == "AIC" | x$input$selModel[1]=="aveAIC")
cat("\n",paste("AIC criterion:"),"\n",sep="")
else cat("\n",paste("BIC criterion:"),"\n",sep="")
print(round(attr(x$fm,"IC"), 2))
} # model selection
cat("\n","Selected for dose estimation:","\n", sep="")
cat(" ", paste(x$model2),"\n\n", sep=" ")
if(is.element(x$input$selModel[1], c("maxT", "AIC", "BIC"))){
cat("Parameter estimates:","\n")
cat(paste(x$model2), "model:\n")
cof <- coef(x$fm[[1]])
namcof <- names(cof)
namcof <- gsub(" ", "", namcof) # remove white spaces for GUI
names(cof) <- gsub("doseM", "dose", namcof) # use more obvious names
print(round(cof, digits))
} else { # model averaging
cat("Model weights:","\n", sep="")
print(round(attr(x$fm, "weights"), digits))
cat("\nParameter estimates:","\n")
for(i in which(!is.na(attr(x$fm, c("IC"))))){
nam <- names(x$fm)[i]
cof <- coef(x$fm[[i]])
namcof <- names(cof)
namcof <- gsub(" ", "", namcof) # remove white spaces for GUI
names(cof) <- gsub("doseM", "dose", namcof) # use more obvious names
cat(paste(nam), "model:\n")
print(round(cof, digits))
}
}
} # information about dose estimate
cat("\nDose estimate","\n")
attr(x$tdose,"tdType") <- NULL # remove attr for output
if(is.element(x$input$selModel, c("AIC", "BIC", "maxT"))) print(x$tdose)
else {
cat("Estimates for models\n")
print(round(attr(x$tdose,"tdModels"), digits))
attr(x$tdose,"tdModels") <- NULL
cat("Model averaged dose estimate\n")
print(round(x$tdose, digits))
}
}
plot.MCPMod <-
function(x, complData = FALSE, CI = FALSE, clinRel = FALSE, doseEst = FALSE,
gamma = NULL, models = "all", nrDoseGam = 1,
colors = c("black","blue","black","gray","blue"),
uGrad = NULL, ...)
{
## complData - logical indicating whether complete
## data set (T) or just group means
## should be plotted (F)
## CI - logical indicating whether confidence intervals
## should be plotted
## clinRel - logical indicating whether clinical relevance
## threshold should be plotted
## doseEst - logical whether dose estimate should be plotted
## in case a vector was specified for dePar in the
## MCPMod function, nrDoseGam determines which is plotted
## gamma - value for the 1-2*gamma pointwise CI around
## the predicted mean. if equal to NULL the value
## determined in the MCPMod call is used. In case
## a vector of gamma values nrDoseGam determines which
## is used (see nrDoseGam)
## models - a character vector determining, which of the fitted
## models should be plotted (only for model averaging)
## nrDoseGam - in case a vector is specified for gamma in
## the MCPMod function (and gamma in the plot.MCPMod
## function is NULL), nrDoseGam determines which of
## these values should be used for the conf. interval
## and the dose estimate (if doseEst = T)
## colors - numeric of length 5 with number of colors for:
## predictions, CI, data, clinical relevance threshold,
## dose estimator
## ... - additional arguments for xyplot
selModel <- x$input$selModel[1]
## model selection or model averaging?
ms <- is.element(selModel, c("maxT", "AIC", "BIC"))
off <- x$input$off
scal <- x$input$scal
if(!is.null(x$input$resp)) resp <- x$input$resp # name of resp column
else resp <- "resp"
if(!is.null(x$input$dose)) dose <- x$input$dose # name of dose column
else dose <- "dose"
ind <- match(c(dose, resp), names(x$data))
Data <- x$data[, ind]
## if no model is significant or no model converged plot raw data
if(!x$signf | all(is.na(x$model2))){
plot(Data[,1], Data[,2], xlab="Dose", ylab="Response")
if(!x$signf) msg <- "No model significant"
else msg <- "None of significant models converged"
title(sub = msg)
return(cat(msg,"\n"))
}
if(is.null(gamma)){ # default for gamma
if(x$input$doseEst != "ED"){
gamma <- x$input$dePar[nrDoseGam] # use the same gamma as
} else { # for dose estimate
gamma <- 0.025 # if "ED" use reas. default
}
}
if(models == "all"){ # select a subset of fitted models
fm <- x$fm # (only possible in case of model averaging)
nam <- x$model2
} else {
if(ms) stop("models argument can only be used with model averaging")
fm <- x$fm[models]
nam <- models
}
lenDose <- x$input$lenDose # fit models and obtain CI
doseSeq <- seq(min(Data[,1]), max(Data[,1]), length=lenDose)
predList <- DoseInd <- list()
z <- 0
for(i in 1:length(fm)){
if(is.list(fm[[i]])){
z <- z + 1
pred <- predict(fm[[i]], doseSeq, se.fit=CI, uGrad = uGrad)
if(CI){
LL <- pred$fit - qt(1-gamma, pred$df)*pred$se.fit
UL <- pred$fit + qt(1-gamma, pred$df)*pred$se.fit
}
if(ms){
tdose <- x$tdose[nrDoseGam] # plot just dose corresp. to
# selected gamma value
} else {
tdose <- attr(x$tdose, "tdModels")
modNum <- which(names(fm)[i]==dimnames(tdose)[[2]])
tdose <- tdose[nrDoseGam, modNum]
}
DoseInd[[z]] <- ifelse(is.na(tdose), NA, which(doseSeq == tdose))
if(CI) predList[[z]] <- c(pred$fit, LL, UL)
else predList[[z]] <- pred
}
}
plotVec <- do.call("c", predList)
DoseInd <- do.call("c", DoseInd)
if(CI) groups <- c("pred","LL","UL")
else groups <- "pred"
lg <- length(groups)
plotMat <- data.frame(pred=plotVec, dose=rep(doseSeq, lg*z),
nam=rep(nam, each=lg*lenDose),
group=rep(rep(groups, rep(lenDose, lg)), z))
if(complData) { # plot all data, not just means
dat <- Data
} else {
me <- tapply(Data[,2], Data[,1], mean)
d <- as.numeric(names(me))
dat <- data.frame(d, me)
}
rg <- range(dat[,2], plotMat$pred)
yl <- c(rg[1] - 0.1*diff(rg), rg[2] + 0.1*diff(rg))
panDat <- list(data=dat, clinRel=x$input$clinRel, complData=complData,
cR=clinRel, doseEst=doseEst, lenDose=lenDose, DoseInd=DoseInd,
colors=colors, lg=lg) # data for panels
## information for key argument
txt <- ifelse(complData, "Responses", "Group Means")
pchar <- ifelse(complData, 1, 16)
keyList <- list(points= list(col = colors[3], pch=pchar),
text = list(lab = txt),
lines = list(col = colors[1]),
text = list(lab = "Model Predictions"))
if(CI){
keyList <- c(keyList, list(lines = list(col = colors[2])),
list(text = list(lab = paste(1-2*gamma, "Pointwise CI"))))
col <- c(colors[2],colors[1],colors[2])
} else col <- c(colors[1])
if(doseEst){
keyList <- c(keyList, list(points=list(col=colors[5], pch=18)),
list(text=list(lab="Estim. Dose")))
}
xyplot(pred~dose|nam, data=plotMat, xlab = "Dose", type="l",
col=col, ylab = "Response",
groups=plotMat$group, pD=panDat, ylim = yl,
panel=function(x, y, subscripts, groups,..., pD) {
panel.superpose(x,y,subscripts,groups, ...)
if(!pD$complData){
panel.xyplot(pD$data[,1],pD$data[,2], pch=16,
col=pD$colors[3]) # plot data/means
} else {
panel.xyplot(pD$data[,1],pD$data[,2], col=pD$colors[3])
}
if(pD$cR) panel.abline(h=y[1]+pD$clinRel, lty=8,
col=pD$colors[4]) # plot clinical
# relevance threshold
if(pD$doseEst) {
ind <- max(subscripts)/(pD$lenDose*pD$lg) # locate dose
# estimator in x
panel.dotplot(x[pD$DoseInd[ind]], y[pD$DoseInd[ind]],
pch=18, col=pD$colors[5], col.line=0) # plot dose
}
},
strip = function(...) strip.default(..., style = 1), as.table=TRUE,
key = keyList, ...)
}
#### example code
###
### mods <- list(linear = NULL, linlog = NULL, emax = 0.3,
### quadratic = -0.001, logistic = c(0.4, 0.09))
###
### mn <- c(0.1, 0.4, 0.55, 0.75, 0.9, 1)
### ds <- c(0, 0.05, 0.2, 0.4, 0.6, 1)
### DRdata <- genDFdata(mu = mn, sigma = 0.8, n = 20, doses = ds)
### MM <- MCPMod(DRdata, mods, clinRel = 0.5*max(mn), selModel="aveAIC")
### MM
### summary(MM)
### plot(MM)
###
### user defined model
### emx2 <- function(dose, a, b, d){
### sigEmax(dose, a, b, d, 1)
### }
### emx2Grad <- function(dose, a, b, d) cbind(1, dose/(dose+d), -b*dose/(dose+d)^2)
### models <- list(emx2=c(0,1,0.1), linear = NULL)
### dats <- genDFdata(mu = mn, doses = ds, sigma=0.8, n=50)
### start <- list(emx2=c(a=0.2, b=0.6, d=0.2))
### MM1 <- MCPMod(dats, models, clinRel = 1, scal = 1.2, selModel="AIC", start = start,
### uGrad = emx2Grad)
## function to generate DF data
genDFdata <-
function(model, argsMod, doses, n, sigma, mu = NULL)
{
## generates data.frame with resp and dose columns
## corresponding to specified design, mean (either
## determined by model or mu) and sigma
##
## model - string with model function name (first argument
## must be dose, must be vectorized)
## argsMod - a named list with the arguments for the model function
##
## doses - set of doses to be used in the trial
## n - sample size (scalar is repeated for each dose)
## sigma - error std. deviation
## mu - vector of mean values can alternatively specified
nD <- length(doses)
dose <- sort(doses)
if (length(n) == 1) n <- rep(n, nD)
dose <- rep(dose, n)
if(!missing(model)){
args <- c(list(dose), argsMod)
mu <- do.call(model, args)
} else if(!is.null(mu)){
if(length(doses) != length(mu)){
stop("'mu' needs to be of the same length as doses")
}
mu <- rep(mu, n)
} else {
stop("either 'model' or 'mu' needs to be specified")
}
data.frame(dose = dose,
resp = mu + rnorm(sum(n), sd = sigma))
}
### examples
### genDFdata("emax", c(e0 = 0.2, eMax = 1, ed50 = 0.05), c(0,0.05,0.2,0.6,1), 20, 1)
### genDFdata(mu = 1:5, doses = 0:4, n = c(20, 20, 10, 5, 1), sigma = 1)
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.