Nothing
R/anova.R
In LMMstar: Repeated Measurement Models for Discrete Times
Defines functions
.simplifyContrast
dfSigma
anova.mlmm
.anova_LRT
.anova_Wald
anova.lmm
Documented in
anova.lmm
### anova.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: mar 5 2021 (21:38)
## Version:
## Last-Updated: aug 1 2023 (15:45)
## By: Brice Ozenne
## Update #: 1346
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * anova.lmm (documentation)
##' @title Multivariate Tests For Linear Mixed Model
##' @description Simultaneous tests of linear combinations of the model paramaters using Wald tests or Likelihood Ratio Test (LRT).
##'
##' @param object a \code{lmm} object. Only relevant for the anova function.
##' @param effects [character or numeric matrix] Should the Wald test be computed for all variables (\code{"all"}),
##' or only variables relative to the mean (\code{"mean"} or \code{"fixed"}),
##' or only variables relative to the variance structure (\code{"variance"}),
##' or only variables relative to the correlation structure (\code{"correlation"}).
##' Can also be use to specify linear combinations of coefficients or a contrast matrix, similarly to the \code{linfct} argument of the \code{multcomp::glht} function.
##' @param robust [logical] Should robust standard errors (aka sandwich estimator) be output instead of the model-based standard errors.
##' @param rhs [numeric vector] the right hand side of the hypothesis. Only used when the argument effects is a matrix.
##' @param df [logical] Should a F-distribution be used to model the distribution of the Wald statistic. Otherwise a chi-squared distribution is used.
##' @param ci [logical] Should an estimate, standard error, confidence interval, and p-value be output for each hypothesis?
##' @param transform.sigma,transform.k,transform.rho,transform.names are passed to the \code{vcov} method. See details section in \code{\link{coef.lmm}}.
##' @param ... Not used. For compatibility with the generic method.
##'
##' @return A data.frame (LRT) or a list of containing the following elements (Wald):\itemize{
##' \item \code{multivariate}: data.frame containing the multivariate Wald test.
##' The column \code{df.num} refers to the degrees of freedom for the numerator (i.e. number of hypotheses)
##' wherease the column \code{df.denum} refers to the degrees of freedom for the denominator (i.e. Satterthwaite approximation).
##' \item \code{univariate}: data.frame containing each univariate Wald test.
##' \item \code{glht}: used internally to call functions from the multcomp package.
##' \item \code{object}: list containing key information about the linear mixed model.
##' \item \code{vcov}: variance-covariance matrix associated to each parameter of interest (i.e. hypothesis).
##' \item \code{iid}: matrix containing the influence function relative to each parameter of interest (i.e. hypothesis).
##' \item \code{args}: list containing argument values from the function call.
##' }
##'
##' @details By default adjustment of confidence intervals and p-values for multiple comparisons is based on the distribution of the maximum-statistic.
##' This is refered to as a single-step Dunnett multiple testing procedures in table II of Dmitrienko et al. (2013).
##' It is performed using the multcomp package with the option \code{test = adjusted("single-step")} with equal degrees of freedom
##' or by simulation using a Student's t copula with unequal degrees of freedom (more in the note of the details section of \code{\link{confint.Wald_lmm}}).
##'
##' @seealso
##' \code{\link{summary.Wald_lmm}} or \code{\link{confint.Wald_lmm}} for a summary of the results. \cr
##' \code{\link{autoplot.Wald_lmm}} for a graphical display of the results. \cr
##' \code{\link{rbind.Wald_lmm}} for combining result across models and adjust for multiple comparisons. \cr
##'
##' @references Dmitrienko, A. and D'Agostino, R., Sr (2013), Traditional multiplicity adjustment methods in clinical trials. Statist. Med., 32: 5172-5218. https://doi.org/10.1002/sim.5990.
##'
##' @keywords htest
##'
##' @examples
##' #### simulate data in the long format ####
##' set.seed(10)
##' dL <- sampleRem(100, n.times = 3, format = "long")
##'
##' #### fit Linear Mixed Model ####
##' eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5,
##' repetition = ~visit|id, structure = "UN", data = dL)
##'
##' #### Multivariate Wald test ####
##' ## F-tests
##' anova(eUN.lmm)
##' anova(eUN.lmm, effects = "all")
##' anova(eUN.lmm, robust = TRUE, df = FALSE)
##' summary(anova(eUN.lmm))
##'
##' ## user defined F-test
##' summary(anova(eUN.lmm, effects = c("X1=0","X2+X5=10")))
##' print(anova(eUN.lmm, effects = "mean_visit"), columns = add("null"))
##'
##' ## chi2-tests
##' anova(eUN.lmm, df = FALSE)
##'
##'
##' ## with standard contrast
##' if(require(multcomp)){
##' amod <- lmm(breaks ~ tension, data = warpbreaks)
##' e.amod <- anova(amod, effect = mcp(tension = "Tukey"))
##' summary(e.amod)
##' }
##'
##' #### Likelihood ratio test ####
##' eUN0.lmm <- lmm(Y ~ X1 + X2, repetition = ~visit|id, structure = "UN", data = dL)
##' anova(eUN.lmm, eUN0.lmm)
##'
##' eCS.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "CS", data = dL)
##' anova(eUN.lmm, eCS.lmm)
## * anova.lmm (code)
##' @export
anova.lmm <- function(object, effects = NULL, robust = FALSE, rhs = NULL, df = !is.null(object$df), ci = TRUE,
transform.sigma = NULL, transform.k = NULL, transform.rho = NULL, transform.names = TRUE, ...){
call <- match.call()
## ** normalized user input
dots <- list(...)
options <- LMMstar.options()
if(length(dots)>0){
stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n")
}
if(is.null(effects)){
effects <- options$effects
}
if(inherits(effects,"lmm")){ ## likelihood ratio test
out <- .anova_LRT(object1 = object, object2 = effects)
attr(out,"call") <- call
if(!inherits(out,"LRT_lmm")){ ## in the case of re-estimation of the model via ML no need to re-assign the class
class(out) <- append("LRT_lmm",class(out))
}
}else{ ## Wald test
attr(robust, "call") <- "robust" %in% names(call)
out <- .anova_Wald(object, effects = effects, robust = robust, rhs = rhs, df = df, ci = ci,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names)
attr(out,"call") <- call
class(out) <- append("Wald_lmm",class(out))
}
## ** export
return(out)
}
## * .anova_Wald
.anova_Wald <- function(object, effects, robust, rhs, df, ci,
transform.sigma, transform.k, transform.rho, transform.names){
options <- LMMstar.options()
type.information <- object$args$type.information
original.transform.sigma <- transform.sigma
original.transform.k <- transform.k
original.transform.rho <- transform.rho
## ** normalized user input
terms.mean <- attr(stats::terms(object$formula$mean.design),"term.labels")
subeffect <- NULL
if(!inherits(effects,"mcp") && length(effects)==1){
if(identical(effects,"all")){
effects <- c("mean","variance","correlation")
}else if(grepl("^mean_",effects)){
iLabels <- attr(stats::terms(object$formula$mean.design),"term.labels")
if(any(effects == paste0("mean_",iLabels))){
subeffect <- iLabels[effects == paste0("mean_",iLabels)]
effects <- "mean"
}
}else if(grepl("^variance_",effects) && !is.null(object$design$vcov$var$X)){
iLabels <- attr(stats::terms(object$formula$var.design),"term.labels")
if(any(effects == paste0("variance_",iLabels))){
subeffect <- iLabels[effects == paste0("variance_",iLabels)]
effects <- "variance"
}
}else if(grepl("^cor_",effects) && !is.null(object$design$vcov$cor$X)){
iLabels <- attr(stats::terms(object$formula$cor.design),"term.labels")
if(any(effects == paste0("correlation_",iLabels))){
subeffect <- iLabels[effects == paste0("correlation_",iLabels)]
effects <- "correlation"
}
}
}
init <- .init_transform(transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho,
x.transform.sigma = object$reparametrize$transform.sigma, x.transform.k = object$reparametrize$transform.k, x.transform.rho = object$reparametrize$transform.rho)
transform.sigma <- init$transform.sigma
transform.k <- init$transform.k
transform.rho <- init$transform.rho
name.coef <- names(stats::coef(object, effects = "all"))
out <- list(multivariate = NULL, univariate = NULL, glht = NULL)
if(inherits(effects,"mcp")){
out.glht <- try(multcomp::glht(object, linfct = effects), ## only used for generating contrast matrix
silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
type <- "all"
if(length(names(effects))==1){
ls.nameTerms <- list(all = names(effects))
}else{
ls.nameTerms <- list(all = NULL)
}
ls.nameTerms.num <- list(all = 1)
ls.contrast <- list(all = matrix(0, nrow = NROW(out.glht$linfct), ncol = length(name.coef), dimnames = list(rownames(out.glht$linfct),name.coef)))
ls.contrast$all[,colnames(out.glht$linfct)] <- out.glht$linfct
ls.null <- list(all = out.glht$rhs)
name.effects <- NULL
}else if(is.matrix(effects)){
## try to re-size the matrix if necessary
if(NCOL(effects)!=length(name.coef)){
if(is.null(colnames(effects))){
stop("Argument \'effect\' should have column names when a matrix. \n")
}
if(any(duplicated(colnames(effects)))){
stop("Argument \'effect\' should not have duplicated column names when a matrix. \n")
}
if(any(colnames(effects) %in% name.coef == FALSE)){
stop("Argument \'effect\' should not have column names matching the coefficient names when a matrix. \n")
}
effects.save <- effects
effects <- matrix(0, nrow = NROW(effects.save), ncol = length(name.coef), dimnames = list(rownames(effects.save),name.coef))
effects[,colnames(effects.save)] <- effects.save
}
if(is.null(rhs)){
rhs <- rep(0, NROW(effects))
}
## run glht
out.glht <- try(multcomp::glht(object, linfct = effects, rhs = rhs, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = FALSE, effects = "all")}),
silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
type <- "all"
ls.nameTerms <- list(all = NULL)
ls.nameTerms.num <- list(all = 1)
ls.contrast <- list(all = out.glht$linfct)
ls.null <- list(all = out.glht$rhs)
name.effects <- rownames(effects)
}else if(all(tolower(effects) %in% c("mean","fixed","variance","correlation"))){
if(transform.k %in% c("sd","var","logsd","logvar")){
stop("Cannot use \'transform.rho\' equal \"sd\", \"var\", \"logsd\", or \"logvar\". \n",
"anova does not handle tests where the null hypothesis is at a boundary of the support of a random variable. \n")
}
if(transform.rho %in% c("cov")){
stop("Cannot use \'transform.rho\' equal \"cov\". \n",
"anova does not handle tests where the null hypothesis is at a boundary of the support of a random variable. \n")
}
effects <- match.arg(effects, c("mean","fixed","variance","correlation"), several.ok = TRUE)
effects[effects=="fixed"] <- "mean"
name.effects <- NULL
type <- NULL
ls.assign <- list()
ls.nameTerms <- list()
ls.contrast <- list()
ls.null <- list()
if("mean" %in% effects){
type <- c(type, "mean")
ls.assign$mean <- attr(object$design$mean,"assign")
ls.nameTerms$mean <- attr(stats::terms(object$formula$mean.design),"term.labels")
ls.contrast <- c(ls.contrast,list(mean = NULL))
null.mean <- 0
ls.null$mean <- rep(null.mean,length(ls.nameTerms$mean))
}
if("variance" %in% effects){
type <- c(type, "variance")
ls.assign$variance <- attr(object$design$vcov$var$X,"assign")
ls.nameTerms$variance <- attr(stats::terms(object$formula$var),"term.labels")
if(!is.na(object$design$vcov$name$strata)){
ls.assign$variance[ls.nameTerms$variance[ls.assign$variance]==object$design$vcov$name$strata] <- 0
}
ls.contrast <- c(ls.contrast,list(variance = NULL))
null.variance <- switch(transform.k,
"none" = 1,
"square" = 1,
"log" = 0,
"logsquare" = 0)
ls.null$variance <- rep(null.variance,length(ls.nameTerms$variance))
}
if("correlation" %in% effects){
type <- c(type, "correlation")
ls.assign$correlation <- rep(1,sum(object$design$param$type=="rho"))
ls.nameTerms$correlation <- if(!is.null(ls.assign$correlation)){object$time$var}else{NULL}
ls.contrast <- c(ls.contrast,list(correlation = NULL))
null.correlation <- switch(transform.rho,
"none" = 0,
"atanh" = 0,
"cov" = 0)
ls.null$correlation <- rep(null.correlation,length(ls.nameTerms$correlation))
}
ls.nameTerms.num <- lapply(ls.nameTerms, function(iName){as.numeric(factor(iName, levels = iName))})
}else if(all(grepl("=",effects)==FALSE)){
stop("Incorrect argument \'effects\': can be \"mean\", \"variance\", \"correlation\", \"all\", \n",
"or something compatible with the argument \'linfct\' of multcomp::glht. \n ")
}else{ ## symbolic definition of effects using equations (characters)
## run glht
out.glht <- try(multcomp::glht(object, linfct = effects, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = FALSE, effects = "all")}),
silent = TRUE)
newname.coef <- names(stats::coef(object, effects = "all"))
if(inherits(out.glht,"try-error")){
## maybe the error is due to white space in the name of the coefficients
if(any(grepl(" ",effects))){
lhs <- sapply(strsplit(effects, split = "=",fixed = TRUE),"[",1)
lhs.term <- unlist(strsplit(unlist(strsplit(lhs, split = "+", fixed = TRUE)), split = "-", fixed = TRUE))
lhs.variable <- trimws(lapply(strsplit(lhs.term, split = "*", fixed = TRUE), utils::tail, 1), which = "both")
if(any(grepl(" ",lhs.variable))){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht,
"There seems to be whitespace(s) in the name of certain coefficients which is likely to be the cause of the error. \n",
"Consider using a contrast matrix to specific the argument \'effects\' \n or renaming the values of categorical variables without whitespace and refiting the lmm object. \n"
)
}
}
test.reparametrize <- any(grepl("log", c(object$reparametrize$transform.sigma,object$reparametrize$transform.k))) || grep("atanh", object$reparametrize$transform.rho)
## restaure untransformed parametrization (glht does not handle log(k). or atanh(cor))
if(test.reparametrize){
object2 <- object
index.var <- which(object$design$param$type %in% c("sigma","k","rho"))
object2$reparametrize <- .reparametrize(p = object$param[index.var],
type = object$design$param$type[index.var],
level = object$design$param$level[index.var],
sigma = object$design$param$sigma[index.var],
k.x = object$design$param$k.x[index.var],
k.y = object$design$param$k.y[index.var],
Jacobian = FALSE, dJacobian = FALSE, inverse = FALSE,
transform.sigma = "none",
transform.k = "none",
transform.rho = "none",
transform.names = TRUE)
out.glht <- try(multcomp::glht(object2, linfct = effects, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = robust, effects = "all")}), silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
oldname.coef <- colnames(out.glht$linfct)
newname.hypo <- rownames(out.glht$linfct)
for(iSub in 1:length(oldname.coef)){
newname.hypo <- gsub(pattern = oldname.coef[iSub], replacement = newname.coef[iSub], x = newname.hypo, fixed = TRUE)
}
dimnames(out.glht$linfct) <- list(newname.hypo,newname.coef)
message("The coefficient names have been transformed but not the null hypotheses. \n")
}else{
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
}
type <- "all"
ls.contrast <- list(all = out.glht$linfct)
ls.null <- list(all = out.glht$rhs)
name.effects <- names(effects)
if(!is.null(name.effects)){
rownames(ls.contrast$all) <- name.effects
}
if(length(effects)==1){
ls.nameTerms <- list(all = name.effects)
}else{
ls.nameTerms <- list(all = NULL)
}
ls.nameTerms.num <- list(all = 1)
}
## ** prepare
if(robust && object$args$method.fit=="REML"){
name.mean <- names(coef(object, effects = "mean"))
if(all(effects %in% c("mean","variance","correlation"))){
if(all(effects %in% c("variance","correlation"))){
stop("Cannot test variance, covariance, or correlation parameters using robust standard errors with REML. \n")
}
}else if(any(names(which(colSums(abs(do.call(rbind,ls.contrast)))>0)) %in% name.mean == FALSE)){
stop("Cannot test variance, covariance, or correlation parameters using robust standard errors with REML. \n")
}
ls.contrast <- lapply(ls.contrast, function(iC){
if(!is.null(iC)){iC[,name.mean,drop=FALSE]}else{iC}
})
effects <- "mean"
}else{
effects <- "all"
}
type.param <- stats::setNames(object$design$param$type,object$design$param$name)
param <- coef(object, effects = effects,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
newparam <- coef(object, effects = effects,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names)
newname <- names(newparam)
name.param <- names(param)
name.paramSigma <- names(type.param)[type.param=="sigma"]
name.paramK <- names(type.param)[type.param=="k"]
name.paramRho <- names(type.param)[type.param=="rho"]
n.param <- length(param)
vcov.param <- vcov(object, df = df*2, effects = effects, robust = robust,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
dVcov.param <- attr(vcov.param,"dVcov")
if(df>0 && object$args$method.fit=="REML" && type.information == "expected"){
warning("when using REML with expected information, the degree of freedom of the F-statistic may depend on the parametrisation of the variance parameters. \n")
}
if(any(sapply(ls.contrast, function(iC){is.null(iC) || identical(colnames(iC), names(param))}) == FALSE)){
warning("Names of the columns of the contrast matrix do not match the names of the model coefficients. \n")
}
## ** F-tests
out$glht <- stats::setNames(vector(mode = "list", length = length(type)), type)
for(iType in type){ ## iType <- "variance"
## skip empty type
if(length(ls.nameTerms.num[[iType]])==0 || (is.null(ls.contrast[[iType]]) && (all(ls.assign[[iType]]==0)))){ next }
iParam <- coef(object, effects = iType,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
name.iParam <- names(iParam)
if(is.null(ls.nameTerms[[iType]])){
ls.nameTerms[[iType]] <- ls.nameTerms.num[[iType]]
}
out$glht[[iType]] <- stats::setNames(vector(mode = "list", length = length(ls.nameTerms.num[[iType]])), ls.nameTerms[[iType]])
for(iTerm in ls.nameTerms.num[[iType]]){ ## iTerm <- 1
iNameTerm <- ls.nameTerms[[iType]][[which(iTerm == ls.nameTerms.num[[iType]])]]
## *** contrast matrix
if(is.null(ls.contrast[[iType]])){
if(all(ls.assign[[iType]]!=iTerm)){next}
if(!is.null(subeffect) && ls.nameTerms$mean[iTerm]!=subeffect){next}
iIndex.param <- which(ls.assign[[iType]]==iTerm)
iN.hypo <- length(iIndex.param)
iNull <- rep(ls.null[[iType]][iTerm],iN.hypo)
iName.hypo <- paste(paste0(name.iParam[iIndex.param],"=",iNull), collapse = ", ")
iC <- matrix(0, nrow = iN.hypo, ncol = n.param, dimnames = list(name.iParam[iIndex.param], name.param))
if(length(iIndex.param)==1){
iC[name.iParam[iIndex.param],name.iParam[iIndex.param]] <- 1
}else{
diag(iC[name.iParam[iIndex.param],name.iParam[iIndex.param]]) <- 1
}
iC.uni <- iC
}else{
iC <- ls.contrast[[iType]]
iC.uni <- iC
iN.hypo <- NROW(iC)
iNull <- ls.null[[iType]]
iName.hypo <- paste(paste0(rownames(iC),"=",iNull),collapse=", ")
}
## *** statistic
outSimp <- .simplifyContrast(iC, iNull) ## remove extra lines
iC.vcov.C_M1 <- try(solve(outSimp$C %*% vcov.param %*% t(outSimp$C)), silent = TRUE)
if(inherits(iC.vcov.C_M1,"try-error")){
iStat <- NA
iDf <- c(outSimp$dim,Inf)
attr(iStat,"error") <- "\n Could not invert the covariance matrix for the proposed contrast."
}else{
iStat <- as.double(t(outSimp$C %*% param - outSimp$rhs) %*% iC.vcov.C_M1 %*% (outSimp$C %*% param - outSimp$rhs))/outSimp$dim
iDf <- c(outSimp$dim,Inf)
if(outSimp$rm>0){
iName.hypo <- paste(paste0(rownames(outSimp$C),"=",outSimp$rhs),collapse=", ")
}
}
## *** degree of freedom
if(df>0 && !inherits(iC.vcov.C_M1,"try-error")){
svd.tempo <- eigen(iC.vcov.C_M1)
D.svd <- diag(svd.tempo$values, nrow = outSimp$dim, ncol = outSimp$dim)
P.svd <- svd.tempo$vectors
contrast.svd <- sqrt(D.svd) %*% t(P.svd) %*% outSimp$C
colnames(contrast.svd) <- name.param
iNu_m <- dfSigma(contrast = contrast.svd,
vcov = vcov.param,
dVcov = dVcov.param,
keep.param = name.param)
iEQ <- sum(iNu_m/(iNu_m - 2))
iDf[2] <- 2 * iEQ/(iEQ - outSimp$dim)
}
## *** confidence interval
if(ci>0){
if(df>0){
ci.df <- .dfX(X.beta = iC.uni, vcov.param = vcov.param, dVcov.param = dVcov.param, return.vcov = ci>0.5)
}else{
ci.df <- Inf
}
CI <- data.frame(estimate = as.double(iC %*% param),
se = sqrt(diag(iC %*% vcov.param %*% t(iC))),
df = ci.df,
statistic = NA,
lower = NA,
upper = NA,
null = iNull,
p.value = NA,
stringsAsFactors = FALSE)
CI$statistic <- (CI$estimate-iNull)/CI$se
rownames(CI) <- rownames(iC)
if(!is.null(names(effects)) && !inherits(effects,"mcp")){
indexName <- intersect(which(names(effects)!=""),which(!is.na(names(effects))))
rownames(CI)[indexName] <- names(effects)[indexName]
}
CI.glht <- multcomp::glht(object, linfct = iC, rhs = iNull, df = ceiling(stats::median(ci.df)),
coef. = function(iX){coef.lmm(iX, effects = effects)},
vcov. = function(iX){vcov.lmm(iX, robust = robust, effects = effects)})
CI.glht$model <- NULL
attr(CI.glht$vcov,"robust") <- robust
}else{
CI <- NULL
CI.glht <- NULL
}
## *** test
iRes <- data.frame("null" = iName.hypo,
"statistic" = iStat,
"df.num" = iDf[1],
"df.denom" = iDf[2],
"p.value" = 1 - stats::pf(iStat, df1 = iDf[1], df2 = iDf[2]),
stringsAsFactors = FALSE)
## R2 calculation from
## "An R2 statistic for fixed effects in the linear mixed model" by Lloyd J. Edwards et al. 2008 (Statistic in medicine)
## Equation 19
## DOI: 10.1002/sim.3429
if(iType=="all"){
iType2 <- apply(iC, MARGIN = 1, FUN = function(iRow){
iType <- unique(type.param[names(iRow)[which(abs(iRow)>0)]])
if(length(iType)>1){return("all")}else{return(iType)}
})
}else{
iType2 <- switch(iType,
"mean"="mu",
"variance"="k",
"correlation"="rho")
}
if(length(unique(iType2))==1){
out$multivariate <- rbind(out$multivariate, cbind(type = unname(iType2[1]), test = iNameTerm, iRes))
}else{
out$multivariate <- rbind(out$multivariate, cbind(type = "all", test = iNameTerm, iRes))
}
if(ci){
out$univariate <- rbind(out$univariate, cbind(type = iType2, test = iNameTerm, CI))
}
out$glht[[iType]][[iNameTerm]] <- CI.glht
}
}
## ** save some of the objects for possible use of rbind.Wald_lmm
if(ci > 0.5){
out$object <- list(outcome = object$outcome$var,
method.fit = object$args$method.fit,
type.information = type.information,
cluster.var = object$cluster$var,
structure = c(type = object$design$vcov$type, class = object$design$vcov$class))
if(!is.na(attr(object$cluster$var,"original"))){
out$object$cluster <- object$cluster$levels
}
globalC <- do.call(rbind,lapply(unlist(out$glht, recursive=FALSE),"[[","linfct"))
## default: non-robust vcov and robust for iid
if(attr(robust,"call")){
out$vcov <- globalC %*% vcov.param %*% t(globalC)
robust2 <- robust
}else{
out$vcov <- globalC %*% vcov(object, df = FALSE, effects = "all", robust = FALSE,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE) %*% t(globalC)
robust2 <- TRUE
}
if(object$args$method.fit == "ML"){
out$iid <- iid(object, effects = effects, robust = robust2,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names) %*% t(globalC)
}else if(object$args$method.fit == "REML" && all(type.param[which(colSums(globalC!=0)>0)]=="mu")){
globalC <- globalC[,name.param[type.param=="mu"],drop=FALSE]
out$iid <- iid(object, effects = "mean", robust = robust2,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names) %*% t(globalC)
}
}
## ** prepare for back-transformation
out$args <- data.frame(type = NA, robust = robust, df = df, ci = ci,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho,
transform.names = transform.names)
if(all(name.paramSigma %in% colnames(globalC) == FALSE) || all(globalC[,name.paramSigma,drop=FALSE]==0)){
out$args$transform.sigma <- NA
}
if(all(name.paramK %in% colnames(globalC) == FALSE) || all(globalC[,name.paramK,drop=FALSE]==0)){
out$args$transform.k <- NA
}
if(all(name.paramRho %in% colnames(globalC) == FALSE) || all(globalC[,name.paramRho,drop=FALSE]==0)){
out$args$transform.rho <- NA
}
test.original <- is.null(original.transform.sigma) && is.null(original.transform.k) && is.null(original.transform.rho)
test.backtransform <- stats::na.omit(c(sigma = out$args$transform.sigma, k = out$args$transform.k, rho = out$args$transform.rho))
## should back transformation be performed when latter calling coef/confint?
## - no if no sigma/k/rho parameters in the linear hypotheses
## - no for the p-value, yes for the estimate, and NA for CI when evaluating a contrast on a tranformed scaled
out$args$backtransform <- length(test.backtransform[test.backtransform != "none"])>0
if(test.original && out$args$backtransform){
## find parameters subject to backtransformation
param.with.backtransform <- object$design$param[object$design$param$type %in% names(test.backtransform)[test.backtransform!="none"],"name"]
## restrict contrast matrix to hypotheses with backtransformation
M.allContrast <- do.call(rbind,lapply(out$glht[lengths(out$glht)>0], function(iLs){do.call(rbind,lapply(iLs,"[[","linfct"))}))
col.with.backtransform <- intersect(colnames(M.allContrast),param.with.backtransform)
row.with.backtransform <- rowSums(M.allContrast[,col.with.backtransform,drop=FALSE]!=0)>0
Mback.allContrast <- M.allContrast[row.with.backtransform,,drop=FALSE]
## sanity check
if(any(rowSums(M.allContrast[row.with.backtransform,setdiff(colnames(M.allContrast),param.with.backtransform),drop=FALSE]!=0)>0)){
warning("Hypothesis with mean and variance-covariance parameter(s): possible misleading results due to transformation of the variance-covariance parameter(s). \n")
}
## evaluate un-transformed value
paramMback.allContrast <- rowSums(Mback.allContrast!=0)
if(any(paramMback.allContrast>1)){
param.untransformed <- coef(object, effects = effects, transform.sigma = "none", transform.k = "none", transform.rho = "none")
attr(out$univariate,"backtransform") <- (Mback.allContrast[paramMback.allContrast>1,,drop=FALSE] %*% param.untransformed)[,1]
}
}
## nice naming
if(transform.names && is.null(names(effects))){
backtransform.names <- names(coef(object,
effects = effects,
transform.sigma = gsub("log","",transform.sigma),
transform.k = gsub("log","",transform.k),
transform.rho = gsub("atanh","",transform.rho),
transform.names = transform.names))
out$args$backtransform.names <- list(stats::setNames(rownames(out$univariate),rownames(out$univariate)))
index <- match(newname,out$args$backtransform.names[[1]])
out$args$backtransform.names[[1]][stats::na.omit(index)] <- backtransform.names[which(!is.na(index))]
}
if(identical(type,"all")){
out$args$type <- list("all")
}else{
out$args$type <- list(c("mu"="mean", "k"="variance", "rho"="correlation")[out$multivariate$type])
}
## ** export
return(out)
}
## * .anova_LRT
.anova_LRT <- function(object1,object2){
tol <- 1e-10
## ** normalize user input
## *** re-order models
logLik1 <- logLik(object1)
logLik2 <- logLik(object2)
if(is.na(logLik1) || is.na(logLik2)){
stop("Cannot perform a likelihood ratio test when the log-likelihood is NA for one of the models.\n")
}else if(logLik2>=logLik1){
type <- "2-1"
objectH0 <- object1
objectH1 <- object2
}else if(logLik1>=logLik2){
type <- "1-2"
objectH0 <- object2
objectH1 <- object1
}
## *** extract coefficients
name.paramH0 <- names(coef(objectH0, effects = "all"))
name.paramH1 <- names(coef(objectH1, effects = "all"))
table.paramH0 <- objectH0$design$param
table.paramH1 <- objectH1$design$param
type.paramH0 <- stats::setNames(table.paramH0[match(name.paramH0, table.paramH0$name),"type"], name.paramH0)
type.paramH1 <- stats::setNames(table.paramH1[match(name.paramH1, table.paramH1$name),"type"], name.paramH1)
mismatchH0 <- stats::setNames(name.paramH0 %in% name.paramH1 == FALSE, name.paramH0)
mismatchH1 <- stats::setNames(name.paramH1 %in% name.paramH0 == FALSE, name.paramH1)
rhs <- stats::setNames(rep("0", sum(mismatchH1)), names(which(mismatchH1)))
current.mismatchH0 <- mismatchH0
current.mismatchH1 <- mismatchH1
## *** objective function
if(object1$args$method.fit!=object2$args$method.fit){
stop("The two models should use the same type of objective function for the likelihood ratio test to be valid. \n")
}
if(objectH1$args$method.fit=="REML" && (any(type.paramH1[mismatchH1]=="mu") || any(type.paramH0[mismatchH0]=="mu"))){
objectH0$call$method.fit <- "ML"
objectH1$call$method.fit <- "ML"
message("Cannot use a likelihood ratio test to compare mean parameters when the objective function is REML. \n",
"Will re-estimate the model via ML and re-run the likelihood ratio test. \n")
out <- anova(eval(objectH0$call),eval(objectH1$call))
attr(out,"type") <- type
return(out)
}
## *** number of observations
nobsH0 <- nobs(objectH0)
nobsH1 <- nobs(objectH1)
if(any(nobsH0 != nobsH1)){
if(nobsH0["missing"]!=nobsH1["missing"]){
stop("Mismatch between the number of observations between the two models - could be due to missing data. \n",
"H0: ",paste(paste(names(nobsH0),"=",nobsH0), collapse = ", "),".\n",
"H1: ",paste(paste(names(nobsH1),"=",nobsH0), collapse = ", "),".\n")
}else{
stop("Mismatch between the number of observations between the two models. \n",
"H0: ",paste(paste(names(nobsH0),"=",nobsH0), collapse = ", "),".\n",
"H1: ",paste(paste(names(nobsH1),"=",nobsH0), collapse = ", "),".\n")
}
}
## *** check nesting
if(any(abs(objectH1$design$Y-objectH0$design$Y)>tol)){
stop("Mismatch in outcome between the two models. \n")
}
diff.strata <- setdiff(objectH0$design$vcov$name$strata,objectH0$design$vcov$name$strata)
if(length(diff.strata)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"Variance-covariance parameters are stratified with respect to \"",diff.strata,"\" in the model with the smallest likelihood but not in the model with the largest likelihood. \n",sep="")
}
if(any(mismatchH0)){
## name may not match even with nested model when using different covariance patterns
## e.g. CS: rho while UN gives rho(1,2), rho(1,3), ...
response <- objectH1$design$Y
if(length(name.paramH0)==length(name.paramH1)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The two models have the same number of parameters. \n")
}
if(length(name.paramH0)>length(name.paramH1)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood had less parameters than the model with the smallest likelihood. \n")
}
## mean
if("mu" %in% type.paramH0[mismatchH0]){
mismatchH0.mu <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("mu")]
mismatchH1.mu <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("mu")]
if(length(mismatchH0.mu)>length(mismatchH1.mu)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less mean parameters than the model with the smallest likelihood. \n")
}
XH0mu <- model.matrix(objectH0, effects = "mean")
XH1mu <- model.matrix(objectH1, effects = "mean")
H.X0mu <- XH0mu %*% solve(crossprod(XH0mu)) %*% t(XH0mu)
H.X1mu <- XH1mu %*% solve(crossprod(XH1mu)) %*% t(XH1mu)
## attempt to check whether the design matrices are equivalent
if(any(abs(H.X0mu-H.X1mu)>tol)){
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had mean parameters \"",paste(mismatchH0.mu, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}else{
current.mismatchH0[mismatchH0.mu] <- FALSE
current.mismatchH1[mismatchH1.mu] <- FALSE
rhs <- rhs[names(current.mismatchH1)[current.mismatchH1]]
}
}
## variance
if("sigma" %in% type.paramH0[mismatchH0] || "k" %in% type.paramH0[mismatchH0]){
mismatchH0.ksigma <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("sigma","k")]
mismatchH1.ksigma <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("sigma","k")]
if(length(mismatchH0.ksigma)>length(mismatchH1.ksigma)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less variance parameters than the model with the smallest likelihood. \n")
}
if(identical(objectH1$design$vcov$name$strata,objectH0$design$vcov$name$strata)){
XH0ksigma <- model.matrix(objectH0, effects = "variance")$var$X
XH1ksigma <- model.matrix(objectH1, effects = "variance")$var$X
H.X0ksigma <- XH0ksigma %*% solve(crossprod(XH0ksigma)) %*% t(XH0ksigma)
H.X1ksigma <- XH1ksigma %*% solve(crossprod(XH1ksigma)) %*% t(XH1ksigma)
## attempt to check whether the design matrices are equivalent
if(any(abs(H.X0ksigma-H.X1ksigma)>tol)){
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}else{
current.mismatchH0[mismatchH0.ksigma] <- FALSE
current.mismatchH1[mismatchH1.ksigma] <- FALSE
rhs <- rhs[names(current.mismatchH1)[current.mismatchH1]]
}
}else if(all(is.na(objectH0$design$vcov$name$strata))){
mismatchH1.ksigma2 <- sapply(strsplit(mismatchH1.ksigma,":",fixed = TRUE),"[[",1)
if(all(mismatchH1.ksigma2 %in% mismatchH0.ksigma)){
current.mismatchH0[mismatchH0.ksigma] <- FALSE
rhs[mismatchH1.ksigma[mismatchH1.ksigma2 %in% mismatchH0.ksigma]] <- mismatchH1.ksigma2[mismatchH1.ksigma2 %in% mismatchH0.ksigma]
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}
}
## correlation
if("rho" %in% type.paramH0[mismatchH0]){
mismatchH0.rho <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("rho")]
mismatchH1.rho <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("rho")]
if(length(mismatchH0.rho)>length(mismatchH1.rho)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less correlation parameters than the model with the smallest likelihood. \n")
}
if(identical(objectH1$design$vcov$name$strata,objectH0$design$vcov$name$strata)){
## no strata or same strata: effect of the type of structure
n.strata <- objectH1$strata$n
test.nested1 <- objectH1$design$vcov$class=="UN" && objectH0$design$vcov$class=="CS" && is.na(objectH0$design$vcov$name$var) && is.na(objectH0$design$vcov$name$cor)
if(test.nested1){
mismatchH0 <- setdiff(mismatchH0,mismatchH0.rho)
for(iS in 1:n.strata){
iRhoH0 <- table.paramH0[which((table.paramH0$type=="rho")*(table.paramH0$index.strata==1)==1),"name"]
iRhoH1 <- table.paramH1[which((table.paramH1$type=="rho")*(table.paramH1$index.strata==1)==1),"name"]
rhs[iRhoH1] <- iRhoH0
}
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had correlation parameters \"",paste(mismatchH0.rho, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n")
}
}else{
## effect of strata
mismatchH1.rho2 <- sapply(strsplit(mismatchH1.rho,":",fixed = TRUE),"[[",1)
if(all(mismatchH1.rho2 %in% mismatchH0.rho)){
current.mismatchH0[mismatchH0.rho] <- FALSE
rhs[mismatchH1.rho[mismatchH1.rho2 %in% mismatchH0.rho]] <- mismatchH1.rho2[mismatchH1.rho2 %in% mismatchH0.rho]
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had correlation parameters \"",paste(mismatchH0.rho, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n")
}
}
}
}
n.paramTest <- length(name.paramH1)-length(name.paramH0)
## ** LRT
out <- data.frame(null = paste(paste0(names(rhs),"==",rhs), collapse = ", "),
logLikH1 = stats::logLik(objectH1),
logLikH0 = stats::logLik(objectH0),
statistic = NA,
df = n.paramTest,
p.value = NA,
stringsAsFactors = FALSE)
out$statistic <- 2*(out$logLikH1 - out$logLikH0)
out$p.value <- 1 - stats::pchisq(out$statistic, df = out$df)
## ** export
attr(out,"type") <- type
return(out)
}
## * anova.mlmm
##' @export
anova.mlmm <- function(object, effects = NULL, rhs = NULL, ...){
## ** normalize argument
if(is.null(effects)){
class(object) <- setdiff(class(object), c("mlmm"))
return(object)
}
if(inherits(effects, "mcp")){
if(length(effects)!=1){
stop("Argument \'effects\' must specify a single hypothesis test when being of class \"mcp\". \n",
"Something like mcp(group = \"Dunnett\") or mcp(group = \"Tukey\") \n")
}
effects.save <- effects
constraint <- effects.save[[1]]
effects <- names(effects.save)
if(!grepl("=",effects)){
effects <- paste0(effects,"=0")
}
}else{
constraint <- NULL
}
## ** test linear combinations
robust <- object$args$robust
df <- object$args$df
ci <- object$args$ci
transform.sigma <- if(is.na(object$args$transform.sigma)){NULL}else{object$args$transform.sigma}
transform.k <- if(is.na(object$args$transform.k)){NULL}else{object$args$transform.k}
transform.rho <- if(is.na(object$args$transform.rho)){NULL}else{object$args$transform.rho}
ls.lmm <- object$model
name.lmm <- names(ls.lmm)
ls.anova <- stats::setNames(lapply(name.lmm, function(iName){ ## iName <- name.lmm[1]
anova(ls.lmm[[iName]], effects = effects, rhs = rhs, df = df, ci = ci, robust = robust,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho)
}), name.lmm)
## ** regenerate a new mlmm object
out <- do.call("rbind.Wald_lmm",
args = c(list(model = ls.anova[[1]], effects = constraint, rhs = rhs, name = names(object$model), sep = object$args$sep), unname(ls.anova[-1]))
)
return(out)
}
## * dfSigma
##' @title Degree of Freedom for the Chi-Square Test
##' @description Computation of the degrees of freedom of the chi-squared distribution
##' relative to the model-based variance. Copied of lavaSearch2:::dfSigmaRobust.
##' @noRd
##'
##' @param contrast [numeric vector] the linear combination of parameters to test
##' @param vcov [numeric matrix] the variance-covariance matrix of the parameters.
##' @param dVcov [numeric array] the first derivative of the variance-covariance matrix of the parameters.
##' @param keep.param [character vector] the name of the parameters with non-zero first derivative of their variance parameter.
##'
dfSigma <- function(contrast, vcov, dVcov, keep.param){
## iLink <- "LogCau~eta"
C.vcov.C <- rowSums(contrast %*% vcov * contrast) ## variance matrix of the linear combination
## C.vcov.C - vcov[iLink,iLink]
C.dVcov.C <- sapply(keep.param, function(x){
rowSums(contrast %*% dVcov[,,x] * contrast)
})
## C.dVcov.C - dVcov[iLink,iLink,]
numerator <- 2 *(C.vcov.C)^2
## numerator - 2*vcov[iLink,iLink]^2
denom <- rowSums(C.dVcov.C %*% vcov[keep.param,keep.param,drop=FALSE] * C.dVcov.C)
## denom - t(dVcov[iLink,iLink,]) %*% vcov[keep.param,keep.param,drop=FALSE] %*% dVcov[iLink,iLink,]
df <- numerator/denom
return(df)
}
## * .simplifyContrast
## remove contrasts making the contrast matrix singular
.simplifyContrast <- function(object, rhs, tol = 1e-10, trace = TRUE){
object.eigen <- eigen(tcrossprod(object))
n.zero <- sum(abs(object.eigen$values) < tol)
if(length(rhs)==1 && NROW(object)>1){
rhs <- rep(rhs, NROW(object))
}
if(n.zero==0){return(list(C = object, rhs = rhs, dim = NROW(object), rm = 0))}
keep.lines <- 1:NROW(object)
for(iLine in NROW(object):1){ ## iLine <- 3
iN.zero <- sum(abs(eigen(tcrossprod(object[setdiff(keep.lines,iLine),,drop=FALSE]))$values) < tol)
if(iN.zero<n.zero){
keep.lines <- setdiff(keep.lines,iLine)
n.zero <- iN.zero
}
if(n.zero==0){
if(trace){
name.rm <- rownames(object)[-keep.lines]
if(length(name.rm)==1){
message("Singular contrast matrix: contrast \"",name.rm,"\" has been removed. \n")
}else{
message("Singular contrast matrix: contrasts \"",paste(name.rm,collapse= "\" \""),"\" have been removed. \n")
}
}
return(list(C = object[keep.lines,,drop=FALSE], rhs = rhs[keep.lines], dim = length(keep.lines), rm = NROW(object)-length(keep.lines)))
}
}
## n.zero>0 so failure
return(NULL)
}
##----------------------------------------------------------------------
### anova.R ends here
Try the LMMstar package in your browser
Any scripts or data that you put into this service are public.
LMMstar documentation built on Nov. 9, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .simplifyContrast dfSigma anova.mlmm .anova_LRT .anova_Wald anova.lmm
Documented in anova.lmm
### anova.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: mar 5 2021 (21:38)
## Version:
## Last-Updated: aug 1 2023 (15:45)
## By: Brice Ozenne
## Update #: 1346
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * anova.lmm (documentation)
##' @title Multivariate Tests For Linear Mixed Model
##' @description Simultaneous tests of linear combinations of the model paramaters using Wald tests or Likelihood Ratio Test (LRT).
##'
##' @param object a \code{lmm} object. Only relevant for the anova function.
##' @param effects [character or numeric matrix] Should the Wald test be computed for all variables (\code{"all"}),
##' or only variables relative to the mean (\code{"mean"} or \code{"fixed"}),
##' or only variables relative to the variance structure (\code{"variance"}),
##' or only variables relative to the correlation structure (\code{"correlation"}).
##' Can also be use to specify linear combinations of coefficients or a contrast matrix, similarly to the \code{linfct} argument of the \code{multcomp::glht} function.
##' @param robust [logical] Should robust standard errors (aka sandwich estimator) be output instead of the model-based standard errors.
##' @param rhs [numeric vector] the right hand side of the hypothesis. Only used when the argument effects is a matrix.
##' @param df [logical] Should a F-distribution be used to model the distribution of the Wald statistic. Otherwise a chi-squared distribution is used.
##' @param ci [logical] Should an estimate, standard error, confidence interval, and p-value be output for each hypothesis?
##' @param transform.sigma,transform.k,transform.rho,transform.names are passed to the \code{vcov} method. See details section in \code{\link{coef.lmm}}.
##' @param ... Not used. For compatibility with the generic method.
##'
##' @return A data.frame (LRT) or a list of containing the following elements (Wald):\itemize{
##' \item \code{multivariate}: data.frame containing the multivariate Wald test.
##' The column \code{df.num} refers to the degrees of freedom for the numerator (i.e. number of hypotheses)
##' wherease the column \code{df.denum} refers to the degrees of freedom for the denominator (i.e. Satterthwaite approximation).
##' \item \code{univariate}: data.frame containing each univariate Wald test.
##' \item \code{glht}: used internally to call functions from the multcomp package.
##' \item \code{object}: list containing key information about the linear mixed model.
##' \item \code{vcov}: variance-covariance matrix associated to each parameter of interest (i.e. hypothesis).
##' \item \code{iid}: matrix containing the influence function relative to each parameter of interest (i.e. hypothesis).
##' \item \code{args}: list containing argument values from the function call.
##' }
##'
##' @details By default adjustment of confidence intervals and p-values for multiple comparisons is based on the distribution of the maximum-statistic.
##' This is refered to as a single-step Dunnett multiple testing procedures in table II of Dmitrienko et al. (2013).
##' It is performed using the multcomp package with the option \code{test = adjusted("single-step")} with equal degrees of freedom
##' or by simulation using a Student's t copula with unequal degrees of freedom (more in the note of the details section of \code{\link{confint.Wald_lmm}}).
##'
##' @seealso
##' \code{\link{summary.Wald_lmm}} or \code{\link{confint.Wald_lmm}} for a summary of the results. \cr
##' \code{\link{autoplot.Wald_lmm}} for a graphical display of the results. \cr
##' \code{\link{rbind.Wald_lmm}} for combining result across models and adjust for multiple comparisons. \cr
##'
##' @references Dmitrienko, A. and D'Agostino, R., Sr (2013), Traditional multiplicity adjustment methods in clinical trials. Statist. Med., 32: 5172-5218. https://doi.org/10.1002/sim.5990.
##'
##' @keywords htest
##'
##' @examples
##' #### simulate data in the long format ####
##' set.seed(10)
##' dL <- sampleRem(100, n.times = 3, format = "long")
##'
##' #### fit Linear Mixed Model ####
##' eUN.lmm <- lmm(Y ~ visit + X1 + X2 + X5,
##' repetition = ~visit|id, structure = "UN", data = dL)
##'
##' #### Multivariate Wald test ####
##' ## F-tests
##' anova(eUN.lmm)
##' anova(eUN.lmm, effects = "all")
##' anova(eUN.lmm, robust = TRUE, df = FALSE)
##' summary(anova(eUN.lmm))
##'
##' ## user defined F-test
##' summary(anova(eUN.lmm, effects = c("X1=0","X2+X5=10")))
##' print(anova(eUN.lmm, effects = "mean_visit"), columns = add("null"))
##'
##' ## chi2-tests
##' anova(eUN.lmm, df = FALSE)
##'
##'
##' ## with standard contrast
##' if(require(multcomp)){
##' amod <- lmm(breaks ~ tension, data = warpbreaks)
##' e.amod <- anova(amod, effect = mcp(tension = "Tukey"))
##' summary(e.amod)
##' }
##'
##' #### Likelihood ratio test ####
##' eUN0.lmm <- lmm(Y ~ X1 + X2, repetition = ~visit|id, structure = "UN", data = dL)
##' anova(eUN.lmm, eUN0.lmm)
##'
##' eCS.lmm <- lmm(Y ~ X1 + X2 + X5, repetition = ~visit|id, structure = "CS", data = dL)
##' anova(eUN.lmm, eCS.lmm)
## * anova.lmm (code)
##' @export
anova.lmm <- function(object, effects = NULL, robust = FALSE, rhs = NULL, df = !is.null(object$df), ci = TRUE,
transform.sigma = NULL, transform.k = NULL, transform.rho = NULL, transform.names = TRUE, ...){
call <- match.call()
## ** normalized user input
dots <- list(...)
options <- LMMstar.options()
if(length(dots)>0){
stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n")
}
if(is.null(effects)){
effects <- options$effects
}
if(inherits(effects,"lmm")){ ## likelihood ratio test
out <- .anova_LRT(object1 = object, object2 = effects)
attr(out,"call") <- call
if(!inherits(out,"LRT_lmm")){ ## in the case of re-estimation of the model via ML no need to re-assign the class
class(out) <- append("LRT_lmm",class(out))
}
}else{ ## Wald test
attr(robust, "call") <- "robust" %in% names(call)
out <- .anova_Wald(object, effects = effects, robust = robust, rhs = rhs, df = df, ci = ci,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names)
attr(out,"call") <- call
class(out) <- append("Wald_lmm",class(out))
}
## ** export
return(out)
}
## * .anova_Wald
.anova_Wald <- function(object, effects, robust, rhs, df, ci,
transform.sigma, transform.k, transform.rho, transform.names){
options <- LMMstar.options()
type.information <- object$args$type.information
original.transform.sigma <- transform.sigma
original.transform.k <- transform.k
original.transform.rho <- transform.rho
## ** normalized user input
terms.mean <- attr(stats::terms(object$formula$mean.design),"term.labels")
subeffect <- NULL
if(!inherits(effects,"mcp") && length(effects)==1){
if(identical(effects,"all")){
effects <- c("mean","variance","correlation")
}else if(grepl("^mean_",effects)){
iLabels <- attr(stats::terms(object$formula$mean.design),"term.labels")
if(any(effects == paste0("mean_",iLabels))){
subeffect <- iLabels[effects == paste0("mean_",iLabels)]
effects <- "mean"
}
}else if(grepl("^variance_",effects) && !is.null(object$design$vcov$var$X)){
iLabels <- attr(stats::terms(object$formula$var.design),"term.labels")
if(any(effects == paste0("variance_",iLabels))){
subeffect <- iLabels[effects == paste0("variance_",iLabels)]
effects <- "variance"
}
}else if(grepl("^cor_",effects) && !is.null(object$design$vcov$cor$X)){
iLabels <- attr(stats::terms(object$formula$cor.design),"term.labels")
if(any(effects == paste0("correlation_",iLabels))){
subeffect <- iLabels[effects == paste0("correlation_",iLabels)]
effects <- "correlation"
}
}
}
init <- .init_transform(transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho,
x.transform.sigma = object$reparametrize$transform.sigma, x.transform.k = object$reparametrize$transform.k, x.transform.rho = object$reparametrize$transform.rho)
transform.sigma <- init$transform.sigma
transform.k <- init$transform.k
transform.rho <- init$transform.rho
name.coef <- names(stats::coef(object, effects = "all"))
out <- list(multivariate = NULL, univariate = NULL, glht = NULL)
if(inherits(effects,"mcp")){
out.glht <- try(multcomp::glht(object, linfct = effects), ## only used for generating contrast matrix
silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
type <- "all"
if(length(names(effects))==1){
ls.nameTerms <- list(all = names(effects))
}else{
ls.nameTerms <- list(all = NULL)
}
ls.nameTerms.num <- list(all = 1)
ls.contrast <- list(all = matrix(0, nrow = NROW(out.glht$linfct), ncol = length(name.coef), dimnames = list(rownames(out.glht$linfct),name.coef)))
ls.contrast$all[,colnames(out.glht$linfct)] <- out.glht$linfct
ls.null <- list(all = out.glht$rhs)
name.effects <- NULL
}else if(is.matrix(effects)){
## try to re-size the matrix if necessary
if(NCOL(effects)!=length(name.coef)){
if(is.null(colnames(effects))){
stop("Argument \'effect\' should have column names when a matrix. \n")
}
if(any(duplicated(colnames(effects)))){
stop("Argument \'effect\' should not have duplicated column names when a matrix. \n")
}
if(any(colnames(effects) %in% name.coef == FALSE)){
stop("Argument \'effect\' should not have column names matching the coefficient names when a matrix. \n")
}
effects.save <- effects
effects <- matrix(0, nrow = NROW(effects.save), ncol = length(name.coef), dimnames = list(rownames(effects.save),name.coef))
effects[,colnames(effects.save)] <- effects.save
}
if(is.null(rhs)){
rhs <- rep(0, NROW(effects))
}
## run glht
out.glht <- try(multcomp::glht(object, linfct = effects, rhs = rhs, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = FALSE, effects = "all")}),
silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
type <- "all"
ls.nameTerms <- list(all = NULL)
ls.nameTerms.num <- list(all = 1)
ls.contrast <- list(all = out.glht$linfct)
ls.null <- list(all = out.glht$rhs)
name.effects <- rownames(effects)
}else if(all(tolower(effects) %in% c("mean","fixed","variance","correlation"))){
if(transform.k %in% c("sd","var","logsd","logvar")){
stop("Cannot use \'transform.rho\' equal \"sd\", \"var\", \"logsd\", or \"logvar\". \n",
"anova does not handle tests where the null hypothesis is at a boundary of the support of a random variable. \n")
}
if(transform.rho %in% c("cov")){
stop("Cannot use \'transform.rho\' equal \"cov\". \n",
"anova does not handle tests where the null hypothesis is at a boundary of the support of a random variable. \n")
}
effects <- match.arg(effects, c("mean","fixed","variance","correlation"), several.ok = TRUE)
effects[effects=="fixed"] <- "mean"
name.effects <- NULL
type <- NULL
ls.assign <- list()
ls.nameTerms <- list()
ls.contrast <- list()
ls.null <- list()
if("mean" %in% effects){
type <- c(type, "mean")
ls.assign$mean <- attr(object$design$mean,"assign")
ls.nameTerms$mean <- attr(stats::terms(object$formula$mean.design),"term.labels")
ls.contrast <- c(ls.contrast,list(mean = NULL))
null.mean <- 0
ls.null$mean <- rep(null.mean,length(ls.nameTerms$mean))
}
if("variance" %in% effects){
type <- c(type, "variance")
ls.assign$variance <- attr(object$design$vcov$var$X,"assign")
ls.nameTerms$variance <- attr(stats::terms(object$formula$var),"term.labels")
if(!is.na(object$design$vcov$name$strata)){
ls.assign$variance[ls.nameTerms$variance[ls.assign$variance]==object$design$vcov$name$strata] <- 0
}
ls.contrast <- c(ls.contrast,list(variance = NULL))
null.variance <- switch(transform.k,
"none" = 1,
"square" = 1,
"log" = 0,
"logsquare" = 0)
ls.null$variance <- rep(null.variance,length(ls.nameTerms$variance))
}
if("correlation" %in% effects){
type <- c(type, "correlation")
ls.assign$correlation <- rep(1,sum(object$design$param$type=="rho"))
ls.nameTerms$correlation <- if(!is.null(ls.assign$correlation)){object$time$var}else{NULL}
ls.contrast <- c(ls.contrast,list(correlation = NULL))
null.correlation <- switch(transform.rho,
"none" = 0,
"atanh" = 0,
"cov" = 0)
ls.null$correlation <- rep(null.correlation,length(ls.nameTerms$correlation))
}
ls.nameTerms.num <- lapply(ls.nameTerms, function(iName){as.numeric(factor(iName, levels = iName))})
}else if(all(grepl("=",effects)==FALSE)){
stop("Incorrect argument \'effects\': can be \"mean\", \"variance\", \"correlation\", \"all\", \n",
"or something compatible with the argument \'linfct\' of multcomp::glht. \n ")
}else{ ## symbolic definition of effects using equations (characters)
## run glht
out.glht <- try(multcomp::glht(object, linfct = effects, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = FALSE, effects = "all")}),
silent = TRUE)
newname.coef <- names(stats::coef(object, effects = "all"))
if(inherits(out.glht,"try-error")){
## maybe the error is due to white space in the name of the coefficients
if(any(grepl(" ",effects))){
lhs <- sapply(strsplit(effects, split = "=",fixed = TRUE),"[",1)
lhs.term <- unlist(strsplit(unlist(strsplit(lhs, split = "+", fixed = TRUE)), split = "-", fixed = TRUE))
lhs.variable <- trimws(lapply(strsplit(lhs.term, split = "*", fixed = TRUE), utils::tail, 1), which = "both")
if(any(grepl(" ",lhs.variable))){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht,
"There seems to be whitespace(s) in the name of certain coefficients which is likely to be the cause of the error. \n",
"Consider using a contrast matrix to specific the argument \'effects\' \n or renaming the values of categorical variables without whitespace and refiting the lmm object. \n"
)
}
}
test.reparametrize <- any(grepl("log", c(object$reparametrize$transform.sigma,object$reparametrize$transform.k))) || grep("atanh", object$reparametrize$transform.rho)
## restaure untransformed parametrization (glht does not handle log(k). or atanh(cor))
if(test.reparametrize){
object2 <- object
index.var <- which(object$design$param$type %in% c("sigma","k","rho"))
object2$reparametrize <- .reparametrize(p = object$param[index.var],
type = object$design$param$type[index.var],
level = object$design$param$level[index.var],
sigma = object$design$param$sigma[index.var],
k.x = object$design$param$k.x[index.var],
k.y = object$design$param$k.y[index.var],
Jacobian = FALSE, dJacobian = FALSE, inverse = FALSE,
transform.sigma = "none",
transform.k = "none",
transform.rho = "none",
transform.names = TRUE)
out.glht <- try(multcomp::glht(object2, linfct = effects, ## only used for generating contrast matrix
coef. = function(iX){coef.lmm(iX, effects = "all")},
vcov. = function(iX){vcov.lmm(iX, robust = robust, effects = "all")}), silent = TRUE)
if(inherits(out.glht,"try-error")){
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
oldname.coef <- colnames(out.glht$linfct)
newname.hypo <- rownames(out.glht$linfct)
for(iSub in 1:length(oldname.coef)){
newname.hypo <- gsub(pattern = oldname.coef[iSub], replacement = newname.coef[iSub], x = newname.hypo, fixed = TRUE)
}
dimnames(out.glht$linfct) <- list(newname.hypo,newname.coef)
message("The coefficient names have been transformed but not the null hypotheses. \n")
}else{
stop("Possible mispecification of the argument \'effects\' as running mulcomp::glht lead to the following error: \n",
out.glht)
}
}
type <- "all"
ls.contrast <- list(all = out.glht$linfct)
ls.null <- list(all = out.glht$rhs)
name.effects <- names(effects)
if(!is.null(name.effects)){
rownames(ls.contrast$all) <- name.effects
}
if(length(effects)==1){
ls.nameTerms <- list(all = name.effects)
}else{
ls.nameTerms <- list(all = NULL)
}
ls.nameTerms.num <- list(all = 1)
}
## ** prepare
if(robust && object$args$method.fit=="REML"){
name.mean <- names(coef(object, effects = "mean"))
if(all(effects %in% c("mean","variance","correlation"))){
if(all(effects %in% c("variance","correlation"))){
stop("Cannot test variance, covariance, or correlation parameters using robust standard errors with REML. \n")
}
}else if(any(names(which(colSums(abs(do.call(rbind,ls.contrast)))>0)) %in% name.mean == FALSE)){
stop("Cannot test variance, covariance, or correlation parameters using robust standard errors with REML. \n")
}
ls.contrast <- lapply(ls.contrast, function(iC){
if(!is.null(iC)){iC[,name.mean,drop=FALSE]}else{iC}
})
effects <- "mean"
}else{
effects <- "all"
}
type.param <- stats::setNames(object$design$param$type,object$design$param$name)
param <- coef(object, effects = effects,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
newparam <- coef(object, effects = effects,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names)
newname <- names(newparam)
name.param <- names(param)
name.paramSigma <- names(type.param)[type.param=="sigma"]
name.paramK <- names(type.param)[type.param=="k"]
name.paramRho <- names(type.param)[type.param=="rho"]
n.param <- length(param)
vcov.param <- vcov(object, df = df*2, effects = effects, robust = robust,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
dVcov.param <- attr(vcov.param,"dVcov")
if(df>0 && object$args$method.fit=="REML" && type.information == "expected"){
warning("when using REML with expected information, the degree of freedom of the F-statistic may depend on the parametrisation of the variance parameters. \n")
}
if(any(sapply(ls.contrast, function(iC){is.null(iC) || identical(colnames(iC), names(param))}) == FALSE)){
warning("Names of the columns of the contrast matrix do not match the names of the model coefficients. \n")
}
## ** F-tests
out$glht <- stats::setNames(vector(mode = "list", length = length(type)), type)
for(iType in type){ ## iType <- "variance"
## skip empty type
if(length(ls.nameTerms.num[[iType]])==0 || (is.null(ls.contrast[[iType]]) && (all(ls.assign[[iType]]==0)))){ next }
iParam <- coef(object, effects = iType,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE)
name.iParam <- names(iParam)
if(is.null(ls.nameTerms[[iType]])){
ls.nameTerms[[iType]] <- ls.nameTerms.num[[iType]]
}
out$glht[[iType]] <- stats::setNames(vector(mode = "list", length = length(ls.nameTerms.num[[iType]])), ls.nameTerms[[iType]])
for(iTerm in ls.nameTerms.num[[iType]]){ ## iTerm <- 1
iNameTerm <- ls.nameTerms[[iType]][[which(iTerm == ls.nameTerms.num[[iType]])]]
## *** contrast matrix
if(is.null(ls.contrast[[iType]])){
if(all(ls.assign[[iType]]!=iTerm)){next}
if(!is.null(subeffect) && ls.nameTerms$mean[iTerm]!=subeffect){next}
iIndex.param <- which(ls.assign[[iType]]==iTerm)
iN.hypo <- length(iIndex.param)
iNull <- rep(ls.null[[iType]][iTerm],iN.hypo)
iName.hypo <- paste(paste0(name.iParam[iIndex.param],"=",iNull), collapse = ", ")
iC <- matrix(0, nrow = iN.hypo, ncol = n.param, dimnames = list(name.iParam[iIndex.param], name.param))
if(length(iIndex.param)==1){
iC[name.iParam[iIndex.param],name.iParam[iIndex.param]] <- 1
}else{
diag(iC[name.iParam[iIndex.param],name.iParam[iIndex.param]]) <- 1
}
iC.uni <- iC
}else{
iC <- ls.contrast[[iType]]
iC.uni <- iC
iN.hypo <- NROW(iC)
iNull <- ls.null[[iType]]
iName.hypo <- paste(paste0(rownames(iC),"=",iNull),collapse=", ")
}
## *** statistic
outSimp <- .simplifyContrast(iC, iNull) ## remove extra lines
iC.vcov.C_M1 <- try(solve(outSimp$C %*% vcov.param %*% t(outSimp$C)), silent = TRUE)
if(inherits(iC.vcov.C_M1,"try-error")){
iStat <- NA
iDf <- c(outSimp$dim,Inf)
attr(iStat,"error") <- "\n Could not invert the covariance matrix for the proposed contrast."
}else{
iStat <- as.double(t(outSimp$C %*% param - outSimp$rhs) %*% iC.vcov.C_M1 %*% (outSimp$C %*% param - outSimp$rhs))/outSimp$dim
iDf <- c(outSimp$dim,Inf)
if(outSimp$rm>0){
iName.hypo <- paste(paste0(rownames(outSimp$C),"=",outSimp$rhs),collapse=", ")
}
}
## *** degree of freedom
if(df>0 && !inherits(iC.vcov.C_M1,"try-error")){
svd.tempo <- eigen(iC.vcov.C_M1)
D.svd <- diag(svd.tempo$values, nrow = outSimp$dim, ncol = outSimp$dim)
P.svd <- svd.tempo$vectors
contrast.svd <- sqrt(D.svd) %*% t(P.svd) %*% outSimp$C
colnames(contrast.svd) <- name.param
iNu_m <- dfSigma(contrast = contrast.svd,
vcov = vcov.param,
dVcov = dVcov.param,
keep.param = name.param)
iEQ <- sum(iNu_m/(iNu_m - 2))
iDf[2] <- 2 * iEQ/(iEQ - outSimp$dim)
}
## *** confidence interval
if(ci>0){
if(df>0){
ci.df <- .dfX(X.beta = iC.uni, vcov.param = vcov.param, dVcov.param = dVcov.param, return.vcov = ci>0.5)
}else{
ci.df <- Inf
}
CI <- data.frame(estimate = as.double(iC %*% param),
se = sqrt(diag(iC %*% vcov.param %*% t(iC))),
df = ci.df,
statistic = NA,
lower = NA,
upper = NA,
null = iNull,
p.value = NA,
stringsAsFactors = FALSE)
CI$statistic <- (CI$estimate-iNull)/CI$se
rownames(CI) <- rownames(iC)
if(!is.null(names(effects)) && !inherits(effects,"mcp")){
indexName <- intersect(which(names(effects)!=""),which(!is.na(names(effects))))
rownames(CI)[indexName] <- names(effects)[indexName]
}
CI.glht <- multcomp::glht(object, linfct = iC, rhs = iNull, df = ceiling(stats::median(ci.df)),
coef. = function(iX){coef.lmm(iX, effects = effects)},
vcov. = function(iX){vcov.lmm(iX, robust = robust, effects = effects)})
CI.glht$model <- NULL
attr(CI.glht$vcov,"robust") <- robust
}else{
CI <- NULL
CI.glht <- NULL
}
## *** test
iRes <- data.frame("null" = iName.hypo,
"statistic" = iStat,
"df.num" = iDf[1],
"df.denom" = iDf[2],
"p.value" = 1 - stats::pf(iStat, df1 = iDf[1], df2 = iDf[2]),
stringsAsFactors = FALSE)
## R2 calculation from
## "An R2 statistic for fixed effects in the linear mixed model" by Lloyd J. Edwards et al. 2008 (Statistic in medicine)
## Equation 19
## DOI: 10.1002/sim.3429
if(iType=="all"){
iType2 <- apply(iC, MARGIN = 1, FUN = function(iRow){
iType <- unique(type.param[names(iRow)[which(abs(iRow)>0)]])
if(length(iType)>1){return("all")}else{return(iType)}
})
}else{
iType2 <- switch(iType,
"mean"="mu",
"variance"="k",
"correlation"="rho")
}
if(length(unique(iType2))==1){
out$multivariate <- rbind(out$multivariate, cbind(type = unname(iType2[1]), test = iNameTerm, iRes))
}else{
out$multivariate <- rbind(out$multivariate, cbind(type = "all", test = iNameTerm, iRes))
}
if(ci){
out$univariate <- rbind(out$univariate, cbind(type = iType2, test = iNameTerm, CI))
}
out$glht[[iType]][[iNameTerm]] <- CI.glht
}
}
## ** save some of the objects for possible use of rbind.Wald_lmm
if(ci > 0.5){
out$object <- list(outcome = object$outcome$var,
method.fit = object$args$method.fit,
type.information = type.information,
cluster.var = object$cluster$var,
structure = c(type = object$design$vcov$type, class = object$design$vcov$class))
if(!is.na(attr(object$cluster$var,"original"))){
out$object$cluster <- object$cluster$levels
}
globalC <- do.call(rbind,lapply(unlist(out$glht, recursive=FALSE),"[[","linfct"))
## default: non-robust vcov and robust for iid
if(attr(robust,"call")){
out$vcov <- globalC %*% vcov.param %*% t(globalC)
robust2 <- robust
}else{
out$vcov <- globalC %*% vcov(object, df = FALSE, effects = "all", robust = FALSE,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = FALSE) %*% t(globalC)
robust2 <- TRUE
}
if(object$args$method.fit == "ML"){
out$iid <- iid(object, effects = effects, robust = robust2,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names) %*% t(globalC)
}else if(object$args$method.fit == "REML" && all(type.param[which(colSums(globalC!=0)>0)]=="mu")){
globalC <- globalC[,name.param[type.param=="mu"],drop=FALSE]
out$iid <- iid(object, effects = "mean", robust = robust2,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, transform.names = transform.names) %*% t(globalC)
}
}
## ** prepare for back-transformation
out$args <- data.frame(type = NA, robust = robust, df = df, ci = ci,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho,
transform.names = transform.names)
if(all(name.paramSigma %in% colnames(globalC) == FALSE) || all(globalC[,name.paramSigma,drop=FALSE]==0)){
out$args$transform.sigma <- NA
}
if(all(name.paramK %in% colnames(globalC) == FALSE) || all(globalC[,name.paramK,drop=FALSE]==0)){
out$args$transform.k <- NA
}
if(all(name.paramRho %in% colnames(globalC) == FALSE) || all(globalC[,name.paramRho,drop=FALSE]==0)){
out$args$transform.rho <- NA
}
test.original <- is.null(original.transform.sigma) && is.null(original.transform.k) && is.null(original.transform.rho)
test.backtransform <- stats::na.omit(c(sigma = out$args$transform.sigma, k = out$args$transform.k, rho = out$args$transform.rho))
## should back transformation be performed when latter calling coef/confint?
## - no if no sigma/k/rho parameters in the linear hypotheses
## - no for the p-value, yes for the estimate, and NA for CI when evaluating a contrast on a tranformed scaled
out$args$backtransform <- length(test.backtransform[test.backtransform != "none"])>0
if(test.original && out$args$backtransform){
## find parameters subject to backtransformation
param.with.backtransform <- object$design$param[object$design$param$type %in% names(test.backtransform)[test.backtransform!="none"],"name"]
## restrict contrast matrix to hypotheses with backtransformation
M.allContrast <- do.call(rbind,lapply(out$glht[lengths(out$glht)>0], function(iLs){do.call(rbind,lapply(iLs,"[[","linfct"))}))
col.with.backtransform <- intersect(colnames(M.allContrast),param.with.backtransform)
row.with.backtransform <- rowSums(M.allContrast[,col.with.backtransform,drop=FALSE]!=0)>0
Mback.allContrast <- M.allContrast[row.with.backtransform,,drop=FALSE]
## sanity check
if(any(rowSums(M.allContrast[row.with.backtransform,setdiff(colnames(M.allContrast),param.with.backtransform),drop=FALSE]!=0)>0)){
warning("Hypothesis with mean and variance-covariance parameter(s): possible misleading results due to transformation of the variance-covariance parameter(s). \n")
}
## evaluate un-transformed value
paramMback.allContrast <- rowSums(Mback.allContrast!=0)
if(any(paramMback.allContrast>1)){
param.untransformed <- coef(object, effects = effects, transform.sigma = "none", transform.k = "none", transform.rho = "none")
attr(out$univariate,"backtransform") <- (Mback.allContrast[paramMback.allContrast>1,,drop=FALSE] %*% param.untransformed)[,1]
}
}
## nice naming
if(transform.names && is.null(names(effects))){
backtransform.names <- names(coef(object,
effects = effects,
transform.sigma = gsub("log","",transform.sigma),
transform.k = gsub("log","",transform.k),
transform.rho = gsub("atanh","",transform.rho),
transform.names = transform.names))
out$args$backtransform.names <- list(stats::setNames(rownames(out$univariate),rownames(out$univariate)))
index <- match(newname,out$args$backtransform.names[[1]])
out$args$backtransform.names[[1]][stats::na.omit(index)] <- backtransform.names[which(!is.na(index))]
}
if(identical(type,"all")){
out$args$type <- list("all")
}else{
out$args$type <- list(c("mu"="mean", "k"="variance", "rho"="correlation")[out$multivariate$type])
}
## ** export
return(out)
}
## * .anova_LRT
.anova_LRT <- function(object1,object2){
tol <- 1e-10
## ** normalize user input
## *** re-order models
logLik1 <- logLik(object1)
logLik2 <- logLik(object2)
if(is.na(logLik1) || is.na(logLik2)){
stop("Cannot perform a likelihood ratio test when the log-likelihood is NA for one of the models.\n")
}else if(logLik2>=logLik1){
type <- "2-1"
objectH0 <- object1
objectH1 <- object2
}else if(logLik1>=logLik2){
type <- "1-2"
objectH0 <- object2
objectH1 <- object1
}
## *** extract coefficients
name.paramH0 <- names(coef(objectH0, effects = "all"))
name.paramH1 <- names(coef(objectH1, effects = "all"))
table.paramH0 <- objectH0$design$param
table.paramH1 <- objectH1$design$param
type.paramH0 <- stats::setNames(table.paramH0[match(name.paramH0, table.paramH0$name),"type"], name.paramH0)
type.paramH1 <- stats::setNames(table.paramH1[match(name.paramH1, table.paramH1$name),"type"], name.paramH1)
mismatchH0 <- stats::setNames(name.paramH0 %in% name.paramH1 == FALSE, name.paramH0)
mismatchH1 <- stats::setNames(name.paramH1 %in% name.paramH0 == FALSE, name.paramH1)
rhs <- stats::setNames(rep("0", sum(mismatchH1)), names(which(mismatchH1)))
current.mismatchH0 <- mismatchH0
current.mismatchH1 <- mismatchH1
## *** objective function
if(object1$args$method.fit!=object2$args$method.fit){
stop("The two models should use the same type of objective function for the likelihood ratio test to be valid. \n")
}
if(objectH1$args$method.fit=="REML" && (any(type.paramH1[mismatchH1]=="mu") || any(type.paramH0[mismatchH0]=="mu"))){
objectH0$call$method.fit <- "ML"
objectH1$call$method.fit <- "ML"
message("Cannot use a likelihood ratio test to compare mean parameters when the objective function is REML. \n",
"Will re-estimate the model via ML and re-run the likelihood ratio test. \n")
out <- anova(eval(objectH0$call),eval(objectH1$call))
attr(out,"type") <- type
return(out)
}
## *** number of observations
nobsH0 <- nobs(objectH0)
nobsH1 <- nobs(objectH1)
if(any(nobsH0 != nobsH1)){
if(nobsH0["missing"]!=nobsH1["missing"]){
stop("Mismatch between the number of observations between the two models - could be due to missing data. \n",
"H0: ",paste(paste(names(nobsH0),"=",nobsH0), collapse = ", "),".\n",
"H1: ",paste(paste(names(nobsH1),"=",nobsH0), collapse = ", "),".\n")
}else{
stop("Mismatch between the number of observations between the two models. \n",
"H0: ",paste(paste(names(nobsH0),"=",nobsH0), collapse = ", "),".\n",
"H1: ",paste(paste(names(nobsH1),"=",nobsH0), collapse = ", "),".\n")
}
}
## *** check nesting
if(any(abs(objectH1$design$Y-objectH0$design$Y)>tol)){
stop("Mismatch in outcome between the two models. \n")
}
diff.strata <- setdiff(objectH0$design$vcov$name$strata,objectH0$design$vcov$name$strata)
if(length(diff.strata)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"Variance-covariance parameters are stratified with respect to \"",diff.strata,"\" in the model with the smallest likelihood but not in the model with the largest likelihood. \n",sep="")
}
if(any(mismatchH0)){
## name may not match even with nested model when using different covariance patterns
## e.g. CS: rho while UN gives rho(1,2), rho(1,3), ...
response <- objectH1$design$Y
if(length(name.paramH0)==length(name.paramH1)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The two models have the same number of parameters. \n")
}
if(length(name.paramH0)>length(name.paramH1)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood had less parameters than the model with the smallest likelihood. \n")
}
## mean
if("mu" %in% type.paramH0[mismatchH0]){
mismatchH0.mu <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("mu")]
mismatchH1.mu <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("mu")]
if(length(mismatchH0.mu)>length(mismatchH1.mu)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less mean parameters than the model with the smallest likelihood. \n")
}
XH0mu <- model.matrix(objectH0, effects = "mean")
XH1mu <- model.matrix(objectH1, effects = "mean")
H.X0mu <- XH0mu %*% solve(crossprod(XH0mu)) %*% t(XH0mu)
H.X1mu <- XH1mu %*% solve(crossprod(XH1mu)) %*% t(XH1mu)
## attempt to check whether the design matrices are equivalent
if(any(abs(H.X0mu-H.X1mu)>tol)){
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had mean parameters \"",paste(mismatchH0.mu, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}else{
current.mismatchH0[mismatchH0.mu] <- FALSE
current.mismatchH1[mismatchH1.mu] <- FALSE
rhs <- rhs[names(current.mismatchH1)[current.mismatchH1]]
}
}
## variance
if("sigma" %in% type.paramH0[mismatchH0] || "k" %in% type.paramH0[mismatchH0]){
mismatchH0.ksigma <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("sigma","k")]
mismatchH1.ksigma <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("sigma","k")]
if(length(mismatchH0.ksigma)>length(mismatchH1.ksigma)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less variance parameters than the model with the smallest likelihood. \n")
}
if(identical(objectH1$design$vcov$name$strata,objectH0$design$vcov$name$strata)){
XH0ksigma <- model.matrix(objectH0, effects = "variance")$var$X
XH1ksigma <- model.matrix(objectH1, effects = "variance")$var$X
H.X0ksigma <- XH0ksigma %*% solve(crossprod(XH0ksigma)) %*% t(XH0ksigma)
H.X1ksigma <- XH1ksigma %*% solve(crossprod(XH1ksigma)) %*% t(XH1ksigma)
## attempt to check whether the design matrices are equivalent
if(any(abs(H.X0ksigma-H.X1ksigma)>tol)){
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}else{
current.mismatchH0[mismatchH0.ksigma] <- FALSE
current.mismatchH1[mismatchH1.ksigma] <- FALSE
rhs <- rhs[names(current.mismatchH1)[current.mismatchH1]]
}
}else if(all(is.na(objectH0$design$vcov$name$strata))){
mismatchH1.ksigma2 <- sapply(strsplit(mismatchH1.ksigma,":",fixed = TRUE),"[[",1)
if(all(mismatchH1.ksigma2 %in% mismatchH0.ksigma)){
current.mismatchH0[mismatchH0.ksigma] <- FALSE
rhs[mismatchH1.ksigma[mismatchH1.ksigma2 %in% mismatchH0.ksigma]] <- mismatchH1.ksigma2[mismatchH1.ksigma2 %in% mismatchH0.ksigma]
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had variance parameters \"",paste(mismatchH0.ksigma, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n",
"Consider reparametrizing the mean structure so the nesting is more obvious.\n")
}
}
## correlation
if("rho" %in% type.paramH0[mismatchH0]){
mismatchH0.rho <- names(type.paramH0[mismatchH0])[type.paramH0[mismatchH0] %in% c("rho")]
mismatchH1.rho <- names(type.paramH1[mismatchH1])[type.paramH1[mismatchH1] %in% c("rho")]
if(length(mismatchH0.rho)>length(mismatchH1.rho)){
stop("Cannot perform a likelihood ratio test when the model are not nested. \n",
"The model with the largest likelihood has less correlation parameters than the model with the smallest likelihood. \n")
}
if(identical(objectH1$design$vcov$name$strata,objectH0$design$vcov$name$strata)){
## no strata or same strata: effect of the type of structure
n.strata <- objectH1$strata$n
test.nested1 <- objectH1$design$vcov$class=="UN" && objectH0$design$vcov$class=="CS" && is.na(objectH0$design$vcov$name$var) && is.na(objectH0$design$vcov$name$cor)
if(test.nested1){
mismatchH0 <- setdiff(mismatchH0,mismatchH0.rho)
for(iS in 1:n.strata){
iRhoH0 <- table.paramH0[which((table.paramH0$type=="rho")*(table.paramH0$index.strata==1)==1),"name"]
iRhoH1 <- table.paramH1[which((table.paramH1$type=="rho")*(table.paramH1$index.strata==1)==1),"name"]
rhs[iRhoH1] <- iRhoH0
}
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had correlation parameters \"",paste(mismatchH0.rho, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n")
}
}else{
## effect of strata
mismatchH1.rho2 <- sapply(strsplit(mismatchH1.rho,":",fixed = TRUE),"[[",1)
if(all(mismatchH1.rho2 %in% mismatchH0.rho)){
current.mismatchH0[mismatchH0.rho] <- FALSE
rhs[mismatchH1.rho[mismatchH1.rho2 %in% mismatchH0.rho]] <- mismatchH1.rho2[mismatchH1.rho2 %in% mismatchH0.rho]
}else{
stop("Do not understand how the two models are nested. \n",
"The model with the lowest likelihood had correlation parameters \"",paste(mismatchH0.rho, collapse="\" \""),"\" that were not present in the model with the largest likelihood. \n")
}
}
}
}
n.paramTest <- length(name.paramH1)-length(name.paramH0)
## ** LRT
out <- data.frame(null = paste(paste0(names(rhs),"==",rhs), collapse = ", "),
logLikH1 = stats::logLik(objectH1),
logLikH0 = stats::logLik(objectH0),
statistic = NA,
df = n.paramTest,
p.value = NA,
stringsAsFactors = FALSE)
out$statistic <- 2*(out$logLikH1 - out$logLikH0)
out$p.value <- 1 - stats::pchisq(out$statistic, df = out$df)
## ** export
attr(out,"type") <- type
return(out)
}
## * anova.mlmm
##' @export
anova.mlmm <- function(object, effects = NULL, rhs = NULL, ...){
## ** normalize argument
if(is.null(effects)){
class(object) <- setdiff(class(object), c("mlmm"))
return(object)
}
if(inherits(effects, "mcp")){
if(length(effects)!=1){
stop("Argument \'effects\' must specify a single hypothesis test when being of class \"mcp\". \n",
"Something like mcp(group = \"Dunnett\") or mcp(group = \"Tukey\") \n")
}
effects.save <- effects
constraint <- effects.save[[1]]
effects <- names(effects.save)
if(!grepl("=",effects)){
effects <- paste0(effects,"=0")
}
}else{
constraint <- NULL
}
## ** test linear combinations
robust <- object$args$robust
df <- object$args$df
ci <- object$args$ci
transform.sigma <- if(is.na(object$args$transform.sigma)){NULL}else{object$args$transform.sigma}
transform.k <- if(is.na(object$args$transform.k)){NULL}else{object$args$transform.k}
transform.rho <- if(is.na(object$args$transform.rho)){NULL}else{object$args$transform.rho}
ls.lmm <- object$model
name.lmm <- names(ls.lmm)
ls.anova <- stats::setNames(lapply(name.lmm, function(iName){ ## iName <- name.lmm[1]
anova(ls.lmm[[iName]], effects = effects, rhs = rhs, df = df, ci = ci, robust = robust,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho)
}), name.lmm)
## ** regenerate a new mlmm object
out <- do.call("rbind.Wald_lmm",
args = c(list(model = ls.anova[[1]], effects = constraint, rhs = rhs, name = names(object$model), sep = object$args$sep), unname(ls.anova[-1]))
)
return(out)
}
## * dfSigma
##' @title Degree of Freedom for the Chi-Square Test
##' @description Computation of the degrees of freedom of the chi-squared distribution
##' relative to the model-based variance. Copied of lavaSearch2:::dfSigmaRobust.
##' @noRd
##'
##' @param contrast [numeric vector] the linear combination of parameters to test
##' @param vcov [numeric matrix] the variance-covariance matrix of the parameters.
##' @param dVcov [numeric array] the first derivative of the variance-covariance matrix of the parameters.
##' @param keep.param [character vector] the name of the parameters with non-zero first derivative of their variance parameter.
##'
dfSigma <- function(contrast, vcov, dVcov, keep.param){
## iLink <- "LogCau~eta"
C.vcov.C <- rowSums(contrast %*% vcov * contrast) ## variance matrix of the linear combination
## C.vcov.C - vcov[iLink,iLink]
C.dVcov.C <- sapply(keep.param, function(x){
rowSums(contrast %*% dVcov[,,x] * contrast)
})
## C.dVcov.C - dVcov[iLink,iLink,]
numerator <- 2 *(C.vcov.C)^2
## numerator - 2*vcov[iLink,iLink]^2
denom <- rowSums(C.dVcov.C %*% vcov[keep.param,keep.param,drop=FALSE] * C.dVcov.C)
## denom - t(dVcov[iLink,iLink,]) %*% vcov[keep.param,keep.param,drop=FALSE] %*% dVcov[iLink,iLink,]
df <- numerator/denom
return(df)
}
## * .simplifyContrast
## remove contrasts making the contrast matrix singular
.simplifyContrast <- function(object, rhs, tol = 1e-10, trace = TRUE){
object.eigen <- eigen(tcrossprod(object))
n.zero <- sum(abs(object.eigen$values) < tol)
if(length(rhs)==1 && NROW(object)>1){
rhs <- rep(rhs, NROW(object))
}
if(n.zero==0){return(list(C = object, rhs = rhs, dim = NROW(object), rm = 0))}
keep.lines <- 1:NROW(object)
for(iLine in NROW(object):1){ ## iLine <- 3
iN.zero <- sum(abs(eigen(tcrossprod(object[setdiff(keep.lines,iLine),,drop=FALSE]))$values) < tol)
if(iN.zero<n.zero){
keep.lines <- setdiff(keep.lines,iLine)
n.zero <- iN.zero
}
if(n.zero==0){
if(trace){
name.rm <- rownames(object)[-keep.lines]
if(length(name.rm)==1){
message("Singular contrast matrix: contrast \"",name.rm,"\" has been removed. \n")
}else{
message("Singular contrast matrix: contrasts \"",paste(name.rm,collapse= "\" \""),"\" have been removed. \n")
}
}
return(list(C = object[keep.lines,,drop=FALSE], rhs = rhs[keep.lines], dim = length(keep.lines), rm = NROW(object)-length(keep.lines)))
}
}
## n.zero>0 so failure
return(NULL)
}
##----------------------------------------------------------------------
### anova.R ends here
Try the LMMstar package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.