R/predict_gMAP.R

In RBesT: R Bayesian Evidence Synthesis Tools

#' Predictions from gMAP analyses
#'
#' @name predict.gMAP
#'
#' @description
#' Produces a sample of the predictive distribution.
#'
#' @param x,object gMAP analysis object for which predictions are performed
#' @param newdata data.frame which must contain the same columns as
#' input into the gMAP analysis. If left out, then a posterior prediction for
#' the fitted data entries from the gMAP object is performed (shrinkage estimates).
#' @param probs defines quantiles to be reported.
#' @param type sets reported scale (\code{response} (default) or \code{link}).
#' @param na.action how to handle missings.
#' @param thin thinning applied is derived from the \code{gMAP} object.
#' @param digits number of displayed significant digits.
#' @param ... ignored.
#'
#' @details Predictions are made using the \eqn{\tau} prediction
#' stratum of the gMAP object. For details on the syntax, please refer
#' to \code{\link{predict.glm}} and the example below.
#'
#' @seealso \code{\link{gMAP}}, \code{\link{predict.glm}}
#'
#' @template example-start
#' @examples
#' # create a fake data set with a covariate
#' trans_cov <- transform(transplant, country=cut(1:11, c(0,5,8,Inf), c("CH", "US", "DE")))
#' set.seed(34246)
#' map <- gMAP(cbind(r, n-r) ~ 1 + country | study,
#'             data=trans_cov,
#'             tau.dist="HalfNormal",
#'             tau.prior=1,
#'             # Note on priors: we make the overall intercept weakly-informative
#'             # and the regression coefficients must have tighter sd as these are
#'             # deviations in the default contrast parametrization
#'             beta.prior=rbind(c(0,2), c(0,1), c(0,1)),
#'             family=binomial,
#'             ## ensure fast example runtime
#'             thin=1, chains=1)
#'
#' # posterior predictive distribution for each input data item (shrinkage estimates)
#' pred_cov <- predict(map)
#' pred_cov
#'
#' # extract sample as matrix
#' samp <- as.matrix(pred_cov)
#'
#' # predictive distribution for each input data item (if the input studies were new ones)
#' pred_cov_pred <- predict(map, trans_cov)
#' pred_cov_pred
#'
#'
#' # a summary function returns the results as matrix
#' summary(pred_cov)
#'
#' # obtain a prediction for new data with specific covariates
#' pred_new <- predict(map, data.frame(country="CH", study=12))
#' pred_new
#'
#' @template example-stop
#' @rdname predict.gMAP
#' @method predict gMAP
#' @export
predict.gMAP <- function(object, newdata, type=c("response", "link"), probs = c(0.025, 0.5, 0.975), na.action = na.pass, thin, ...) {
    f <- object$formula
    mf <- object$model
    tt <- terms(f, data=mf, lhs=1, rhs=1)
    type <- match.arg(type)
    if (missing(newdata)) {
        posterior_predict <- TRUE
        X <- model.matrix(f, mf, rhs=1)
        log_offset <- object$log_offset
        group.factor <- model.part(f, data = mf, rhs = 2)
    } else {
        posterior_predict <- FALSE
        Terms <- delete.response(tt)
        ## replace model frame with newdata context
        m <- model.frame(Terms, newdata, na.action = na.action,
                         xlev = .getXlevels(tt, mf))
        if (!is.null(cl <- attr(Terms, "dataClasses")))
            .checkMFClasses(cl, m)
        X <- model.matrix(Terms, m, contrasts.arg = attr(model.matrix(f, mf, rhs=1), "contrasts"))
        log_offset <- rep(0, nrow(X))
        if (!is.null(off.num <- attr(tt, "offset")))
            for (i in off.num) log_offset <- log_offset + eval(attr(tt,
                                                            "variables")[[i + 1]], newdata)
        if (!is.null(object$call$offset))
            log_offset <- log_offset + eval(object$call$offset, newdata)

        group.factor <- model.part(f, data = newdata, rhs = 2)
    }

    if(ncol(group.factor) != 1)
        stop("Grouping factor must be a single term (study).")
    group.factor <- group.factor[,1]

    if (!is.factor(group.factor)) {
        group.factor <- factor(group.factor)
    }
    labels <- as.character(group.factor)
    group.index <- array(as.integer(group.factor))

    ## nubmer of groups coded by the factor
    n.groups <- nlevels(group.factor)
    ## number of groups actually observed in the data
    n.groups.obs <- length(unique(group.index))

    if(missing(thin)) {
        thin <- object$thin
    }

    beta <- rstan::extract(object$fit, inc_warmup=FALSE, permuted=FALSE, pars="beta")
    n.pred <- nrow(X)
    n.iter <- dim(beta)[1]
    n.chains <- dim(beta)[2]

    if(posterior_predict) {
        pred <- aperm(rstan::extract(object$fit, inc_warmup=FALSE, permuted=FALSE, pars="theta"), c(3,1,2))
    } else {
        pred <- apply(beta, c(1,2), function(x) X %*% x)
        if(n.pred==1)
            pred <- array(pred, dim=c(1,dim(pred)))
    }

    sub_ind <- seq(1,n.iter,by=thin)

    pred <- t(matrix(pred[,sub_ind,], nrow=n.pred))

    if(!posterior_predict) {
        ## in case we make a prediction unconditional on the fitted
        ## data, we draw random effects here (one for each study per
        ## iteration)

        S <- nrow(pred)
        ## sample random effects for as many groups defined, which can
        ## be more than the ones in the data set, since we sample for
        ## all defined factor levels
        tau <- as.vector(rstan::extract(object$fit, inc_warmup=FALSE, permuted=FALSE, pars=paste0("tau[", object$tau.strata.pred, "]"))[sub_ind,,])
        if(object$REdist == "normal") {
            re <- tau * matrix(rnorm(n.groups * S, 0, 1), nrow=S)
        }
        if(object$REdist == "t") {
            re <- tau * matrix(rt(n.groups * S, df=object$t.df), nrow=S)
        }

        ## ... and add it to predictions
        pred <- pred + re[,group.index]
    }

    if(type == "response")
        pred <- object$family$linkinv(pred)

    predNames <- NULL
    if(!is.null(rownames(X)))
        predNames <- rownames(X)
    dimnames(pred) <- list(NULL, predNames)

    stat <- SimSum(pred, probs=probs, margin=2)
    attr(pred, "summary") <- stat
    attr(pred, "type") <- type
    attr(pred, "family") <- object$family
    attr(pred, "sigma_ref") <- object$sigma_ref
    invisible(structure(pred, class=c("gMAPpred")))
}

#' @rdname predict.gMAP
#' @method print gMAPpred
#' @export
print.gMAPpred <- function(x, digits=3, ...) {
    cat("Meta-Analytic-Predictive Prior Predictions\n")
    cat("Scale:", attr(x, "type"),"\n")
    cat("\n")
    cat("Summary:\n")
    print(signif(attr(x, "summary"), digits=digits))
}

#' @rdname predict.gMAP
#' @method summary gMAPpred
#' @export
summary.gMAPpred <- function(object, ...) {
    attr(object, "summary")
}

#' @rdname predict.gMAP
#' @method as.matrix gMAPpred
#' @export
as.matrix.gMAPpred <- function(x, ...) {
    class(x) <- "matrix"
    x
}