R/drmEMls.R

In DoseResponse/drc: Analysis of Dose-Response Data

"drmEMls" <- 
function(dose, resp, multCurves, startVec, robustFct, weights, rmNA, dmf = NULL, 
doseScaling = 1, respScaling = 1, varcov = NULL)
{
    ## Defining the objective function and its derivative
    opfct <- function(parm)
    {
        sum(robustFct(((resp / respScaling) - multCurves((dose / doseScaling), parm)) / weights), 
            na.rm = rmNA)          
        # inverse weights enter multiplicatively before being squared 
        # (they have to be on same scale as the response)
    }
    
    if (!is.null(varcov))  # Note: robustFct() not used; derivatives don't work with this version
    {
        opfct <- function(parm)
        {  
            resVec <- (resp / respScaling) - multCurves((dose / doseScaling), parm)  
            resVec %*% solve(varcov) %*% resVec  # maybe a generalized inverse
        }    
    }  
    
    if (!is.null(dmf))
    {
        opdfct1 <- function(parm)
        {
#            apply(-2*(resp - multCurves(dose, parm))*dmf(dose, parm), 2, sum)
#            apply(-2*(resp - multCurves(dose, parm))*dmf(dose, parm), 2, appFct, cid)

            -2*((resp / respScaling) - multCurves((dose / doseScaling), parm)) * dmf((dose / doseScaling), parm)
        }
    } else {
        opdfct1 <- NULL
    }
       
    ## Defining additional self starter function (none needed)
    ssfct <- NULL

    ## Defining the log likelihood function
    llfct <- function(object)
    {
        degfre <- object$"sumList"$"lenData"  # "df.residual"  # object$summary[6]
        c( -(degfre/2)*(log(2*pi)+log(object$"fit"$"value")-log(degfre)+1), 
        object$"sumList"$"lenData" - object$"sumList"$"df.residual" + 1)
#        length(object$"fit"$"par") + 1 )
    }   
    
    ## Defining functions returning the residual variance, the variance-covariance matrix and the fixed effects estimates
    rvfct <- function(object)
    {
        object$"fit"$"value" / df.residual(object)  # object$"sumList"$"df.residual"
    }

    vcovfct <- function(object)
    {
#        scaledH <- (object$"fit"$"hessian")*(1/(2*object$"fit"$"ovalue"/object$"sumList"$"df.residual"))  # /2
        scaledH <- (object$"fit"$"hessian") / (2 * rvfct(object))
        invMat <- try(solve(scaledH), silent = TRUE)
    
        if (inherits(invMat, "try-error"))
        {
            ## More stable than 'solve' (suggested by Nicholas Lewin-Koh - 2007-02-12)
            ch <- try(chol(scaledH))
            if(inherits(ch, "try-error")) 
            {
                ch <- try(chol(0.99 * object$fit$hessian + 0.01 * diag(dim(object$fit$hessian)[1])))
            }
            ## Try regularizing if the varcov is unstable
            if(!inherits(ch, "try-error")) return(chol2inv(ch))
        } else {
            return(invMat)
        }
#        solve((object$"fit"$"hessian")*(1/rvfct(object))/2)  
#        solve((object$"fit"$"hessian")*(1/(object$"fit"$"ovalue"/object$"sumList"$"df.residual"))/2)   
    }
    
    parmfct <- function(fit, fixed = TRUE)
    {
        fit$par
    }

    rstanfct <- function(object)
    {
        rep(1, object$"sumList"$"lenData")*sqrt(summary(object)$"resVar")
    }


    return(list(llfct = llfct, opfct = opfct, opdfct1 = opdfct1, ssfct = ssfct, 
    rvfct = rvfct, vcovfct = vcovfct, parmfct = parmfct, rstanfct = rstanfct))
}