#'
#' @title Average predictions from several (non)linear models based on IC weights
#' @name predict_nlme
#' @rdname predict_nlme
#' @description Computes weights based on AIC, AICc, or BIC and it generates weighted predictions by
#' the relative value of the IC values
#' @param ... nlme, lme, gls or gnls objects.
#' @param criteria either \sQuote{AIC}, \sQuote{AICc} or \sQuote{BIC}.
#' @param interval either \sQuote{none}, \sQuote{confidence} or \sQuote{prediction}.
#' It is also possible to choose \sQuote{new-prediction}, which is a prediction that
#' resamples the random effects (only relevant for \sQuote{lme} or \sQuote{nlme} objects.)
#' @param level probability level for the interval (default 0.95)
#' @param nsim number of simulations to perform for intervals. Default 1000.
#' @param plevel parameter level prediction to be passed to prediciton functions.
#' @param newdata new data frame for predictions
#' @param weights vector of weights of the same length as the number of models. It should sum up to one and
#' it will override the information-criteria based weights. The weights should match the order of the models.
#' @return numeric vector of the same length as the fitted object.
#' @note all the objects should be fitted to the same data. The weights are
#' based on the IC value.
#' @seealso \code{\link{predict.nlme}} \code{\link{predict.lme}} \code{\link{predict.gnls}}
#' @export
#' @examples
#' \donttest{
#' ## Example
#' require(ggplot2)
#' require(nlme)
#' data(Orange)
#'
#' ## All models should be fitted using Maximum Likelihood
#' fm.L <- nlme(circumference ~ SSlogis(age, Asym, xmid, scal),
#' random = pdDiag(Asym + xmid + scal ~ 1),
#' method = "ML", data = Orange)
#' fm.G <- nlme(circumference ~ SSgompertz(age, Asym, b2, b3),
#' random = pdDiag(Asym + b2 + b3 ~ 1),
#' method = "ML", data = Orange)
#' fm.F <- nlme(circumference ~ SSfpl(age, A, B, xmid, scal),
#' random = pdDiag(A + B + xmid + scal ~ 1),
#' method = "ML", data = Orange)
#' fm.B <- nlme(circumference ~ SSbg4rp(age, w.max, lt.e, ldtm, ldtb),
#' random = pdDiag(w.max + lt.e + ldtm + ldtb ~ 1),
#' method = "ML", data = Orange)
#'
#' ## Print the table with weights
#' IC_tab(fm.L, fm.G, fm.F, fm.B)
#'
#' ## Each model prediction is weighted according to their AIC values
#' prd <- predict_nlme(fm.L, fm.G, fm.F, fm.B)
#'
#' ggplot(data = Orange, aes(x = age, y = circumference)) +
#' geom_point() +
#' geom_line(aes(y = predict(fm.L, level = 0), color = "Logistic")) +
#' geom_line(aes(y = predict(fm.G, level = 0), color = "Gompertz")) +
#' geom_line(aes(y = predict(fm.F, level = 0), color = "4P-Logistic")) +
#' geom_line(aes(y = predict(fm.B, level = 0), color = "Beta")) +
#' geom_line(aes(y = prd, color = "Avg. Model"), linewidth = 1.2)
#' }
predict_nlme <- function(..., criteria = c("AIC", "AICc", "BIC"),
interval = c("none", "confidence", "prediction", "new-prediction"),
level = 0.95, nsim = 1e3, plevel = 0,
newdata = NULL, weights){
objs <- list(...)
criteria <- match.arg(criteria)
interval <- match.arg(interval)
nms <- get_mnames(match.call())
lobjs <- length(objs)
wtab <- data.frame(model = character(lobjs), IC = NA)
nr <- stats::nobs(objs[[1]])
if(!is.null(newdata)) nr <- nrow(newdata)
if(interval == "none"){
prd.mat <- matrix(nrow = nr, ncol = lobjs)
}else{
prd.mat <- matrix(nrow = nr, ncol = lobjs)
prd.mat.se <- matrix(nrow = nr, ncol = lobjs)
prd.mat.lwr <- matrix(nrow = nr, ncol = lobjs)
prd.mat.upr <- matrix(nrow = nr, ncol = lobjs)
}
data.name <- as.character(deparse(objs[[1]]$call$data))
for(i in seq_len(lobjs)){
obj <- objs[[i]]
if(!inherits(obj, c("nlme","lme","gnls","gls"))) stop("All objects should be of class 'nlme', 'lme', 'gnls' or 'gls'")
if(data.name != as.character(deparse(obj$call$data))) stop("All models should be fitted to the same data")
wtab$model[i] <- nms[i]
if(criteria == "AIC") wtab$IC[i] <- stats::AIC(obj)
if(criteria == "AICc") wtab$IC[i] <- AICc(obj)
if(criteria == "BIC") wtab$IC[i] <- stats::BIC(obj)
}
## Check weights
if(!missing(weights)){
if(length(weights) != lobjs)
stop("'weights' should be a vector of length equal to the number of models", call. = FALSE)
if(isFALSE(all.equal(sum(weights), 1)))
stop("'sum of 'weights' should be equal to 1.", call. = FALSE)
if(any(weights < 0))
stop("all weights should be greater than zero", call. = FALSE)
if(any(weights > 1))
stop("all weights should be less than one", call. = FALSE)
}
## Predictions
if(interval == "none"){
for(i in seq_len(lobjs)){
obj <- objs[[i]]
if(inherits(obj, "lme")){
if(is.null(newdata)){
prd.mat[,i] <- predict(obj, level = plevel)
}else{
prd.mat[,i] <- predict(obj, newdata = newdata, level = plevel)
}
}
if(inherits(obj, "gls")){
if(is.null(newdata)){
prd.mat[,i] <- predict(obj)
}else{
prd.mat[,i] <- predict(obj, newdata = newdata)
}
}
}
}
if(interval == "confidence" || interval == "prediction" || interval == "new-prediction"){
if(interval == "confidence") psim <- 1
if(interval == "prediction") psim <- 2
if(interval == "new-prediction") psim <- 3
lb <- (1 - level)/2
ub <- 1 - lb
for(i in seq_len(lobjs)){
obj <- objs[[i]]
if(inherits(obj, "lme")){
if(inherits(obj, "nlme")){
tmp.sim <- simulate_nlme(obj, psim = psim, nsim = nsim,
level = plevel, newdata = newdata)
}else{
tmp.sim <- simulate_lme(obj, psim = psim, nsim = nsim,
level = plevel, newdata = newdata)
}
}
if(inherits(obj, "gls")){
if(psim == 3) stop("new-prediction is only possible for lme or nlme objects")
if(inherits(obj, "gnls")){
tmp.sim <- simulate_nlme(obj, psim = psim, nsim = nsim, newdata = newdata)
}else{
tmp.sim <- simulate_lme(obj, psim = psim, nsim = nsim, newdata = newdata)
}
}
prd.mat[,i] <- apply(tmp.sim, 1, quantile, probs = 0.5)
prd.mat.se[,i] <- apply(tmp.sim, 1, sd)
prd.mat.lwr[,i] <- apply(tmp.sim, 1, quantile, probs = lb)
prd.mat.upr[,i] <- apply(tmp.sim, 1, quantile, probs = ub)
}
}
wtab$dIC <- wtab$IC - min(wtab$IC)
if(missing(weights)){
wtab$weight <- exp(-0.5 * wtab$dIC) / sum(exp(-0.5 * wtab$dIC))
}else{
wtab$weight <- weights
}
if(interval == "none"){
ans <- rowSums(sweep(prd.mat, 2, wtab$weight, "*"))
}else{
prd <- rowSums(sweep(prd.mat, 2, wtab$weight, "*"))
se <- rowSums(sweep(prd.mat.se, 2, wtab$weight, "*"))
lwr <- rowSums(sweep(prd.mat.lwr, 2, wtab$weight, "*"))
upr <- rowSums(sweep(prd.mat.upr, 2, wtab$weight, "*"))
ans <- cbind(prd, se, lwr, upr)
colnames(ans) <- c("Estimate", "Est.Error",
paste0("Q", (1 - level)/2 * 100),
paste0("Q", (1 - (1 - level)/2)*100))
}
return(ans)
}
#' @rdname predict_nlme
#' @description predict function for objects of class \code{\link[nlme]{lme}}
#' @export
predict_lme <- predict_nlme
#' @rdname predict_nlme
#' @description predict function for objects of class \code{\link[nlme]{gnls}}
#' @export
predict_gnls <- predict_nlme
#' @rdname predict_nlme
#' @description predict function for objects of class \code{\link[nlme]{gls}}
#' @export
predict_gls <- predict_nlme
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.