R/pairwise.R

In walmes/wzRfun: Walmes Zeviani's collection of functions

#' @author Walmes Zeviani, \email{walmes@@ufpr.br}.
#' @export
#' @title Generate Matrix of All Pairwise Comparisons (Tukey contrasts)
#' @description This function takes a matrix where each line defines a
#'     linear function of the parameters to estimate a marginal mean
#'     (aka least squares mean) and return the matrix that define the
#'     contrasts among these means. All pairwise contrasts are returned
#'     (aka Tukey contrasts). The matrix with these contrasts can be
#'     passed to \code{\link[multcomp]{glht}()} to estimate them or used
#'     in explicit matricial calculus.
#' @param lfm a \eqn{k \times p} matrix where each line defines a linear
#'     function to estimate a lsmean. In general, these matrices are
#'     obtained by using \code{\link[doBy]{LE_matrix}()}.
#' @param lev a character vector with length equals to the numbers of
#'     lines of \code{lfm} matrix, (\eqn{k}). Default is \code{NULL} and
#'     the row names of code{lfm} is used. If row names is also
#'     \code{NULL}, incremental integer values are used to identify the
#'     comparisons.
#' @return a \eqn{K\times p} matrix with the linear functions that
#'     define all pairwise contrasts. \eqn{K} is \eqn{{k}\choose{2}}.
#' @seealso \code{\link{apmc}()}, \code{\link[doBy]{LE_matrix}()}.
#' @examples
#'
#' X <- diag(3)
#' rownames(X)
#' apc(X)
#'
#' rownames(X) <- letters[nrow(X):1]
#' apc(X)
#'
#' apc(X, lev = LETTERS[1:nrow(X)])
#'
#' # Objects from doBy::LE_matrix() have an "grid" attribute.
#' attr(X, "grid") <- data.frame(n = LETTERS[1:nrow(X)])
#' rownames(X) <- NULL
#' apc(X)
#'
apc <- function(lfm, lev = NULL) {
    nlev <- nrow(lfm)
    rn <- rownames(lfm)
    a <- attr(lfm, "grid")
    if (is.null(lev)) {
        if (!is.null(a)) {
            lev <- apply(a, 1, paste, collapse = ":")
        } else if (!is.null(rn)) {
            lev <- rn
        } else {
            lev <- as.character(1:nlev)
        }
    }
    cbn <- utils::combn(seq_along(lev), 2)
    M <- lfm[cbn[1, ], ] - lfm[cbn[2, ], ]
    if (is.vector(M)) {
        dim(M) <- c(1, length(M))
    }
    rownames(M) <- paste(lev[cbn[1, ]], lev[cbn[2, ]], sep = "-")
    return(M)
}

