##' Apply function for all generated datasets.
##'
##' @title Lapply for generated datasets
##' @param datasets Datasets list to be used to generate datasets.
##' @param FUN the function to be applied to each generated dataset. The
##' function will be called like \code{FUN(data, ...)}.
##' @param ... optional arguments to \code{FUN}.
##' @param label optional parameter, if present, each result is added an
##' attribute named \emph{label} with the value of \code{label}.
##' @param POST_FUN function to be applied to the result of \code{FUN}. While
##' one could just modify \code{FUN} instead, this additional argument makes
##' it a bit easier to combine different kinds of methods together.
##' @param datasetIndices optional vector of dataset indices to fit, useful to
##' try only a few datasets instead of all of them. Use \code{"all"} to
##' process all datasets (default).
##' @author Manuel Koller
##' @return list of results. The items in the resulting list will have two
##' additional attributes: \code{datasetIndex} and \code{proc.time}. If
##' \code{FUN} failed for an item, then the item will be the error as
##' returned by try, i.e., it ill be of class \code{try-error}.
##' @examples
##' oneWay <- generateAnovaDatasets(2, 1, 5, 4)
##' lapplyDatasets(oneWay, function(data) sum(data$y))
##' lapplyDatasets(oneWay, function(data) sum(data$y), POST_FUN = function(x) x^2)
##' @export
lapplyDatasets <- function(datasets,
FUN,
...,
label,
POST_FUN,
datasetIndices = "all") {
haveLabel <- !missing(label)
if (missing(POST_FUN)) {
POST_FUN <- identity
}
datasetIndices <- checkDatasetIndices(datasetIndices, datasets)
return(lapply(datasetIndices, function(datasetIndex) {
data <- datasets$generateData(datasetIndex)
warnings <- NULL
time <-
system.time(result <- withCallingHandlers(
try(POST_FUN(FUN(data, ...)))
,
warning = function(w) {
warnings <<- append(warnings, list(w))
invokeRestart("muffleWarning")
}
))
attr(result, "datasetIndex") <- attr(data, "datasetIndex")
attr(result, "proc.time") <- time
if (!is.null(warnings)) {
attr(result, "warnings") <- warnings
}
if (haveLabel) {
attr(result, "label") <- label
}
return(result)
}))
}
checkDatasetIndices <- function(datasetIndices, datasets) {
if (length(datasetIndices) == 1L) {
if (datasetIndices == "all") {
datasetIndices <- 1:datasets[["numberOfDatasets"]]
} else {
datasetIndices <- as.integer(datasetIndices)
}
} else if (!is.integer(datasetIndices)) {
stop("Expected a vector of integers for 'datasetIndices' but got ",
datasetIndices)
}
return(datasetIndices)
}
##' Methods to fit various mixed effects estimators to all generated datasets.
##'
##' @details Existing fitting functions are:
##' @title Fitting Functions
##' @param datasets Datasets list to be used to generate datasets.
##' @param control a list (of correct class for the respective fitting function)
##' containing control parameters to be passed through.
##' @param label a string used to identify which fits have been created by which
##' function.
##' @param postFit a function, taking one argument, the resulting fit. This
##' makes it easy to add an additional step after fitting.
##' @param datasetIndices optional vector of dataset indices to fit, useful to
##' try only a few datasets instead of all of them.
##' @author Manuel Koller
##' @return list of fitted models. See also \code{\link{lapplyDatasets}} which
##' is called internally.
##' @examples
##' set.seed(1)
##' oneWay <- generateAnovaDatasets(1, 1, 10, 4,
##' lmeFormula = y ~ 1,
##' heavyLmeRandom = ~ 1,
##' heavyLmeGroups = ~ Var2,
##' lqmmRandom = ~ 1,
##' lqmmGroup = "Var2",
##' groups = cbind(rep(1:4, each = 10), rep(1:10, 4)),
##' varcov = matrix(1, 4, 4),
##' lower = 0)
##' fitDatasets_lmer(oneWay)
##' @aliases fitDatasets
##' @details \code{fitDatasets_lmer}: Fits datasets using \code{\link{lmer}}
##' using its default options.
##' @rdname fitDatasets
##' @export
##' @importFrom lmer lmerControl
fitDatasets_lmer <-
function(datasets,
control,
label,
postFit,
datasetIndices = "all") {
if (missing(control)) {
control <- lmerControl()
}
if (missing(label)) {
label <- "fitDatasets_lmer"
}
fun <- function(data) {
return(eval(substitute(
lmer(formula, data = data, control = control),
list(formula = datasets[["formula"]])
)))
}
return(
lapplyDatasets(
datasets,
fun,
label = label,
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_lmer_bobyqa}: Fits datasets using \code{\link{lmer}} using
##' the \code{bobyqa} optimizer.
##' @rdname fitDatasets
##' @export
fitDatasets_lmer_bobyqa <-
function(datasets, postFit, datasetIndices = "all") {
return(
fitDatasets_lmer(
datasets,
lmerControl(optimizer = "bobyqa"),
"fitDatasets_lmer_bobyqa",
postFit,
datasetIndices
)
)
}
##' @details \code{fitDatasets_lmer_Nelder_Mead}: Fits datasets using
##' \code{\link{lmer}} using the \code{Nelder Mead} optimizer.
##' @rdname fitDatasets
##' @export
fitDatasets_lmer_Nelder_Mead <-
function(datasets, postFit, datasetIndices = "all") {
return(
fitDatasets_lmer(
datasets,
lmerControl(optimizer = "Nelder_Mead"),
"fitDatasets_lmer_Nelder_Mead",
postFit,
datasetIndices
)
)
}
##' @importFrom Matrix isDiagonal
isDiagonalDataset <- function(datasets) {
data <- datasets$generateData(1)
Lambda <- datasets$createLambdaMatrix(data)
return(isDiagonal(Lambda))
}
##' @param method argument passed on to \code{\link{rlmer}}.
##' @param tuningParameter argument passed on to
##' \code{\link{extractTuningParameter}}.
##' @param init optional argument passed on to \code{\link{rlmer}}.
##' @param ... argument passed on to \code{\link{createRhoFunction}}.
##' @details \code{fitDatasets_rlmer}: Fits datasets using \code{\link{rlmer}}
##' using a custom configuration. The argument 'tuningParameter' is passed to
##' \code{\link{extractTuningParameter}}, details are documented there.
##' @examples
##' ## call rlmer with custom arguments
##' fitDatasets_rlmer_custom <- function(datasets) {
##' return(fitDatasets_rlmer(datasets,
##' method = "DASvar",
##' tuningParameter = c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14),
##' label = "fitDatasets_rlmer_custom"))
##' }
##' fitDatasets_rlmer_custom(oneWay)
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer <-
function(datasets,
method,
tuningParameter,
label,
postFit,
datasetIndices = "all",
...,
init) {
if (missing(method)) {
method <- "DAStau"
}
if (missing(tuningParameter)) {
tuningParameter <- c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14)
}
if (missing(label)) {
label <- "fitDatasets_rlmer_DAStau"
}
rho.e <- createRhoFunction(tuningParameter, "rho.e", ...)
rho.sigma.e <-
createRhoFunction(tuningParameter, "rho.sigma.e", ...)
if (isDiagonalDataset(datasets)) {
rho.b <- createRhoFunction(tuningParameter, "rho.b.diagonal", ...)
rho.sigma.b <-
createRhoFunction(tuningParameter, "rho.sigma.b.diagonal", ...)
} else {
rho.b <-
createRhoFunction(tuningParameter, "rho.b.blockDiagonal", ...)
rho.sigma.b <-
createRhoFunction(tuningParameter, "rho.sigma.b.blockDiagonal", ...)
}
if (missing(init)) {
fun <- function(data) {
eval(substitute(
rlmer(
formula,
data = data,
method = method,
rho.e = rho.e,
rho.b = rho.b,
rho.sigma.e = rho.sigma.e,
rho.sigma.b = rho.sigma.b
),
list(formula = datasets[["formula"]])
))
}
} else {
fun <- function(data) {
eval(substitute(
rlmer(
formula,
data = data,
method = method,
rho.e = rho.e,
rho.b = rho.b,
rho.sigma.e = rho.sigma.e,
rho.sigma.b = rho.sigma.b,
init = init
),
list(formula = datasets[["formula"]])
))
}
}
return(
lapplyDatasets(
datasets,
fun,
label = label,
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
##' Convenience function to create rho-functions with custom tuning parameter.
##'
##' 'rho.b.diagonal' denotes the tuning parameter to be used for 'rho.b' for
##' models with diagonal random effects covariance matrix. 'rho.b.blockDiagonal'
##' is the tuning parameter to be used in the block diagonal case, respectively.
##'
##' For arguments \code{rho.sigma.e} (and \code{rho.sigma.b.diagonal}), the
##' Proposal 2 variant of the function specified for \code{rho.e} (and
##' \code{rho.b}) is used.
##' @title Create Rho-Functions With Custom Tuning Parameter
##' @param tuningParameter argument passed on to
##' \code{\link{extractTuningParameter}}. See its documentation for details.
##' @param which string specifiying which tuning parameter should be extracted.
##' @param rho.e \code{\link{PsiFunction}} to be used for \code{rho.e}.
##' @param rho.sigma.e \code{\link{PsiFunction}} to be used for
##' \code{rho.sigma.e}.
##' @param rho.b.diagonal \code{\link{PsiFunction}} to be used for \code{rho.b}
##' for models with diagonal random effects covariance matrix.
##' @param rho.sigma.b.diagonal \code{\link{PsiFunction}} to be used for
##' \code{rho.sigma.b} for models with diagonal random effects covariance
##' matrix.
##' @param rho.b.blockDiagonal \code{\link{PsiFunction}} to be used for
##' \code{rho.b} for models with block-diagonal random effects covariance
##' matrix.
##' @param rho.sigma.b.blockDiagonal \code{\link{PsiFunction}} to be used for
##' \code{rho.sigma.b} for models with block-diagonal random effects
##' covariance matrix.
##' @param ... passed on to \code{\link{chgDefaults}}.
##' @author Manuel Koller
##' @examples
##' createRhoFunction(c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14), "rho.sigma.e")
##' @export
createRhoFunction <-
function(tuningParameter,
which = c(
"rho.e",
"rho.sigma.e",
"rho.b.diagonal",
"rho.sigma.b.diagonal",
"rho.b.blockDiagonal",
"rho.sigma.b.blockDiagonal"
),
rho.e = smoothPsi,
rho.sigma.e = psi2propII(rho.e),
rho.b.diagonal = rho.e,
rho.sigma.b.diagonal = psi2propII(rho.b.diagonal),
rho.b.blockDiagonal = rho.e,
rho.sigma.b.blockDiagonal = rho.b.blockDiagonal,
...) {
which <- match.arg(which)
k <- extractTuningParameter(tuningParameter, which)
psiFunction <- switch (
which,
rho.e = rho.e,
rho.sigma.e = rho.sigma.e,
rho.b.diagonal = rho.b.diagonal,
rho.sigma.b.diagonal = rho.sigma.b.diagonal,
rho.b.blockDiagonal = rho.b.blockDiagonal,
rho.sigma.b.blockDiagonal = rho.sigma.b.blockDiagonal
)
return(chgDefaults(psiFunction, k = k, ...))
}
##' Methods to extract which tuning parameters have been used for fitting
##' models. Use \code{extractTuningParameter} for custom configurations and
##' \code{extractPredefinedTuningParameter} for predefined configurations
##' provided in this package.
##' @title Extract Tuning Parameters Used In Fitting
##' @param tuningParameter vector of tuning parameters. The vector is expected
##' to be of length 6, containing the tuning parameters for rho.e,
##' rho.sigma.e, rho.b.diagonal, rho.sigma.b.diagonal, rho.b.blockDiagonal and
##' rho.sigma.b.blockDiagonal. 'rho.b.diagonal' denotes the tuning parameter
##' to be used for 'rho.b' for models with diagonal random effects covariance
##' matrix. Names are optional.
##' @param which string specifiying which tuning parameter should be extracted.
##' @author Manuel Koller
##' @return scalar tuning parameter
##' @export
extractTuningParameter <-
function(tuningParameter,
which = c(
"rho.e",
"rho.sigma.e",
"rho.b.diagonal",
"rho.sigma.b.diagonal",
"rho.b.blockDiagonal",
"rho.sigma.b.blockDiagonal"
)) {
which <- match.arg(which)
expectedNames <- c(
"rho.e",
"rho.sigma.e",
"rho.b.diagonal",
"rho.sigma.b.diagonal",
"rho.b.blockDiagonal",
"rho.sigma.b.blockDiagonal"
)
if (length(tuningParameter) != 6) {
stop(
"Argument 'tuningParameter' needs to be of length 6. ",
"The contents are assumed to be the tuning parameters for ",
"the rho functions in the order of [",
paste(expectedNames, collapse = ", "),
"]."
)
}
if (is.null(names(tuningParameter))) {
names(tuningParameter) <- expectedNames
} else {
diff <- setdiff(expectedNames, names(tuningParameter))
if (length(diff) > 0) {
stop(
"Expected names of argument 'tuningParameter' to be [",
paste(expectedNames, collapse = ", "),
"] but got [",
paste(names(tuningParameter), collapse = ", "),
"] instead."
)
}
}
return(tuningParameter[[which]])
}
predefinedTuningParameter <- list()
##' @param label label or vector of labels in results. Only predefined labels of
##' the form 'fitDatasets_rlmer_...' are supported (for others NA is
##' returned).
##' @examples
##' extractPredefinedTuningParameter("fitDatasets_rlmer_DAStau", "rho.e")
##' @rdname extractTuningParameter
##' @export
extractPredefinedTuningParameter <- function(label, which) {
extractForOneLabel <- function(oneLabel) {
tuningParameter <- predefinedTuningParameter[[oneLabel]]
if (is.null(tuningParameter)) {
return(NA_real_)
} else {
return(extractTuningParameter(tuningParameter, which))
}
}
return(sapply(label, extractForOneLabel))
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau"]] <-
c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14)
##' @details \code{fitDatasets_rlmer_DAStau}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau and \code{\link{smoothPsi}} for
##' the rho functions. The tuning parameters are k = 1.345 for \code{rho.e}.
##' For \code{rho.sigma.e}, the Proposal 2 variant is used using k = 2.28.
##' The choices for \code{rho.b} and \code{rho.sigma.b} depend on whether the
##' model uses a diagonal or a block diagonal matrix for Lambda. In the former
##' case, the same psi functions and tuning parameters are use as for
##' \code{rho.e} and \code{rho.sigma.b}. In the block diagonal case,
##' \code{rho.b} and \code{rho.sigma.b} both use \code{\link{smoothPsi}} using
##' a tuning parameter k = 5.14 (assuming blocks of dimension 2).
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_rlmer_DAStau_lmerNoFit}: Fits datasets using
##' \code{\link{rlmer}} using the same configuration as
##' \code{fitDatasets_rlmer_DAStau} except for that it is using
##' \code{\link{lmerNoFit}} as initial estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_lmerNoFit <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_lmerNoFit"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[["fitDatasets_rlmer_DAStau"]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices,
init = lmerNoFit
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DASvar"]] <-
c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14)
##' @details \code{fitDatasets_rlmer_DASvar}: Fits datasets using
##' \code{\link{rlmer}} using method DASvar. The same rho functions and tuning
##' parameters are used as for \code{fitDatasets_rlmer_DAStau}.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DASvar <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DASvar"
return(
fitDatasets_rlmer(
datasets,
method = "DASvar",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_noAdj"]] <-
c(1.345, 1.345, 1.345, 1.345, 5.14, 5.14)
##' @details \code{fitDatasets_rlmer_DAStau_noAdj}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. The same rho functions and tuning
##' parameters are used as for \code{fitDatasets_rlmer_DAStau}, except for
##' \code{rho.sigma.e} (and \code{rho.sigma.b} in the diagonal case) for which
##' the Proposal 2 variant of \code{smoothPsi} using k = 1.345 is used.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_noAdj <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_noAdj"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices,
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_0_5"]] <-
c(0.5, 1.47, 0.5, 1.47, 2.17, 2.17)
##' @details \code{fitDatasets_rlmer_DAStau_k_0_5}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 0.5} for \code{rho.e} and
##' \code{k = 1.47} for \code{rho.sigma.e}, the latter adjusted to reach the
##' same asymptotic efficiency. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 2.17} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_0_5 <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_0_5"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_0_5_noAdj"]] <-
c(0.5, 0.5, 0.5, 0.5, 2.17, 2.17)
##' @details \code{fitDatasets_rlmer_DAStau_k_0_5_noAdj}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 0.5} for \code{rho.e} and
##' \code{rho.sigma.e}. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 2.17} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_0_5_noAdj <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_0_5_noAdj"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_2"]] <-
c(2, 2.9, 2, 2.9, 8.44, 8.44)
##' @details \code{fitDatasets_rlmer_DAStau_k_2}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 2} for \code{rho.e} and
##' \code{k = 2.9} \code{rho.sigma.e}, the latter adjusted to reach the same
##' asymptotic efficiency. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 8.44} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_2 <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_2"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_2_noAdj"]] <-
c(2, 2, 2, 2, 8.44, 8.44)
##' @details \code{fitDatasets_rlmer_DAStau_k_2_noAdj}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 2} for \code{rho.e} and
##' \code{rho.sigma.e}. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 8.44} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_2_noAdj <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_2_noAdj"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_5"]] <-
c(5, 5.03, 5, 5.03, 8.44, 8.44)
##' @details \code{fitDatasets_rlmer_DAStau_k_5}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 5} for \code{rho.e} and
##' \code{k = 5.03} \code{rho.sigma.e}, the latter adjusted to reach the same
##' asymptotic efficiency. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 34.21} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_5 <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_5"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
predefinedTuningParameter[["fitDatasets_rlmer_DAStau_k_5_noAdj"]] <-
c(5, 5, 5, 5, 8.44, 8.44)
##' @details \code{fitDatasets_rlmer_DAStau_k_5_noAdj}: Fits datasets using
##' \code{\link{rlmer}} using method DAStau. Use \code{\link{smoothPsi}}
##' psi-function with tuning parameter \code{k = 5} for \code{rho.e} and
##' \code{rho.sigma.e}. In the diagonal case, the same are used for
##' \code{rho.b} and \code{rho.sigma.b} as well. In the block-diagonal case,
##' the tuning parameter \code{k = 34.21} is used for \code{rho.b} and
##' \code{rho.sigma.b}. The tuning parameter is chosen to reach about the same
##' asymptotic efficiency for theta as for the fixed effects.
##' @rdname fitDatasets
##' @export
fitDatasets_rlmer_DAStau_k_5_noAdj <-
function(datasets, postFit, datasetIndices = "all") {
label <- "fitDatasets_rlmer_DAStau_k_5_noAdj"
return(
fitDatasets_rlmer(
datasets,
method = "DAStau",
tuningParameter = predefinedTuningParameter[[label]],
label = label,
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_heavyLme}: Fits datasets using
##' \code{heavyLme} from package \code{heavy}. Additional
##' required arguments are: \code{lmeFormula}, \code{heavyLmeRandom} and
##' \code{heavyLmeGroups}. They are passed to the \code{formula},
##' \code{random} and \code{groups} arguments of \code{heavyLme}.
##' @rdname fitDatasets
##' @export
fitDatasets_heavyLme <-
function(datasets, postFit, datasetIndices = "all") {
packageInstalled <- isPackageInstalled("heavy")
fun <- function(data) {
if (!packageInstalled) {
return(createPackageMissingReturnValue("heavy", "fitDatasets_heavyLme"))
}
fun <- eval(parse(text="heavy::heavyLme"))
eval(substitute(
fun(
formula,
random = random,
groups = groups,
data = data
),
list(
formula = datasets[["lmeFormula"]],
random = datasets[["heavyLmeRandom"]],
groups = datasets[["heavyLmeGroups"]]
)
))
}
return(
lapplyDatasets(
datasets,
fun,
label = "fitDatasets_heavyLme",
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
isPackageInstalled <- function(package) {
return(length(find.package(package, quiet = TRUE)) > 0)
}
createPackageMissingReturnValue <- function(package, functionName) {
value <-
paste0(
"Package '",
package,
"' is not installed. Please install it to run '",
functionName,
"'."
)
attr(value, "class") <- "package-not-installed"
return(value)
}
##' @details \code{fitDatasets_lqmm}: Fits datasets using
##' \code{\link[lqmm]{lqmm}} from package \code{lqmm}. Additional required
##' arguments are: \code{lmeFormula}, \code{lqmmRandom}, \code{lqmmGroup} and
##' \code{lqmmCovariance}. They are passed to the \code{formula},
##' \code{random}, \code{groups} and \code{covariance} arguments of
##' \code{lqmm}. \code{lqmmCovariance} is optional, if omitted \code{pdDiag}
##' is used.
##' @rdname fitDatasets
##' @export
fitDatasets_lqmm <-
function(datasets, postFit, datasetIndices = "all") {
packageInstalled <- isPackageInstalled("lqmm")
covariance <- datasets[["lqmmCovariance"]]
if (is.null(covariance))
covariance <- "pdDiag"
fun <- function(data) {
if (!packageInstalled) {
return(createPackageMissingReturnValue("lqmm", "fitDatasets_lqmm"))
}
eval(substitute(
lqmm::lqmm(
formula,
random = random,
group = group,
data = data,
covariance = covariance
),
list(
formula = datasets[["lmeFormula"]],
random = datasets[["lqmmRandom"]],
group = datasets[["lqmmGroup"]],
covariance = covariance
)
))
}
return(
lapplyDatasets(
datasets,
fun,
label = "fitDatasets_lqmm",
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_rlme}: Fits datasets using
##' \code{\link[rlme]{rlme}} from package \code{rlme}.
##' @rdname fitDatasets
##' @export
fitDatasets_rlme <-
function(datasets, postFit, datasetIndices = "all") {
packageInstalled <- isPackageInstalled("rlme")
fun <- function(data) {
if (!packageInstalled) {
return(createPackageMissingReturnValue("rlme", "fitDatasets_rlme"))
}
return(rlme::rlme(datasets[["formula"]], data))
}
return(
lapplyDatasets(
datasets,
fun,
label = "fitDatasets_rlme",
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob}: Prototype method to fit datasets
##' using \code{varComprob} from package
##' \code{robustvarComp}. Additional required items in \code{datasets} are:
##' \code{lmeFormula}, \code{groups}, \code{varcov} and \code{lower}. They are
##' passed to the \code{fixed}, \code{groups}, \code{varcov} and \code{lower}
##' arguments of \code{varComprob}. The running of this method produces many
##' warnings of the form "passing a char vector to .Fortran is not portable"
##' which are suppressed.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob <-
function(datasets,
control,
label,
postFit,
datasetIndices = "all") {
packageInstalled <- isPackageInstalled("robustvarComp")
fun <- function(data) {
if (!packageInstalled) {
return(createPackageMissingReturnValue("robustvarComp", label))
}
warnings <- NULL
set.seed(attr(data, "datasetIndex"))
result <- withCallingHandlers(
eval(substitute(
robustvarComp::varComprob(
fixed = formula,
groups = groups,
data = data,
varcov = varcov,
control = control,
model = FALSE,
X = FALSE,
Y = FALSE,
K = FALSE
),
list(
formula = datasets[["lmeFormula"]],
groups = datasets[["groups"]],
varcov = datasets[["varcov"]]
)
)),
warning = function(w) {
warnings <<- append(warnings, list(w))
invokeRestart("muffleWarning")
}
)
for (condition in warnings) {
if (!grepl(
"passing a char vector to .Fortran is not portable",
conditionMessage(condition),
fixed = TRUE
)) {
warning(condition)
}
}
return(result)
}
return(
lapplyDatasets(
datasets,
fun,
label = label,
POST_FUN = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_compositeTau}: Fits datasets with the
##' composite Tau method using \code{varComprob} from
##' package \code{robustvarComp}. See \code{fitDatasets_varComprob} for
##' additional details.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeTau <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]])
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeTau",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
varComrob.control.or.null <- function(...) {
if (isPackageInstalled("robustvarComp")) {
return(robustvarComp::varComprob.control(...))
} else {
return(NULL)
}
}
##' @details \code{fitDatasets_varComprob_compositeTau_OGK}: Similar to
##' \code{fitDatasets_varComprob_compositeTau} but using \code{covOGK} as initial
##' covariance matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeTau_OGK <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
cov.init = "covOGK")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeTau_OGK",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_compositeTau_2SGS}: Similar to
##' \code{fitDatasets_varComprob_compositeTau} but using \code{2SGS} as initial covariance
##' matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeTau_2SGS <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
cov.init = "2SGS")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeTau_2SGS",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_compositeS}: Similar to
##' \code{fitDatasets_varComprob_compositeTau} but using method composite S.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeS <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
method = "compositeS")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeS",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_compositeS_OGK}: Similar to
##' \code{fitDatasets_varComprob_compositeS} but using \code{covOGK} as
##' initial covariance matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeS_OGK <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
method = "compositeS",
cov.init = "covOGK")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeS_OGK",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_compositeS_2SGS}: Similar to
##' \code{fitDatasets_varComprob_compositeS} but using \code{2SGS} as initial
##' covariance matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_compositeS_2SGS <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
method = "compositeS",
cov.init = "2SGS")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_compositeS_2SGS",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_S}: Similar to
##' \code{fitDatasets_varComprob_compositeTau} but using method S and the
##' Rocke psi-function.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_S <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(lower = datasets[["lower"]],
method = "S",
psi = "rocke")
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_S",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_S_OGK}: Similar to
##' \code{fitDatasets_varComprob_S} but using \code{covOGK} as initial
##' covariance matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_S_OGK <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(
lower = datasets[["lower"]],
method = "S",
psi = "rocke",
cov.init = "covOGK"
)
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_S_OGK",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' @details \code{fitDatasets_varComprob_S_2SGS}: Similar to
##' \code{fitDatasets_varComprob_S} but using \code{2SGS} as initial
##' covariance matrix estimator.
##' @rdname fitDatasets
##' @export
fitDatasets_varComprob_S_2SGS <-
function(datasets, postFit, datasetIndices = "all") {
lcontrol <-
varComrob.control.or.null(
lower = datasets[["lower"]],
method = "S",
psi = "rocke",
cov.init = "2SGS"
)
return(
fitDatasets_varComprob(
datasets,
lcontrol,
"fitDatasets_varComprob_S_2SGS",
postFit = postFit,
datasetIndices = datasetIndices
)
)
}
##' Shorten labels created by the various \code{\link{fitDatasets}} functions,
##' for use in plotting, etc.
##'
##' The labels are shortened as they are in the simulation study published in
##' Koller and Stahel (2022).
##' @title Shorten Labels
##' @param labels vector of labels as assigned by \code{\link{fitDatasets}}
##' @author Manuel Koller
##' @return Vector of shortened labels
##' @references Koller M, Stahel WA (2022). "Robust Estimation of General Linear
##' Mixed Effects Models.” In PM Yi, PK Nordhausen (eds.), Robust and Multivariate
##' Statistical Methods, Springer Nature Switzerland AG.
##' @examples
##' labels <- c("fitDatasets_lmer", "fitDatasets_rlmer_DAStau",
##' "fitDatasets_rlmer_DAStau_noAdj",
##' "fitDatasets_varComprob_compositeTau_OGK",
##' "fitDatasets_varComprob_S_OGK",
##' "fitDatasets_heavyLme",
##' "fitDatasets_lqmm")
##' shortenLabelsKS2022(labels)
##' @export
shortenLabelsKS2022 <- function(labels) {
idx <- labels == "fitDatasets_lmer"
labels[idx] <- "lme"
idx <- labels == "fitDatasets_rlmer_DAStau"
labels[idx] <- "RSEa"
idx <- labels == "fitDatasets_rlmer_DAStau_noAdj"
labels[idx] <- "RSEn"
idx <- grepl("fitDatasets_varComprob_compositeTau", labels)
labels[idx] <- "cTau"
idx <- grepl("fitDatasets_varComprob_S", labels)
labels[idx] <- "S"
idx <- grepl("fitDatasets_heavyLme", labels)
labels[idx] <- "t4"
idx <- grepl("fitDatasets_lqmm", labels)
labels[idx] <- "med"
return(labels)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.