R/modified_drm.R

"drm" <- function(
formula, curveid, pmodels, weights, data = NULL, subset, fct, 
type = c("continuous", "binomial", "Poisson", "negbin1", "negbin2", "event", "ssd"), bcVal = NULL, bcAdd = 0, 
start, na.action = na.omit, robust = "mean", logDose = NULL, 
control = drmc(), lowerl = NULL, upperl = NULL, separate = FALSE,
pshifts = NULL, varcov = NULL)
{
#    ## Matching 'adjust' argument
#    adjust <- match.arg(adjust)

    ## Matching argument values
    type <- match.arg(type)
    
    ## Loading MASS
#    require(MASS, quietly = TRUE)  # used for boxcox and ginv

    ## Setting na.action option
    options(na.action = deparse(substitute(na.action)))

    ## Setting control parameters
    useD <- control$"useD"
    constrained <- control$"constr"
#    maxDose <- control$"maxDose"
    maxIt <- control$"maxIt"
    optMethod <- control$"method"
    relTol <- control$"relTol"
    warnVal <- control$"warnVal"
#    zeroTol <- control$"zeroTol"
#    bcConstant <- bcAdd  
    rmNA <- control$"rmNA"  # in drmEM...
    errorMessage <- control$"errorm"  # in drmOpt
    noMessage <- control$"noMessage"  # reporting finding control measurements? 
#    trace <- control$"trace"
#    otrace <- control$"otrace"
    dscaleThres <- control$"dscaleThres"
    rscaleThres <- control$"rscaleThres"
    conCheck <- control$"conCheck"
            
    ## Setting warnings policy
    options(warn = warnVal)

#    ## Setting adjustment
#    if (adjust == "none") {boxcox <- FALSE; varPower <- FALSE}
#    if (adjust == "bc1") {boxcox <- TRUE; varPower <- FALSE}
#    if (adjust == "vp") {boxcox <- FALSE; varPower <- TRUE}

#    if ( (!is.null(bcVal)) && (is.numeric(bcVal))) {boxcox <- bc}
        
#    if (!(robust == "mean"))
#    {
#        boxcox <- FALSE
#        varPower <- FALSE
#    }

    ## Handling 'start' argument
    if (missing(start)) {selfStart <- TRUE} else {selfStart <- FALSE}

    ## Handling 'fct' argument
    if ( (!is.list(fct)) && (!is.function(fct)) ) {stop("No function or list given in argument 'fct'")}
    if (is.function(fct)) 
    {
        fct <- fct2list(fct, 2)
    }
    
    ## Converting a user specified list
    if (is.null(names(fct))) {fct$"fct" <- fct[[1]]; fct$"ssfct" <- fct[[2]]; fct$"names" <- fct[[3]]}
    
    if (!is.function(fct$"fct")) 
    {
        stop("First entry in list to 'fct' NOT a function")
    } else {
        drcFct <- fct$"fct"
    }
    
    if (is.null(fct$"ssfct")) {noSSfct <- TRUE} else {noSSfct <- FALSE}
    if ((!is.function(fct$"ssfct")) && selfStart)
    {
        stop("Neither self starter function nor starting values provided")
    } else {
        ssfct <- fct$"ssfct"
    }
    
    if (is.null(fct$"names") || (!is.character(fct$"names"))) 
    {
        stop("Parameter names (as vector a strings) are NOT supplied")
    } else {
        parNames <- fct$"names" 
        numNames <- length(parNames)
    }
    
#    ## Coercing two arguments in 'ssfct' into one argument
#    lenASS <- length(formals(ssfct))
#    if (lenASS > 1)
#    {
#        stop("Self starter function should only have one argument, which takes a data frame")
##        ssTemp <- ssfct
##        ssfct <- function(dataset) {ssTemp(dataset[, head(1:lenASS, -1)], dataset[, lenASS])}
#    }

    ## Checking whether or not first derivates are supplied    
    isDF <- is.function(fct$"deriv1")
    if ( (useD) && (isDF) )
    {
        dfct1 <- fct$"deriv1"  # deriv1  # [[4]]
#        drcDer2 <- fct$deriv2  # [[5]]
    } else {
        dfct1 <- NULL
    }

    ## Checking whether or not second derivates are supplied    
    if ( (useD) && (is.function(fct$"deriv2")) )
    {
        dfct2 <- fct$"deriv2"
    } else {
        dfct2 <- NULL
    }    
#    fct$"anovaYes"$"bin" <- NULL
#    fct$"anovaYes"$"cont" <- TRUE

    ## Storing call details
    callDetail <- match.call()

    ## Handling the 'formula', 'curveid' and 'data' arguments
    anName <- deparse(substitute(curveid))  # storing name for later use
    if (length(anName) > 1) {anName <- anName[1]}  # to circumvent the behaviour of 'substitute' in do.call("multdrc", ...)
    if (nchar(anName) < 1) {anName <- "1"}  # in case only one curve is analysed


    mf <- match.call(expand.dots = FALSE)   
    nmf <- names(mf) 
    mnmf <- match(c("formula", "curveid", "data", "subset", "na.action", "weights"), nmf, 0) 

    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf[c(1,mnmf)], parent.frame())  #, globalenv())
    mt <- attr(mf, "terms")

    varNames <- names(mf)[c(2, 1)]  
    varNames0 <- names(mf) 
    # only used once, but mf is overwritten later on
    
    dose <- model.matrix(mt, mf)[,-c(1)]  # with no intercept
    xDim <- ncol(as.matrix(dose))
    resp <- model.response(mf, "numeric")
    if (is.null(resp)) 
    {
        if (xDim > 1) {doseForResp <- dose[, 1]} else {doseForResp <- dose}
        resp <- ppoints(doseForResp, 0.5)[order(doseForResp)]  # just one option
        varNames[1] <- varNames[2] 
        varNames[2] <- "proportion"
    }
    origDose <- dose
    origResp <- resp  # in case of transformation of the response    
    lenData <- length(resp)
    numObs <- length(resp)

#    xDim <- ncol(as.matrix(dose))

#    if (xDim > 1)
#    {
#        stop("drm() is only designed for 1-dim. dose vectors")
#    }

#    dimData <- xDim + 1  # dimension of dose plus 1 dimensional response
    
#    varNames <- names(mf)
#    varNames <- varNames[c(2:dimData,1)]

#    print(names(mf))
#    print(model.extract(mf, "weights"))
#    print(model.weights(mf))

    ## Retrieving weights
    wVec <- model.weights(mf)
    if (is.null(wVec))
    {
        wVec <- rep(1, numObs)
    }

#    ## Extracting variable for heterogeneous variances
#    vvar <- model.extract(mf, "hetvar")    
    
    ## Finding indices for missing values
    missingIndices <- attr(mf, "na.action")
    if (is.null(missingIndices)) {removeMI <- function(x){x}} else {removeMI <- function(x){x[-missingIndices,]}}

    ## Handling "curveid" argument
    assayNo <- model.extract(mf, "curveid") 
    if (is.null(assayNo))  # in case not supplied
    {
        assayNo <- rep(1, numObs)
    }
    uniqueNames <- unique(assayNo)
    colOrder <- order(uniqueNames)
#    print(colOrder)
    uniqueNames <- as.character(uniqueNames)

    ## Re-enumerating the levels in 'assayNo' and 'pmodels'
    assayNoOld <- assayNo
#    ciOrigIndex <- uniqueNames  # unique(assayNo)
#    ciOrigLength <- length(unique(assayNoOld))

    ## Detecting control measurements 
    
    ## Defining helper function     
    colConvert <- function(vec)
    {
        len <- length(vec)
        assLev <- unique(vec)

        retVec <- rep(0,len)
        j <- 1
        for (i in 1:length(assLev)) {retVec[vec == assLev[i]] <- j; j <- j + 1}

        return(retVec)
    } 
    assayNo <- colConvert(assayNoOld)
    assayNames <- as.character(unique(assayNoOld))
    numAss <- length(assayNames)
      
#    lenDose <- unlist(lapply(tapply(dose, assayNoOld, unique), length))
#    conDose <- names(lenDose)[lenDose == 1]
#    nconDose <- names(lenDose)[lenDose > 1]
#    if (length(conDose) > 0) 
#    {        
#        if (!noMessage) 
#        {
#            cat(paste("Control measurements detected for level: ", conDose, "\n", sep = ""))
#        }
#            
#        assayNo[assayNoOld %in% conDose] <- nconDose[1]
#        ciOrigIndex <- unique(assayNo)
##        ciOrigLength <- length(unique(assayNoOld))  # numAss
#      
#      
#        ## Updating names, number of curves and the enumeration (starting from 1)
#        assayNames <- nconDose
##        numAss <- length(assayNames) 
#        assayNo <- colConvert(assayNo)
#                       
#        cm <- NULL
##        }
#         
#        uniqueDose <- lapply(tapply(dose, assayNoOld, unique), length)
#        udNames <- names(uniqueDose[uniqueDose == 1])
#        if (length(udNames) > 0) 
#        {
#            cm <- udNames
#            if (!noMessage) {cat(paste("Control measurements detected for level: ", udNames, "\n", sep = ""))}
#            ## add a check to see if at least one component in pmodels results in a single column
            