#' @author Walmes Zeviani, \email{walmes@@ufpr.br}.
#' @export
#' @title A Wraper of glht to Get All Pairwise Mean Comparisons
#' @description This function performs all pairwise compararisons among
#'     means returning pontual and intervalar estimates followed by
#'     letters to easy discriminate values. It is in fact a wraper of
#'     \code{\link[multcomp]{glht}()}.
#' @param X a matrix where each line is a linear function of the model
#'     parameters to estimate a least squares mean. In most pratical
#'     cases, it is an object from the \code{\link[doBy]{LE_matrix}()}.
#' @param model a model with class recognized by
#'     \code{\link[multcomp]{glht}()}.
#' @param focus a string with the name of the factor which levels are
#'     being compared.
#' @param test a p-value correction method. See
#'     \code{\link[stats]{p.adjust.methods}()}.
#' @param level the experimentwise significance level for the multiple
#'     comparisons. The individual coverage of the confidence interval
#'     is \code{1-level}. Default is 0.05.
#' @param cld2 Logical, if \code{TRUE} uses the \code{\link{cld2}()}
#'     functions, otherwise uses the \code{\link[multcomp]{cld}()}
#'     function.
#' @return a \code{data.frame} with interval estimates and compact
#'     letter display for the means comparisons.
#' @seealso \code{\link{apc}()}, \code{\link[doBy]{LE_matrix}()},
#'     \code{\link[multcomp]{glht}()}.
#' @examples
#'
#' library(doBy)
#' library(multcomp)
#'
#' # Single factor.
#' m0 <- lm(weight ~ feed, data = chickwts)
#' anova(m0)
#'
#' # Prepare the matrix to estimate lsmeans.
#' L <- LE_matrix(m0, effect = "feed")
#' rownames(L) <- levels(chickwts$feed)
#' apmc(L, model = m0, focus = "feed", test = "fdr")
#'
#' data(warpbreaks)
#'
#' # Two factors (complete factorial).
#' m1 <- lm(breaks ~ wool * tension, data = warpbreaks)
#' anova(m1)
#'
#' L <- LE_matrix(m1, effect = c("wool", "tension"))
#' attributes(L)
#'
#' Ls <- by(L, INDICES = attr(L, "grid")$tension, FUN = as.matrix)
#' Ls <- lapply(Ls, "rownames<-", levels(warpbreaks$wool))
#'
#' # Comparing means of wool in each tension.
#' lapply(Ls, apmc, model = m1, focus = "wool",
#'        test = "single-step", level = 0.1)
#'
#' # Two factors (incomplete factorial).
#' warpbreaks <- subset(warpbreaks, !(tension == "H" & wool == "A"))
#' xtabs(~tension + wool, data = warpbreaks)
#'
#' # There is NA in the estimated parameters.
#' m2 <- lm(breaks ~ wool * tension, data = warpbreaks)
#' coef(m2)
#'
#' X <- model.matrix(m2)
#' b <- coef(m2)
#'
#' X <- X[, !is.na(b)]
#'
#' # unique(X)
#'
#' # Uses the full estimable model matriz.
#' m3 <- update(m2, . ~ 0 + X)
#'
#' # These models are in fact the same.
#' anova(m2, m3)
#'
#' # LS matrix has all cells.
#' L <- LE_matrix(m2, effect = c("wool", "tension"))
#' g <- attr(L, "grid")
#' L <- L[, !is.na(b)]
#' i <- 5
#' L <- L[-i, ]
#' g <- g[-i, ]
#'
#' rownames(L) <- make.names(g$tension, unique = FALSE)
#' Ls <- split.data.frame(L, g$wool)
#'
#' # LSmeans with MCP test.
#' lapply(Ls, apmc, model = m3, focus = "tension",
#'        test = "single-step", level = 0.1, cld2 = TRUE)
#'
#' # Sample means.
#' aggregate(breaks ~ tension + wool, data = warpbreaks, FUN = mean)
#'
apmc <- function(X,
                 model,
                 focus,
                 test = "single-step",
                 level = 0.05,
                 cld2 = FALSE) {
    if (is.null(rownames(X))) {
        stop("The X matrix must have row names.")
    }
    Xc <- apc(X)
    g <- multcomp::glht(model, linfct = X)
    ci <- stats::confint(g,
                         level = 1 - level,
                         calpha = multcomp::univariate_calpha())$confint
    ci <- as.data.frame(ci)
    names(ci) <- tolower(names(ci))
    names(ci)[1] <- "fit"
    h <- summary(multcomp::glht(model, linfct = Xc),
                 test = multcomp::adjusted(type = test))
    h$type <- "Tukey"
    h$focus <- focus
    if (cld2) {
        ci$cld <- cld2(h, level = level)$mcletters$Letters
    } else {
        ci$cld <- multcomp::cld(h, level = level,
                                decreasing = TRUE)$mcletters$Letters
    }
    ci <- cbind(rownames(ci), ci)
    names(ci)[1] <- focus
    rownames(ci) <- NULL
    return(ci)
}

