bmd: Benchmark dose estimation for dose-response data

mjust<-function (modelList,expressions, dataused, level=0.95, seType="san") {
  require(sandwich, quietly = TRUE)
  require(car, quietly = TRUE)
  require(RLRsim, quietly = TRUE)
  require(Matrix, quietly = TRUE)
  deltab<-function (object, g, func = g, ...)
  {
    if (!is.character(g))
      stop("The argument 'g' must be a character string")
    if(inherits(object,c("lmerMod","lme"))){
      para<-fixef(object)
      coefVec<-fixef(object)
      para.names<-sapply(strsplit(names(coefVec), ":"), "[[", 1)
      para.names[1] <- gsub("\\(Intercept\\)", "Intercept", para.names[1])
      names(para)<-para.names
    }else{
      para <- coef(object)
      if(inherits(object,"lm")){
        para.names <- names(coef(object))
        para.names[1] <- gsub("\\(Intercept\\)", "Intercept", para.names[1])
        names(para) <- para.names
      }else{
        if(inherits(object,"drc")){
          coefVec<-coef(object)
          #para.names<-sapply(strsplit(names(coefVec), ":"), "[[", 1)
          para.names.tmp<-gsub(":", "_",names(coefVec))
          para.names<-sapply(strsplit(para.names.tmp, "_\\(Intercept\\)"), "[[", 1)
          names(para)<-para.names
        }else{
          para.names <- names(para)
        }}}
    g <- parse(text = g)
    q <- length(para)
    for (i in 1:q) {
      assign(names(para)[i], para[i])
    }
    gd <- rep(0,q)
    for (i in 1:q) {
      gd[i] <- eval(D(g, names(para)[i]))
    }
    gd
  }
  makeIIDdecomp <- function(modelObject,g)
  {
    if(inherits(modelObject, "lmerMod")){
      numObsUsed <- length(predict(modelObject))
      #data.tmp<-dataused
      #data.tmp[is.na(data.tmp)]<-10
      allRepUnits<-unique(dataused[,names(getME(modelObject,"flist"))])
      repUnitsUsed<-unique(unlist(getME(modelObject,"flist")))
      naRepUnits<-as.numeric(setdiff(allRepUnits, repUnitsUsed))
      numInd<-length(repUnitsUsed)
      beta<-getME(modelObject,"beta")
      X<-getME(modelObject,"X")
      Y<-getME(modelObject,"y")
      Z<-getME(modelObject,"Z")
      A<-getME(modelObject,"A")
      Sigma<-sigma(modelObject)
      R<-diag(Sigma^2,numObsUsed)
      V<-sigma(modelObject)^2*t(A)%*%A+R
      Vminus<-solve(V)
      f<-function(i){
        if(allRepUnits[i]%in%repUnitsUsed){
          tmpList<-which(getME(modelObject,"flist")[[1]]==allRepUnits[i])
          Xi<-matrix(X[tmpList,],nrow=length(tmpList))
          Zi<-matrix(Z[tmpList,],nrow=length(tmpList))
          Vi<-V[tmpList,tmpList]
          Yi<-matrix(Y[tmpList],nrow=length(tmpList))
          as.matrix(-t(Xi)%*%solve(Vi)%*%(Yi-Xi%*%beta))} else{matrix(rep(0,ncol(X)),ncol=1)}}
      EstFun<-matrix(unlist(lapply(seq_along(allRepUnits),f)),nrow=length(allRepUnits),byrow=T)
      db<-deltab(modelObject,g)
      iidVec0<-as.matrix(-db%*%solve(t(X)%*%Vminus%*%X)*length(repUnitsUsed))%*%t(EstFun)
      if (!is.null(naRepUnits)) {
        iidVec <- sqrt(length(allRepUnits)/numInd) * iidVec0
      } else {
        iidVec <- iidVec0
      }
    }else{
      if(inherits(modelObject, "lme")){
        numObsUsed <- length(predict(modelObject))
        allRepUnits<-unique(modelObject$data[,attr(getGroups(modelObject),"label")])
        repUnitsUsed<-unique(getGroups(modelObject))
        naRepUnits<-as.numeric(setdiff(allRepUnits, repUnitsUsed))
        numInd<-length(repUnitsUsed)
        beta<-fixed.effects(modelObject)
        X<-extract.lmeDesign(modelObject)$X
        Y<-as.matrix(extract.lmeDesign(modelObject)$y)
        Z<-extract.lmeDesign(modelObject)$Z
        GB<-getVarCov(modelObject)
        G.Block<-matrix(as.numeric(GB),nrow=sqrt(length(GB)))
        G<-bdiag(rep(list(G.Block),numInd))
        Sigma<-modelObject$sigma
        if(is.null(modelObject$modelStruct$corStruct)){
          R<-diag(Sigma^2,numObsUsed)
        }else{
          R<-Sigma^2*bdiag(corMatrix(modelObject$modelStruct$corStruct))
        }
        V<-Z%*%G%*%t(Z)+R
        Vminus<-solve(V)
        f<-function(i){
          if(allRepUnits[i]%in%repUnitsUsed){
            tmpList<-which(getGroups(modelObject)==allRepUnits[i])
            Xi<-matrix(X[tmpList,],nrow=length(tmpList))
            Zi<-matrix(Z[tmpList,],nrow=length(tmpList))
            Vi<-V[tmpList,tmpList]
            Yi<-matrix(Y[tmpList],nrow=length(tmpList))
            as.matrix(-t(Xi)%*%solve(Vi)%*%(Yi-Xi%*%beta))} else{matrix(rep(0,ncol(X)),ncol=1)}}
        EstFun<-matrix(unlist(lapply(seq_along(allRepUnits),f)),nrow=length(allRepUnits),byrow=T)
        db<-deltab(modelObject,g)
        iidVec0<-as.matrix(-db%*%solve(t(X)%*%Vminus%*%X)*length(repUnitsUsed))%*%t(EstFun)
        if (!is.null(naRepUnits)) {
          iidVec <- sqrt(length(allRepUnits)/numInd) * iidVec0
        } else {
          iidVec <- iidVec0
        }
      }else{
        numObsUsed <- ifelse(inherits(modelObject, "coxph"),
                             modelObject$n, ifelse(inherits(modelObject,"nls"),
                                                   length(predict(modelObject)),ifelse(inherits(modelObject,"drc"),
                                                                                       length(predict(modelObject)), nrow(modelObject$model))))
        db<-deltab(modelObject,g)
        iidVec0 <- db %*%bread(modelObject) %*% t(estfun(modelObject))
        moNAac <- modelObject$na.action
        numObs <- numObsUsed + length(moNAac)
        numInd<-numObs
        iidVec <- rep(0, numObs)
        if (!is.null(moNAac)) {
          iidVec[-moNAac] <- sqrt(numObs/numObsUsed) * iidVec0[!is.na(iidVec0)]
        } else {
          iidVec <- iidVec0
        }}}
    list(iidVec = iidVec, numObsUsed = numObsUsed, numInd = numInd)
  }
  iidList <- list()
  numModels <- length(modelList)
  for(i in 1:numModels)
  {
    iidList[[i]]<- makeIIDdecomp(modelList[[i]], expressions[[i]])
  }
  iidresp <- matrix(as.vector(unlist(lapply(iidList, function(listElt){listElt[[1]]}))), nrow = numModels, byrow = TRUE)
  numObsUsed <- as.vector(unlist(lapply(iidList, function(listElt) {
    listElt[[2]]})))
  thetaEst <- rep(NA, numModels)
  thetaSe <- rep(NA, numModels)
  for(i in 1:numModels)
  {
    if (inherits(modelList[[i]], "drc")){
      coefVec <- coef(modelList[[i]])
      #names(coefVec) <- sapply(strsplit(names(coefVec), ":"), "[[", 1)
      para.names.tmp<-gsub(":", "_",names(coefVec))
      names(coefVec)<-sapply(strsplit(para.names.tmp, "_\\(Intercept\\)"), "[[", 1)
      deltaRes <- deltaMethod(coefVec,expressions[[i]],vcov(modelList[[i]]))
    } else {
      if (inherits(modelList[[i]], "lmerMod")){
        coefVec <- fixef(modelList[[i]])
        names(coefVec) <- sapply(strsplit(names(coefVec), ":"), "[[", 1)
        deltaRes <- deltaMethod(coefVec,expressions[[i]],vcov(modelList[[i]]))
      }else{
        deltaRes <- deltaMethod(modelList[[i]],expressions[[i]])
      }}
    thetaEst[i] <- deltaRes[1]
    thetaSe[i] <- deltaRes[2]
  }
  thetaEst <- unlist(thetaEst)
  thetaSe <- unlist(thetaSe)
  ## Calculating the estimated variance-covariance matrix of the parameter estimates
  numInd <- iidList[[1]]$numInd
  covar <- (iidresp %*% t(iidresp)) / numInd
  vcMat <- covar / numInd # Defining the finite-sample variance-covariance matrix
  ## Replacing sandwich estimates by model-based standard errors
  modbas <- seType == "mod"
  if (any(modbas))
  {
    corMat <- cov2cor(vcMat)
    ## Retrieving standard errors for the specified estimate from the individual fits
    modSE <- thetaSe
    sanSE <- sqrt(diag(vcMat))
    sanSE[modbas] <- modSE[modbas]
    vcMat <- diag(sanSE,nrow=length(sanSE)) %*% corMat %*% diag(sanSE,nrow=length(sanSE))
  }
  quant <- qnorm(1 - (1 - level)/2)
  numInd <- iidList[[1]]$numInd
  varMA <- vcMat
  seMA <- sqrt(diag(varMA))
  quantVal <- quant * seMA
  zVec <- thetaEst*(1/seMA)
  pvals <- 1 - pchisq(zVec * zVec, 1)
  retMat <- as.matrix(cbind(thetaEst, seMA, thetaEst - quantVal, thetaEst + quantVal,pvals))
  colnames(retMat) <- c("Estimate", "Std. Error", "Lower", "Upper", "Pr(>|z|)")
  output<-list(retMat,varMA)
  names(output)<-c("coef","covar")
  return(invisible(output))
}