metafor: Meta-Analysis Package for R

############################################################################

.selmodel.pval <- function(yi, vi, alternative) {
   zi <- yi / sqrt(vi)
   if (alternative == "two.sided") {
      pval <- 2 * pnorm(abs(zi), lower.tail=FALSE)
   } else {
      pval <- pnorm(zi, lower.tail=alternative=="less")
   }
   return(pval)
}

.selmodel.verbose <- function(ll, beta, tau2, delta, mstyle, digits) {
   cat(mstyle$verbose(paste0("ll = ",         fmtx(ll,    digits[["fit"]], flag=" "),                "  ")))
   cat(mstyle$verbose(paste0("beta =",  paste(fmtx(beta,  digits[["est"]], flag=" "), collapse=" "), "  ")))
   cat(mstyle$verbose(paste0("tau2 =",        fmtx(tau2,  digits[["var"]], flag=" "),                "  ")))
   cat(mstyle$verbose(paste0("delta =", paste(fmtx(delta, digits[["est"]], flag=" "), collapse=" "))))
   cat("\n")
}

.mapfun <- function(x, lb, ub, fun=NA) {
   if (is.na(fun)) {
      if (lb==0 && ub==1) {
         plogis(x)
      } else {
         lb + (ub-lb) / (1 + exp(-x)) # map (-inf,inf) to (lb,ub)
      }
   } else {
      x <- sapply(x, fun)
      pmin(pmax(x, lb), ub)
   }
}

.mapinvfun <- function(x, lb, ub, fun=NA) {
   if (is.na(fun)) {
      if (lb==0 && ub==1) {
         qlogis(x)
      } else {
         log((x-lb)/(ub-x)) # map (lb,ub) to (-inf,inf)
      }
   } else {
      sapply(x, fun)
   }
}

############################################################################

.selmodel.int <- function(yvals, yi, vi, preci, yhat, wi.fun, delta, tau2, alternative, pval.min, steps) {
   pval <- .selmodel.pval(yvals, vi, alternative)
   pval[pval < pval.min]     <- pval.min
   pval[pval > (1-pval.min)] <- 1-pval.min
   wi.fun(pval, delta, yi, vi, preci, alternative, steps) * dnorm(yvals, yhat, sqrt(vi+tau2))
}

.selmodel.ll.cont <- function(par, yi, vi, X, preci, k, pX, pvals,
                              deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                              tau2.arg, tau2.transf, tau2.max, beta.arg,
                              wi.fun, steps, pgrp,
                              alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {

   mstyle <- .get.mstyle()

   beta  <- par[1:pX]
   tau2  <- par[pX+1]
   delta <- par[(pX+2):(pX+1+deltas)]

   beta <- ifelse(is.na(beta.arg), beta, beta.arg)

   if (tau2.transf)
      tau2 <- exp(tau2)

   tau2[!is.na(tau2.arg)] <- tau2.arg

   tau2[tau2 < .Machine$double.eps*10] <- 0
   tau2[tau2 > tau2.max] <- tau2.max

   if (delta.transf)
      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)

   delta <- ifelse(is.na(delta.arg), delta, delta.arg)

   yhat <- c(X %*% beta)

   Ai <- rep(NA_real_, k)

   for (i in seq_len(k)) {
      tmp <- try(integrate(.selmodel.int, lower=intCtrl$lower, upper=intCtrl$upper,
                           yi=yi[i], vi=vi[i], preci=preci[i], yhat=yhat[i], wi.fun=wi.fun,
                           delta=delta, tau2=tau2, alternative=alternative, pval.min=pval.min, steps=steps,
                           subdivisions=intCtrl$subdivisions, rel.tol=intCtrl$rel.tol)$value, silent=TRUE)
      if (inherits(tmp, "try-error"))
         stop(mstyle$stop(paste0("Could not integrate over density in study ", i, ".")), call.=FALSE)
      Ai[i] <- tmp
   }

   llval <- sum(log(wi.fun(pvals, delta, yi, vi, preci, alternative, steps)) + dnorm(yi, yhat, sqrt(vi+tau2), log=TRUE) - log(Ai))

   if (dofit) {
      res <- list(ll=llval, beta=beta, tau2=tau2, delta=delta)
      return(res)
   }

   if (verbose)
      .selmodel.verbose(ll=llval, beta=beta, tau2=tau2, delta=delta, mstyle=mstyle, digits=digits)

   if (verbose > 2) {
      xs <- seq(pval.min, 1-pval.min, length=1001)
      ys <- wi.fun(xs, delta, yi, vi, preci=1, alternative, steps)
      plot(xs, ys, type="l", lwd=2, xlab="p-value", ylab="Relative Likelihood of Selection")
   }

   return(-1*llval)

}

