diff --git a/R/glance.mlm.R b/R/glance.mlm.R index 558055f..fdd0394 100644 --- a/R/glance.mlm.R +++ b/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")) diff --git a/dev/uniStats.R b/dev/uniStats.R new file mode 100644 index 0000000..b7d322c --- /dev/null +++ b/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 +} \ No newline at end of file diff --git a/extra/school-step.R b/extra/school-step.R new file mode 100644 index 0000000..196b9eb --- /dev/null +++ b/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)