Feature/refactor summary

.github/install_deps.R

@@ -1,7 +1,7 @@
 # R script
 github_deps <- c(
     "iNZightVIT/[email protected]",
-    "iNZightVIT/iNZightPlots@2.12",
+    "iNZightVIT/iNZightPlots@dev",
     "iNZightVIT/[email protected]"     # dependency of iNZightPlots
 )
 

DESCRIPTION

@@ -22,7 +22,7 @@ Imports:
     graphics,
     grDevices,
     grid,
-    iNZightPlots (>= 2.12),
+    iNZightPlots (>= 2.12.6.9000),
     multcomp,
     patchwork,
     stats,

NAMESPACE

@@ -3,6 +3,7 @@
 S3method(compare_models,default)
 S3method(compare_models,svyglm)
 S3method(inzplot,lm)
+S3method(inzsummary,lm)
 S3method(print,inzfactorcomp)
 S3method(print,inzmodelcomp)
 export(Poly)
@@ -13,6 +14,7 @@ export(histogramArray)
 export(iNZightQQplot)
 export(iNZightSummary)
 export(inzplot)
+export(inzsummary)
 export(partialResPlot)
 export(plotlm6)
 import(ggplot2)
@@ -24,3 +26,4 @@ import(grid)
 import(patchwork)
 import(stats)
 importFrom(iNZightPlots,inzplot)
+importFrom(iNZightPlots,inzsummary)

R/inzsummary.R

