diff --git a/R/MCPMod.R b/R/MCPMod.R
index 4c19de4..dab255a 100644
--- a/R/MCPMod.R
+++ b/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)
diff --git a/R/MCTtest.R b/R/MCTtest.R
index 7dfadfe..f6aa2de 100644
--- a/R/MCTtest.R
+++ b/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
diff --git a/R/Mods.R b/R/Mods.R
index 196958e..4ee3001 100644
--- a/R/Mods.R
+++ b/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)
diff --git a/R/Mods_helpers.R b/R/Mods_helpers.R
index 5f6f7be..ed20e14 100644
--- a/R/Mods_helpers.R
+++ b/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 {
diff --git a/R/bFitMod.R b/R/bFitMod.R
index 2887647..9dea8f3 100644
--- a/R/bFitMod.R
+++ b/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)
       }
diff --git a/R/bFitMod_helpers.R b/R/bFitMod_helpers.R
index 58c9d07..8f19c68 100644
--- a/R/bFitMod_helpers.R
+++ b/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")
       }
diff --git a/R/bMCTtest.R b/R/bMCTtest.R
index 8197755..93aba24 100644
--- a/R/bMCTtest.R
+++ b/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,7 +304,7 @@ 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)
@@ -312,7 +312,7 @@ mvpostmix <- function(priormix, mu_hat, S_hat)
   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
 }
