R/reg.2Relations.R

In lessR: Less Code with More Comprehensive Results

.reg2Relations <-
function(lm.out, n.keep,
         collinear, subsets, best_sub, max.sublns,
         in.data.frame, sterrs, MSW) {

  nm <- all.vars(lm.out$terms)  # names of vars in the model
  n.vars <- length(nm)
  n.pred <- n.vars - 1L
  n.obs <- nrow(lm.out$model)

  # Correlations among the model variables are no longer reported as a
  # text matrix.  They are shown visually in the scatterplot matrix, with
  # the coefficients in the upper triangle (see the  scatter_coef  option).

  # -------------------------
  # collinearity

  txcln <- ""
  if (collinear) {
    tx <- character(length = 0)

    if (is.null(options()$knitr.in.progress)) {
      tx[length(tx)+1] <- "-- Collinearity"
      tx[length(tx)+1] <- ""
    }

    vif <- numeric(length = 0)
    tol <- numeric(length = 0)
    for (i in seq_len(n.pred)) {
      v <- var(lm.out$model[i+1])
      vif[i] <- (v * (n.keep-1) * sterrs[i+1]^2) / MSW
      tol[i] <- 1 / vif[i]
    }

    out <- cbind(tol, vif)
    colnames(out) <- c("Tolerance", "      VIF")
    rownames(out) <- nm[seq(2L, length(nm))]

    txcln <- .prntbl(out, 3, cc=NULL)
    for (i in seq_along(txcln)) tx[length(tx)+1] <- txcln[i]

    txcln <- tx
  } 

    else {  # !collinear
      vif <- NA
      tol <- NA
    }

  # -------------------------
  # all possible subsets of predictor variables    
  txsbs <- ""
  if (subsets) {
    tx <- character(length = 0)

    if (is.null(options()$knitr.in.progress)) {
      tx[length(tx)+1] <- "-- Best Subset Regression Models"
      if (n.pred > 5)
        tx[length(tx)+1] <- "up to 10 subsets of each number of predictors"
      tx[length(tx)+1] <- ""
    }

    X <- data.frame(lm.out$model[nm[seq(2,n.vars)]])
    Y <- lm.out$model[[nm[1]]]  # atomic vector for leaps
    lp.out <- leaps(X, Y, method=best_sub)  # leaps function
    md <- lp.out$which  # md is logical
    rownames(md) <- seq_len(nrow(md))  # matrix md does not have proper row names
    # gives 0, 1
    if (best_sub == "adjr2") {
      models <- data.frame(md, lp.out$adjr2, lp.out$size-1)
      names(models) <- c(names(X), "R2adj", "X's")
      mod.srt <- models[order(models$R2adj, decreasing=TRUE),]

    }
    else {
      models <- data.frame(md, lp.out$Cp, lp.out$size-1)
      names(models) <- c(names(X), "   Cp", "X's")
      mod.srt <- models[order(models$'   Cp', decreasing=FALSE),]
    }
    best_sub <- if (best_sub == "adjr2") "R2adj" else best_sub
    names(mod.srt)[ncol(mod.srt)-1L] <- paste("   ", best_sub, sep="")
    lines <- min(max.sublns, nrow(mod.srt))

    # width of data columns
    max.ln <- integer(length=ncol(X))
    for (i in seq_len(ncol(X))) {
      ln.nm <- nchar(colnames(X)[i])
      max.ln[i] <- ln.nm + 1
      if (max.ln[i] < 4) max.ln[i] <- 4L
    }

    tx[length(tx)+1] <- ""  # build the line of variable names
    for (i in seq_len(n.pred)) {
      if (i <= n.pred) ww <- max.ln[i]
      tx[length(tx)] <- paste(tx[length(tx)],
          .fmtc(names(mod.srt)[i], w=ww), sep="")
      nms <- tx[length(tx)]
    }
    for (i in seq_len(2L)) {
      if (i == 1) ww <- 9L
      if (i == 2) ww <- 7L
      tx[length(tx)] <- paste(tx[length(tx)], .fmtc(names(mod.srt)[n.pred+i],
                              w=ww), sep="")
    }
    for (i in seq_len(lines)) {
      if (lines > 40) if (i %% 30 == 0) tx[length(tx)+1] <- nms
      tx[length(tx)+1] <- ""
      for (j in seq_len(n.pred))
        tx[length(tx)] <- paste(tx[length(tx)], .fmti(mod.srt[i,j],
          w=max.ln[j]), sep="")
      tx[length(tx)] <- paste(tx[length(tx)], .fmt(mod.srt[i,n.pred+1],
          d=3, w=8)) # R2adj
      tx[length(tx)] <- paste(tx[length(tx)], .fmti(mod.srt[i,n.pred+2],
          w=6))  # num X's
    }

    if (nrow(mod.srt) > max.sublns)
      tx[length(tx)+1] <- paste(
        "\n>>> Only first", lines, "of", nrow(mod.srt), "rows printed\n",
        "   To indicate more, add subset=n, where n is the number of lines")

    tx[length(tx)+1] <- ""
    tx[length(tx)+1] <- paste("[based on Thomas Lumley's leaps function",
                        "from the leaps package]")

    txsbs <- tx
  }  # end subsets

  return(list(out_collinear=txcln, out_subsets=txsbs, tol=tol, vif=vif))

}