Nothing
### rbind.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: feb 9 2022 (14:51)
## Version:
## Last-Updated: nov 8 2023 (15:09)
## By: Brice Ozenne
## Update #: 489
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * rbind.Wald_lmm (documentation)
##' @title Linear Hypothesis Testing Across Linear Mixed Models
##' @description Linear hypothesis testing accross linear mixed model.
##'
##' @param model a \code{Wald_lmm} object (output of \code{anova} applied to a \code{lmm} object)
##' @param ... possibly other \code{Wald_lmm} objects
##' @param effects [character or numeric matrix] how to combine the left-hand side of the hypotheses.
##' By default identity matrix but can also be \code{"Dunnett"}, \code{"Tukey"}, or \code{"Sequen"} (see function \code{multcomp::contrMat} from the multcomp package).
##' @param rhs [numeric vector] the right hand side of the hypothesis. Should have the same length as the number of row of argument \code{effects}.
##' @param name [character vector or NULL] character used to identify each model in the output.
##' By default, use the name of the outcome of the model.
##' @param sep [character] character used to separate the outcome and the covariate when naming the tests.
##'
##' @details WARNING: in presence of measurements from the same cluster across several models,
##' the influence function is used to estimate the covariance between the model parameters.
##' This is why the (robust) standard errors may not match the (model-based) standard error from the linear mixed
##' Setting the argument \code{robust} to \code{FALSE} when calling \code{anova.lmm} will rescale the (robust) standard errors to mimic the original model-based standard errors.
##' @keywords methods
##'
##' @examples
##' ## simulate data
##' set.seed(10)
##' dL <- sampleRem(1e2, n.times = 3, format = "long")
##'
##' ## estimate mixed models
##' e.lmm1 <- lmm(Y ~ X1+X2+X3, repetition = ~visit|id, data = dL,
##' structure = "CS", df = FALSE)
##' e.lmm2 <- lmm(Y ~ X1+X8+X9, repetition = ~visit|id, data = dL,
##' structure = "CS", df = FALSE)
##'
##' model.tables(e.lmm1) ## model-based standard errors
##' model.tables(e.lmm1, robust = TRUE) ## robust standard errors
##'
##' ## select null hypotheses & combine (robust standard errors)
##' AAA <- anova(e.lmm1, ci = TRUE, effect = c("X1|X2,X3"="X1=0","X2|X1,X3"="X2=0"))
##' BBB <- anova(e.lmm2, ci = TRUE, effect = c("X1|X8,X9"="X1=0"))
##'
##' ZZZ <- rbind(AAA,BBB)
##'
##' ## select null hypotheses & combine (model-based like standard errors)
##' AA <- anova(e.lmm1, ci = TRUE, effect = c("X1|X2,X3"="X1=0","X2|X1,X3"="X2=0"),
##' robust = FALSE)
##' BB <- anova(e.lmm2, ci = TRUE, effect = c("X1|X8,X9"="X1=0"),
##' robust = FALSE)
##' ZZ <- rbind(AA,BB)
## * rbind.Wald_lmm (code)
##' @export
rbind.Wald_lmm <- function(model, ..., effects = NULL, rhs = NULL, name = NULL, sep = ": "){
default <- LMMstar.options()
call <- match.call()
## ** Check user input
dots <- list(...)
if(length(dots)==0){
return(model) ## nothing to combine
}else if(any(sapply(dots,inherits,"Wald_lmm")==FALSE)){
stop("Extra arguments should inherit from Wald_lmm. \n")
}
ls.object <- c(list(model),dots)
n.object <- length(ls.object)
if(!is.null(name) && length(name)!=n.object){
stop("Argument \'name\' should have length ",n.object,", i.e. as many elements as anova objects. \n")
}
if(any(sapply(ls.object, function(iO){NROW(iO$multivariate)})!=1)){
stop("Cannot handle multiple Multivariate Wald test in a single object. \n")
}
table.args <- cbind(do.call(rbind,lapply(ls.object,"[[","args")),
do.call(rbind,lapply(ls.object, function(iO){data.frame(n.test = NROW(iO$univariate),
iO$object[c("outcome","method.fit","type.information")])
})))
Utype <- unique(unlist(table.args$type))
newtable.args <- data.frame(type = ifelse(length(Utype)>1,"all",Utype), sep = sep,
table.args[1,c("robust","df","ci","transform.sigma","transform.k","transform.rho","backtransform")])
newobject <- list()
if(length(unique(table.args$outcome))==1){
newobject$outcome <- table.args$outcome[1]
}
if(length(unique(table.args$method.fit))==1){
newobject$method.fit <- table.args$method.fit[1]
}
if(length(unique(table.args$type.information))==1){
newobject$type.information <- table.args$type.information[1]
}
if(any(table.args$ci==FALSE)){
stop("All argument should contain a \"glht\" object, i.e. call anova with argument ci=TRUE. \n")
}
if(any(newtable.args$robust!=table.args$robust[-1])){
stop("Element \'robust\' should take the same value for all objects. \n")
}
if(any(newtable.args$df!=table.args$df[-1])){
stop("Element \'df\' should take the same value for all objects. \n")
}
if(!all(is.na(table.args$transform.sigma)) && any(newtable.args$transform.sigma!=table.args$transform.sigma[-1])){
stop("Element \'transform.sigma\' should take the same value for all objects. \n")
}
if(!all(is.na(table.args$transform.k)) && any(newtable.args$transform.k!=table.args$transform.k[-1])){
stop("Element \'transform.k\' should take the same value for all objects. \n")
}
if(!all(is.na(table.args$transform.rho)) && any(newtable.args$transform.rho!=table.args$transform.rho[-1])){
stop("Element \'transform.rho\' should take the same value for all objects. \n")
}
if(any(newtable.args$backtransform!=table.args$backtransform[-1])){
stop("Element \'backtransform\' should take the same value for all objects. \n")
}
ls.glht <- unname(unlist(unlist(lapply(ls.object,"[[","glht"),recursive = FALSE),recursive = FALSE))
alternative <- ls.glht[[1]]$alternative
if(any(alternative != sapply(ls.glht,"[[","alternative")[-1])){
stop("Element \'alternative\' should take the same value for all glht objects. \n")
}
if(is.null(name)){
name.modelparam <- unlist(lapply(ls.object, function(iO){rownames(iO$univariate)}))
if(any(duplicated(name.modelparam))){
name.modelparam <- unlist(lapply(ls.object, function(iO){paste0(iO$object$outcome, sep, rownames(iO$univariate))}))
if(is.na(sep) || any(duplicated(name.modelparam))){
stop("Univariate test should have distinct names. \n",
"Consider naming them when calling anova, e.g. anova(object, effect = c(\"myname1\"=\"X1=0\",\"myname2\"=\"X2=0\"))")
}else{
for(iO in 1:n.object){
rownames(ls.object[[iO]]$univariate) <- paste0(ls.object[[iO]]$object$outcome, sep, rownames(ls.object[[iO]]$univariate))
}
}
}
}else{
name.modelparam <- unlist(lapply(1:n.object, function(iO){paste0(name[iO],sep,rownames(ls.object[[iO]]$univariate))}))
for(iO in 1:n.object){ ## iO <- 1
rownames(ls.object[[iO]]$univariate) <- paste0(name[iO],sep,rownames(ls.object[[iO]]$univariate))
if(!is.null(attr(ls.object[[iO]]$univariate, "backtransform"))){
names(attr(ls.object[[iO]]$univariate, "backtransform")) <- paste0(name[iO],sep,attr(ls.object[[iO]]$univariate, "backtransform"))
}
}
if(any(duplicated(name.modelparam))){
stop("Univariate test should have distinct names. \n",
"Consider naming them when calling anova, e.g. anova(object, effect = c(\"myname1\"=\"X1=0\",\"myname2\"=\"X2=0\"))")
}
}
n.modelparam <- length(name.modelparam)
if(!is.null(effects)){
if(is.matrix(effects)){
contrast <- effects
if(NCOL(contrast)!=n.modelparam){
stop("Incorrect contrast matrix: should have ",n.modelparam," columns.\n",
"(one for each univariate test) \n")
}
if(is.null(colnames(contrast))){
colnames(contrast) <- name.modelparam
}else if(!identical(sort(colnames(contrast)), sort(name.modelparam))){
stop("Incorrect column names for argument contrast.\n",
"Should match \"",paste(sort(name.modelparam), collapse="\" \""),"\".\n")
}else{
contrast <- contrast[,name.modelparam,drop=FALSE]
}
}else if(all(is.character(effects))){
if(length(effects)>1){
stop("When a character, argument \'effects\' should have length 1. \n")
}
valid.contrast <- c("Dunnett","Tukey","Sequen")
if(effects %in% valid.contrast == FALSE){
stop("When a character, argument \'effects\' should be one of \"",paste(valid.contrast, collapse = "\" \""),"\". \n")
}
contrast <- multcomp::contrMat(rep(1,length(name.modelparam)), type = effects)
colnames(contrast) <- name.modelparam
try(rownames(contrast) <- unlist(lapply(strsplit(split = "-",rownames(contrast),fixed=TRUE), function(iVec){
paste(name.modelparam[as.numeric(trimws(iVec))], collapse = " - ")
})), silent = TRUE)
}
}else{
contrast <- diag(1, nrow = n.modelparam, ncol = n.modelparam)
dimnames(contrast) <- list(name.modelparam, name.modelparam)
}
n.test <- NROW(contrast)
name.test <- rownames(contrast)
if(!is.null(rhs)){
if(length(rhs)!=n.test){
stop("Incorrect rhs: should have ",n.test," values.\n",
"(one for each univariate test) \n")
}
}else if(is.null(effects)){
rhs <- stats::setNames(unlist(lapply(ls.object,function(iO){iO$univariate$null})), name.test)
}else{
rhs <- rep(0, NROW(effects))
}
cluster.var <- ls.object[[1]]$object$cluster.var
if(any(cluster.var!=unlist(lapply(ls.object[-1], function(iO){iO$object$cluster.var})))){
stop("Cluster variable differs among objects. \n")
}
newobject$cluster.var <- cluster.var
## ** Try to find unique names
## *** for the model parameters
## find all names
ls.testname <- lapply(ls.object, function(iO){iO$multivariate$test})
vec.testname <- unlist(ls.testname)
if(!is.null(name)){
outcome <- name
}else if(all(sapply(ls.testname, function(iO){all(is.character(iO))})) && all(duplicated(vec.testname)==FALSE)){
outcome <- vec.testname
}else{
outcome <- sapply(ls.object, function(iO){iO$object$outcome})
if(any(duplicated(outcome))){
outcome <- c(deparse(call$model),
sapply(setdiff(which(names(call)==""),1), function(iIndex){deparse(call[[iIndex]])}))
}
}
## *** for the hypotheses
gridTest <- do.call(rbind,lapply(1:n.object, function(iO){
cbind(outcome = outcome[iO], ls.object[[iO]]$multivariate[,c("type","test"),drop=FALSE])
}))
if(is.na(sep) && all(duplicated(gridTest$test)==FALSE)){
gridTest <- gridTest[,"test",drop=FALSE]
}else if(!is.na(sep) && all(duplicated(outcome)==FALSE)){
gridTest <- gridTest[,"outcome",drop=FALSE]
}else if(all(duplicated(gridTest$test)==FALSE)){
gridTest <- gridTest[,"test",drop=FALSE]
}else if(all(duplicated(gridTest$type)==FALSE)){
gridTest <- gridTest[,"type",drop=FALSE]
}else{
if(all(duplicated(gridTest$outcome))){
gridTest$outcome <- NULL
}
if(all(duplicated(gridTest$type))){
gridTest$type <- NULL
}
if(all(duplicated(gridTest$test))){
gridTest$test <- NULL
}
}
col.nametest <- colnames(gridTest)
name.test <- unique(nlme::collapse(gridTest[,col.nametest,drop=FALSE], as.factor = FALSE))
## ** Extract elements from anova object
## *** univariate Wald test
ls.univariate <- lapply(1:n.object,function(iO){ ## iO <- 1
iTable <- cbind(outcome = unname(outcome[iO]), ls.object[[iO]]$univariate)
iTable$name.test <- factor(nlme::collapse(iTable[,col.nametest,drop=FALSE], as.factor = FALSE),levels = name.test)
rownames(iTable) <- rownames(ls.object[[iO]]$univariate)
attr(iTable,"backtransform") <- attr(ls.object[[iO]]$univariate, "backtransform")
return(iTable)
})
table.univariate <- do.call(rbind,ls.univariate)
attr(table.univariate, "backtransform") <- unlist(lapply(ls.univariate,attr,"backtransform"))
## *** cluster
ls.cluster <- lapply(ls.object, function(iO){iO$object$cluster})
seq.cluster <- unique(unlist(ls.cluster))
n.cluster <- length(seq.cluster)
newobject$cluster <- seq.cluster
independence <- (cluster.var=="XXclusterXX" || all(duplicated(unlist(ls.cluster))==FALSE)) ## each individual only appear once
if(all(sapply(ls.object, function(iO){attr(iO$args$robust,"call")})==FALSE)){
newtable.args$robust <- independence==FALSE
}
## *** estimate
beta.estimate <- stats::setNames(table.univariate$estimate, name.modelparam)
if(independence){
beta.vcov <- as.matrix(do.call(Matrix::bdiag,lapply(ls.object, "[[", "vcov")))
if(all(diag(contrast)==1) && all(contrast[upper.tri(contrast, diag = FALSE)] == 0) && all(contrast[lower.tri(contrast, diag = FALSE)] == 0)){
beta.df <- table.univariate$df
}else{
beta.df <- rep(Inf, length(beta.estimate))
newtable.args$df <- FALSE
}
}else{
if(any(sapply(ls.object, function(iO){is.null(iO$iid)}))){
stop("Cannot evaluate the joint influence function of the variance-covariance parameters. \n",
"Consider using method.fit=\"ML\" when calling lmm. \n")
}
M.iid <- matrix(0, nrow = n.cluster, ncol = length(beta.estimate),
dimnames = list(seq.cluster, names(beta.estimate)))
for(iO in 1:n.object){
M.iid[ls.cluster[[iO]],rownames(ls.object[[iO]]$univariate)] <- ls.object[[iO]]$iid
}
if(any(is.na(M.iid))){
beta.vcov <- tcrossprod(sqrt(apply(M.iid^2, 2, sum, na.rm = TRUE))) * stats::cor(M.iid, use = "pairwise")
## usually better compared to formula 11.43 from chapter 11.4 of the book High-dimensional statistics by WAINWRIGHT
## iIDD0 <- M.iid/(1-mean(is.na(M.iid)))
## iIDD0[is.na(M.iid)] <- 0
## out$vcov <- crossprod(iIDD0) - mean(is.na(M.iid))*diag(diag(crossprod(iIDD0)))
## out$vcov - crossprod(M.iid)
}else{
beta.vcov <- crossprod(M.iid)
}
if(any(!is.infinite(table.univariate$df)) && any(!is.na(table.univariate$df)) && abs(diff(range(table.univariate$df, na.rm = TRUE)))<0.1){
beta.df <- rep(mean(table.univariate$df), length(beta.estimate))
newtable.args$df <- TRUE
}else{
beta.df <- rep(Inf, length(beta.estimate))
newtable.args$df <- FALSE
}
}
dimnames(beta.vcov) <- list(name.modelparam,name.modelparam)
names(beta.df) <- name.modelparam
## ** Combine elements
## *** multivariate tests
outSimp <- .simplifyContrast(contrast, rhs = rhs) ## remove extra lines
C.vcov.C <- outSimp$C %*% beta.vcov %*% t(outSimp$C)
C.vcov.C_M1 <- try(solve(C.vcov.C), silent = TRUE)
if(inherits(C.vcov.C_M1,"try-error")){
multistat <- NA
attr(multistat,"error") <- "\n Could not invert the covariance matrix for the proposed contrast."
}else{
multistat <- as.double(t(outSimp$C %*% beta.estimate - outSimp$rhs) %*% C.vcov.C_M1 %*% (outSimp$C %*% beta.estimate - outSimp$rhs))/outSimp$dim
}
Utype2 <- unique(unlist(lapply(ls.object, function(iO){iO$multivariate$type})))
vec.null <- paste0(rownames(outSimp$C),"=",outSimp$rhs)
newtable.multivariate <- data.frame(type = ifelse(length(Utype2)>1,"all",Utype2),
test = 1,
null = paste(vec.null, collapse = ", "),
statistic = multistat,
df.num = outSimp$dim,
df.denom = Inf,
p.value = 1 - stats::pf(multistat, df1 = n.test, df2 = Inf))
## *** univariate
if(!is.null(effects)){
C.vcov.C <- contrast %*% beta.vcov %*% t(contrast)
newtable.univariate <- data.frame(outcome = NA,
type = ifelse(length(Utype2)>1,"all",Utype2),
test = 1,
estimate = (contrast %*% beta.estimate)[,1],
se = sqrt(diag(C.vcov.C)),
df = Inf,
statistic = NA,
lower = NA,
upper = NA,
null = rhs,
p.value = NA)
rownames(newtable.univariate) <- rownames(contrast)
newtable.univariate$outcome <- lapply(1:NROW(contrast), function(iRow){table.univariate[colnames(contrast)[contrast[iRow,]!=0],"outcome"]})
newtable.univariate$statistic <- newtable.univariate$estimate/newtable.univariate$se
newtable.args$df <- FALSE
## take care of parameters subject to transformation for statistical inference
if(!is.null(attr(table.univariate,"backtransform"))){ ## previously already a linear combination
param.untransformed <- stats::setNames(attr(table.univariate,"backtransform"),name.modelparam)
attr(newtable.univariate,"backtransform") <- (contrast %*% param.untransformed[colnames(contrast)])[,1]
}else if(newtable.args$backtransform && (any(contrast %in% 0:1 == FALSE) || any(rowSums(contrast != 0)>1))){ ## only now a linear combinations (backtransform no more ok)
index.nobacktransform <- apply(contrast, MARGIN = 1, function(iRow){sum(iRow %in% 0:1 == FALSE) + sum(iRow != 0)>1})
param.untransformed <- stats::setNames(unlist(lapply(ls.object, coef.Wald_lmm, backtransform = TRUE)), name.modelparam)
attr(newtable.univariate,"backtransform") <- (contrast %*% param.untransformed[colnames(contrast)])[index.nobacktransform,1]
}
}else{
C.vcov.C <- beta.vcov
newtable.univariate <- table.univariate
## update se and statistic in case robust standard errors have been used
newtable.univariate$se <- sqrt(diag(beta.vcov))
newtable.univariate$statistic <- newtable.univariate$estimate/newtable.univariate$se
}
if(newtable.args$df){
e.glht <- list(linfct = contrast, rhs = rhs,
coef = table.univariate$estimate, vcov = beta.vcov, df = floor(stats::median(newtable.univariate$df)), alternative = alternative)
if(!is.na(default$df)){
e.glht$df <- pmax(e.glht$df, default$df)
}
}else{
e.glht <- list(linfct = contrast, rhs = rhs,
coef = table.univariate$estimate, vcov = beta.vcov, df = FALSE, alternative = alternative)
}
class(e.glht) <- "glht"
## ** retrieve parameter associated to each test
nParam.hypo <- unlist(lapply(ls.glht,function(iO){rowSums(iO$linfct!=0)}))
if(all(contrast %in% 0:1) && all(rowSums(contrast)==1) && all(nParam.hypo==1)){ ## single parameter per test
index.param <- which(contrast==1, arr.ind = TRUE)[,"col"]
newtable.univariate[,"parameter"] <- names(nParam.hypo[index.param])
}else{
ls.names <- do.call(c, lapply(ls.glht, FUN = function(iO){
lapply(1:NROW(iO$linfct), FUN = function(iRow){colnames(iO$linfct)[iO$linfct[iRow,]!=0]})
}))
if(all(contrast %in% 0:1) && all(rowSums(contrast)==1)){ ## several parameters per test but all from the same model
newtable.univariate$parameter <- ls.names
}else if(length(unique(ls.names))==1){
if(length(ls.names[[1]])==1){ ## contrasting one parameter accross models
newtable.univariate$parameter <- ls.names[[1]]
}else{ ## contrasting one set of parameters accross models
newtable.univariate$parameter <- ls.names[1]
}
}
}
## ** export
e.glht$linfct.original <- lapply(ls.glht, "[[", "linfct")
out <- list(multivariate = newtable.multivariate[,names(model$multivariate),drop=FALSE],
univariate = newtable.univariate[,c("outcome","parameter",names(model$univariate)),drop=FALSE],
glht = list(all = list("1" = e.glht)),
object = newobject,
args = newtable.args,
vcov = C.vcov.C)
attr(out$univariate,"backtransform") <- attr(newtable.univariate,"backtransform")
attr(out$object,"independence") <- independence
attr(out,"call") <- list(anova = lapply(ls.object,attr,"call"),
rbind = match.call())
class(out) <- append(c("rbindWald_lmm","Wald_lmm"),class(out))
return(out)
}
## * rbind.rbindWald_lmm (code)
##' @export
rbind.rbindWald_lmm <- function(...){
stop("Cannot use rbind on the output of rbind.Wald_lmm. \n")
}
##----------------------------------------------------------------------
### rbind.R ends here
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.