##            conInd <- assayNoOld%in%udNames
##            assayNo[conInd] <- (assayNo[!conInd])[1]
##            cm <- NULL
##assayNew <- assayNo
##assayNew[conInd] <- (assayNo[!conInd])[1]
##print(assayNew)
##
#            conInd <- assayNoOld%in%udNames
#            assayNo[conInd] <- (assayNo[!conInd])[1]
#            ciOrigIndex <- unique(assayNo)
#            ciOrigLength <- numAss
#            
#            ## Updating names, number of curves and the enumeration (starting from 1)
#            assayNames <- as.character(unique(assayNoOld[!conInd]))
#            numAss <- length(assayNames) 
#            assayNo <- colConvert(assayNo)
#                       
#            cm <- NULL
#

## New -commented out
#    } else {
#        cm <- NULL
#        ciOrigIndex <- unique(assayNo)
##        ciOrigLength <- numAss 
#        
#        assayNames <- as.character(unique(assayNoOld))          
#        assayNo <- colConvert(assayNoOld)  # re-enumerating from 1 to numAss       
#    }
#    numAss <- length(assayNames)   
#    print(ciOrigIndex)
#    print(ciOrigLength)

#    print(xDim)
#    print(cbind(dose, assayNoOld))
    
    if (xDim > 1) {tempDoseVec <- dose[, 1]} else {tempDoseVec <- dose} 
#    print(tempDoseVec)
#    print(assayNoOld)
#    uniqueDose <- lapply(tapply(dose, assayNoOld, unique), length)
    uniqueDose <- lapply(tapply(tempDoseVec, assayNoOld, unique), length)
    udNames <- names(uniqueDose[uniqueDose == 1])
    if ( (conCheck) && (length(udNames) > 0) ) 
    {
        cm <- udNames
        if (!noMessage) 
        {
            cat(paste("Control measurements detected for level: ", udNames, "\n", sep = ""))
            
            if (separate)
            {
                stop("Having a common control when fitting separate models does not make sense!\n")
            }
        }
        conInd <- assayNoOld %in% udNames
        assayNo[conInd] <- (assayNo[!conInd])[1]
        ciOrigIndex <- unique(assayNo)
        ciOrigLength <- numAss
        
        ## Updating names, number of curves and the enumeration (starting from 1)
        assayNames <- as.character(unique(assayNoOld[!conInd]))
        numAss <- length(assayNames)
        assayNo <- colConvert(assayNo)
        cm <- NULL
    } else {
        cm <- NULL
        ciOrigIndex <- unique(assayNo)
        ciOrigLength <- numAss
    }
#    print(assayNo)
    
    ## Pooling data from different curves
    if ((separate) && (numAss < 2))
    {
#        warning("Nothing to pool", call. = FALSE)
        warning("Only one level: separate = TRUE has no effect", call. = FALSE)
        separate <- FALSE 
    }    
    if ((separate) && (!missing(pmodels)))
    {
        warning("Separate fitting switched off", call. = FALSE)
        separate <- FALSE
    }
    if (separate)
    {
#        return(idrm(dose, resp, assayNo, wVec, fct, type))
        return(idrm(dose, resp, assayNoOld, wVec, fct, type, control))
    }    
    
    ## Handling "pmodels" argument
    pmodelsList <- list()
    if (missing(pmodels)) 
    {
#        pmodels <- as.data.frame(matrix(assayNo, numObs, numNames))
#        
        if (length(unique(assayNo)) == 1) 
        {
            for (i in 1:numNames) 
            {
                pmodelsList[[i]] <- matrix(1, numObs, 1)
            }
        } else {
            modelMat <- model.matrix(~ factor(assayNo) - 1, level = unique(assayNo))  # no intercept term
            colnames(modelMat) <- assayNames
            for (i in 1:numNames) 
            {
                pmodelsList[[i]] <- modelMat
#                print(head(modelMat))
#                pmodelsList[[i]] <- pmodelsList[[i]][, colOrder]
            }
        }   
    } else {
        ## Handling a list or data.frame argument of "pmodels"
        if (is.null(data)) 
        {
            pmodels <- eval(substitute(pmodels), envir = .GlobalEnv)
        } else {
            pmodels <- eval(substitute(pmodels), envir = data, enclos = parent.frame())
        }
   
        if (is.data.frame(pmodels))
        {
            lenCol <- ncol(pmodels)
            pmodelsMat <- matrix(0, numObs, lenCol)    
    
            for (i in 1:lenCol) 
            {
                if (length(unique(pmodels[,i])) == 1) 
                {
                    pmodelsList[[i]] <- matrix(1, numObs, 1)
                    pmodelsMat[,i] <- rep(1, numObs)    
                } else {
                    mf <- eval(model.frame(~factor(pmodels[,i]) - 1), parent.frame())  # converting to factors
                    mt <- attr(mf, "terms")    
    
                    mf2 <- model.matrix(mt, mf)
                    ncmf2 <- ncol(mf2)

                    mf3 <- removeMI(mf2)
                    pmodelsList[[i]] <- mf3
                    pmodelsMat[, i] <- mf3 %*% c(1:ncmf2)
                }
            }
        } else {

            if (is.list(pmodels))
            {   
                lenCol <- length(pmodels)
                pmodelsMat <- matrix(0, length(resp), lenCol)
    
                for (i in 1:lenCol) 
                {
                    if (paste(as.character(pmodels[[i]]), collapse = "") == "~1") 
                    {
                        pmodelsList[[i]] <- matrix(1, numObs, 1)
                        pmodelsMat[,i] <- rep(1, numObs)
                    } else {
                        mf <- eval(model.frame(pmodels[[i]], data=data), parent.frame())   
                        mt <- attr(mf, "terms")    
                        
                        mf2 <- model.matrix(mt, mf)
                        ncmf2 <- ncol(mf2)

                        mf3 <- removeMI(mf2)                    
                        pmodelsList[[i]] <- mf3  
                    
                        pmodelsMat[,i] <- mf3%*%c(1:ncmf2)                    
                    }
                }
            }
        }     
#        pmodelsOld <- pmodels
#        pmodels <- as.data.frame(pmodelsMat)  # pmodelsMat not used any more
    }
#    for (i in 1:numNames) {pmodels[, i] <- colConvert(pmodels[, i])}

    
    ## Re-setting na.action
    options(na.action = "na.omit")  # the default

    ## Transforming dose value if they are provided as log dose
    if ( !is.null(logDose) && is.numeric(logDose) ) 
    {
       origDose <- dose
       dose <- logDose^dose
    }

#    ## Handling one-dimensional x     
#    if (xDim == 1)
#    {
        ## Defining ANOVA model
#        bcc <- rep(bcAdd, numObs)    
#        if (numAss > 1) 
#        {
#            anovaFormula <- (resp + bcc) ~ offset(bcc) + doseFactor*factor(assayNo)
#            alternative <- 2
#        } else {
#            anovaFormula <- (resp + bcc) ~ offset(bcc) + doseFactor
#            alternative <- 1
#        }
         

#        ## Checking whether there is enough df to perform Box-Cox transformation
#        if ( boxcox && ( (numObs - numAss*length(unique(dose))) < numObs/10) )
#        {
#            if (boxcox) {warning("Box-Cox transformation based on clustering of dose values", call. = FALSE)}
#            doseFactor <- factor(cutree(hclust(dist(dose), method = "average"), numObs/3))  
#            # constructing groups containing roughly 3 observations 
#
#            ## Re-defining ANOVA model
#            if (numAss > 1) 
#            {
#                anovaFormula <- (resp + bcc) ~ offset(bcc) + doseFactor*factor(assayNo)
#                dset <- data.frame(doseFactor, resp, assayNo, bcc)
#                alternative <- 2
#            } else {
#                anovaFormula <- (resp + bcc) ~ offset(bcc) + doseFactor
#                dset <- data.frame(doseFactor, resp, bcc)
#                alternative <- 1
#            }
#        } else {
#            doseFactor <- factor(dose)
#        }
         
#        ## Fitting ANOVA model
#        if (type == "continuous")
#        {
#            testList <- drmLOFls()
##            if (varPower) {testList <- drmLOFvp()}
##            if (!is.null(vvar)) {testList <- drmLOFhv()}
#        }
#        if (type == "binomial")
#        {
#            testList <- drmLOFbinomial()
#        }
#        if (type == "Poisson")
#        {
#            testList <- drmLOFPoisson()
#        }
#        
##        if (varPower) {testList <- mdrcVp(anovaYes = TRUE)} else {testList <- drmEMls(anovaYes = TRUE)}
##        if (!is.null(vvar)) {testList <- mdrcHetVar(anovaYes = TRUE)}
#
#        gofTest <- testList$"gofTest"            
#        lofTest <- testList$"anovaTest"
#        if (!is.null(lofTest))
#        {
#            dset <- data.frame(dose, factor(dose), resp, assayNo, bcc)
#            anovaModel0 <- lofTest(anovaFormula, dset)
#        } else {
#            anovaModel0 <- NULL
#            alternative <- 0
#        }
        
