Fix guesst() sigEmax label and refactor internal model handling

R/MCPMod.R

@@ -147,8 +147,7 @@ MCPMod <- function(dose, resp, data = NULL, models = NULL, S=NULL,
   ## fit models and calculate model selection criteria
   addArgs <- list(off=attr(models, "off"), scal=attr(models, "scal"))
   selModel <- match.arg(selModel)
-  builtIn <- c("linlog", "linear", "quadratic", "linInt", "emax",
-               "exponential", "logistic", "betaMod", "sigEmax")
+  builtIn <- builtInMods
   nams <- gsub("[0-9]", "", names(tstat)) ## remove numbers from model-names
   namsU <- unique(nams)
 
@@ -175,7 +174,7 @@ MCPMod <- function(dose, resp, data = NULL, models = NULL, S=NULL,
     modcrit <- function(x)
       max(tstat[attr(x, "model") == nams])
   }
-  for(i in 1:length(namsU)){
+  for(i in seq_along(namsU)){
     if(!is.null(data)){
       callMod <- list(deparse(substitute(dose)), deparse(substitute(resp)), data,
                       namsU[i], S, type, addCovars, placAdj, bnds[[namsU[i]]],
@@ -286,7 +285,7 @@ print.MCPMod <- function(x, digits=3, eps=1e-03, ...){
   cat("\n")
 
   cat("Estimated Dose Response Models:")
-  for(i in 1:length(x$mods)){
+  for(i in seq_along(x$mods)){
     cat("\n")
     cat(names(x$mods)[i], "model\n")
     cofList <- coef(x$mods[[i]], sep = TRUE)
@@ -336,7 +335,7 @@ print.summary.MCPMod <- function(x, ...){
   cat(rep("*", 39), "\n", sep="")
   cat("Mod part \n")
   cat(rep("*", 39), "\n", sep="")
-  for(i in 1:length(x$mods)){
+  for(i in seq_along(x$mods)){
     if(i > 1)
       cat("\n")
     if(length(x$mods) > 1)

R/MCTtest.R

@@ -326,7 +326,7 @@ MCTpval <- function(contMat, corMat, df, tStat,
                   one.sided = matrix(rep(tStat, each = nMod), nrow = nMod),
                   two.sided = matrix(rep(tStat, each = nMod), nrow = nMod))
   pVals <- numeric(nMod)
-  for(i in 1:nMod){
+  for(i in seq_len(nMod)){
     tmp <- 1 - mvtnorm::pmvt(lower[,i], upper[,i], df = df,
                     corr = corMat, algorithm = ctrl)
     pVals[i] <- tmp

R/Mods.R

@@ -392,7 +392,7 @@ TD <- function(object, Delta, TDtype = c("continuous", "discrete"),
     for(nam in names(object)){
       par <- object[[nam]]
       if(is.matrix(par)){
-        for(i in 1:nrow(par)){
+        for(i in seq_len(nrow(par))){
           td <- calcTD(nam, par[i,], Delta, TDtype, direction, doses, off, scal, nodes)
           modNams <- c(modNams, paste(nam, i, sep=""))
           tds <- c(tds, td)
@@ -501,7 +501,7 @@ ED <- function(object, p, EDtype = c("continuous", "discrete"),
     for(nam in names(object)){
       par <- object[[nam]]
       if(is.matrix(par)){
-        for(i in 1:nrow(par)){
+        for(i in seq_len(nrow(par))){
           ed <- calcED(nam, par[i,], p, maxD, EDtype, doses, off, scal, nodes)
           modNams <- c(modNams, paste(nam, i, sep=""))
           eds <- c(eds, ed)

R/Mods_helpers.R

@@ -49,9 +49,7 @@ getAddArgs <- function(addArgs, doses = NULL){
 }
 
 checkEntries <- function(modL, doses, fullMod){
-  biModels <- c("emax", "linlog", "linear", "quadratic",
-                "exponential", "logistic", "betaMod", "sigEmax",
-                "linInt")
+  biModels <- builtInMods
   checkNam <- function(nam){
     if(is.na(match(nam, biModels)))
       stop("Invalid model specified: ", nam)
@@ -131,13 +129,13 @@ fullMod <-  function(models, doses, placEff, maxEff, scal, off){
       nmod <- length(pars)
       if(nmod > 1){
         Pars <- matrix(ncol=3, nrow=nmod)
-        for(j in 1:length(pars)){
+        for(j in seq_along(pars)){
           tmp <- getLinPars(nm, doses, as.vector(pars[j]), placEff[z], maxEff[z])
           Pars[j,] <- tmp
           z <- z+1
         }
         colnames(Pars) <- names(tmp)
-        rownames(Pars) <- 1:length(pars)
+        rownames(Pars) <- seq_along(pars)
         i <- i+1
       } else {
         Pars <-  getLinPars(nm, doses, as.vector(pars), placEff[z], maxEff[z])
@@ -147,13 +145,13 @@ fullMod <-  function(models, doses, placEff, maxEff, scal, off){
     if(is.element(nm,c("logistic", "betaMod", "sigEmax"))){
       if(is.matrix(pars)){
         Pars <- matrix(ncol=4, nrow=nrow(pars))
-        for(j in 1:nrow(pars)){
+        for(j in seq_len(nrow(pars))){
           tmp <- getLinPars(nm, doses, as.vector(pars[j,]), placEff[z], maxEff[z])
           Pars[j,] <- tmp
           z <- z+1
         }
         colnames(Pars) <- names(tmp)
-        rownames(Pars) <- 1:nrow(pars)
+        rownames(Pars) <- seq_len(nrow(pars))
         i <- i+1
       } else {
         Pars <-  getLinPars(nm, doses, as.vector(pars), placEff[z], maxEff[z]); i <- i+1; z <- z+1
@@ -162,12 +160,12 @@ fullMod <-  function(models, doses, placEff, maxEff, scal, off){
     if(nm == "linInt"){
       if(is.matrix(pars)){
         Pars <- matrix(ncol=length(nodes), nrow=nrow(pars))
-        for(j in 1:nrow(pars)){
+        for(j in seq_len(nrow(pars))){
           Pars[j,] <-  getLinPars(nm, doses, as.vector(pars[j,]), placEff[z], maxEff[z])
           z <- z+1
         }
         colnames(Pars) <- paste("d", doses, sep="")
-        rownames(Pars) <- 1:nrow(pars)
+        rownames(Pars) <- seq_len(nrow(pars))
         i <- i+1
       } else {
         Pars <- getLinPars(nm, doses, as.vector(pars), placEff[z], maxEff[z]); i <- i+1; z <- z+1
@@ -250,7 +248,7 @@ plotModels <- function(models, nPoints = 200, superpose = FALSE,
                                 lattice::panel.superpose(x[ind], y[ind], subscripts[ind], 
                                                          groups, ...)
                                 if(plotTD){
-                                  for(z in 1:length(pdos)){
+                                  for(z in seq_along(pdos)){
                                     lattice::panel.lines(c(0, pdos[z]), c(yax[z], yax[z]),lty=2, col=2)
                                     lattice::panel.lines(c(pdos[z], pdos[z]), c(0, yax[z]),lty=2, col=2)
                                   }
@@ -384,7 +382,7 @@ calcTD <- function(model, pars, Delta, TDtype = c("continuous", "discrete"),
       inds <- cf < cf[1] + Delta
       if(all(inds))
         return(NA)
-      ind <- min((1:length(cf))[!inds])-1
+      ind <- min(seq_along(cf)[!inds])-1
       tmp <- (cf[1]+Delta-cf[ind])/(cf[ind+1]-cf[ind])
       td <- nodes[ind] + tmp*(nodes[ind+1]-nodes[ind])
       if(td > 0)
@@ -669,9 +667,9 @@ calcResp <- function(models, doses, off, scal, nodes){
         pars <- cbind(pars, off)
       if(nm == "betaMod")
         pars <- cbind(pars, scal)
-      ind <- 1:nmod
+      ind <- seq_len(nmod)
       nams <- c(nams, paste(nm, ind, sep = ""))
-      for(j in 1:nmod) {
+      for(j in seq_len(nmod)) {
         if(nm != "linInt"){
           val[[k]] <- do.call(nm, c(list(doses), as.list(pars[j,])))
         } else {

R/bFitMod.R

@@ -152,6 +152,8 @@ bFitMod <- function(dose, resp, model, S, placAdj = FALSE,
     nodes <- dose
 
   ## model number
+  ## NB: this ordering is coupled to the C code (modNr is passed to
+  ## bFitMod.Bayes) and therefore differs from builtInMods on purpose.
   builtIn <- c("linear", "linlog", "quadratic", "linInt", "emax",
                "logistic", "exponential", "sigEmax", "betaMod")
   modNr <- match(model, builtIn)
@@ -329,7 +331,7 @@ plot.bFitMod <- function (x, plotType = c("dr-curve", "effect-curve"),
         sdev <- sqrt(diag(attr(x, "data")$S))
         q <- qnorm(1 - (1 - level)/2)
         LBm <- UBm <- numeric(length(dose))
-        for (i in 1:length(dose)) {
+        for (i in seq_along(dose)) {
           LBm[i] <- resp[i] - q * sdev[i]
           UBm[i] <- resp[i] + q * sdev[i]
         }
@@ -350,7 +352,7 @@ plot.bFitMod <- function (x, plotType = c("dr-curve", "effect-curve"),
       sdev <- sqrt(diag(attr(x, "data")$S))
       q <- qnorm(1 - (1 - level)/2)
       LBm <- UBm <- numeric(length(dose))
-      for (i in 1:length(dose)) {
+      for (i in seq_along(dose)) {
         LBm[i] <- resp[i] - q * sdev[i]
         UBm[i] <- resp[i] + q * sdev[i]
       }
@@ -373,7 +375,7 @@ plot.bFitMod <- function (x, plotType = c("dr-curve", "effect-curve"),
       points(dose, resp, pch = 19, cex = 0.75)
     if (plotData == "meansCI") {
       points(dose, resp, pch = 19, cex = 0.75)
-      for (i in 1:length(dose)) {
+      for (i in seq_along(dose)) {
         lines(c(dose[i], dose[i]), c(LBm[i], UBm[i]), 
               lty = 2)
       }

R/bFitMod_helpers.R

@@ -2,7 +2,7 @@
 
 checkPrior <- function(prior){
   z <- 1
-  for(z in 1:length(prior)){
+  for(z in seq_along(prior)){
     prvec <- prior[[z]]
     nam <- names(prior)[z]
     if(!all(is.numeric(prvec)))
@@ -32,7 +32,7 @@ checkPrior <- function(prior){
 
 getPrBnds <- function(prior){
   prbnds <- matrix(ncol = 2, nrow = length(prior))
-  for(z in 1:length(prior)){
+  for(z in seq_along(prior)){
     prvec <- prior[[z]]
     nam <- names(prior)[z]
     if(nam %in% c("norm", "t"))
@@ -104,11 +104,11 @@ bFitMod.Bayes <- function(dose, resp, S, model, placAdj,
                    placAdj = placAdj, addArgs=list(off = off, scal = scal))
     if(is.null(start)){
       start <- coef(gfit)
-      for(i in 1:length(start)){
+      for(i in seq_along(start)){
         start[i] <- projPrBnds(start[i], prBnds[i,1], prBnds[i,2])
       }
     } else {
-      for(i in 1:length(start)){
+      for(i in seq_along(start)){
         if((start[i] < prBnds[i,1]) | (start[i] > prBnds[i,2]))
           stop("specified start value not consistent with bounds on prior distribution")
       }

R/bMCTtest.R

@@ -242,14 +242,14 @@ print.bMCTtest <- function(x, digits = 3, eps = 1e-3, ...){
   print(round(x$contMat, digits))
   cat("\n","Posterior Mixture Weights:","\n",sep="")
   w <- round(unlist(x$posterior[[1]]), digits = digits)
-  names(w) <- paste("Comp.", 1:length(w))
+  names(w) <- paste("Comp.", seq_along(w))
   print(w)
   ord <- rev(order(attr(x$tStat, "pVal")))
   pval <- format.pval(attr(x$tStat, "pVal"),
                       digits = digits, eps = eps)
   dfrm <- data.frame(round(x$tStat, digits)[ord, , drop = FALSE],
                      pval[ord])
-  names(dfrm) <- c(paste0("Comp. ", 1:ncol(x$tStat)), "posterior probability")
+  names(dfrm) <- c(paste0("Comp. ", seq_len(ncol(x$tStat))), "posterior probability")
   cat("\n","Bayesian t-statistics:","\n",sep="")
   print(dfrm)
   if(!is.null(x$critVal)){
@@ -304,15 +304,15 @@ mvpostmix <- function(priormix, mu_hat, S_hat)
   ## prior predictive distributions are MVN distribution with mean vectors equal to the prior components' mean vectors 
   ## and covariance matrices which are the sum of the prior components' covariance matrices and the "known" covariance 
   ## matrix of the data (for which S_hat is plugged in here)
-  for(i in 1:length(lw)){
+  for(i in seq_along(lw)){
     lw[i] <- log(priormix[[1]][[i]]) + mvtnorm::dmvnorm(mu_hat, priormix[[2]][[i]], SigmaPred[[i]], log = TRUE)
     postmix[[2]][[i]] <- solve(priorPrec[[i]] + dataPrec) %*% (priorPrec[[i]] %*% priormix[[2]][[i]] + dataPrec %*% mu_hat)
     postmix[[3]][[i]] <- solve(priorPrec[[i]] + dataPrec)
   }
   postmix[[1]] <- as.list(exp(lw - logSumExp(lw)))
 
   for(i in 1:3)
-    names(postmix[[i]]) <- paste0("Comp", 1:length(lw))
+    names(postmix[[i]]) <- paste0("Comp", seq_along(lw))
 
    postmix
 }

R/drmodels.R

@@ -1,4 +1,29 @@
 
+## internal metadata and helpers shared across the package
+
+## built-in dose-response model names, ordered so that the first four
+## entries are the linear models. This ordering is relied upon in fitMod,
+## maFitMod and MCPMod, where models 1-4 are treated as linear models.
+## Note: bFitMod uses a *different* ordering that is coupled to the C code.
+builtInMods <- c("linlog", "linear", "quadratic", "linInt", "emax",
+                 "exponential", "logistic", "betaMod", "sigEmax")
+
+## number of parameters for each built-in model (linInt depends on the
+## data and is handled separately)
+nParMod <- c(linlog = 2, linear = 2, quadratic = 3, emax = 3,
+             exponential = 3, logistic = 4, betaMod = 4, sigEmax = 4)
+
+## log(x) with the convention log(0) := 0 (vectorised)
+lg2 <- function(x){
+  l <- x
+  l[x == 0] <- 0
+  l[x != 0] <- log(x[x != 0])
+  l
+}
+
+## x*log(x) with the convention 0*log(0) := 0 (scalar x only)
+xlogx <- function(x) if(x == 0) 0 else x * log(x)
+
 ## model functions
 #' @rdname drmodels
 #' @param dose Dose variable
@@ -36,7 +61,6 @@ sigEmax <- function(dose, e0, eMax, ed50, h){
 #' @usage NULL
 #' @export 
 sigEmaxGrad <- function(dose, eMax, ed50, h, ...){
-  lg2 <- function(x) {l<-x; l[x==0] <- 0; l[x!=0] <- log(x[x!=0]); l}
   a <-  1 / (1 + (dose/ed50)^h)
   g1 <- 1 / (1 + (ed50/dose)^h)
   g2 <- -(h * eMax / ed50) * g1 * a
@@ -86,7 +110,6 @@ quadraticGrad <- function(dose, ...){
 #' @usage NULL
 #' @export 
 betaMod <- function(dose, e0, eMax, delta1, delta2, scal){
-  xlogx <- function(x) if(x == 0) 0 else x * log(x) # will not be called with vector x
   logMaxDens <- xlogx(delta1) + xlogx(delta2) - xlogx(delta1 + delta2)
   dose <- dose/scal
   e0 + eMax/exp(logMaxDens) * (dose^delta1) * (1 - dose)^delta2
@@ -96,8 +119,6 @@ betaMod <- function(dose, e0, eMax, delta1, delta2, scal){
 #' @usage NULL
 #' @export 
 betaModGrad <- function(dose, eMax, delta1, delta2, scal, ...){
-  lg2 <- function(x) {l<-x; l[x==0] <- 0; l[x!=0] <- log(x[x!=0]); l}
-  xlogx <- function(x) if(x == 0) 0 else x * log(x) # will not be called with vector x
   dose <- dose/scal
   if(any(dose > 1)) {
     stop("doses cannot be larger than scal in betaModel")

R/fitMod.R

@@ -210,8 +210,7 @@ fitMod <- function(dose, resp, data = NULL, model = NULL, S = NULL,
       df <- Inf
   }
   ## check whether model has been specified correctly
-  builtIn <- c("linlog", "linear", "quadratic", "linInt", "emax",
-               "exponential", "logistic", "betaMod", "sigEmax")
+  builtIn <- builtInMods
   if(missing(model))
     stop("Need to specify the model that should be fitted")
   modelNum <- match(model, builtIn)
@@ -617,7 +616,7 @@ predict.DRMod <- function(object, predType = c("full-model", "ls-means", "effect
     } else { ## calculate st. error (no need to calculate full covMat here)
       covMat <- vcov(object)
       if(addCovars != ~1) ## remove columns corresponding to covariates
-        covMat <- covMat[1:length(DRpars), 1:length(DRpars)]
+        covMat <- covMat[seq_along(DRpars), seq_along(DRpars)]
       if(!attr(object, "placAdj")){ ## remove intercept from cov-matrix
         if(model != "linInt"){
           covMat <- covMat[-1,-1]

R/fitMod_helpers.R

@@ -64,8 +64,7 @@ fitMod.raw <- function(dose, resp, data, model, S, type,
   ## differences to fitMod:
   ## - dose, resp need to be vectors containing the data
   ## - additional args: doseNam, respNam, off, scal
-  builtIn <- c("linlog", "linear", "quadratic", "linInt", "emax",
-               "exponential", "logistic", "betaMod", "sigEmax")
+  builtIn <- builtInMods
   modelNum <- match(model, builtIn)
 
   weights <- NULL;clinS <- NULL
@@ -509,7 +508,7 @@ plotFunc <- function(x, CI = FALSE, level = 0.95,
   if(inherits(x, "MCPMod")){
     nmods <- length(x$mods)
     lst <- vector(mode = "list", nmods)
-    for(i in 1:nmods){
+    for(i in seq_len(nmods)){
       pred <- predict(x$mods[[i]], predType = predtype, doseSeq = doseSeq, se.fit = CI)
       lbnd <- ubnd <- rep(NA, length(doseSeq))
       if(CI){
@@ -562,7 +561,7 @@ plotFunc <- function(x, CI = FALSE, level = 0.95,
                  lattice::lpoints(pList$dos, pList$mns, pch=19, col = colMn)
                  if(plotData == "meansCI"){
                    quant <- qnorm(1 - (1 - level)/2)
-                   for(i in 1:length(pList$dos)){
+                   for(i in seq_along(pList$dos)){
                      lattice::llines(rep(pList$dos[i], 2),
                             c(pList$lbndm[i], pList$ubndm[i]),
                             lty=2, col = colMn, ...)

R/guesst.R

@@ -190,7 +190,7 @@ guesst <- function(d, p, model = c("emax", "exponential", "logistic", "quadratic
       x <- log(d)
       par <- coef(lm(y~x))
       names(par) <- NULL
-      res <- c(ed50 = exp(par[1]/-par[2]), delta = -par[2])
+      res <- c(ed50 = exp(par[1]/-par[2]), h = -par[2])
     }
     if(local) {
       fooSE <- function(par, d, p, Maxd) {