#' @name cld2
#' @author Walmes Zeviani, \email{walmes@@ufpr.r}.
#' @export
#' @title Modified Compact Letter Display to Irregular Designs
#' @description This functions get the compact letter display for
#'     objects of class \code{"glht"}. Modification was done to get the
#'     letters to design with missing cells, non completelly crossed
#'     factorial designs and nested factorial designs. These models are
#'     usually declared by a model matrix to have all effects
#'     estimated. It is assumed that Tukey contrasts was used.
#' @param object an object returned by \code{\link[multcomp]{glht}()}. It
#'     is assumed that the matrix used as the \code{linfct} argument in
#'     \code{glht} corresponds to a matrix to get Tukey contrasts of
#'     least squares means.
#' @param level the nominal significance level.
#' @return an object of class \code{"cld"} with letters to resume mean
#'     comparisons.
#' @seealso \code{\link{apc}()}, \code{\link[doBy]{LE_matrix}()},
#'     \code{\link[multcomp]{glht}()}.
# @import multcomp
#' @examples
#'
#' # Toy data 1: experiment with cultivars in several locations.
#' td1 <- expand.grid(loc = gl(5, 1),
#'                    block = gl(3, 1),
#'                    cult = LETTERS[1:6])
#' td1 <- subset(td1, !(loc == 1 & cult == "A"))
#' td1 <- subset(td1, !(loc == 2 & cult == "B"))
#' xtabs(~loc + cult, td1)
#' td1$y <- seq_len(nrow(td1))
#'
#' library(lme4)
#'
#' # Fit the mixed model.
#' m0 <- lmer(y ~ loc * cult + (1 | loc:block), data = td1)
#' logLik(m0)
#'
#' # The same model but without rank deficience.
#' td1$loccult <- with(td1, interaction(loc, cult, drop = TRUE))
#' m1 <- lmer(y ~ loccult + (1 | loc:block), data = td1)
#' logLik(m1)
#'
#' library(doBy)
#'
#' X <- LE_matrix(lm(nobars(formula(m1)), data = td1), effect = "loccult")
#' rownames(X) <- levels(td1$loccult)
#' dim(X)
#'
#' Xs <- X[grepl(x = rownames(X), "^1\\."),]
#' Xc <- apc(Xs)
#'
#' library(multcomp)
#'
#' g <- summary(glht(m1, linfct = Xc), test = adjusted(type = "fdr"))
#'
#' cld2(g)
#'
#' confint(glht(m1, linfct = Xs), calpha = univariate_calpha())
#'
cld2 <- function(object,
                 level = 0.05) {
    lvl_order <- unique(unlist(
        strsplit(rownames(object$linfct), "-")))
    signif <- (object$test$pvalues < level)
    ret <- list()
    ret$signif <- signif
    ret$comps <- do.call(rbind,
                         strsplit(rownames(object$linfct), "-"))
    ret$mcletters <-
        insert_absorb(x = signif,
                      decreasing = TRUE,
                      comps = ret$comps,
                      lvl_order = lvl_order)
    class(ret) <- "cld"
    return(ret)
}