#        ## Fitting ANOVA model
#        testList <- switch(type, 
#        "continuous" = drmLOFls(), 
#        "binomial" = drmLOFbinomial(), 
#        "Poisson" = drmLOFPoisson())
#
#        gofTest <- testList$"gofTest"            
#        lofTest <- testList$"anovaTest"
#        if (!is.null(lofTest))
#        {
#            afList <- anovaFormula(dose, resp, assayNo, bcAdd)
#            anovaForm <- afList$"anovaFormula"
#            anovaData <- afList$"anovaData"        
#
#            anovaModel0 <- lofTest(anovaForm, anovaData)
#        } else {
#            anovaModel0 <- NULL
#        }


#        ## Applying the Box-Cox transformation (lambda is defined here!)
#        bcResult <- drmBoxcox(boxcox, anovaFormula, dset)  
#        lambda <- bcResult[[1]]
#        boxcoxci <- bcResult[[2]] 
#        boxcox <- bcResult[[3]]      

#        lambda <- 0
#        isNumeric <- is.numeric(boxcox)
#        if ( (isNumeric) || (is.logical(boxcox) && boxcox)  ) 
#        {
#            if (!isNumeric)
#            {
#                profLik <- boxcox(anovaFormula, lambda = seq(-2.6, 2.6, 1/10), plotit = FALSE, data = dset)  
#                # boxcox in MASS
#                
#                maxIndex <- which.max(profLik$y)
#                lambda <- (profLik$x)[maxIndex]
#                boxcoxci <- drmBoxcoxCI(profLik)
#            }
#            if (isNumeric)
#            {
#                lambda <- boxcox
#                boxcoxci <- c(NA, NA)                
#            }
#        } else {
#            lambda <- NA
#            boxcoxci <- c(NA, NA)
#        }
 

#        ## Using self starter 
#        if (!noSSfct)
#        {
#            ## Calculating initial estimates for the parameters using the self starter
#            startMat <- matrix(0, numAss, numNames)
#            doseresp <- data.frame(dose, origResp)
#    
#            for (i in 1:numAss)
#            {
#                indexT1 <- (assayNo == i)
#                if (any(indexT1)) 
#                {
#                    isfi <- is.finite(dose)  # removing infinite dose values
#                    logVec <- indexT1 & isfi
#                    startMat[i, ] <- ssfct(doseresp[logVec, ])  # ssfct(dose[logVec], origResp[logVec] )
#                } else {
#                    startMat[i, ] <- rep(NA, numNames)
#                }
#
#                ## Identifying a dose response curve only consisting of control measurements
#                if (sum(!is.na(startMat[i, ])) == 1) {upperPos <- (1:numNames)[!is.na(startMat[i, ])]}
#            }
##            colMat <- matrix(0, numNames, numAss)
##            maxParm <- rep(0, numNames)  # storing the max number of parameters
#        }

        
#    ## Handling multi-dimensional x   
#    } else { 
#        stop("Currently multi-dimensional dose values are not supported") 
#        alternative <- NULL
#        anovaModel0 <- NULL
#        anovaModel <- NULL
#        gofTest <- NULL
        
#        if (!is.null(bcVal))
#        {
#            lambda <- boxcox                  
#            boxcoxci <- c(NA, NA)
#        } else {
#            lambda <- NA                  
#            boxcoxci <- NULL
#        }       

#        ## Using self starter
#        if (!noSSfct)
#        {
#            ## Calculating initial estimates for the parameters using the self starter
#            startMat <- matrix(0, numAss, numNames)
#            doseresp <- data.frame(dose, origResp)
#    
#            for (i in 1:numAss)
#            {
#                indexT1 <- (assayNo == i)
#                if (any(indexT1)) 
#                {
#                    startMat[i, ] <- ssfct(doseresp[indexT1, ])  # ssfct(dose[indexT1], origResp[indexT1])
#                } else {
#                    startMat[i, ] <- rep(NA, numNames)
#                }
#                
#                ## Identifying a dose response curve only consisting of control measurements
#                if (sum(!is.na(startMat[i,]))==1) {upperPos <- (1:numNames)[!is.na(startMat[i,])]}
#            }
#            colMat <- matrix(0, numNames, numAss)
#            maxParm <- rep(0, numNames)  # storing the max number of parameters
#        }
#    }

    ## Finding parameters for the control measurements which will not be estimated
    pmodelsList2 <- list()
    for (i in 1:numNames)
    { 
        colNames <- colnames(pmodelsList[[i]])

        if ( (!is.null(cm)) && (!is.null(colNames)) ) 
        {
            accm <- as.character(cm)
            pos <- grep(accm, colNames)
            if (length(pos) == 0) 
            {
                candCol <- pmodelsList[[i]][, 1]
                if ( !(length(assayNoOld[candCol==1])==0) && (all(assayNoOld[candCol==1] == accm)) )
                {
                    pos <- 1  # the control measurements correspond to the "Intercept" term
                }
            }  
        } else {pos <- numeric(0)}


        ## Defining 'pmodelsList2' from 'pmodelsList'
        if ((length(pos) > 0) && !(upperPos == i) )
        {
            pmodelsList2[[i]] <- as.matrix(pmodelsList[[i]][, -pos])  # column is removed
        } else {
            pmodelsList2[[i]] <- as.matrix(pmodelsList[[i]])  # column is kept
        } 
    }
    
    for (i in 1:numNames)
    {
        if (ncol(pmodelsList[[i]]) > numAss) 
        {
            pmodelsList2[[i]] <- model.matrix(~factor(assayNo) - 1)
            colnames(pmodelsList2[[i]]) <- assayNames 
        } else {
            pmodelsList2[[i]] <- as.matrix(pmodelsList[[i]])  # columns are kept
        }
    }
    
    ## Constructing vectors 'ncclVec' and 'parmPos' used below
    ncclVec <- rep(0, numNames)
    for (i in 1:numNames)
    {
        ncclVec[i] <- ncol(pmodelsList2[[i]])  # ncol(as.matrix(pmodelsList2[[i]]))
    }
    parmPos <- c(0, cumsum(ncclVec)[-numNames])

    ## Constructing parameter names                       
    pnList <- drmParNames(numNames, parNames, pmodelsList2)
    parmVec <- pnList[[1]]
    parmVecA <- pnList[[2]]
    parmVecB <- pnList[[3]]   

    ## Defining with indices for the individual parameters in the model
    parmIndex <- list()
    for (i in 1:numNames)
    {
        parmIndex[[i]] <- parmPos[i] + 1:ncclVec[i]
    }

    ## Scaling of dose and response values 
    scaleFct <- fct$"scaleFct"
    if (!is.null(scaleFct))  # && (is.null(lowerl)) && (is.null(upperl)) )
    # currently the scaling interferes with constraining optimization
    {
        # Defining scaling for dose and response values 
        doseScaling <- 10^(floor(log10(median(dose))))   
#        if ( (is.na(doseScaling)) || (doseScaling < 1e-10) )  # changed May 16 2012
        if ( (is.na(doseScaling)) || (doseScaling < dscaleThres) )
        {
            doseScaling <- 1
        }

        respScaling <- 10^(floor(log10(median(resp)))) 
#        if ( (is.na(respScaling)) || (respScaling < 1e-10) || (!identical(type, "continuous")) || (!is.null(bcVal)) )  # changed May 16 2012
        if ( (is.na(respScaling)) || (respScaling < rscaleThres) || (!identical(type, "continuous")) || (!is.null(bcVal)) )
        {
            respScaling <- 1
        }   
#        print(resp)
#        print(median(resp))

#        doseScaling <- 1
#        respScaling <- 1

## Starting values need to be calculated after BC transformation!!!        

        # Retrieving scaling vector
        longScaleVec <- rep(scaleFct(doseScaling, respScaling), 
                            as.vector(unlist(lapply(parmIndex, length))))
        
    } else {
        doseScaling <- 1
        respScaling <- 1
        
        longScaleVec <- 1
#        
#        startVecSc <- startVec
    }      
#    print(c(doseScaling, respScaling, longScaleVec))    

    ## Constructing vector of initial parameter values
    startVecList <- list()
    
    ## Calculating initial estimates for the parameters using the self starter