@@ -0,0 +1,262 @@
+#' inzsummary method
+#'
+#' @importFrom iNZightPlots inzsummary
+#' @name inzsummary
+#' @rdname inzsummary.lm
+#' @export
+NULL
+
+#' Informative Summary Information for Regression Models
+#'
+#' Displays a summary of a regression model fitted using `lm`, `glm`, or `survey::svyglm`.
+#'
+#' @param x an object of class \code{"lm"}, \code{"glm"} or \code{"svyglm"}, usually the result of a call to the corresponding function.
+#' @param method the type of inference to perform: "standard" uses Normal theory, "bootstrap" uses bootstrap resampling.
+#' @param reorder.factors if `FALSE` (default), factor order is left unchanged, if `TRUE`, the most common level is set as the baseline.
+#' @param digits the number of significant digits to display in results
+#' @param symbolic.cor logical, if \code{TRUE}, print the correlations in a symbolic form (see \code{\link{symnum}}), rather than as numbers.
+#' @param signif.stars logical, if \code{TRUE}, \sQuote{significance stars} are printed for each coefficient.
+#' @param exclude a character vector of names of variables to be excluded from the summary output (i.e., confounding variables).
+#' @param exponentiate.ci logical, if \code{TRUE}, the exponential of the confidence intervals will be printed if appropriate (log/logit link or log transformed response)
+#' @param env the environment for evaluating things (e.g., bootstraps)
+#'
+#' @param ... additional arguments, ignored.
+#'
+#' @return Summary information is printed to the console.
+#'
+#' @author Tom Elliott
+#' @md
+#' @export
+inzsummary.lm <- function(x,
+                          method = c("standard", "bootstrap"),
+                          reorder.factors = FALSE,
+                          digits = max(3, getOption("digits") - 3),
+                          symbolic.cor = x$symbolic.cor,
+                          signif.stars= getOption("show.signif.stars"),
+                          exclude = NULL,
+                          exponentiate.ci = FALSE,
+                          ...,
+                          env = parent.frame()) {
+
+    method <- match.arg(method)
+
+    ## section 1: model information
+
+    # response:
+    cat(glue::glue("Linear Model for '{attr(x$model, 'names')[1]}'"), "\n\n")
+
+    # confounding variables/adjusted for:
+    if (!is.null(exclude)) {
+        cat('The model has been adjusted for the',
+            'following confounder(s):\n', sep = ' ')
+        cat('\t')
+        cat(exclude, sep = ', ')
+        cat('\n\n')
+    }
+
+    ## section 2: coefficient matrix
+    coef_matrix(x, method = method, signif.stars = signif.stars, exclude = exclude)
+
+    ## section 3: errors, df, R-squared, etc (model dependent)
+    x.lm <- x
+    x.data <- x.lm$model
+    x <- summary(x)
+    if (isGlm(x.lm)) {
+        cat("\n(Dispersion parameter for ",
+            x.lm$family$family, " family taken to be ",
+            format(x$dispersion), ")\n\n",
+            apply(cbind(paste(format(c("Null", "Residual"),
+                                     justify = "right"), "deviance:"),
+                        format(unlist(x[c("null.deviance", "deviance")]),
+                               digits = max(5, digits + 1)), " on",
+                        format(unlist(x[c("df.null", "df.residual")])),
+                        " degrees of freedom\n"),
+                  1L, paste, collapse = " "), sep = "")
+        if (nzchar(mess <- naprint(x$na.action)))
+            cat("  (", mess, ")\n", sep = "")
+        if (!is.na(x$aic))
+            cat("AIC: ", format(x$aic, digits = max(4, digits + 1)),
+                '\n', sep = '')
+        cat("\n", "Number of Fisher Scoring iterations: ", x$iter,
+            "\n", sep = "")
+        correl <- x$correlation
+        if (!is.null(correl)) {
+            p <- NCOL(correl)
+            if (p > 1) {
+                cat("\nCorrelation of Coefficients:\n")
+                if (is.logical(symbolic.cor) && symbolic.cor) {
+                    print(symnum(correl, abbr.colnames = NULL))
+                }
+                else {
+                    correl <- format(round(correl, 2), nsmall = 2,
+                                     digits = digits)
+                    correl[!lower.tri(correl)] <- ""
+                    print(correl[-1, -p, drop = FALSE], quote = FALSE)
+                }
+            }
+        }
+    } else if (!isCox(x.lm)) {
+        cat("\nResidual standard error:",
+            format(signif(x$sigma, digits)), "on", rdf,
+            "degrees of freedom\n")
+        if(nzchar(mess <- naprint(x$na.action))) cat("  (",mess, ")\n",
+                                                     sep="")
+        if (!is.null(x$fstatistic)) {
+            cat("Multiple R-squared:", formatC(x$r.squared, digits=digits))
+            cat(",\tAdjusted R-squared:",formatC(x$adj.r.squared,
+                                                 digits=digits), "\n")
+        }
+        correl <- x$correlation
+        if (!is.null(correl)) {
+            p <- NCOL(correl)
+            if (p > 1L) {
+                cat("\nCorrelation of Coefficients:\n")
+                if(is.logical(symbolic.cor) && symbolic.cor) {
+                    print(symnum(correl, abbr.colnames = NULL))
+                } else {
+                    correl <-
+                        format(round(correl, 2), nsmall = 2,
+                               digits = digits)
+                    correl[!lower.tri(correl)] <- ""
+                    print(correl[-1, -p, drop=FALSE], quote = FALSE)
+                }
+            }
+        }
+    } else if (isCox(x.lm)) {
+      ## For Cox PH models, just print the last few lines of summary output
+      other.stats <- utils::capture.output(x)
+      s.len <- length(other.stats)
+
+      other.stats <- other.stats[(s.len - 4):(s.len - 1)]
+      cat("\n", other.stats, sep = "\n")
+    }
+    cat("\n")
+    invisible(x)
+}
+
+coef_matrix <- function(x, method, signif.stars, exclude) {
+    z <- summary(x)
+
+    coef <- z$coefficients
+    if (any(aliased <- z$aliased)) {
+        cn <- names(aliased)
+        coef <- matrix(NA, length(aliased), 4, dimnames = list(cn, colnames(coef)))
+        coef[!aliased, ] <- z$coefficients
+    }
+    coefs <- coef[,1,drop=FALSE]
+    err <- coef[,2,drop=FALSE]
+    tval <- coef[,3,drop=FALSE]
+    pval <- coef[,4,drop=FALSE]
+    ci <- confint(x)
+
+    mat <- cbind(
+        est = coefs,
+        err = err,
+        tval = tval,
+        pval = pval,
+        low = ci[,1],
+        upp = ci[,2]
+    )
+
+    terms <- terms(z)
+
+    # Now terms: should be one group of row(s) per column of `terms`
+    dc <- attr(terms, "dataClasses")
+    tf <- attr(terms, "factors")
+    aov <- anova(x)
+    coef_mat <- lapply(colnames(tf),
+        function(term) {
+            if (sum(tf[,term]) == 1) {
+                # just a single variable
+                if (dc[term] == "numeric")
+                    return(list(names = term, rows = mat[term,,drop=FALSE]))
+                # else is a factor; factor row + baseline + other rows from `mat`
+                lvls <- x$xlevels[[term]]
+                rows <- rbind(
+                    c(NA, NA, NA, aov[term,"Pr(>F)"], NA, NA),
+                    c(0, NA, NA, NA, NA, NA),
+                    mat[paste0(term, lvls[-1]), , drop = FALSE]
+                )
+                return(
+                    list(
+                        names = c(term, paste0("  ", lvls)),
+                        rows = rows
+                    )
+                )
+            } else {
+                # an interaction
+                vars <- strsplit(term, ":")[[1]]
+                lvls <- lapply(vars,
+                    function(v) if (dc[v] == "numeric") "" else x$xlevels[[v]][-1]
+                )
+                names(lvls) <- vars
+                labels <- lapply(vars, function(v) paste0(v, lvls[[v]]))
+                labels <- do.call(expand.grid, labels)
+                labels <- apply(as.matrix(labels), 1, paste, collapse = ":")
+                rows <- mat[labels, , drop = FALSE]
+
+                names <- do.call(expand.grid, lvls)
+                names <- paste0("  ", apply(as.matrix(names), 1, paste, collapse = ":"))
+
+                names <- c(term, names)
+                rows <- rbind(c(NA, NA, NA, aov[term,"Pr(>F)"], NA, NA), rows)
+                return(list(names = names, rows = rows))
+            }
+        }
+    )
+    names <- do.call(c, lapply(coef_mat, function(x) x$names))
+    coef_mat <- do.call(rbind, lapply(coef_mat, function(x) x$rows))
+
+    # Intercept:
+    if (attr(terms, "intercept") == 1) {
+        coef_mat <- rbind(mat["(Intercept)", , drop=FALSE], coef_mat)
+        names <- c("Intercept", names)
+    }
+    rownames(coef_mat) <- names
+    colnames(coef_mat) <- c("Estimate", "Std. error", "t value", "p value", "lower", "upper")
+
+    class(coef_mat) <- "inzcoefmat"
+    cat("Coefficients:\n")
+    print(coef_mat, signif.stars = signif.stars)
+}
+
+print.inzcoefmat <- function(x, digits = 3, signif.stars = TRUE, ...) {
+    # assume (for now) that everything is there, we just need to format the column values
+    # cols 1-3,5-6 formatted normally; col 4 formatted as p-value
+
+    if (signif.stars) {
+        pval <- x[, 4]
+        Signif <- symnum(pval, corr = FALSE, na = FALSE,
+            cutpoints = c(0,  .001,.01,.05, .1, 1),
+            symbols   =  c("***","**","*","."," "))
+    }
+
+    mat <- apply(x, 2, function(col) {
+        format(col, digits = digits, scientific = FALSE)
+    })
+    mat <- matrix(mat, ncol = ncol(x))
+    mat[,4] <- format.pval(x[,4], digits = digits)
+
+    mat <- rbind(colnames(x), mat)
+    mat <- matrix(
+        apply(mat, 2, function(col) format(col, justify = "right")),
+        nrow = nrow(mat)
+    )
+    mat <- gsub("NA", "-", mat)
+
+    mat <- cbind(mat[,1:4], c("", as.character(Signif)), mat[,5:6])
+
+    mat <- cbind(c("", rownames(x)), mat)
+    mat <- matrix(
+        apply(mat, 2, function(col) format(col, justify = "left")),
+        nrow = nrow(mat)
+    )
+
+    mat <- apply(mat, 1,
+        function(x) paste0("   ", paste(x, collapse = "   "))
+    )
+
+    cat(mat, sep = "\n")
+    if (signif.stars)
+        cat("---\nSignif. codes: ", attr(Signif, "legend"), "\n")
+}

R/summary.R

@@ -397,8 +397,7 @@ iNZightSummary <- function (x, method = "standard", reorder.factors = FALSE,
                       # substring matching an existing variable
                       # name, and there is a level present
 
-                        if (current.var == summary.row.subs &&
-                            is.level.of.cvar) {
+                        if (current.var == summary.row.subs && is.level.of.cvar) {
                             levels.of.cvar <- levels(x.data[, current.var])
                             base.level <- levels.of.cvar[1]
                             nlines.to.add <- length(levels.of.cvar) + 1
@@ -425,6 +424,10 @@ iNZightSummary <- function (x, method = "standard", reorder.factors = FALSE,
                                             }
                                         pvalue <- type3pval
                                     }
+                                    # print(coefs.copy)
+                                    # print(name.k)
+                                    # print(k)
+                                    # print(i)
                                     coefs.copy <-
                                         insert.lines(
                                             name.k, k,