diff --git a/.travis.yml b/.travis.yml index bfb7745..3b7beea 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,7 +9,7 @@ r: - devel r_github_packages: - - DoseResponse/drcData +# - DoseResponse/drcData - r-lib/covr after_success: diff --git a/DESCRIPTION b/DESCRIPTION index c61fc23..9fa82f7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,12 @@ Package: drc Version: 3.2-0 -Date: 2021-01-13 +Date: 2023-04-25 Title: Analysis of Dose-Response Data Authors@R: c(person("Christian" , "Ritz", role = c("aut", "cre"), email="ritz@bioassay.dk"), - person("Jens", "Streibig", middle="C.", email="streibig@bioassay.dk", role = "aut")) + person("Jens", "Streibig", middle="C.", email="streibig@bioassay.dk", role = "aut"), + person("Ryan", "Ward", role = "ctb")) Maintainer: Christian Ritz -Depends: R (>= 2.0.0), MASS, stats, drcData +Depends: R (>= 2.0.0), MASS, stats Imports: car, gtools, multcomp, plotrix, sandwich, scales LazyLoad: yes LazyData: yes @@ -14,3 +15,4 @@ License: GPL-2 URL: http://www.r-project.org, http://www.bioassay.dk RoxygenNote: 6.1.1 BugReports: https://github.com/DoseResponse/drc/issues/ + diff --git a/NAMESPACE b/NAMESPACE index b841b48..3386207 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -39,7 +39,7 @@ importFrom(scales, alpha) #importFrom(stats, cooks.distance, hatvalues) #exportMethods(BIC) -import(drcData) +# import(drcData) importFrom("graphics", "abline", "axTicks", "axis", "lines", "par", "plot", "points", "polygon", "segments", "text") diff --git a/R/drmOpt.R b/R/drmOpt.R index a8bb0d2..85a6563 100644 --- a/R/drmOpt.R +++ b/R/drmOpt.R @@ -1,3 +1,12 @@ +generateOptimControl <- function(maxIt, relTol, parscale, traceVal, method) { + if (method == "L-BFGS-B") { + control <- list(maxit = maxIt, factr = relTol, parscale = parscale, trace = traceVal) + } else { + control <- list(maxit = maxIt, reltol = relTol, parscale = parscale, trace = traceVal) + } + return(control) +} + "drmOpt" <- function(opfct, opdfct1, startVec, optMethod, constrained, warnVal, upperLimits, lowerLimits, errorMessage, maxIt, relTol, opdfct2 = NULL, parmVec, traceVal, silentVal = TRUE, @@ -21,7 +30,7 @@ matchCall) { nlsObj <- try(optim(startVec, opfct, opdfct1, hessian = hes, method = "L-BFGS-B", lower = lowerLimits, upper = upperLimits, - control = list(maxit = maxIt, reltol = relTol, parscale = psVec)), silent = silentVal) + control = generateOptimControl(maxIt, relTol, psVec, traceVal, "L-BFGS-B")), silent = silentVal) } else { nlsObj <- try(optim(startVec, opfct, opdfct1, hessian = hes, method = optMethod, control = list(maxit = maxIt, reltol = relTol, parscale = psVec)), silent = silentVal) @@ -54,7 +63,7 @@ matchCall) # print(opfct(startVec)) nlsObj <- try(optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits, upper = upperLimits, - control = list(maxit = maxIt, parscale = psVec, reltol = relTol, trace = traceVal)), silent = silentVal) + control = generateOptimControl(maxIt, relTol, psVec, traceVal, "L-BFGS-B")), silent = silentVal) # parscale is needed for the example in methionine.Rd } else { # psVec <- abs(startVec) diff --git a/R/predict.drc.R b/R/predict.drc.R index d9de777..b856794 100644 --- a/R/predict.drc.R +++ b/R/predict.drc.R @@ -1,18 +1,18 @@ -"predict.drc" <- function(object, newdata, se.fit = FALSE, - interval = c("none", "confidence", "prediction", "ssd"), - level = 0.95, na.action = na.pass, od = FALSE, vcov. = vcov, +"predict.drc" <- function(object, newdata, se.fit = FALSE, + interval = c("none", "confidence", "prediction", "ssd"), + level = 0.95, na.action = na.pass, od = FALSE, vcov. = vcov, ssdSEfct = NULL, constrain = TRUE, checkND = TRUE, ...) { ## Checking arguments interval <- match.arg(interval) respType <- object[["type"]] - dataList <- object[["dataList"]] + dataList <- object[["dataList"]] doseDim <- ncol(dataList[["dose"]]) - if (is.null(doseDim)) {doseDim <- 1} + if (is.null(doseDim)) {doseDim <- 1} ## Assigning dataset from object if no data frame is provided - if (missing(newdata)) + if (missing(newdata)) { # predValues <- fitted(object) # not used # newdata <- data.frame(object$data[, 1], object$data[, 3]) @@ -25,8 +25,8 @@ groupLevels <- as.character(dataList[["plotid"]]) } else { groupLevels <- as.character(dataList[["curveid"]]) - } -# + } +# # if (identical(respType, "event")) # { # newdata <- data.frame(dataList[["dose"]], dataList[["plotid"]]) @@ -34,13 +34,13 @@ # newdata <- data.frame(dataList[["dose"]], dataList[["curveid"]]) # } } else { - + if (checkND) - { + { dName <- dataList[["names"]][["dName"]] if (any(names(newdata) %in% dName)) { - doseVec <- newdata[, dName] + doseVec <- newdata[, dName] } else { doseVec <- newdata[, 1] # warning("Dose variable not in 'newdata'") @@ -48,33 +48,33 @@ } else { doseVec <- newdata } - + cName <- dataList[["names"]][["cNames"]] if (any(names(newdata) %in% cName)) { - groupLevels <- as.character(newdata[, cName]) - # as.character() removes factor encoding - + groupLevels <- as.character(newdata[, cName]) + # as.character() removes factor encoding + } else { groupLevels <- rep(1, nrow(newdata)) } -# -# +# +# # if (ncol(newdata) < (doseDim + 1)) {newdata <- data.frame(newdata, rep(1, nrow(newdata)))} # # ndncol <- ncol(newdata) # # doseVec <- newdata[, 1:(ndncol-1)] # doseVec <- newdata[, 1:doseDim] -# # groupLevels <- as.character(newdata[, ndncol]) # 'as.character()' used to suppress factor levels -# groupLevels <- as.character(newdata[, doseDim + 1]) # 'as.character()' used to suppress factor levels +# # groupLevels <- as.character(newdata[, ndncol]) # 'as.character()' used to suppress factor levels +# groupLevels <- as.character(newdata[, doseDim + 1]) # 'as.character()' used to suppress factor levels } noNewData <- length(groupLevels) - + # if (ncol(newdata) < 2) {newdata <- data.frame(newdata, rep(1, nrow(newdata)))} # if (ncol(newdata) > 2) {stop("More than 2 variables in 'newdata' argument")} - - ## Defining dose values -- dose in the first column! + + ## Defining dose values -- dose in the first column! # doseVec <- newdata[, 1] -# groupLevels <- as.character(newdata[, 2]) # 'as.character()' used to suppress factor levels +# groupLevels <- as.character(newdata[, 2]) # 'as.character()' used to suppress factor levels # noNewData <- length(doseVec) @@ -86,9 +86,9 @@ } ## Retrieving matrix of parameter estimates - parmMat <- object[["parmMat"]] + parmMat <- object[["parmMat"]] pm <- t(parmMat[, groupLevels, drop = FALSE]) - + # parmNames <- colnames(parmMat) # lenCN <- length(parmNames) # indVec <- 1:lenCN @@ -96,22 +96,22 @@ # if (lenCN > 1) # { # indVec <- indVec[as.character(newdata[, 2])] -# +# ## groupLevels <- newdata[, 2] # if (!all(is.numeric(groupLevels))) # { -## pm <- parmMat[, as.character(groupLevels)] # 'as.character()' used to suppress factor levels -# pm <- parmMat[, groupLevels] +## pm <- parmMat[, as.character(groupLevels)] # 'as.character()' used to suppress factor levels +# pm <- parmMat[, groupLevels] # } else { # pm <- parmMat[, groupLevels] # } # pm <- parmMat[, groupLevels] -# +# # } else { # lenDV <- length(doseVec) ## indVec <- rep(1, lenDV) # pm <- matrix(parmMat[, 1], length(parmMat[, 1]), lenDV) -# } +# } # ## Checking for NAs in matrix of parameter estimates # naVec <- rep(FALSE, lenCN) @@ -119,23 +119,23 @@ # { # naVec[i] <- any(is.na(parmMat[, i])) # } -# parmMat <- parmMat[, !naVec, drop = FALSE] +# parmMat <- parmMat[, !naVec, drop = FALSE] ## Retrieving variance-covariance matrix sumObj <- summary(object, od = od) -# varMat <- sumObj[["varMat"]] - vcovMat <- vcov.(object) +# varMat <- sumObj[["varMat"]] + vcovMat <- vcov.(object) ## Defining index matrix for parameter estimates indexMat <- object[["indexMat"]] - - ## Calculating predicted values -# indexVec <- as.vector(indVec) -# print(indexVec) -# lenIV <- length(indexVec) - - + + ## Calculating predicted values +# indexVec <- as.vector(indVec) +# print(indexVec) +# lenIV <- length(indexVec) + + # retMat <- matrix(0, lenIV, 4) retMat <- matrix(0, noNewData, 4) colnames(retMat) <- c("Prediction", "SE", "Lower", "Upper") @@ -144,25 +144,25 @@ # print(doseVec) retMat[, 1] <- objFct$"fct"(doseVec, pm) # print(pm) - + ## Checking if derivatives are available deriv1 <- objFct$"deriv1" if (is.null(deriv1)) { - return(retMat[, 1]) - } + return(retMat[, 1]) + } ## Calculating the quantile to be used in the confidence intervals if (!identical(interval, "none")) - { + { if (identical(respType, "continuous")) { - tquan <- qt(1 - (1 - level)/2, df.residual(object)) + tquan <- as.vector(qt(1 - (1 - level)/2, df.residual(object))) } else { - tquan <- qnorm(1 - (1 - level)/2) + tquan <- as.vector(qnorm(1 - (1 - level)/2)) } - } - + } + ## Calculating standard errors and/or confidence intervals if (se.fit || (!identical(interval, "none"))) { @@ -171,8 +171,8 @@ { estVec <- object[["dataList"]][["dose"]] seVec <- object[["dataList"]][["weights"]] - - if (is.null(ssdSEfct)) + + if (is.null(ssdSEfct)) { # lmObj <- lm(seVec ~ estVec) # linear, not great lmObj <- lm(log(seVec) ~ log(estVec)) @@ -186,24 +186,24 @@ # print(derivxRes) # if (is.finite(derivxRes)) # { - # sumObjRV <- (derivxRes * sePred)^2 + # sumObjRV <- (derivxRes * sePred)^2 # } else { - # sumObjRV <- 0 # setting Inf * 0 = 0 (would be NaN otherwise) - # } + # sumObjRV <- 0 # setting Inf * 0 = 0 (would be NaN otherwise) + # } sumObjRV <- rep(0, length(derivxRes)) - isFinDR <- is.finite(derivxRes) + isFinDR <- is.finite(derivxRes) sumObjRV[isFinDR] <- ((derivxRes * sePred)^2)[isFinDR] - + # print(sumObjRV) - } + } if (identical(interval, "prediction")) { sumObjRV <- rep(sumObj$"resVar", noNewData) - } + } #else { # sumObjRV <- 0 # } -# rowIndex <- 1 +# rowIndex <- 1 # for (i in indexVec) # for (i in 1:ncol(indexMat)) @@ -214,17 +214,17 @@ for (rowIndex in 1:noNewData) { # parmInd <- indexMat[, i] -# print(indexVec) +# print(indexVec) # print(varMat) -# print(parmInd) - +# print(parmInd) + # varCov <- varMat[parmInd, parmInd] # print(varCov) # groupLevels <- newdata[, 2] # parmInd <- indexMat[, groupLevels[rowIndex]] # varCov <- varMat[parmInd, parmInd] - parmInd <- piMat[, rowIndex] + parmInd <- piMat[, rowIndex] varCov <- vcovMat[parmInd, parmInd] # parmChosen <- t(parmMat[, i, drop = FALSE]) @@ -233,28 +233,28 @@ dfEval <- deriv1(doseVec[rowIndex], pm[rowIndex, , drop = FALSE]) varVal <- dfEval %*% varCov %*% dfEval - retMat[rowIndex, 2] <- sqrt(varVal) -# retMat[rowIndex, 2] <- sqrt(dfEval %*% varCov %*% dfEval) + retMat[rowIndex, 2] <- sqrt(varVal) +# retMat[rowIndex, 2] <- sqrt(dfEval %*% varCov %*% dfEval) if (!se.fit) { - #retMat[rowIndex, 3:4] <- rep(retMat[rowIndex, 1], 2) + + #retMat[rowIndex, 3:4] <- rep(retMat[rowIndex, 1], 2) + # (tquan * sqrt(varVal + sumObjRV[rowIndex])) * c(-1, 1) retMat[rowIndex, 3] <- retMat[rowIndex, 1] - tquan * sqrt(varVal + sumObjRV[rowIndex]) - retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal + sumObjRV[rowIndex]) - } + retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal + sumObjRV[rowIndex]) + } # if (identical(interval, "confidence")) # { # retMat[rowIndex, 3] <- retMat[rowIndex, 1] - tquan * sqrt(varVal) -# retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal) +# retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal) # } # if (identical(interval, "prediction")) # { # sumObjRV <- sumObj$"resVar" # retMat[rowIndex, 3] <- retMat[rowIndex, 1] - tquan * sqrt(varVal + sumObjRV) -# retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal + sumObjRV) -# } -# rowIndex <- rowIndex + 1 +# retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal + sumObjRV) +# } +# rowIndex <- rowIndex + 1 } } ## Imposing constraints on predicted values @@ -270,13 +270,11 @@ retMat[, 3] <- pmax(retMat[, 3], 0) } } - + ## Keeping relevant indices keepInd <- 1 if (se.fit) {keepInd <- c(keepInd, 2)} if (!identical(interval, "none")) {keepInd <- c(keepInd, 3, 4)} - + return(retMat[, keepInd]) # , drop = FALSE]) } - - diff --git a/README.Rmd b/README.Rmd index a051d23..8528a45 100644 --- a/README.Rmd +++ b/README.Rmd @@ -29,7 +29,7 @@ Analysis of dose-response data is made available through a suite of flexible and ## You can install drc from GitHub # install.packages("devtools") ## first installing drcData -devtools::install_github("DoseResponse/drcData") +# devtools::install_github("DoseResponse/drcData") ## then installing the development version of drc devtools::install_github("DoseResponse/drc") ``` diff --git a/README.md b/README.md index ebea6e5..ca0fd69 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ Installation ## You can install drc from GitHub # install.packages("devtools") ## first installing drcData -devtools::install_github("DoseResponse/drcData") +# devtools::install_github("DoseResponse/drcData") ## then installing the development version of drc -devtools::install_github("DoseResponse/drc") +devtools::install_github("ryandward/drc") ``` diff --git a/docs/index.html b/docs/index.html index c1eea05..12f4967 100644 --- a/docs/index.html +++ b/docs/index.html @@ -75,7 +75,7 @@

## You can install drc from GitHub
 # install.packages("devtools")
 ## first installing drcData
-devtools::install_github("DoseResponse/drcData")
+# devtools::install_github("DoseResponse/drcData")
 ## then installing the development version of drc
 devtools::install_github("DoseResponse/drc")