#    if (identical(type, "ssd") & inherits(fct, "llogistic"))
#    {
#        ssFct <- function(dframe) {c(-1, 1, 1, 1, 1) * fct[["ssfct"]](dframe)}
#    }  
    
    if(!noSSfct)
    {
        startMat <- matrix(0, numAss, numNames)
        lenASS <- length(formals(ssfct))
        if (lenASS > 1)  
        # in case doseScaling and respScaling arguments are available
        # scaling is done inside ssfct()
        {
            doseresp <- data.frame(x = dose, y = origResp)
            ssFct <- function(dframe){ssfct(dframe, doseScaling, respScaling)}
        } else {
        # scaling is explicitly applied to the dose and response values
            doseresp <- data.frame(x = dose / doseScaling, y = origResp / respScaling)
            ssFct <- ssfct
        }
#        doseresp <- data.frame(x = dose / doseScaling, y = origResp / respScaling)
#        doseresp <- data.frame(dose, origResp)

# Not sure this indicator is needed?! Only used once below!
# Note is.finite() only works with vectors!
# Commented out 2010-12-13   
        isfi <- is.finite(dose)  # removing infinite dose values
        
        if (identical(type, "event"))
        {   
            # 9/4/19 - Andrea Onofri. It is necessary to take the rightmost
            # column as the response and not, rigidly, the 3rd. Furthermore,
            # it is necessary to allow for the possibility that, when
            # 'curveid' is used, there are replicates for each curveid level
            #dr2 <- doseresp[, 3]
            dr2 <- doseresp[, length(doseresp[1,])]
            isFinite <- is.finite(doseresp[, 2])
            respVec <- rep(NA, length(dr2))
            # 9/4/19 - Andrea Onofri.
            # The following command works well only without replicates.
            # With replicates, the cumulative proportion needs to be built
            # for each replicate (Petri Dish). The following code assumes
            # that data are ordered by Time and the last observation is 
            # always Inf per each dish
            #respVec[isFinite] <- cumsum(dr2[isFinite]) / sum(dr2)
            
            Dish <- c(); cont <- 1
            for(i in 1:length(dr2)) {Dish[i] <- cont; if(is.finite(doseresp[i,2]) == F ) cont <- cont+1}
            
            total <- ( tapply(dr2, list(as.factor(Dish)), sum ) )
            number <- ( tapply(dr2, list(as.factor(Dish)), length ) )
            divideBy <- rep(total, number)
            
            #divideBy <- rep(total, number)
            respVec[isFinite] <- unlist( tapply(dr2, Dish, cumsum ) )[isFinite]/divideBy[isFinite]
            #respVec <- unlist( tapply(dr2, Dish, cumsum ) )/divideBy
#            doseresp[, 3] <- cumsum(dr2[isFinite]) / sum(dr2)
##            doseresp[!is.finite(doseresp[, 2]), 1] <- NA              
#            doseresp <- doseresp[isFinite, c(1, 3)]
#            names(doseresp) <- c("x", "y")
            nCols <- c(2:(length(doseresp[1,])-1))
            doseresp <- (data.frame(x = doseresp[, nCols], y = respVec))[isFinite, ]
        } else { 
            isFinite <- is.finite(doseresp[, 2])
        }
        
        ## Finding starting values for each curve
        for (i in 1:numAss)
        {
            indexT1 <- (assayNo[isFinite] == i)
            if (any(indexT1)) 
            {
# Commented out 2010-12-13            
#                logVec <- indexT1 & isfi
                logVec <- indexT1
                                
#                startMat[i, ] <- ssfct(doseresp[logVec, ])  # ssfct(dose[logVec], origResp[logVec] )
#                startMat[i, ] <- ssfct(doseresp[logVec, ], doseScaling, respScaling)
                startMat[i, ] <- ssFct(doseresp[logVec, ])
            } else {
                 startMat[i, ] <- rep(NA, numNames)
            }
    
            ## Identifying a dose response curve only consisting of control measurements
            if (sum(!is.na(startMat[i, ])) == 1) 
            {
                upperPos <- (1:numNames)[!is.na(startMat[i, ])]
#                print(upperPos)
            }
        }
#        print(startMat)
          
#        startMat2 <- matrix(unlist(lapply(split(doseresp, assayNo[isFinite]), ssFct)), nrow = numAss, byrow = TRUE)
#        upperPos2 <- c(rep(1:numNames, numAss))[t(is.na(startMat))]
#        print(upperPos2)
#        print(startMat2)
  

#       New approach?
#        timeUsed <- 0 
#        ssFctWrapper <- function(dframeSubset)
#        {
#            ssFct(dframeSubset[is.finite(dframeSubset[1, ]), ])
#            timeUsed <- timeUsed + system.time(ssFct(dframeSubset[is.finite(dframeSubset[1, ]), ]))[3]
#        }  
#        startMat2 <- matrix(as.vector(unlist(lapply(split(doseresp, assayNo), ssFctWrapper))), 
#        numAss, numNames, byrow = TRUE)
#        print(startMat2)  # for comparison  
  
    
        ## Transforming matrix of starting values into a vector
        nrsm <- nrow(startMat)
        for (i in 1:numNames)
        {
            sv <- rep(0, max(nrsm, ncclVec[i]))
            indVec <- 1:ncclVec[i]
            sv[1:nrsm] <- startMat[, i]
            sv <- sv[!is.na(sv)]
            
            isZero <- (sv == 0)
            sv[isZero] <- mean(sv)
            
            startVecList[[i]] <- sv[indVec]
#            print(startVecList[[i]])
        }
        startVec <- unlist(startVecList)        
    } else {
        startVec <- start  # no checking if no self starter function is provided!!!
    }
    
    ## Checking the number of start values provided
    if (!selfStart && !noSSfct) 
    {
        lenReq <- length(startVec)  # generated from self starter
        if (length(start) == lenReq) 
        {
            startVec <- start / longScaleVec
        } else {
            stop(paste("Wrong number of initial parameter values. ", lenReq, " values should be supplied", sep = ""))
        }
    }

    ## Converting parameters
    if (selfStart)
    {
        startVec <- drmConvertParm(startVec, startMat, assayNo, pmodelsList2) 
    }

    # Scaling starting values (currently not done in drmEMls)
#    startVecSc <- startVec / longScaleVec 
    startVecSc <- startVec
#    print(startVecSc)    

    ## Defining function which converts parameter vector to parameter matrix            
    parmMatrix <- matrix(0, numObs, numNames)
    parm2mat <- function(parm)
    {
#        parmMatrix <- matrix(0, lenData, numNames)
        for (i in 1:numNames)
        {
#           print(as.matrix(pmodelsList2[[i]]))
#           print(parmPos[i] + 1:ncclVec[i])
#           print(parm[parmPos[i] + 1:ncclVec[i]])          
#           parmMatrix[, i] <- pmodelsList2[[i]] %*% parm[parmPos[i] + 1:ncclVec[i]]
           parmMatrix[, i] <- pmodelsList2[[i]] %*% parm[parmIndex[[i]]]
        }
        return(parmMatrix)
    }        

    ## Defining non-linear function
#    if (!is.null(fctList))
#    {
#        ivList <- list()
#        ivList2 <- list()        
#        matList <- list()
#        svList <- list()
#        for (i in 1:numAss)
#        {
#            indexT1 <- (assayNo == i)
#            isfi <- is.finite(dose)  # removing infinite dose values
#
#            ivList[[i]] <- indexT1
##            svList[[i]] <- fctList[[i]]$"ssfct"( doseresp[(indexT1 & isfi), ] )
#            logVec <- indexT1 & isfi
#            svList[[i]] <- fctList[[i]]$"ssfct"(doseresp[logVec, ])  # dose[logVec], origResp[logVec])
#            matList[[i]] <- c( sum(indexT1), length(svList[[i]]) )
#            
#            ivList2[[i]] <- match(fctList[[i]]$names, fct$names)             
#        }
#
#    
#        posVec <- rep(0, numAss)
#        for (i in 1:numAss)
#        {
#            posVec[i] <- matList[[i]][2]
#        }
#        posVec <- cumsum(posVec)
#        posVec <- c(0, posVec)
##        print(posVec)
#    
#        drcFct1 <- function(dose, parm)
#        {
#            retVec <- rep(0, numObs)
#            for (i in 1:numAss)
#            {
#                iVec <- ivList[[i]]
#                pMat <- matrix(parm[(posVec[i]+1):posVec[i+1]], matList[[i]][1], matList[[i]][2], byrow = TRUE) 
#                retVec[iVec] <- fctList[[i]]$"fct"( dose[iVec], pMat )
#            }
#            return(retVec)
#        }
#        
#        startVec <- as.vector(unlist(svList))
#    } else {

    ## Defining model function 
    multCurves <- modelFunction(dose, parm2mat, drcFct, cm, assayNoOld, upperPos, fct$"retFct", 
                                doseScaling, respScaling, isFinite = rep(TRUE, lenData), pshifts)