############################################################################

.selmodel.ll.stepfun <- function(par, yi, vi, X, preci, k, pX, pvals,
                                 deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                                 tau2.arg, tau2.transf, tau2.max, beta.arg,
                                 wi.fun, steps, pgrp,
                                 alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {

   mstyle <- .get.mstyle()

   beta  <- par[1:pX]
   tau2  <- par[pX+1]
   delta <- par[(pX+2):(pX+1+deltas)]

   beta <- ifelse(is.na(beta.arg), beta, beta.arg)

   if (tau2.transf)
      tau2 <- exp(tau2)

   tau2[!is.na(tau2.arg)] <- tau2.arg

   tau2[tau2 < .Machine$double.eps*10] <- 0
   tau2[tau2 > tau2.max] <- tau2.max

   if (decreasing) {

      if (delta.transf) {
         delta <- exp(delta)
         delta <- cumsum(c(0, -delta[-1]))
         delta <- exp(delta)
      }

   } else {

      if (delta.transf)
         delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)

   }

   delta <- ifelse(is.na(delta.arg), delta, delta.arg)

   if (decreasing && any(!is.na(delta.arg[-1])))
      delta <- rev(cummax(rev(delta)))

   yhat <- c(X %*% beta)

   sei <- sqrt(vi+tau2)

   N <- length(steps)

   Ai <- rep(NA_real_, k)

   if (alternative == "greater") {
      for (i in seq_len(k)) {
         Ai[i] <- pnorm(qnorm(steps[1], 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE)
         for (j in 2:N) {
            if (j < N) {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * (pnorm(qnorm(steps[j], 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE) - pnorm(qnorm(steps[j-1], 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE))
            } else {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * pnorm(qnorm(steps[j-1], 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=TRUE)
            }
         }
      }
   }

   if (alternative == "less") {
      for (i in seq_len(k)) {
         Ai[i] <- pnorm(qnorm(steps[1], 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE)
         for (j in 2:N) {
            if (j < N) {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * (pnorm(qnorm(steps[j], 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE) - pnorm(qnorm(steps[j-1], 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE))
            } else {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * pnorm(qnorm(steps[j-1], 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=FALSE)
            }
         }
      }
   }

   if (alternative == "two.sided") {
      for (i in seq_len(k)) {
         Ai[i] <- pnorm(qnorm(steps[1]/2, 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE) + pnorm(qnorm(steps[1]/2, 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE)
         for (j in 2:N) {
            if (j < N) {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * ((pnorm(qnorm(steps[j]/2, 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE) - pnorm(qnorm(steps[j-1]/2, 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=FALSE)) + (pnorm(qnorm(steps[j]/2, 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE) - pnorm(qnorm(steps[j-1]/2, 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE)))
            } else {
               Ai[i] <- Ai[i] + delta[j] / preci[i] * (pnorm(qnorm(steps[j-1]/2, 0, sqrt(vi[i]), lower.tail=FALSE), yhat[i], sei[i], lower.tail=TRUE) - pnorm(qnorm(steps[j-1]/2, 0, sqrt(vi[i]), lower.tail=TRUE), yhat[i], sei[i], lower.tail=TRUE))
            }
         }
      }
   }

   llval <- sum(log(delta[pgrp] / preci) + dnorm(yi, yhat, sei, log=TRUE) - log(Ai))

   if (dofit) {

      res <- list(ll=llval, beta=beta, tau2=tau2, delta=delta)
      return(res)

   }

   if (verbose)
      .selmodel.verbose(ll=llval, beta=beta, tau2=tau2, delta=delta, mstyle=mstyle, digits=digits)

   if (verbose > 2) {
      xs <- seq(0, 1, length=1001)
      ys <- wi.fun(xs, delta, yi, vi, preci=1, alternative, steps)
      plot(xs, ys, type="l", lwd=2, xlab="p-value", ylab="Relative Likelihood of Selection")
   }

   return(-1*llval)

}

############################################################################

.selmodel.ll.trunc <- function(par, yi, vi, X, preci, k, pX, pvals,
                               deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                               tau2.arg, tau2.transf, tau2.max, beta.arg,
                               wi.fun, steps, pgrp,
                               alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {

   mstyle <- .get.mstyle()

   beta  <- par[1:pX]
   tau2  <- par[pX+1]
   delta <- par[(pX+2):(pX+1+deltas)]

   beta <- ifelse(is.na(beta.arg), beta, beta.arg)

   if (tau2.transf)
      tau2 <- exp(tau2)

   tau2[!is.na(tau2.arg)] <- tau2.arg

   tau2[tau2 < .Machine$double.eps*10] <- 0
   tau2[tau2 > tau2.max] <- tau2.max

   if (delta.transf)
      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)

   delta <- ifelse(is.na(delta.arg), delta, delta.arg)

   yhat <- c(X %*% beta)

   sei <- sqrt(vi+tau2)

   if (is.na(steps)) {

      if (alternative == "greater")
         lli <- ifelse(yi > delta[2], 0, log(delta[1])) + dnorm(yi, yhat, sei, log=TRUE) - log(1 - (1-delta[1]) * pnorm(delta[2], yhat, sei, lower.tail=TRUE))

      if (alternative == "less")
         lli <- ifelse(yi < delta[2], 0, log(delta[1])) + dnorm(yi, yhat, sei, log=TRUE) - log(1 - (1-delta[1]) * pnorm(delta[2], yhat, sei, lower.tail=FALSE))

   } else {

      if (alternative == "greater")
         lli <- ifelse(yi > steps, 0, log(delta[1])) + dnorm(yi, yhat, sei, log=TRUE) - log(1 - (1-delta[1]) * pnorm(steps, yhat, sei, lower.tail=TRUE))

      if (alternative == "less")
         lli <- ifelse(yi < steps, 0, log(delta[1])) + dnorm(yi, yhat, sei, log=TRUE) - log(1 - (1-delta[1]) * pnorm(steps, yhat, sei, lower.tail=FALSE))

   }

   llval <- sum(lli)

   if (dofit) {

      res <- list(ll=llval, beta=beta, tau2=tau2, delta=delta)
      return(res)

   }

   if (verbose)
      .selmodel.verbose(ll=llval, beta=beta, tau2=tau2, delta=delta, mstyle=mstyle, digits=digits)

   return(-1*llval)

}

############################################################################

.rma.selmodel.ineqfun.pos <- function(par, yi, vi, X, preci, k, pX, pvals,
                                      deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                                      tau2.arg, tau2.transf, tau2.max, beta.arg,
                                      wi.fun, steps, pgrp, alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {

   delta <- par[-seq_len(pX+1)]

   if (delta.transf)
      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)

   delta <- ifelse(is.na(delta.arg), delta, delta.arg)

   diffs <- -diff(delta) # -1 * differences (delta1-delta2, delta2-delta3, ...) must be positive

   return(diffs)

}

.rma.selmodel.ineqfun.neg <- function(par, yi, vi, X, preci, k, pX, pvals,
                                      deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                                      tau2.arg, tau2.transf, tau2.max, beta.arg,
                                      wi.fun, steps, pgrp, alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {

   delta <- par[-seq_len(pX+1)]

   if (delta.transf)
      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)

   delta <- ifelse(is.na(delta.arg), delta, delta.arg)

   diffs <- diff(delta) # differences (delta1-delta2, delta2-delta3, ...) must be negative

   return(diffs)

}

############################################################################