demo/emaxfit.R

In MSToolkit: The MSToolkit library for clinical trial design

emax.fit<- function (y, dose, iparm = c("ed50", "e0", "emax"), ed50cutoff = 2.5 * max(doselev)) 
{
    dm <- tapply(y, dose, mean)
    dsd <- sqrt(tapply(y, dose, var))
    intercept <- rep(1, length(dose))
    doselev <- sort(unique(dose))
    est3 <- rep(NA, 3)
    estA <- rep(NA, 2)

        e0 <- min(dm)
        emax <- max(dm) - e0
		ed50 <- max(doselev)/2
        Sparm <- c(ed50 = ed50, e0 = e0, emax = emax)

        fit <- try(nls(y ~ e0 + (emax * dose)/(dose + ed50), 
            start = Sparm, control = nls.control(maxiter = 100), 
            trace = F), silent = T)

    AcceptConv <- F
    if (class(fit) == "nls") {
        est3 <- coef(fit)
        if (coef(fit)[1] > 0) {
            if (coef(fit)[1] <= ed50cutoff) {
                vc <- as.vector(vcov(fit))
                seout <- SeEmax3(fit, doselev)
                fitpred <- seout$fitpred
                sdpred <- seout$sdpred
				sddif <- seout$sddif
                AcceptConv <- T
            }
        }
    }
    if (!AcceptConv) {
        fitL <- lm(y ~ dose)
        sigL <- summary(fitL)$sigma
        predobj <- predict(fitL, data.frame(dose = doselev), se.fit = T)
        fitpred <- predobj$fit
        sdpred <- predobj$se.fit
		sddif <- doselev * sqrt(vcov(fitL)[2, 2])
    }
    return(list(fitpred = fitpred, 
        sdpred = sdpred, sddif = sddif))
}

SeEmax3<-function (fit, doselev) 
{
    E0<-coef(fit)[2]
    ED50<-coef(fit)[1]
    EMAX<-coef(fit)[3]
    fitpred <- E0+(EMAX*doselev)/(ED50+doselev)
    fitdif <- fitpred - E0
    vfcov <- vcov(fit)
    sddif <- sqrt((fitdif^2) * (vfcov[3, 3]/coef(fit)[3]^2 + 
        vfcov[1, 1]/(coef(fit)[1] + doselev)^2 - (2 * vfcov[1, 
        3])/(coef(fit)[3] * (coef(fit)[1] + doselev))))
    covint <- -((doselev * coef(fit)[3])/(coef(fit)[1] + doselev)^2) * 
        vfcov[2, 1] + (doselev/(coef(fit)[1] + doselev)) * vfcov[2, 
        3]
    sdpred <- sqrt(vfcov[2, 2] + sddif^2 + 2 * covint)
    return(list(fitpred = fitpred, sdpred = sdpred, fitdif = fitdif, 
        sddif = sddif))
}