#    drcFct1 <- function(dose, parm)
#    {
#        drcFct(dose, parm2mat(parm))
#    }
##    }
#
#
#    ## Defining model function        
#    if (!is.null(fct$"retFct"))
#    {
#        drcFct <- fct$"retFct"(doseScaling, respScaling)  #, numObs)
#        drcFct1 <- function(dose, parm)
#        {
#            drcFct(dose, parm2mat(parm))
#        }
#    }
#
#    if (is.null(cm)) 
#    {
#        multCurves <- function(dose, parm)
#        {          
#           drcFct1(dose, parm)  # fctList           
#        }
#    } else {
#        iv <- assayNoOld == cm
#        niv <- !iv
#        fctEval <- rep(0, numObs)
#
#        multCurves <- function(dose, parm) 
#        {
#            parmVal <- parm2mat(parm)           
#            fctEval[iv] <- parmVal[iv, upperPos, drop = FALSE]
#            fctEval[niv] <- drcFct(dose[niv], parmVal[niv, , drop = FALSE])
#
#            fctEval
#        }
#    }
##    print(startVec)
##    print(multCurves(dose, startVec))

    
    ## Defining first derivative (if available) ... used once in drmEMls()
    if (!is.null(dfct1))
    {
        dmatfct <- function(dose, parm)
        {
            dfct1(dose, parm2mat(parm))
        }
    } else {
        dmatfct <-NULL
    } 

    ## Box-Cox transformation is applied
    if (!is.null(bcVal))  # (boxcox)
    {
#        varPower <- FALSE  # not both boxcox and varPower at the same time

        ## Defining Box-Cox transformation function
        bcfct <- function(x, lambda, bctol, add = bcAdd)
        {
            if (abs(lambda) > bctol)
            {
                return(((x + add)^lambda - 1)/lambda)
            } else {
                return(log(x + add))    
            }
        }
        
        ## Setting the tolerance for Box-Cox transformation being the logarithm transformation 
        ##  (same as in boxcox.default in MASS package)
        bcTol <- 0.02 
        
#        resp <- bcfct(resp, lambda, bcTol)
        resp <- bcfct(resp, bcVal, bcTol)

        multCurves2 <- function(dose, parm)
        {
            bcfct(multCurves(dose, parm), bcVal, bcTol)
        } 
    } else {multCurves2 <- multCurves}
#    print(startVec)
#    print(multCurves2(dose, startVec))    

    ## Defining estimation method -- perhaps working for continuous data
#    robustFct <- drmRobust(robust, match.call(), numObs, length(startVec))  
    robustFct <- drmRobust(robust, callDetail, numObs, length(startVec))  

    if (type == "continuous")
    {
        ## Ordinary least squares estimation
        estMethod <- drmEMls(dose, resp, multCurves2, startVecSc, robustFct, wVec, rmNA, dmf = dmatfct, 
        doseScaling = doseScaling, respScaling = respScaling, varcov = varcov)
        
#        if (adjust == "vp")  #(varPower)
#        {        
#            estMethod <- drmEMvp(dose, resp, multCurves2)  # mdrcVp(dose, resp, multCurves2)
#            lenStartVec <- length(startVec)
#            
#            start2ss <- estMethod$"ssfct"(cbind(dose, resp))
#            if (missing(start2))
#            {
#                startVec <- c(startVec, start2ss)            
#            } else {
#                if (length(start2) == 2)  # canonical 2?
#                {
#                    startVec <- c(startVec, start2)            
#                }
#            }
##            startVec <- c(startVec, estMethod$"ssfct"(cbind(dose, resp)))
#            parmVec <- c(parmVec, "Sigma", "Power")
#            
#            startVecSc <- startVec
#        }
                  
#        if (!is.null(vvar))
#        {
#            estMethod <- mdrcHetVar(dose, resp, multCurves2, vvar)
#            lenStartVec <- length(startVec)            
#            startVec <- c(startVec, estMethod$"ssfct"(cbind(dose, resp)))
#            parmVec <- c(parmVec, as.character(unique(vvar)))        
#        }
    }
    if (identical(type, "binomial"))
    {
        estMethod <- drmEMbinomial(dose, resp, multCurves2, startVecSc, robustFct, wVec, rmNA, 
        doseScaling = doseScaling)        
    } 
    if (identical(type, "Poisson"))
    {
        estMethod <- drmEMPoisson(dose, resp, multCurves2, startVecSc, weightsVec = wVec, 
                                  doseScaling = doseScaling)
    }
    if (identical(type, "negbin1") || identical(type, "negbin2"))
    {
        estMethod <- drmEMnegbin(dose, resp, multCurves2, startVecSc, weightsVec = wVec, 
                                 doseScaling = doseScaling, 
                                 dist.type = ifelse(type == "negbin1", 1, 2))
    }    
    if (identical(type, "event"))
    {
        estMethod <- drmEMeventtime(dose, resp, multCurves2, doseScaling = doseScaling)
    }
    if (identical(type, "ssd"))
    {
#        if (is.null(fct[["retFctDx"]])) {fct[["retFct"]] <- NULL} else {fct[["retFct"]] <- fct[["retFctDx"]]}
#        fct[["retFct"]] <- NULL
#        print(doseScaling)
#        print(respScaling)
        doseScaling <- 1  # dose is the response!
        respScaling <- 1  # no response variable
        longScaleVec <- rep(1, length(longScaleVec))
        multCurves2loc <- modelFunction(dose, parm2mat, fct$"derivx", cm, assayNoOld, upperPos, 
                                        retFct = fct[["retFctDx"]],  #NULL,
                                        doseScaling = doseScaling, respScaling = respScaling, 
                                        isFinite = rep(TRUE, lenData), pshifts)   
        estMethod <- drmEMssd(dose, resp, multCurves2loc, doseScaling = doseScaling, multCurves2 = multCurves2)
    }  
#    if (identical(type, "standard"))
#    {
#        estMethod <- drmEMstandard(dose, resp, multCurves2, doseScaling = doseScaling)
#    }      
            
    opfct <- estMethod$opfct            


    ## Re-fitting the ANOVA model to incorporate Box-Cox transformation (if necessary)
#    if (type == "continuous")
#    {    
#        if (!is.na(lambda))
#        {
#            dset <- data.frame(dose, doseFactor, resp, assayNo, bcc)  # dataset with new resp values        
#            anovaModel0 <- (testList$"anovaTest")(anovaFormula, dset)            
##            anovaModel <- anovaModel0$"anovaFit"
#        }
#    }

    ## Defining lower and upper limits of parameters
#    if (constrained)
#    {
    if (!is.null(lowerl)) 
    {
        if (!is.numeric(lowerl) || !((length(lowerl) == sum(ncclVec)) || (length(lowerl) == numNames)))
        {
            stop("Not correct 'lowerl' argument")
        } else {
            if (length(lowerl) == numNames) 
            {
                lowerLimits <- rep(lowerl, ncclVec)
            } else {
                lowerLimits <- lowerl
            }
        }
        constrained <- TRUE 
        
    } else {  ## In case lower limits are not specified
        lowerLimits <- rep(-Inf, length(startVec))
    }

    if (!is.null(upperl)) 
    {
        if (!is.numeric(upperl) || !((length(upperl) == sum(ncclVec)) || (length(upperl) == numNames)))
        {
            stop("Not correct 'upperl' argument")
        } else {
            if (length(upperl) == numNames) 
            {
                upperLimits <- rep(upperl, ncclVec)
            } else {
                upperLimits <- upperl
            }
        } 
        constrained <- TRUE
                
    } else {  ## In case upper limits are not specified
        upperLimits <- rep(Inf, length(startVec))
    }
    
    lowerLimits <- lowerLimits  / longScaleVec
    upperLimits <- upperLimits  / longScaleVec
        
#    if (all(!is.finite(lowerLimits)) && all(!is.finite(upperLimits))) 
#    {
#        stop("No constraints are imposed via 'lowerl' and 'upperl' arguments")
#    }
#    }

    ## Optimisation
    
    ## Setting derivatives
    opdfctTemp <- estMethod$"opdfct1"
    appFct <- function(x, y){tapply(x, y, sum)}   
    
    if (!is.null(opdfctTemp))
    {
        opdfct1 <- function(parm)
        {
#            print(as.vector(apply(opdfctTemp(parm), 2, appFct, assayNo)))
            as.vector(apply(opdfctTemp(parm), 2, appFct, assayNo))
        }
    } else {
        opdfct1 <- NULL
    }  

    ## Manipulating before optimisation
        