#' @name ordered_cld
#' @author Walmes Zeviani, \email{walmes@@ufpr.r}.
#' @export
#' @title Order Letters According to Numeric Vector
#' @description This function order the letters in the compact letter
#'     display to the highest estimate receive the letter \code{a}. This
#'     is a convetion in most software for analysis of experiments.
#' @param let Character vector with the letters returned by
#'     \code{\link[multcomp]{cld}()} or \code{\link{cld2}()}.
#' @param means Numeric vector with the corresponding estimates in which
#'     the highest value will have the letter \code{a}.
#' @return A character vector with the letters rearranged.
#' @seealso \code{\link{cld2}()}.
#' @examples
#'
#' # Toy data.
#' set.seed(4321)
#' td <- data.frame(trt = rep(sample(1:8), each = 5))
#' td$y <- rnorm(nrow(td), mean = sort(td$trt), sd = 2)
#'
#' plot(y ~ trt, data = td)
#'
#' # Fit the model.
#' td$trt <- factor(td$trt)
#' m0 <- lm(y ~ trt, data = td)
#' anova(m0)
#' summary(m0)
#'
#' library(multcomp)
#' library(doBy)
#'
#' X <- LE_matrix(m0, effect = "trt")
#' rownames(X) <- levels(td$trt)
#' Xc <- apc(X)
#'
#' g <- summary(glht(m0, linfct = Xc),
#'              test = adjusted(type = "fdr"))
#'
#' res <- data.frame(trt = levels(td$trt),
#'                   mean = X %*% coef(m0))
#'
#' let <- cld2(g)
#' res$cld2 <- let$mcletters$Letters
#' res[order(res$mean, decreasing = TRUE), ]
#'
#' res$let2 <- ordered_cld(res$cld2, res$mean)
#' res[order(res$mean, decreasing = TRUE), ]
#'
#' \dontrun{
#'
#' library(latticeExtra)
#' library(grid)
#'
#' ci <- as.data.frame(
#'     confint(glht(m0, linfct = X),
#'             calpha = univariate_calpha())$confint)
#' ci <- cbind(res, ci)
#'
#' segplot(reorder(trt, Estimate) ~ lwr + upr,
#'         centers = Estimate,
#'         data = ci,
#'         draw = FALSE,
#'         cld = ci$let2,
#'         par.settings = list(layout.widths = list(right.padding = 7))) +
#'     layer(panel.text(x = centers,
#'                      y = z,
#'                      labels = sprintf("%0.2f %s",
#'                                       centers,
#'                                       cld),
#'                      pos = 3))
#'
#' ocld <- with(ci[order(ci$Estimate), ],
#'      ordered_cld(cld2, Estimate))
#' x <- attr(ocld, "ind")
#' index <- which(x, arr.ind = TRUE)
#' trellis.focus("panel", column = 1, row = 1, clip.off = TRUE)
#' xcor <- 1.03 + (index[, 2] - 1)/50
#' grid.segments(x0 = unit(xcor, "npc"),
#'               x1 = unit(xcor, "npc"),
#'               y0 = unit(index[, 1] + 0.5, units = "native"),
#'               y1 = unit(index[, 1] - 0.5, units = "native"),
#'               gp = gpar(lwd = 2, col = "blue"))
#' trellis.unfocus()
#'
#' }
#'
ordered_cld <- function(let, means = let) {
    or <- order(means, decreasing = TRUE)
    let <- as.character(let[or])
    s <- strsplit(let, "")
    u <- unlist(s)
    ul <- unique(u)
    UL <- LETTERS[seq_along(ul)]
    l <- sapply(s, FUN = function(i) {
        paste(sort(UL[match(i, table = ul)]), collapse = "")
    })
    x <- tolower(l[order(or)])
    UL <- tolower(UL)
    attr(x, "match") <- cbind("before" = ul,
                              "after" = UL)
    attr(x, "ind") <- sapply(UL, FUN = grepl, x = x)
    return(x)
}

