#' EMD-R2
#'
#' Performs EMD-R2 by using the given regression function in each MIMF.
#'
#' @param formula A formula describing the models to be fit for each MIMF.
#' @param mimf \code{mimf} object containing both the response and predictor
#' IMFs.
#' @param covariates A list or data.frame containing eventual covariates
#' which are not IMFs.
#' @param tt A vector of time indices for \code{mimf}. The default is to use
#' the \code{tt} attribute.
#' @param reg.fun The regression function to perform either as a function or
#' as a character. Should include \code{formula} and \code{data} arguments.
#' @param reg.args A list containing additional arguments to be passed to
#' \code{reg.fun}.
#' @param pimf.args A list of arguments to be passed to the function
#' \code{\link{pimf}} to prepare the IMFs for the regression.
#' @param xy.args A character vector of length two indicating the argument
#' names corresponding to the predictors and response in the case
#' \code{reg.fun} does not include arguments \code{formula} and \code{data}.
#' The default of \code{xy.args = c("x", "y")} should cover the majority
#' of cases.
#'
#' @details \code{emdr2} wraps the chosen regression function (\code{reg.fun}),
#' and applies it once per MIMF in \code{mimf}. Variables given in
#' \code{covariates} are used with each IMF.
#'
#' \code{reg.fun} can include either a couple \code{formula} and \code{data}
#' arguments, either \code{x} and \code{y} arguments. Particular
#' argument names for the response or predictor variables can be given
#' in \code{xy.args}.
#'
#' @return A \code{emdr2} object which is a list of length nimf containing all
#' results of \code{reg.fun} on individual MIMFs.
#'
#' @seealso \code{\link{extract.emdr2}} to extract a particular element of the
#' result from \code{reg.fun} for each submodel. \code{\link{coef.emdr2}}
#' to specifically extract the coefficients. \code{\link{predict.emdr2}}
#' to obtain predictions from the whole EMD-R2 model.
#'
#' @references
#' Masselot, P., Chebana, F., Belanger, D., St-Hilaire, A., Abdous, B.,
#' Gosselin, P., Ouarda, T.B.M.J., 2018. EMD-regression for modelling
#' multi-scale relationships, and application to weather-related
#' cardiovascular mortality. \emph{Science of The Total Environment}
#' 612, 1018-1029.
#'
#' @examples
#'
#' library(dlnm)
#' library(glmnet)
#'
#' dat <- chicagoNMMAPS[,c("death", "temp", "rhum")]
#'
#' mimfs <- memd(dat)
#' cmimfs <- combine.mimf(mimfs, list(12:13, 14:17, 18:19),
#' new.names = c("C12", "C13", "r"))
#' amps <- mean_amplitude(cmimfs)
#'
#' # EMD-R2 with glm
#' lm.R2 <- emdr2(death ~ temp + rhum, mimf = cmimfs)
#' betas.R2 <- coef(lm.R2)
#' sensitivity.R2 <- sensitivity(amps[,-1], coefs = betas.R2[,-1])
#' plot_emdr(sensitivity.R2, periods = period(cmimfs)[,-1],
#' col = c("red", "blue"), pch = 16:17)
#' abline(h = 0, lty = 2)
#'
#' # EMD-R2 with lasso
#' lasso.R2 <- emdr2(death ~ temp + rhum, mimf = cmimfs,
#' reg.fun = "cv.glmnet")
#' betas.R2 <- coef(lasso.R2, s = "lambda.1se")
#' sensitivity.R2 <- sensitivity(amps[,-1], coefs = betas.R2[,-1])
#' plot_emdr(sensitivity.R2, periods = period(cmimfs)[,-1],
#' show.coef = "nonzero", pch = 16:17)
#' abline(h = 0, lty = 2)
#'
#'
#' @export
emdr2 <- function(formula, mimf, covariates = NULL, tt = attr(mimf,"tt"), reg.fun = "glm", reg.args = list(), pimf.args = list(), xy.args = c("x", "y"))
{
nimf <- dim(mimf)[2]
vars <- all.vars(formula)
if ("." %in% vars){
vars <- vars[-which(vars==".")]
toadd <- dimnames(mimf)[[3]][!dimnames(mimf)[[3]] %in% vars]
vars <- c(vars, toadd)
}
reg.formals <- formals(reg.fun)
result <- vector("list",nimf)
for (i in 1:nimf){
def.pimf <- list(x = mimf[,i,vars[-1]], y = mimf[,i,vars[1]], tt = tt, covariates = covariates)
pimf.args <- c(def.pimf, pimf.args[!names(pimf.args) %in% names(def.pimf)])
iprepared <- do.call(pimf, pimf.args)
if ("formula" %in% names(reg.formals)){
iform <- stats::as.formula(sprintf("y ~ %s", paste(names(iprepared)[-1], collapse = " + ")))
def.args <- list(formula = iform, data = iprepared)
} else {
if (length(xy.args) != 2) stop(sprintf("'xy.args' must be provided since there is no 'formula' argument in '%s'", reg.fun))
def.args <- list()
def.args[[xy.args[1]]] <- as.matrix(iprepared[,-1])
def.args[[xy.args[2]]] <- iprepared[,1]
}
reg.args <- c(def.args, reg.args[!names(reg.args) %in% names(def.args)])
result[[i]] <- do.call(reg.fun,reg.args)
}
names(result) <- dimnames(mimf)[[2]]
class(result) <- "emdr2"
return(result)
}
#' Extract results from an emdr2 object
#'
#' \code{extract.emdr2} extracts particular elements from an \code{emdr2} object
#' such as performance criteria or residuals.
#' \code{coef.emdr2} specifically extract regression coefficients.
#'
#' @param object An \code{emdr2} object.
#' @param what A character giving the name of the element to extract.
#' @param select A character of numeric vector giving a subset of MIMFs
#' for which to extract the element.
#'
#' @return \code{extract.emdr2}: a vector, matrix or list of the extracted
#' elements. The class of the output depend on the type of element.
#'
#' @seealso \code{\link{emdr2}} to produce an \code{emdr2} object.
#'
#' @examples
#' library(dlnm)
#'
#' dat <- chicagoNMMAPS[,c("death", "temp", "rhum")]
#'
#' mimfs <- memd(dat)
#' cmimfs <- combine.mimf(mimfs, list(12:13, 14:17, 18:19),
#' new.names = c("C12", "C13", "r"))
#'
#' lm.R2 <- emdr2(death ~ temp + rhum, mimf = cmimfs)
#' betas.R2 <- coef(lm.R2)
#' aic.R2 <- extract.emdr2(lm.R2, what = "aic")
#'
#' @export
extract.emdr2 <- function(object, what, select = NULL)
{
if (missing(what)) stop("Argument 'what' missing")
if (!(what %in% ls(object[[1]]))) stop("Invalid 'what' argument")
if (is.null(select)) select <- 1:length(object)
vals.list <- sapply(object[select],"[[",what)
return(vals.list)
}
#' @rdname extract.emdr2
#'
#' @param method The coefficient method linked to the regression function
#' used. The default should cover most cases.
#' @param ... Additional arguments to be passed to the \code{method}.
#'
#' @return \code{coef.emdr2}: a nimfs x nvariable matrix of coefficients.
#'
#' @export
coef.emdr2 <- function(object, method = "coef", ...)
{
dots <- list(...)
K <- length(object)
coef.list <- vector("list",K)
for (i in 1:K){
met.args <- c(list(object[[i]]), dots)
coef.list[[i]] <- drop(as.matrix(do.call(method, met.args)))
}
coef.mat <- do.call(rbind,coef.list)
rownames(coef.mat) <- names(object)
colnames(coef.mat) <- names(coef.list[[1]])
return(coef.mat)
}
#' Predictions from emdr2
#'
#' Method predict for \code{emdr2} objects.
#'
#' @param object An \code{emdr2} object.
#' @param newmimfs Array of new MIMFs values for prediction. Number of IMFs and
#' variables must match those of the \code{mimf} argument given in the
#' \code{\link{emdr2}} function. If NULL, the data used to fit the model
#' are used.
#' @param newcovariates Data.frame of new covariates value for prediction.
#' @param method The name of a prediction method linked to the regression
#' function used for \code{emdr2}.
#' @param ... Additional arguments for the predict method.
#'
#' @details If new data are given for prediction, both \code{newmimfs} and
#' \code{newcovariates} must be given (provided covariates have been used in
#' the model).
#'
#' @return A vector of predicted values.
#'
#' @examples
#' library(dlnm)
#'
#' dat <- chicagoNMMAPS[,c("death", "temp", "rhum")]
#'
#' mimfs <- memd(dat)
#' cmimfs <- combine.mimf(mimfs, list(12:13, 14:17, 18:19),
#' new.names = c("C12", "C13", "r"))
#'
#' lm.R2 <- emdr2(death ~ temp + rhum, mimf = cmimfs)
#' deathHat <- predict(lm.R2)
#'
#' @export
predict.emdr2 <- function(object, newmimfs = NULL, newcovariates = NULL, method = "predict", ...)
{
dots <- list(...)
K <- length(object)
pred.list <- vector("list",K)
if (length(dim(newmimfs)) == 2) dim(newmimfs) <- c(dim(newmimfs),1)
for (i in 1:K){
met.args <- list(object[[i]])
if (!is.null(newmimfs)) met.args <- c(met.args, list(cbind(newmimfs[,i,], newcovariates)))
met.args <- c(met.args, dots)
pred.list[[i]] <- drop(do.call(method, met.args))
}
yhat <- Reduce("+",pred.list)
return(yhat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.