#    ## Scaling x values
#if (FALSE)
#{
#    sxInd <- fct$"sxInd"
#    sxYN <- !is.null(sxInd) && ((max(dose)<1e-2) || (min(dose)>1e2) || (diff(range(dose))>1e2) )
#    if ( sxYN && (is.null(fctList)) )
#    {
##        if (!is.null(fctList))
##        {
##            parmIndX <- rep(0, numAss)
##            for (i in 1:numAss)
##            {
##                parmIndX[i] <- fctList[[i]]$"sxInd"
##            }
##            parmIndX <- cumsum(parmIndX)
##        } else {
#            parmIndX <- parmPos[sxInd] + 1:ncclVec[sxInd]        
##        }
#    
#        scaleXConstant <- median(dose)
#        sxFct <- scaleX(scaleXConstant)  # , scaleX(dose, maxDose)        
#        if (adjust == "vp")
#        {
#            dose <- sxFct(dose)    
#            opfct <- drmEMvp(dose, resp, multCurves2)$"opfct"
#        }
#        
#        startVec[parmIndX] <- sxFct(startVec[parmIndX])
#    }
##    print(startVec)  # 2
#}

#    ## Scaling y values
#    ##  based on the original response value
#    ##  not the transformed values
#    syInd <- fct$"syInd"
#    lensy <- length(syInd)
#    parmIndY <- list()
#    
#    lyLim <- 1e-2
#    uyLim <- 1e2
#    syYN <- !is.null(syInd) && ((max(origResp)<lyLim) || (min(origResp)>uyLim) || (diff(range(origResp))>uyLim)) 
#    if ( syYN && (is.null(fctList)) )
#    {
##        if (!is.null(fctList))
##        {
##            parmIndY <- rep(0, numAss)
##            for (i in 1:numAss)
##            {
##                parmIndY[[i]] <- fctList[[i]]$"syInd"
##            }
##            parmIndY <- cumsum(as.vector(unlist(parmIndY)))
##        } else {
#            for (i in 1:lensy)
#            {
#                parmIndY[[i]] <- parmPos[syInd[i]] + c(1:ncclVec[syInd[i]])
#            }
#            tempPIY <- as.vector(unlist(parmIndY))
#            parmIndY <- tempPIY
##        }
#        if (adjust == "bc1")
#        { 
#            scaleYConstant <- bcfct(median(origResp), lambda, bcTol)  # median(origResp)
#        } else {
#            scaleYConstant <- median(origResp)  
#        }        
#        syFct <- scaleY(median(origResp))  # scaleY(scaleYConstant)
#        startVec[parmIndY] <- syFct(startVec[parmIndY])
#    }
#    # scaling of y values through 'opfct' definition
##    print(startVec)  # 3


    ## Testing nonlinear function
#    print(startVecSc)
#    print(multCurves2(dose, startVecSc))
#    print(opfct(startVecSc))
##    print(dose)
##    print(resp)
     
    ## Scaling objective function
#    if (type == "continuous")
#    {
#        ofVal <- opfct(startVec)
#        if ( !is.nan(ofVal) && ( (ofVal < 1e-2) || (ofVal >1e2) ) )
#        {
#            opfct2 <- function(c){opfct(c)/opfct(startVec)}
#        } else {
#            opfct2 <- opfct
#        }
#    } else {
#        opfct2 <- opfct
#    }
#    opfct2 <- opfct  # only used once below
    

    ## Updating starting values
    startVecSc <- as.vector(startVecSc)  # removing names
    if (identical(type, "negbin1") || identical(type, "negbin2"))
    {
      startVecSc <- c(startVecSc, 1)
      parmVec <- c(parmVec, "O:(Intercept)")
      parmVecA <- c(parmVecA, "O")
      parmVecB <- c(parmVecB, "(Intercept)")
      
#      print(opfct(startVecSc))  
    }  
    

    ## Optimising the objective function previously defined
    nlsFit <- drmOpt(opfct, opdfct1, startVecSc, optMethod, constrained, warnVal, 
    upperLimits, lowerLimits, errorMessage, maxIt, relTol, parmVec = parmVec, traceVal = control$"trace",
    matchCall = callDetail, silentVal = control$"otrace") 
#    matchCall = match.call()) 
        
    if (!nlsFit$convergence) {return(nlsFit)}
    
    ## Manipulating after optimisation   
    if (identical(type, "negbin1") || identical(type, "negbin2"))
    {
      longScaleVec <- c(longScaleVec, 1)
        
    }  
    
    if (identical(type, "event"))
    {
    #9/4/19 - Correction by Andrea #####
#        dose <- dose[isFinite, 2]
#        resp <- (as.vector(unlist(tapply(resp, assayNo, function(x){cumsum(x) / sum(x)}))))[isFinite]

#        orderDose <- order(dose0)
#        dose1 <- dose0[orderDose]

      
      assayNo0 <- assayNo[isFinite]
      
        #dose0 <- dose[, 2]
        dose0 <- dose[, -1] #Original dose
        
        #Ok se dose0 is vector
        #Also Create an id for experimental units (Petri dishes or other)
        if(is.vector(dose0)==T){
          dose <- dose0[is.finite(dose0)==T]
          dose1 <- dose0[is.finite(dose0)==T] 
          dose <- as.vector(unlist(tapply(dose1, assayNo0, function(x){unique(sort(x))})))
          Dish <- c(); cont <- 1
          for(i in 1:length(dose0)) {Dish[i] <- cont; if(is.finite(dose0[i]) == F ) cont <- cont+1}
        }else{
          dose <- dose0[is.finite(dose0[,1])==T,]
          dose1 <- dose0[is.finite(dose0[,1])==T,] 
          dose <- dose[order(dose[,1]),]
          Dish <- c(); cont <- 1
          for(i in 1:length(dose0[,1])) {Dish[i] <- cont; if(is.finite(dose0[i,1]) == F ) cont <- cont+1}
        }
        #dose <- dose[is.finite(dose[,2]),] 
        
        #dose0: original dose; dose1: only data without Inf; dose: dose for means
        #dose <- as.vector(unlist(tapply(dose1[,1], assayNo0, function(x){unique(sort(x))})))
        
        ## Rescaling per curve id
        idList <- split(data.frame(dose0, resp), assayNo) #For assay
        idList2 <- split(data.frame(dose0, resp), Dish) #For Dish
        respFct <- function(idListElt)
        {
            doseVec <- idListElt[, 1]
            dose2 <- unique(sort(doseVec))
            orderDose <- order(doseVec)
            resp1 <- tapply(idListElt[orderDose, 2], doseVec[orderDose], sum)  # obtaining one count per time interval
            resp2 <- cumsum(resp1) / sum(resp1)
            
            cbind(dose2, resp2)[is.finite(dose2), , drop = FALSE]
        }
#        drList <- lapply(idList, respFct)
#        lapList <- lapply(drList, function(x){x[, 1]})
#        dose <- as.vector(unlist(lapList))
#        resp <- as.vector(unlist(lapply(drList, function(x){x[, 2]})))
        
        #These functions here do not work properly with replicates.
        drList <- lapply(idList, respFct) #dose/resp per assay
        drList2 <- lapply(idList2, respFct) #dose/resp per dish
        lapList <- lapply(drList, function(x){x[, 1]}) 
         
        #dose <- as.vector(unlist(lapList))
        #Here resp and resp2 represent respectively the predictions 
        # for the means (if there are replicates per each AssayNo) 
        # and the individual responses
        resp <- as.vector(unlist(lapply(drList, function(x){x[, 2]}))) #This are the means
        resp2 <- as.vector(unlist(lapply(drList2, function(x){x[, 2]}))) #This are the data
        
#        listCI <- split(assayNoOld[isFinite], assayNoOld[isFinite])
#        splitFactor <- factor(assayNoOld[isFinite], exclude = NULL)
        splitFactor <- factor(assayNo, exclude = NULL)        
        listCI <- split(splitFactor, splitFactor)
        lenVec <- as.vector(unlist(lapply(lapList, length)))
#        print(listCI)
#        print(lenVec)
        plotid <- as.factor(as.vector(unlist(mapply(function(x,y){x[1:y]}, listCI, lenVec))))
        if(is.vector(dose0)==T){
          plotid2 <- as.factor(as.vector(unlist(listCI))[is.finite(dose0)])
        }else{
          plotid2 <- as.factor(as.vector(unlist(listCI))[is.finite(dose0[,1])])
        }
#        plotid <- plotid[complete.cases(plotid)]
        levels(plotid) <- unique(assayNoOld)
    } else {
        plotid <- NULL
    }

    
    if (identical(type, "ssd"))
    {
        if (ncol(as.matrix(dose)) > 1)
        {
            dose2 <- dose[, 2]
            ifDose2 <- is.finite(dose2)
            dose <- dose2[ifDose2]
            resp <- resp[ifDose2]
        } else {
            ifDose2 <- is.finite(dose)  # in case of no censoring  
          
        }
    }
    
    