#' @name radial_cld
#' @importFrom graphics lines plot points segments text
#' @author Walmes Zeviani, \email{walmes@@ufpr.r}.
#' @export
#' @title Radial Plot for a Compact Letter Display Vector
#' @description This function does a radial plot based on the vector of
#'     letters resulted from pairwise comparisons.
#' @param cld Character vector with strings of letters that indicates
#'     which pair of treatment cells are not different.
#' @param labels Vector of text to be annotated next each point.
#' @param col Vector of colors to be used in the segments that joint
#'     points.
#' @param means Numeric vector with the estimated means of treatment
#'     cells. It is used to place points at distances proportional to
#'     the differences on means.
#' @param perim Logical value (default is \code{FALSE}) that indicates
#'     weather draw or not a circle in the perimeter passing by the
#'     points.
#' @param legend Logical value (default is \code{TRUE}) that indicates
#'     weather daraw or not the legend.
#' @return None is returned, only the plot is done.
#' @seealso \code{\link{cld2}()}.
#' @examples
#'
#' set.seed(4321)
#' td <- data.frame(trt = rep(sample(1:20), each = 5))
#' td$y <- rnorm(nrow(td), mean = 0.15 * sort(td$trt), sd = 1)
#'
#' plot(y ~ trt, data = td)
#'
#' # Fit the model.
#' td$trt <- factor(td$trt)
#' m0 <- lm(y ~ trt, data = td)
#' anova(m0)
#' summary(m0)
#'
#' library(multcomp)
#' library(doBy)
#'
#' X <- LE_matrix(m0, effect = "trt")
#' rownames(X) <- levels(td$trt)
#'
#' ci <- apmc(X, m0, focus = "trt", test = "fdr")
#' ci$cld <- with(ci, ordered_cld(cld, fit))
#' ci <- ci[order(ci$fit, decreasing = TRUE), ]
#'
#' library(latticeExtra)
#'
#' segplot(reorder(trt, fit) ~ lwr + upr,
#'         centers = fit,
#'         data = ci,
#'         draw = FALSE,
#'         cld = ci$cld) +
#'     layer(panel.text(x = centers,
#'                      y = z,
#'                      labels = sprintf("%0.2f %s",
#'                                       centers,
#'                                       cld),
#'                      pos = 3))
#'
#' radial_cld(cld = ci$cld)
#' radial_cld(cld = ci$cld, means = ci$fit, perim = TRUE)
#' radial_cld(cld = ci$cld, col = 1:3)
#' radial_cld(cld = ci$cld, col = 1:3)
#' radial_cld(cld = ci$cld, labels = sprintf("%0.2f %s", ci$fit, ci$cld))
#'
radial_cld <- function(cld,
                       labels = cld,
                       col = NULL,
                       means = NULL,
                       perim = FALSE,
                       legend = TRUE) {
    if (is.null(means)) {
        s <- seq(from = 0,
                 to = 2 * pi,
                 length.out = length(cld) + 1)[-1]
    } else {
        ext <- (2 * pi)/c(length(cld), 1)
        m <- means - min(means)
        m <- m/max(m)
        s <- ext[1] + diff(ext) * m
        s <- rev(s)
    }
    sincos <- cbind(sin = sin(s), cos = cos(s))
    # Quais as letras únicas formadoras das strings?
    u <- unique(unlist(strsplit(cld, split = "")))
    if (is.null(col)) {
        col <- grDevices::palette()
    }
    if (length(col) != length(u)) {
        warning(paste("Length of vector `col` is different",
                      "of the number of unique letters.",
                      "Colors will be recycled."))
    }
    col <- col[seq_along(u) %% length(col) + 1]
    # Membros da mesma família compartilham a mesma letra.
    fam <- sapply(u, grepl, x = cld)
    plot(x = NULL,
         y = NULL,
         xlim = 1.2 * c(-1, 1),
         ylim = 1.2 * c(-1, 1),
         asp = 1,
         axes = FALSE,
         ann = FALSE)
    if (perim) {
        circ <- seq(0, 2 * pi, length.out = 60)
        lines(x = sin(circ), y = cos(circ), col = "gray", lty = 3)
    }
    for (i in 1:ncol(fam)) {
        cb <- utils::combn(x = which(fam[, i]), m = 2)
        apply(cb,
              MARGIN = 2,
              FUN = function(index) {
                  segments(x0 = sincos[index[1], 1],
                           x1 = sincos[index[2], 1],
                           y0 = sincos[index[1], 2],
                           y1 = sincos[index[2], 2],
                           col = col[i],
                           lwd = 2,
                           lty = 2)
              })
    }
    points(x = sincos[, 1], y = sincos[, 2])
    if (legend) {
        legend("topright",
               legend = u,
               col = col,
               lty = 2,
               lwd = 2,
               bty = "n")
    }
    text(sincos[, 1],
         sincos[, 2],
         labels = labels,
         pos = ifelse(sincos[, 1] > 0, 4, 2))
}