Merge pull request #120 from sfcheung/devel

0.2.9.1: Can print R-squared increase

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: stdmod
 Title: Standardized Moderation Effect and Its Confidence Interval
-Version: 0.2.9
+Version: 0.2.9.1
 Authors@R:
     c(person(given = "Shu Fai",
            family = "Cheung",

NAMESPACE

@@ -22,3 +22,4 @@ export(stdmod)
 export(stdmod_boot)
 export(stdmod_lavaan)
 importFrom(rlang,.data)
+importFrom(stats,anova)

NEWS.md

@@ -1,3 +1,18 @@
+# stdmod 0.2.9.1
+
+- Revised `update.std_selected()`. Though
+  still not recommended, it should now
+  work more reliably if it needs to be
+  called. (0.2.9.1)
+- Improved the printout of the `summary()` of `std_selected()`
+  and `std_selected_boot()` outputs.
+  It now prints the R-squared increase
+  of the highest order term, as well as
+  the *F* test for the increase, if the
+  model has one and only one highest
+  order term (e.g., an interaction
+  term). (0.2.9.1)
+
 # stdmod 0.2.9
 
 - Fixed the issue with `stdmod-package`. (0.2.8.9001)

R/print_summary_stdmod.R

@@ -171,6 +171,19 @@ print.summary.std_selected <- function(x, ...,
                        p_digits = ceiling(-log10(pvalue_less_than)))
       cat("\n")
     }
+  if (!is.na(x$highest_order) && !identical(x$f_highest, NA)) {
+      rsq_highest <- formatC(x$f_highest[2, "R.sq.change"],
+                             digits = est_digits,
+                             format = "f")
+      cat("= Test the highest order term =",
+          paste0("The highest order term             : ", x$highest_order),
+          paste0("R-squared increase adding this term: ", rsq_highest),
+          sep = "\n")
+      print_fstatistic_change(x$f_highest,
+                              f_digits = t_digits,
+                              p_digits = ceiling(-log10(pvalue_less_than)))
+      cat("\n")
+    }
   tmp <- character(0)
   if (scaled_or_centered) {
       tmp1 <- paste("- Estimates and their statistics are based on the data after",
@@ -241,12 +254,36 @@ print_fstatistic <- function(fstatistic,
      p <- stats::pf(f, df1, df2, lower.tail = FALSE)
      p_txt <- format_pvalue(p,
                             eps = 10^(-p_digits))
+     if (!grepl("^<", p_txt)) {
+        p_txt <- paste0("= ", p_txt)
+       }
      cat("ANOVA test of R-squared  : ",
          f_txt, ", p ", p_txt, "\n", sep = "")
   }
 
 #' @noRd
 
+print_fstatistic_change <- function(fstatistic,
+                                    f_digits = 4,
+                                    p_digits = 3) {
+     f <- fstatistic[2, "F"]
+     df1 <- fstatistic[2, "Df"]
+     df2 <- fstatistic[2, "Res.Df"]
+     f_txt <- paste0("F(",
+                     df1, ", ", df2, ") = ",
+                     round(f, f_digits))
+     p <- fstatistic[2, "Pr(>F)"]
+     p_txt <- format_pvalue(p,
+                            eps = 10^(-p_digits))
+     if (!grepl("^<", p_txt)) {
+        p_txt <- paste0("= ", p_txt)
+       }
+     cat("F test of R-squared increase       : ",
+         f_txt, ", p ", p_txt, "\n", sep = "")
+  }
+
+#' @noRd
+
 format_pvalue <- function(p,
                           eps = 1e-3) {
     p_digits <- ceiling(-log10(eps))

R/summary_stdmod.R

@@ -42,6 +42,7 @@
 #' summary(lm_std_boot)
 #'
 #' @export
+#' @importFrom stats anova
 
 summary.std_selected <- function(object, ...) {
     out <- stats::summary.lm(object, ...)
@@ -57,6 +58,43 @@ summary.std_selected <- function(object, ...) {
                                 object$boot_ci,
                                 out$coefficients[, -1])
       }
+    out$highest_order <- tryCatch(highest_order(object),
+                                  error = function(e) NA)
+    if (!is.na(out$highest_order)) {
+        lm_out <- eval(object$lm_out_call,
+                       envir = parent.frame())
+        lm_call0 <- stats::update(lm_out,
+                                  paste0("~ .-", out$highest_order),
+                                  evaluate = FALSE)
+        lm_out0 <- eval(lm_call0,
+                        envir = parent.frame())
+        names(lm_out0)
+        anova_out <- anova(lm_out0, lm_out)
+        rsq_change <- summary(lm_out)$r.squared - summary(lm_out0)$r.squared
+        anova_out1 <- cbind(R.sq.change = c(NA, rsq_change), anova_out)
+        class(anova_out1) <- class(anova_out)
+        attr(anova_out1, "heading") <- attr(anova_out, "heading")
+        out$f_highest <- anova_out1
+      } else {
+        out$f_highest <- NA
+      }
     class(out) <- c("summary.std_selected", class(out))
     out
   }
+
+#' @noRd
+# Adapted from lmhelprs
+
+highest_order <- function(lm_out) {
+    terms_x <- stats::terms(lm_out)
+    labels_x <- labels(terms_x)
+    order_x <- attr(terms_x, "order")
+    order_max <- which.max(order_x)
+    order_min <- which.min(order_x)
+    max_n <- sum(order_x == max(order_x))
+    if ((order_max == order_min) ||
+        (max_n != 1)) {
+        stop("No unique highest order term.")
+      }
+    labels_x[order_max]
+  }

R/update_stdmod.R

@@ -43,9 +43,10 @@
 
 update.std_selected <- function(object, formula., ..., evaluate = TRUE) {
     # Adapted from update.default in stats
-    if (is.null(call <- stats::getCall(object))) {
-        stop("need an object with call component")
-      }
+    # if (is.null(call <- stats::getCall(object))) {
+    #     stop("need an object with call component")
+    #   }
+    call <- object$lm_out_call
     extras <- match.call(expand.dots = FALSE)$...
     if (!missing(formula.)) {
         call$formula <- stats::update.formula(stats::formula(object), formula.)
@@ -63,20 +64,18 @@ update.std_selected <- function(object, formula., ..., evaluate = TRUE) {
     if (is.null(extras$data)) {
         lm_out_call$data <- object$lm_out_call$data
       }
+    if (!is.null(object$std_selected_boot_call)) {
+        new_call <- object$std_selected_boot_call
+        new_call$lm_out <- lm_out_call
+      } else {
+        new_call <- object$std_selected_call
+        new_call$lm_out <- lm_out_call
+      }
     if (evaluate) {
         # lm_out <- eval(lm_out_call, parent.frame())
-        if (!is.null(object$std_selected_boot_call)) {
-            new_call <- object$std_selected_boot_call
-            new_call$lm_out <- lm_out_call
-            out <- eval(new_call, parent.frame())
-            return(out)
-          } else {
-            new_call <- object$std_selected_call
-            new_call$lm_out <- lm_out_call
-            out <- eval(new_call, parent.frame())
-            return(out)
-          }
+        out <- eval(new_call, parent.frame())
+        return(out)
       } else {
-        return(call)
+        return(new_call)
       }
   }

README.md

@@ -13,7 +13,7 @@
 
 # stdmod: Standardized Moderation <img src="man/figures/logo.png" align="right" height="150" />
 
-(Version 0.2.9, updated on 2023-09-09, [release history](https://sfcheung.github.io/stdmod/news/index.html))
+(Version 0.2.9.1, updated on 2023-11-14, [release history](https://sfcheung.github.io/stdmod/news/index.html))
 
 (Important changes since 0.2.0.0: Bootstrap confidence intervals and
 variance-covariance matrix of estimates are the defaults of `confint()`

tests/testthat/test_stdmod_selected_print_summary_rsq.R

@@ -0,0 +1,51 @@
+library(testthat)
+library(stdmod)
+
+set.seed(64252)
+dat <- test_x_1_w_1_v_1_cat1_n_500
+
+lm0_out <- lm(dv ~ iv + mod + v1 + cat1, dat)
+lm1_out <- lm(dv ~ iv * mod + v1 + cat1, dat)
+lm2_out <- lm(dv ~ iv * cat1 + mod + v1, dat)
+lm3_out <- lm(dv ~ mod * cat1 + iv + v1, dat)
+lm4_out <- lm(dv ~ mod * cat1 + iv*v1, dat)
+lm5_out <- lm(dv ~ v1, dat)
+lm0_std <- std_selected(lm0_out, to_standardize = ~ .)
+lm1_std <- std_selected(lm1_out, to_standardize = ~ .)
+lm2_std <- std_selected(lm2_out, to_standardize = ~ .)
+lm3_std <- std_selected(lm3_out, to_standardize = ~ .)
+lm4_std <- std_selected(lm4_out, to_standardize = ~ .)
+lm5_std <- std_selected(lm5_out, to_standardize = ~ .)
+lm0_std_boot <- std_selected_boot(lm0_out, to_standardize = ~ .)
+lm1_std_boot <- std_selected_boot(lm1_out, to_standardize = ~ .)
+lm2_std_boot <- std_selected_boot(lm2_out, to_standardize = ~ .)
+lm3_std_boot <- std_selected_boot(lm3_out, to_standardize = ~ .)
+lm4_std_boot <- std_selected_boot(lm4_out, to_standardize = ~ .)
+lm5_std_boot <- std_selected_boot(lm5_out, to_standardize = ~ .)
+
+test_that("summary.std_select with Rsq", {
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm0_std))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm1_std))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm2_std))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm3_std))))))
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm4_std))))))
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm5_std))))))
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm0_std_boot))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm1_std_boot))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm2_std_boot))))))
+    expect_true(any(grepl("adding this term",
+                    capture.output(print(summary(lm3_std_boot))))))
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm4_std_boot))))))
+    expect_false(any(grepl("adding this term",
+                     capture.output(print(summary(lm5_std_boot))))))
+  })