#    print(nlsFit$par)    
    ## Adjusting for pre-fit scaling 
    if (!is.null(scaleFct))
    {
        # Scaling the sums of squares value back
        nlsFit$value <- nlsFit$value * (respScaling^2)
   
        # Scaling estimates and Hessian back
#        print(longScaleVec)
        nlsFit$par <- nlsFit$par * longScaleVec
        nlsFit$hessian <- nlsFit$hessian * (1/outer(longScaleVec/respScaling, longScaleVec/respScaling))
    }
    
    if (!is.null(fct$"retFct"))  #&& (!identical(type, "ssd")))
    {
        drcFct <- fct$"retFct"(1, 1)  #, numObs)  # resetting the scaling
        drcFct1 <- function(dose, parm)
        {
            drcFct(dose, parm2mat(parm)[isFinite, , drop = FALSE])
        }
    }
    
    
#    print(nlsFit$par)
#    nlsFit$value <- opfct(nlsFit$par)  # used in the residual variance

    ## Manipulating after optimisation
    
#    ## Adjusting for scaling of y values
#    if ( syYN && (is.null(fctList)) )
#    {
#        nlsFit$value <- syFct(syFct(nlsFit$value, down = FALSE), down = FALSE)
#        startVec[parmIndY] <- syFct(startVec[parmIndY], down = FALSE)
#        nlsFit$par[parmIndY] <- syFct(nlsFit$par[parmIndY], down = FALSE)
#
#        scaleFct1 <- function(hessian) 
#                     {
#                         newHessian <- hessian
#                         newHessian[, parmIndY] <- syFct(newHessian[, parmIndY], down = FALSE)
#                         newHessian[parmIndY, ] <- syFct(newHessian[parmIndY, ], down = FALSE)
#                         return(newHessian)
#                     }                
#    } else {
#        scaleFct1 <- function(x) {x}    
#    }

    
#    ## Adjusting for scaling of x values
#if (FALSE)
#{    
#    if ( sxYN && (is.null(fctList)) )  # (!is.null(sxInd))
#    {
#        if (adjust == "vp")
#        {
#            dose <- sxFct(dose, down = FALSE)
#        }
#        startVec[parmIndX] <- sxFct(startVec[parmIndX], down = FALSE)
#        nlsFit$par[parmIndX] <- sxFct(nlsFit$par[parmIndX], down = FALSE)
#
#        scaleFct2 <- function(hessian) 
#                     {
#                         newHessian <- scaleFct1(hessian)
#                         newHessian[, parmIndX] <- sxFct(newHessian[, parmIndX], down = FALSE)
#                         newHessian[parmIndX, ] <- sxFct(newHessian[parmIndX, ], down = FALSE)
#                         return(newHessian)
#                     }                
#    } else {
#        scaleFct2 <- function(hessian) 
#        {
#            scaleFct1(hessian)
#        }
#    }
#}

#    ## Handling variance parameters
#    varParm <- NULL
#        
#    if (varPower)
#    {
#        varParm <- list(type = "varPower", index = 1:lenStartVec)        
#    }
#    if (!is.null(vvar))
#    {
#        varParm <- list(type = "hetvar", index = 1:lenStartVec)
#    }

    # Testing against the ANOVA (F-test)
    nlsSS <- nlsFit$value
    nlsDF <- numObs - length(startVec)

    ## Constructing a plot function
        
    ## Picking parameter estimates for each curve. Does only work for factors not changing within a curve!
    if (!is.null(cm)) {iVec <- (1:numAss)[!(uniqueNames==cm)]} else {iVec <- 1:numAss}
    
    pickCurve <- rep(0, length(iVec))
    for (i in iVec)
    {
       pickCurve[i] <- (1:numObs)[assayNo == i][1]
    }
    parmMat <- matrix(NA, numAss, numNames)

    fixedParm <- (estMethod$"parmfct")(nlsFit)
#    print(nlsFit$par)
#    print(fixedParm)
    parmMat[iVec, ] <- (parm2mat(fixedParm))[pickCurve, ]

    indexMat2 <- parm2mat(1:length(fixedParm))
    indexMat2 <- indexMat2[!duplicated(indexMat2), ]

#    if(!is.null(fctList))
#    {
#         parmMat <- matrix(NA, numAss, numNames)
#         for (i in 1:numAss)
#         {
#             parmMat[i, ivList2[[i]]] <- fixedParm[(posVec[i]+1):posVec[i+1]]
#         }
#    }    
        
    if (!is.null(cm))
    {
#        conPos <- upperPos
#        print(conPos)
        parmMat[-iVec, upperPos] <- (parm2mat(fixedParm))[assayNoOld == cm, , drop = FALSE][1, upperPos]  
        # 1: simply picking the first row
    }
    rownames(parmMat) <- assayNames


    pmFct <- function(fixedParm)
    {
        if (!is.null(cm)) {iVec <- (1:numAss)[!(uniqueNames == cm)]} else {iVec <- 1:numAss}
    
        if (!is.null(cm))
        {
#            conPos <- conList$"pos"
            parmMat[-iVec, upperPos] <- (parm2mat(fixedParm))[assayNoOld == cm, , drop = FALSE][1, upperPos]  
            # 1: simply picking the first row
        }
        rownames(parmMat) <- assayNames
    
        return(parmMat)
    }
    parmMat <- pmFct(fixedParm)  # (estMethod$"parmfct")(nlsFit) )
#    print(parmMat)
#    print(pmFct(1:length(fixedParm)))

#    ## Scaling parameters
#    if (!is.null(fct$scaleFct))
#    {
#        scaleFct <- function(parm)
#        {
#            fct$scaleFct(parm, xScaling, yScaling)
#        }
#    
#        parmMat <- apply(parmMat, 1, scaleFct)
#    }
#

    ## Constructing design matrix allowing calculations for each curve
#    colPos <- 1
#    rowPos <- 1
#    Xmat <- matrix(0, numAss*numNames, length(nlsFit$par))
#    Xmat <- matrix(0, numAss*numNames, length(fixedParm))


#    if (!is.null(fctList)) {omitList <- list()}
#    for (i in 1:numNames)
#    {
#        indVec <- iVec
#        lenIV <- length(indVec)
#
#        nccl <- ncol(pmodelsList2[[i]])  # min(maxParm[i], ncol(pmodelsList2[[i]]))        
#
#        XmatPart <- matrix(0, lenIV, nccl)
#        k <- 1
#        if (!is.null(fctList)) {omitVec <- rep(TRUE, lenIV)}
#        for (j in indVec)
#        {
#            if (!is.null(fctList))
#            {
#                parPresent <- !is.na(match(i, ivList2[[j]]))
#                omitVec[k] <- parPresent
#            }
#            
#            XmatPart[k, ] <- (pmodelsList2[[i]])[(1:lenData)[assayNo == j][1], 1:nccl]
#            k <- k + 1
#        }
#        if (!is.null(fctList))
#        {
#            XmatPart <- XmatPart[omitVec, , drop = FALSE]
#            nccl <- nccl - sum(!omitVec)
#            omitList[[i]] <- omitVec
#        }
#
#        Xmat[rowPos:(rowPos+lenIV-1), colPos:(colPos+nccl-1)] <- XmatPart
#        colPos <- colPos + nccl
#        rowPos <- rowPos + lenIV
#    }
#    Xmat <- Xmat[1:(rowPos-1), 1:(colPos-1)]


    ## Defining the plot function    
    pfFct <- function(parmMat)
    {
        plotFct <- function(dose)
        {
#            if (xDim == 1) {lenPts <- length(dose)} else {lenPts <- nrow(dose)}
            if (is.vector(dose)) 
            {
                lenPts <- length(dose)
            } else {
                lenPts <- nrow(dose)
            }
#            print(lenPts)
#            print(ciOrigLength)

            curvePts <- matrix(NA, lenPts, ciOrigLength)  # numAss)
            for (i in 1:numAss)
            {
#                if (!is.null(fctList)) 
#                {
#                    drcFct <- fctList[[i]]$"fct"
#                    numNames <- matList[[i]][2]    
#                }
            
                if (i %in% iVec)
                {
#                    parmChosen <- parmMat[i, ]
                    parmChosen <- parmMat[i, complete.cases(parmMat[i, ])]  # removing NAs 
#                    print(parmChosen)                   
                    
                    parmMat2 <- matrix(parmChosen, lenPts, numNames, byrow = TRUE)
#                    print(parmMat2)
                    curvePts[, ciOrigIndex[i]] <- drcFct(dose, parmMat2)
                } else { curvePts[, i] <- rep(NA, lenPts)}
            }
            return(curvePts)
        }
    
        return(plotFct)    
    }
#    print(parmMat)
    plotFct <- pfFct(parmMat)
