R/predictSG.drc.R

In OnofriAndreaPG/drcSeedGerm: Utilities for data analyses in seed germination/emergence assays

"predictSG" <- function(object, newdata, se.fit = FALSE, interval = c("none", "confidence", "prediction"),
level = 0.95, na.action = na.pass, od = FALSE, vcov. = vcov, ...)
{
    ## Checking arguments
    interval <- match.arg(interval)
    respType <- object[["type"]]

    dataList <- object[["dataList"]]
    doseDim <- ncol(dataList[["dose"]])
    if (is.null(doseDim)) {doseDim <- 1}

    ## Assigning dataset from object if no data frame is provided
    if (missing(newdata))
    {
#        predValues <- fitted(object)  # not used
#        newdata <- data.frame(object$data[, 1], object$data[, 3])
#        dataList <- object[["dataList"]]

        ## New part (25/6-2014)
        doseVec <- dataList[["dose"]]
        if (identical(respType, "event"))
        {
            groupLevels <- as.character(dataList[["plotid"]])
        } else {
            groupLevels <- as.character(dataList[["curveid"]])
        }
#
#        if (identical(respType, "event"))
#        {
#            newdata <- data.frame(dataList[["dose"]], dataList[["plotid"]])
#        } else {
#            newdata <- data.frame(dataList[["dose"]], dataList[["curveid"]])
#        }
    } else {
      #I added these lines in 2018. Not yet on gitHub#################
      if(is.vector(dataList[["dose"]])==T){
        dName <- dataList[["names"]][["dName"]]
        if (any(names(newdata) %in% dName))
        {
          doseVec <- newdata[, dName]
        } else {
            doseVec <- newdata[, 1]
#            warning("Dose variable not in 'newdata'")
        }
        }else{
        #numPred <- length(dataList[["dose"]][1, ])
        #doseVec <- newdata[, 1:numPred]
          dName <- colnames(dataList$dose)
          doseVec <- newdata
        }

        #Non ho capito il ruolo di cName
        cName <- dataList[["names"]][["cNames"]]
        if (any(names(newdata) %in% cName))
        {
            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
    }
    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!
#    doseVec <- newdata[, 1]
#    groupLevels <- as.character(newdata[, 2])  # 'as.character()' used to suppress factor levels
#    noNewData <- length(doseVec)


    ## Transforming to dose scale if necessary
    powerExp <- (object$"curve")[[2]]
    if (!is.null(powerExp))
    {
        doseVec <- powerExp ^ doseVec
    }

    ## Retrieving matrix of parameter estimates
    parmMat <- object[["parmMat"]]
    pm <- t(parmMat[, groupLevels, drop = FALSE])

#    parmNames <- colnames(parmMat)
#    lenCN <- length(parmNames)
#    indVec <- 1:lenCN
#    names(indVec) <- parmNames
#    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]
#        } 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)
#    for (i in 1:lenCN)
#    {
#        naVec[i] <- any(is.na(parmMat[, i]))
#    }
#    parmMat <- parmMat[, !naVec, drop = FALSE]


    ## Retrieving variance-covariance matrix
    sumObj <- summary(object, od = od)
#    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)


#    retMat <- matrix(0, lenIV, 4)
    retMat <- matrix(0, noNewData, 4)
    colnames(retMat) <- c("Prediction", "SE", "Lower", "Upper")
    objFct <- object[["fct"]]
    retMat[, 1] <- objFct$"fct"(doseVec, pm) #Predictions

    ## Checking if derivatives are available
    deriv1 <- objFct$"deriv1"
    if (is.null(deriv1))
    {
        return(retMat[, 1])
    }

    ## Calculating the quantile to be used in the confidence intervals
    if (!identical(interval, "none"))
    { #Se interval != "none" oppure se assente
        if (identical(respType, "continuous"))
          #usa t o z. T è solo per continuo
        {
            tquan <- qt(1 - (1 - level)/2, df.residual(object))
        } else {
            tquan <- qnorm(1 - (1 - level)/2)
        }
    }

    ## Calculating standard errors and/or confidence intervals
    if (se.fit || (!identical(interval, "none")))
    {
         if (identical(interval, "prediction"))
         {
             sumObjRV <- sumObj$"resVar"
         } else {
             sumObjRV <- 0
         }
#        rowIndex <- 1
#        for (i in indexVec)
#        for (i in 1:ncol(indexMat))

#        groupLevels <- newdata[, 2]
        piMat <- indexMat[, groupLevels, drop = FALSE]
#        print(piMat)
#        print(groupLevels)
        for (rowIndex in 1:noNewData)
        {
#            parmInd <- indexMat[, i]
#            print(indexVec)
#            print(varMat)
#            print(parmInd)

#            varCov <- varMat[parmInd, parmInd]
#            print(varCov)
#            groupLevels <- newdata[, 2]
#            parmInd <- indexMat[, groupLevels[rowIndex]]
#            varCov <- varMat[parmInd, parmInd]

            parmInd <- piMat[, rowIndex]
            varCov <- vcovMat[parmInd, parmInd]

#            parmChosen <- t(parmMat[, i, drop = FALSE])
#            parmChosen <- t(pm[, rowIndex, drop = FALSE])
#            dfEval <- deriv1(doseVec[rowIndex], parmChosen)

            #Edited Andrea Onofri 7/12/18
            if(is.vector(dataList[["dose"]])==T){
            dfEval <- deriv1(doseVec[rowIndex], pm[rowIndex, , drop = FALSE])
            varVal <- dfEval %*% varCov %*% dfEval
            }else{
            dfEval <- deriv1(doseVec[rowIndex,], pm[rowIndex, , drop = FALSE])
            varVal <- dfEval %*% varCov %*% t(dfEval)
            }

            #varVal <- dfEval %*% varCov %*% dfEval
            retMat[rowIndex, 2] <- sqrt(varVal)
#            retMat[rowIndex, 2] <- sqrt(dfEval %*% varCov %*% dfEval)

            if (!se.fit)
            {
                retMat[rowIndex, 3:4] <- retMat[rowIndex, 1] + (tquan * sqrt(varVal + sumObjRV)) * c(-1, 1)
#                retMat[rowIndex, 3] <- retMat[rowIndex, 1] - tquan * sqrt(varVal + sumObjRV)
#                retMat[rowIndex, 4] <- retMat[rowIndex, 1] + tquan * sqrt(varVal + sumObjRV)
            }
#            if (identical(interval, "confidence"))
#            {
#                retMat[rowIndex, 3] <- 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
        }
    }
    ## 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])
}