school-styep

R/glance.mlm.R

@@ -1,3 +1,5 @@
+# Removed adj.r.squared as unnecessary in multivariate context
+
 #' Glance at an mlm object
 #' 
 #' This function takes an "mlm" object, fit by \code{\link[stats]{lm}} with a multivariate response.
@@ -6,13 +8,12 @@
 #' In the multivariate case, it returns a \code{\link[tibble]{tibble}} with one row for each
 #' response variable, containing goodness of fit measures, F-tests and p-values.
 #' 
-#' @param x An \code{mlm} object created by \code{\link[stats]{lm}}, i.e., with a multivariate response.
+#' @param x An \code{"mlm"} object created by \code{\link[stats]{lm}}, i.e., with a multivariate response.
 #' @param ... Additional arguments. Not used.
 #' @method glance mlm
 #' @return A \code{\link[tibble]{tibble}} with one row for each response variable and the columns:
 #' \describe{
 #'   \item{\code{r.squared}}{R squared statistic, or the percent of variation explained by the model.}
-#'   \item{\code{adj.r.squared}}{Adjusted R squared statistic, which is like the R squared statistic except taking degrees of freedom into account.}
 #'   \item{\code{sigma}}{Estimated standard error of the residuals}
 #'   \item{\code{fstatitic}}{Overall F statistic for the model}
 #'   \item{\code{numdf}}{Numerator degrees of freedom for the overall test}
@@ -41,7 +42,7 @@ glance.mlm <- function(x, ...) {
   sumry <- summary(x)
   response <- sub("Response ", "", names(sumry))
   # basic model summary statistics
-  stats <- purrr::map_dfr(sumry, magrittr::extract, c("r.squared", "adj.r.squared", "sigma"))
+  stats <- purrr::map_dfr(sumry, magrittr::extract, c("r.squared", "sigma"))
 
   # get the F statistic & df for each response
   fstats <- purrr::map(sumry, magrittr::extract, c("fstatistic"))

dev/uniStats.R

@@ -0,0 +1,46 @@
+# Univariate test statistics for a mlm
+
+uniStats <- function(x, ...) {
+  if (!inherits(x, "mlm")) stop("This function is only for 'mlm' objects")
+  SS <- summary(x)
+  nr <- length(SS)
+  result <- as.data.frame(matrix(0, nrow=nr, 5))
+  for (i in 1:nr) {
+    result[i,1] <- SS[[i]]$r.squared
+    f <- SS[[i]]$fstatistic
+    result[i,2:4] <- f
+    result[i,5] <- pf(f[1], f[2], f[3], lower.tail=FALSE)
+  }
+  result$stars <- c(gtools:::stars.pval(result[,5]))
+  result[,5] <- format.pval(result[,5], eps=0.001)
+  colnames(result) <- c("R^2", "F", "df1", "df2", "Pr (>F)", " ")
+
+  resp <- sub("Response ", "", names(SS))  
+  result <- cbind(response = resp, result)
+  result
+
+
+}
+
+if (FALSE) {
+
+  data(Wine, package="candisc")
+  Wine.mod <- lm(as.matrix(Wine[, -1]) ~ Cultivar, data=Wine)
+  uniStats(Wine.mod)
+
+
+ss <- summary(Wine.mod)
+UniStats <- as.data.frame(matrix(0, nrow=length(ss), 5))
+for (i in 1:length(ss)) {
+  UniStats[i,1] <- ss[[i]]$r.squared
+  f <- ss[[i]]$fstatistic
+  UniStats[i,2:4] <- f
+  UniStats[i,5] <- pf(f[1], f[2], f[3], lower.tail=FALSE)
+}
+
+rownames(UniStats) <- sub("Response ", "", names(ss))
+UniStats$stars <- c(gtools:::stars.pval(UniStats[,5]))
+UniStats[,5] <- format.pval(UniStats[,5], eps=0.001)
+colnames(UniStats) <- c("R^2", "F", "df1", "df2", "Pr (>F)", "")
+UniStats
+}

extra/school-step.R

@@ -0,0 +1,31 @@
+# Stepdown analysis / Choleski
+
+library(heplots)
+library(car)
+
+data(schooldata)
+
+school.mod <- lm(cbind(reading, mathematics, selfesteem) ~ 
+                   education + occupation + visit + counseling + teacher, data=schooldata)
+Anova(school.mod)
+
+pairs(school.mod, fill=TRUE, fill.alpha=0.1)
+
+Y <- schooldata[, c("reading", "mathematics", "selfesteem")] |> as.matrix()
+X <- schooldata[, c("education", "occupation", "visit", "counseling", "teacher")] |> as.matrix()
+
+S <- cov(Y)
+L <- chol(S)
+Linv <- solve(L)
+
+YL <- Y %*% Linv
+colnames(YL) <- c("reading", "math.read", "self.read_math")
+
+school.mat.mod <- lm(Y ~ X)
+Anova(school.mat.mod)
+
+school.step.mod <- lm(YL ~ X[,1] + X[,2] + X[,3] + X[,4] + X[,5])
+Anova(school.mat.mod)
+
+heplot(school.step.mod)
+pairs(school.step.mod)