#    plotFct(0:10)


    ## Computation of fitted values and residuals
    ## 9/4/2019 - Andrea Onofri. The original routine did
    ## not appear to work with type = "event". Therefore I
    ## parted the two routines
    if (identical(type, "event"))
    {   dose <- dose1
        multCurves2 <- modelFunction(dose, parm2mat, drcFct, cm, assayNoOld, upperPos, fct$"retFct", doseScaling, respScaling, isFinite)
        #predVec <- multCurves2(dose[,-1], fixedParm)
        predVec <- multCurves2(dose, fixedParm)
        resVec <- resp2 - predVec # Corrected here #############################
        resVec[is.nan(predVec)] <- 0    
        #diagMat <- matrix(c(predVec, resVec), length(dose[,2]), 2)
        diagMat <- matrix(c(predVec, resVec), length(dose), 2)
        colnames(diagMat) <- c("Predicted values", "Residuals")
        #print(diagMat)
        }
    else{
    #Everything else, but time-to-event models    
    if (identical(type, "ssd"))
    {
        multCurves2 <- modelFunction(dose, parm2mat, drcFct, cm, assayNoOld, upperPos, fct$"retFct", doseScaling, respScaling, isFinite)
    }
    
    predVec <- multCurves2(dose, fixedParm)        
    resVec <- resp - predVec
    resVec[is.nan(predVec)] <- 0    

    diagMat <- matrix(c(predVec, resVec), length(dose), 2)
    colnames(diagMat) <- c("Predicted values", "Residuals")
    }

    ## Adjusting for robust estimation: MAD based on residuals, centered at 0, is used as scale estimate    
    if (robust%in%c("median", "trimmed", "tukey", "winsor"))
    {
        nlsFit$value <- (mad(resVec, 0)^2)*nlsDF 
    }
#    if (robust=="winsor")
#    {
#        K <- 1 + length(startVec)*var(psi.huber(resVec/s, deriv=1))
#    }
    if (robust%in%c("lms", "lts"))  # p. 202 i Rousseeuw and Leroy: Robust Regression and Outlier Detection
    {  
        scaleEst <- 1.4826*(1+5/(numObs-length(nlsFit$par)))*sqrt(median(resVec^2))                                 
        w <- (resVec/scaleEst < 2.5)
        nlsFit$value <- sum(w*resVec^2)/(sum(w)-length(nlsFit$par))    
    }

    
    ## Adding meaningful names for robust methods
    robust <- switch(robust, median="median", trimmed="metric trimming", tukey="Tukey's biweight", 
                             winsor="metric Winsorizing", lms="least median of squares",
                             lts="least trimmed squares")


    ## Collecting summary output
    sumVec <- c(NA, NA, NA, nlsSS, nlsDF, numObs)  # , alternative)
    sumList <- list(lenData = numObs, 
    alternative = NULL,  # alternative, 
    df.residual = numObs - length(startVec))


    ## The function call
#    callDetail <- match.call()
#    if (is.null(callDetail$fct)) {callDetail$fct <- substitute(l4())}


    ## The data set
    if (!is.null(logDose)) 
    {
        dose <- origDose
    }
    dataSet <- data.frame(origDose, origResp, assayNo, assayNoOld, wVec)
    
#    print(varNames0)
    if (identical(type, "event"))
    {
        names(dataSet) <- c(varNames0[c(2, 3, 1)], anName, paste("orig.", anName, sep = ""), "weights")
    } else {
        names(dataSet) <- c(varNames0[c(2, 1)], anName, paste("orig.", anName, sep = ""), "weights")
    }


#    ## Box-Cox information
#    bcVec <- c(lambda, boxcoxci)
#    if (all(is.na(bcVec))) {bcVec <- NULL}
#    if (!is.null(bcVec)) {bcVec <- c(bcVec, bcAdd)}


    ## Evaluating goodness-of-fit test
#    if (!is.null(gofTest)) {gofTest <- gofTest(resp, weights, predVec, sumList$"df.residual")}


#    ## Adjusting in case 'fctList' is specified
#    if (!is.null(fctList))
#    {
#        omitAllVec <- as.vector(unlist(omitList))
#
#        parmVec <- parmVec[omitAllVec] 
#        parmVecA <- parmVecA[omitAllVec] 
#        parmVecB <- parmVecB[omitAllVec] 
#        
#        orderVec <- match(as.vector(parmMat), nlsFit$par)
#        orderVec <- orderVec[complete.cases(orderVec)]       
#        
#        nlsFit$par <- nlsFit$par[orderVec]
#        nlsFit$hessian <- nlsFit$hessian[orderVec, orderVec]
#    }


    ## Constructing an index matrix for use in ED and SI
# (commented out Dec 7 2011, replaced by definition below of the index matrix)
#    hfct1 <- function(x)  # helper function
#    {
#        uniVec <- unique(x[!is.na(x)])
#        rv <- rep(NA, length(x))
#        for (i in 1:length(uniVec))
#        {
#            rv[abs(x-uniVec[i]) < 1e-12] <- i
#        }
#        rv
#    }
#    hfct2 <- function(x)
#    {
#        length(unique(x))
#    }
##    parmMat <- t(parmMat)
#    mat1 <- t(apply(t(parmMat), 1, hfct1))  # , 1:ncol(parmMat)))
#    cnccl <- head(cumsum(ncclVec), -1)
##    mat2 <- mat1
#    if (nrow(mat1) == 1) {mat1 <- t(mat1)}  # in case of only one curve
#    mat1[-1, ] <- mat1[-1, ] + cnccl

    ## Matrix of first derivatives evaluated at the parameter estimates
    if (isDF)
    {
#        print((parmMat[assayNo, , drop = FALSE])[isFinite, , drop = FALSE])
        deriv1Mat <- fct$"deriv1"(dose, (parmMat[assayNo, , drop = FALSE])[isFinite, , drop = FALSE])
    } else {
        deriv1Mat <- NULL
    }
#    deriv1Mat <- NULL

    ## Box-Cox information
    if (!is.null(bcVal))
    {
        bcVec <- list(lambda = bcVal, ci = c(NA, NA), bcAdd = bcAdd)
    } else {
        bcVec <- NULL
    }

    ## Parameter estimates
    coefVec <- nlsFit$par
    names(coefVec) <- parmVec
    
    ## Constructing the index matrix
#    parmMat <- t(parmMat)
    indexMat <- apply(t(parmMat), 2, function(x){match(x, coefVec)})

    ## Constructing data list ... where is it used?
    wName <- callDetail[["weights"]]
    if (is.null(wName)) 
    {
        wName <- "weights"
    } else {
        wName <- deparse(wName)
    }
#    dataList <- list(dose = as.vector(origDose), origResp = as.vector(origResp), weights = wVec, 
    dataList <- list(dose = origDose, origResp = as.vector(origResp), weights = wVec, 
    curveid = assayNoOld, resp = as.vector(resp),
    names = list(dName = varNames[1], orName = varNames[2], wName = wName, cNames = anName, rName = ""))
    if (identical(type, "event"))
    {
        dataList <- list(dose = dose, origResp = resp2, weights = wVec[isFinite], 
        curveid = assayNoOld[isFinite], plotid = plotid2, resp = resp2,
        names = list(dName = varNames[1], orName = varNames[2], wName = wName, cNames = anName, rName = ""))
    }
    if (identical(type, "ssd"))
    {
        dataList <- list(dose = dose, origResp = resp, weights = wVec[ifDose2], 
                         curveid = assayNoOld[ifDose2], resp = resp,
                         names = list(dName = varNames[1], orName = varNames[2], 
                                      wName = wName, cNames = anName, rName = ""))
    }


    ## What about naming the vector of weights?

    ## Returning the fit
#    returnList <- list(varParm, nlsFit, list(plotFct, logDose), sumVec, startVec, list(parmVec, parmVecA, parmVecB), 
    returnList <- list(NULL, nlsFit, list(plotFct, logDose), sumVec, startVecSc * longScaleVec, 
#    returnList <- list(nlsFit, list(plotFct, logDose), sumVec, startVecSc * longScaleVec, 
    list(parmVec, parmVecA, parmVecB), 
    diagMat, callDetail, dataSet, t(parmMat), fct, robust, estMethod, numObs - length(startVec), 
#    anovaModel0, gofTest, 
#    sumList, NULL, pmFct, pfFct, type, mat1, logDose, cm, deriv1Mat, 
    sumList, NULL, pmFct, pfFct, type, indexMat, logDose, cm, deriv1Mat, 
    anName, data, wVec, 
    dataList,
    coefVec, bcVec,
    indexMat2)
    
    names(returnList) <- c("varParm", "fit", "curve", "summary", "start", "parNames", 
                           "predres", "call", "data", 
#    names(returnList) <- c("fit", "curve", "summary", "start", "parNames", "predres", "call", "data", 
    "parmMat", "fct", "robust", "estMethod", "df.residual", 
#    "anova", "gofTest", 
    "sumList", "scaleFct", "pmFct", "pfFct", "type", "indexMat", "logDose", "cm", "deriv1",
    "curveVarNam", "origData", "weights",
    "dataList", "coefficients", "boxcox", "indexMat2")
    ## Argument "scaleFct" not used anymore
    class(returnList) <- c("drc")  # , class(fct))

    return(returnList)
}
OnofriAndreaPG/drcSeedGerm documentation built on Oct. 9, 2019, 3:45 p.m.