diff --git a/R/drmodels.R b/R/drmodels.R
index 0e8bde3..50b225f 100644
--- a/R/drmodels.R
+++ b/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")
diff --git a/R/fitMod.R b/R/fitMod.R
index 8cfec4e..84778e1 100644
--- a/R/fitMod.R
+++ b/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]
diff --git a/R/fitMod_helpers.R b/R/fitMod_helpers.R
index c5507c8..33891cb 100644
--- a/R/fitMod_helpers.R
+++ b/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, ...)
diff --git a/R/guesst.R b/R/guesst.R
index e25225f..00c5356 100644
--- a/R/guesst.R
+++ b/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) {
diff --git a/R/maFitMod.R b/R/maFitMod.R
index 54f6792..6bb8ce1 100644
--- a/R/maFitMod.R
+++ b/R/maFitMod.R
@@ -51,8 +51,7 @@ maFitMod <- function(dose, resp, S, models,
                      nSim = 1000,
                      control, bnds, addArgs = NULL){
 
-  builtIn <- c("linlog", "linear", "quadratic", "linInt", "emax",
-               "exponential", "logistic", "betaMod", "sigEmax")
+  builtIn <- builtInMods
   if(missing(models))
     stop("Need to specify the models that should be fitted")
   modelNum <- match(models, builtIn)
@@ -83,7 +82,7 @@ maFitMod <- function(dose, resp, S, models,
   sims <- mvtnorm::rmvnorm(nSim, resp, S)
   fits <- vector("list", nSim)
   selModel <- character(nSim)
-  for(i in 1:nSim){
+  for(i in seq_len(nSim)){
     mod_fits <- lapply(models, function(mod){
       fitMod(dose, sims[i,], model = mod, S = S,
              type = "general", bnds = bnds[[mod]],
@@ -121,8 +120,8 @@ predict.maFit <- function(object,
   nSim <- length(object$selModel)
   pred <- matrix(nrow = nSim, ncol = length(doseSeq))
   colnames(pred) <- doseSeq
-  rownames(pred) <- 1:nSim
-  for(i in 1:nSim){
+  rownames(pred) <- seq_len(nSim)
+  for(i in seq_len(nSim)){
     pred[i,] <- predict(object$fits[[i]], doseSeq = doseSeq, predType = "ls-means")
   }
   if(!is.null(summaryFct)){
@@ -199,7 +198,7 @@ plot.maFit <- function(x,
     sdev <- sqrt(diag(x$args$S))
     crit <- qnorm(1 - tail_prob)
     LBm <- UBm <- numeric(length(x$args$dose))
-    for (i in 1:length(x$args$dose)) {
+    for (i in seq_along(x$args$dose)) {
       LBm[i] <- trafo(x$args$resp[i] - crit * sdev[i])
       UBm[i] <- trafo(x$args$resp[i] + crit * sdev[i])
     }
diff --git a/R/optDesign.R b/R/optDesign.R
index a19c692..20a1433 100644
--- a/R/optDesign.R
+++ b/R/optDesign.R
@@ -455,7 +455,7 @@ calcCrit <- function(design, models, probs, doses,
   intdesignCrit <- match(designCrit, c("TD", "Dopt", "Dopt&TD"))
   res <- numeric(nrow(design))
   ## check for sufficient number of design points
-  iter <- 1:nrow(design)
+  iter <- seq_len(nrow(design))
   design0 <- sweep(design, 2, nold, "+")
   count <- apply(design0, 1, function(x) sum(x > 0.0001))
   ind <- count < max(p[probs > 0])
@@ -557,7 +557,7 @@ plot.DRdesign <- function(x, models, lwdDes = 10, colDes = rgb(0,0,0,0.3), ...){
       miny <- min(args$y)
       maxy <- max(args$y)
       dy <- maxy-miny
-      for(k in 1:length(x$doses)){
+      for(k in seq_along(x$doses)){
         yy <- c(0,(x$design*dy)[k])+miny
         xx <- rep(x$doses[k],2)
         lattice::panel.xyplot(xx, yy, type="l", col = colDes, lwd = lwdDes)
diff --git a/R/optDesign_helpers.R b/R/optDesign_helpers.R
index 72cf3dc..6e5d852 100644
--- a/R/optDesign_helpers.R
+++ b/R/optDesign_helpers.R
@@ -9,7 +9,7 @@ calcGrads <- function(fmodels, doses, weights,
   for(nam in names(fmodels)){
     pars <- fmodels[[nam]]
     if(is.matrix(pars)){
-      for(i in 1:nrow(pars)){
+      for(i in seq_len(nrow(pars))){
         modgrad[[z]] <- t(gradCalc(nam, pars[i,], doses, off=off, scal=scal)*sqrt(weights))
         if(designCrit != "Dopt")
           TDgrad[[z]] <- calcTDgrad(nam, pars[i,], Delta, direction, off, scal)
@@ -34,14 +34,12 @@ calcGrads <- function(fmodels, doses, weights,
 
 ## returns the number of parameters (needed for C call)
 nPars <- function(mods){
-  builtIn <- c("linlog", "linear", "quadratic", 
-               "emax", "exponential", "logistic", 
-               "betaMod", "sigEmax")
-  ind <- match(mods, builtIn)
+  ## linInt is deliberately excluded (not allowed in optDesign)
+  ind <- match(mods, names(nParMod))
   if(any(is.na(ind))){
     stop(mods[which(is.na(ind))], " model not allowed in optDesign")
   }
-  c(2,2,3,3,3,4,4,4)[ind]
+  unname(nParMod[ind])
 }
 
 ## function which calls different optimizers
diff --git a/R/planMod.R b/R/planMod.R
index 1573395..b7da3ee 100644
--- a/R/planMod.R
+++ b/R/planMod.R
@@ -154,7 +154,7 @@ planMod <- function(model, altModels, n, sigma, S, doses,
                    linInt = length(doses),
                    nPars(model))
     bestPar <- matrix(nrow = ncol(muMat), ncol = npar) ## best fit by model to models in altModels
-    for(i in 1:ncol(muMat)){
+    for(i in seq_len(ncol(muMat))){
       ## if other model-class approximate best fit
       nam <- gsub("[0-9]", "", nams[i]) # model name (number removed)
       if(nam == model){
@@ -205,7 +205,7 @@ planMod <- function(model, altModels, n, sigma, S, doses,
   if(simulation){
     cat("Running simulations\n")
     requireNamespace("parallel", quietly = TRUE)
-    sim <- parallel::mclapply(1:ncol(muMat), function(i){
+    sim <- parallel::mclapply(seq_len(ncol(muMat)), function(i){
       if(showSimProgress){
         if(cores == 1){
           cat(sprintf("Scenario %d/%d\n", i, ncol(muMat)))
@@ -219,7 +219,7 @@ planMod <- function(model, altModels, n, sigma, S, doses,
       mse <- LBmn <- edpred <- resp <- numeric(nSim)
       coefs <- vector("list", length = nSim)
       modelSel <- character(nSim)
-      for(j in 1:nSim){
+      for(j in seq_len(nSim)){
         if(showSimProgress & cores == 1)
           setTxtProgressBar(pb, j/nSim)
         fit <- vector("list", length = length(model))
@@ -279,7 +279,7 @@ planMod <- function(model, altModels, n, sigma, S, doses,
     coefs <- lapply(sim, function(x) attr(x, "coefs"))
     modelSel <- sapply(sim, function(x) attr(x, "model"))
     names(NAind) <- colnames(modelSel) <- names(coefs) <- nams
-    rownames(modelSel) <- 1:nSim
+    rownames(modelSel) <- seq_len(nSim)
     sim <- do.call("rbind", sim)
     colnames(sim) <- c("dRMSE", "Pow(maxDose)", "P(EDp)")
     rownames(sim) <- nams
@@ -390,7 +390,7 @@ print.summary.planMod <- function(x, digits = 3, len = 101,
   if(!is.null(p)){
     eds <- getSimEst(x, "ED", p=p)
   }
-  for(i in 1:ncol(muMat)){
+  for(i in seq_len(ncol(muMat))){
     out[i,1] <- mseANOVA/x$sim[i,1]^2
     ## calculate mse of estimating the plac-adj dose-response at fine grid
     ## first calculate placebo-adjusted predictions
diff --git a/R/planMod_helpers.R b/R/planMod_helpers.R
index 924f4d3..2221162 100644
--- a/R/planMod_helpers.R
+++ b/R/planMod_helpers.R
@@ -53,7 +53,7 @@ aprCov <- function(doses, model, cf, S, off, scal){
 tableMatch <- function(x, match){
   ## like "table", but also returns categories with 0 counts
   out <- numeric(length(match))
-  for(i in 1:length(match)){
+  for(i in seq_along(match)){
     out[i] <- sum(x == match[i], na.rm=TRUE)
   }
   names(out) <- match
@@ -81,13 +81,13 @@ getSimEst <- function(x, type = c("dose-response", "ED", "TD"),
       stop("\"Delta\" needs to be > 0")
   }
   out <- vector("list", nAlt)
-  for(i in 1:nAlt){
+  for(i in seq_len(nAlt)){
     ind <- matrix(ncol = length(model), nrow = nSim)
     if(type == "dose-response"){
       resMat <- matrix(nrow = nSim, ncol = length(doseSeq))
       colnames(resMat) <- doseSeq
-      rownames(resMat) <- 1:nSim
-      for(j in 1:length(model)){
+      rownames(resMat) <- seq_len(nSim)
+      for(j in seq_along(model)){
         ind[,j] <- modelSel[,i] == model[j]
         if(any(ind[,j])){
           cf <- do.call("rbind", (coefs[[i]])[ind[,j]])
@@ -101,7 +101,7 @@ getSimEst <- function(x, type = c("dose-response", "ED", "TD"),
     }
     if(is.element(type, c("TD", "ED"))){
       resVec <- numeric(nSim)
-      for(j in 1:length(model)){
+      for(j in seq_along(model)){
         ind[,j] <- modelSel[,i] == model[j]
         if(any(ind[,j])){
           cf <- do.call("rbind", (coefs[[i]])[ind[,j]])
@@ -139,7 +139,7 @@ plotDoseSims <- function(x, type = c("ED", "TD"), p, Delta, xlab){
   }
   ## write plotting data frame
   nams <- names(out)
-  group <- factor(rep(1:length(nams), each=length(out[[1]])), labels=nams)
+  group <- factor(rep(seq_along(nams), each=length(out[[1]])), labels=nams)
   pdat <- data.frame(est = do.call("c", out),
                      group = group)
   ## determine limits for x-axis
diff --git a/R/powMCT_helpers.R b/R/powMCT_helpers.R
index 7f90f99..c94c5cb 100644
--- a/R/powMCT_helpers.R
+++ b/R/powMCT_helpers.R
@@ -19,7 +19,7 @@ powCalc <- function(alternative, critV, df, corMat, deltaMat, control){
   ctrl$interval <- NULL      # not used with pmvt
   nScen <- ncol(deltaMat)
   res <- numeric(nScen)
-  for(i in 1:nScen){
+  for(i in seq_len(nScen)){
     pmvtCall <- c(list(lower, upper, df = df, corr = corMat, delta = deltaMat[,i],
                        algorithm = ctrl))
     res[i] <- as.vector(1 - do.call(mvtnorm::pmvt, pmvtCall))
@@ -124,7 +124,7 @@ powMCTBinCount <- function(n, doses, candModList = NULL, respModList,
   resp <- getResp(resp_mods)
   nMod <- ncol(resp)
   pow <- numeric(nMod)
-  for(i in 1:nMod){
+  for(i in seq_len(nMod)){
     mu_vec <- resp[,i, drop=FALSE] # column i contains true response vector
     ## calculate covariance matrix
     v <- getVarBinCount(mu_vec, type, theta)
diff --git a/R/sampSize.R b/R/sampSize.R
index dc274be..b85cc2f 100644
--- a/R/sampSize.R
+++ b/R/sampSize.R
@@ -308,7 +308,7 @@ targN <- function(upperN, lowerN, step, targFunc,
     colnames(out) <- ""
   }
   out2 <- out
-  for(i in 1:length(sumFct)){
+  for(i in seq_along(sumFct)){
     out2 <- cbind(out2, apply(out, 1, sumFct[i]))
   }
   dimnames(out2) <- list(nseq, c(colnames(out), sumFct))