Merge pull request #580 from crsh/apa_print.VGAM

Add support for analysis of deviance from car

DESCRIPTION

@@ -46,13 +46,16 @@ Suggests:
   latex2exp,
   lme4,
   lmerTest,
+  MASS,
   MBESS,
   multcomp,
   nlme,
+  nnet,
   R.rsp,
   skimr,
   spelling,
-  testthat
+  testthat,
+  VGAM
 SystemRequirements: Rendering the document template requires
   pandoc (>= 2.0; https://pandoc.org) and for PDFs a TeX distribution,
   such as TinyTeX (>= 0.12; https://yihui.org/tinytex/)

R/apa_print_anova.R

@@ -1,9 +1,10 @@
-#' Typeset Statistical Results from ANOVA
+#' Typeset Statistical Results from Analysis of Variance (or Deviance)
 #'
-#' These methods take objects from various R functions that calculate ANOVA to
-#' create formatted character strings to report the results in accordance with
-#' APA manuscript guidelines. For `anova`-objects from model comparisons see
-#' \code{\link{apa_print.list}}.
+#' These methods take objects from various R functions that calculate analysis
+#' of variance (i.e., ANOVA) or analysis of deviance. They create formatted
+#' character strings to report the results in accordance with APA manuscript
+#' guidelines. For `anova`-objects from model comparisons see
+#' [apa_print.list()].
 #'
 #' @param x An object containing the results from an analysis of variance ANOVA
 #' @param correction Character. For repeated-measures ANOVA, the type of
@@ -415,6 +416,30 @@ apa_print.anova <- function(
         , simplify = TRUE
       )
     )
+  # stats::anova.glm() and car::Anova.glm
+  } else if(any(grepl("Deviance", object_heading))) {
+    x$Term <-rownames(x)
+    y <- canonize(x)
+    y <- remove_residuals_row(y)
+    if(all(colnames(x) != "F values")) y$df.residual <- NULL
+    if(any(colnames(x) == "Cp")) y$df <- NULL
+
+    if(is.null(y$statistic)) {
+      y$statistic <- y$deviance
+      variable_label(y$statistic) <- "$\\chi^2$"
+      y$deviance  <- NULL
+    }
+    y <- beautify(y, ...)
+    return(
+      glue_apa_results(
+        y
+        , est_glue = construct_glue(y, "estimate")
+        , stat_glue = construct_glue(y, "statistic")
+        , in_paren = in_paren
+        , simplify = TRUE
+      )
+    )
+
   } else if(any(grepl("Satterthwaite|Kenward", object_heading))) {
     # lmerTest::anova.merModLmerTest -------------------------------------------
 
@@ -481,7 +506,7 @@ apa_print.anova <- function(
   } else if(identical(object_heading[1], "ANOVA-like table for random-effects: Single term deletions")) {
     stop("Single-term deletions are not supported, yet.\nVisit https://github.com/crsh/papaja/issues to request support.")
   }
-  # anova::lm (single model) ----
+  # anova.lm() (single model) ----
   # Canonize, beautify, glue ----
   y <- as.data.frame(x, stringsAsFactors = FALSE)
   y$Effect <- trimws(rownames(y))

R/lookup_tables.R

@@ -23,6 +23,8 @@ lookup_names <- c(
   , "BIC"     = "BIC"
   , "npar"    = "n.parameters"
   , "alternative" = "alternative"
+  , "Deviance"    = "deviance"
+  , "Resid..Dev"  = "residual.deviance" # from broom
   # term
   , "Effect"     = "term"
   , "Term"       = "term"
@@ -81,7 +83,9 @@ lookup_names <- c(
   , "F.value"   = "statistic"
   , "F"         = "statistic"
   , "F.ratio"   = "statistic"
+  , "F.values"  = "statistic"
   , "LRT"       = "statistic"
+  , "LR.Chisq"  = "statistic"
   , "Chisq"     = "statistic"
   , "chisq"     = "statistic"
   , "X.squared" = "statistic"
@@ -99,6 +103,8 @@ lookup_names <- c(
   , "logbf01"      = "statistic"
   , "Bartlett.s.K.2"          = "statistic"
   , "Bartlett.s.K.squared"    = "statistic"
+  , "Rao" = "statistic"
+  , "Cp"  = "statistic"
   # df, df.residual
   , "multivariate.df1" = "multivariate.df"
   , "multivariate.df2" = "multivariate.df.residual"
@@ -121,9 +127,11 @@ lookup_names <- c(
   , "denom.df"   = "df.residual"
   , "den.df"     = "df.residual"
   , "Res.Df"     = "df.residual"
+  , "Resid..Df"  = "df.residual"
   # p.value
   , "p.value"    = "p.value"
   , "Pr..Chisq." = "p.value"
+  , "Pr..Chi."   = "p.value"
   , "Pr..F."     = "p.value"
   , "Pr..PB."    = "p.value"
   , "Pr...t.."   = "p.value"
@@ -148,6 +156,8 @@ lookup_labels <- c(
   , "BIC"     = "$\\mathit{BIC}$"
   , "npar"    = "$k$"
   , "alternative" = "$\\mathcal{H}_1$"
+  , "Deviance"   = "$\\mathit{Dev}$"
+  , "Resid..Dev" = "$\\mathit{Dev}_{\\mathrm{res}}$"
   # term
   , "Effect"     = "Effect"
   , "Term"       = "Term"
@@ -199,8 +209,10 @@ lookup_labels <- c(
   , "F.value"   = "$F$"
   , "F"         = "$F$"
   , "F.ratio"   = "$F$"
+  , "F.values"  = "$F$"
   , "approx.F"  = "$F$"
   , "LRT"       = "$\\chi^2$"
+  , "LR.Chisq"  = "$\\chi^2$"
   , "chisq"     = "$\\chi^2$"
   , "Chisq"     = "$\\chi^2$"
   , "X.squared" = "$\\chi^2$"
@@ -218,6 +230,8 @@ lookup_labels <- c(
   , "logbf01"      = "$\\log \\mathrm{BF}_{\\textrm{01}}$"
   , "Bartlett.s.K.2"          = "$K^2$"
   , "Bartlett.s.K.squared"    = "$K^2$"
+  , "Rao" = "$\\mathit{RS}$"
+  , "Cp"  = "$C_p$"
   # df, df.residual
   , "multivariate.df" = "$\\mathit{df}$"
   , "multivariate.df.residual" = "$\\mathit{df}_{\\mathrm{res}}$"
@@ -246,6 +260,7 @@ lookup_labels <- c(
   , "Pr...z.."    = "$p$"
   , "pvalues"     = "$p$"
   , "Pr..Chisq."  = "$p$"
+  , "Pr..Chi."    = "$p$"
   , "Pr..F."      = "$p$"
   , "Pr..PB."     = "$p$"
   , "adj.p.value" = "$p_\\mathrm{adj}$"

inst/NEWS.md

@@ -1,6 +1,7 @@
 # Upcoming release
 
 - For ANOVA methods, *MSE*s are again returned if requested by the user (reported by @Sashpta, [#562](https://github.com/crsh/papaja/issues/562)). The global default for reporting *MSE*s now depends on the [**effectsize**](https://CRAN.r-project.org/package=effectsize) package: If **effectsize** is installed, the default for reporting *MSE*s is `FALSE`, if **effectsize** is not installed, the default is `TRUE`.
+- Added `apa_print()` support for analysis of deviance from the **car** package.
 
 ### Existing functions
 

tests/testthat/test_apa_print_anova.R

@@ -511,5 +511,110 @@ test_that(
   }
 )
 
+test_that(
+  "Analysis of deviance from car"
+  , {
+    fit <- bwt.mu <- nnet::multinom(low ~ ., MASS::birthwt)
+    car_out <- papaja::apa_print(car::Anova(fit))
+    expect_apa_results(
+      car_out
+      , labels = list(
+        term          = "Term"
+        , statistic   = "$\\chi^2$"
+        , df          = "$\\mathit{df}$"
+        , p.value     = "$p$"
+      )
+    )
+    # Example 1: a proportional odds model fitted to pneumo.
+    set.seed(1)
+    pneumo <- transform(VGAM::pneumo, let = log(exposure.time), x3 = runif(8))
+    fit1 <- VGAM::vglm(cbind(normal, mild, severe) ~ let     , VGAM::propodds, pneumo)
+    fit2 <- VGAM::vglm(cbind(normal, mild, severe) ~ let + x3, VGAM::propodds, pneumo)
+    fit3 <- VGAM::vglm(cbind(normal, mild, severe) ~ let + x3, VGAM::cumulative, pneumo)
+    car_out <- apa_print(car::Anova(fit1, type = 3))
+    expect_apa_results(
+      car_out
+      , labels = list(
+        term          = "Term"
+        , statistic   = "$\\chi^2$"
+        , df          = "$\\mathit{df}$"
+        , p.value     = "$p$"
+      )
+    )
+  }
+)
+
+test_that(
+  "Analysis of deviance from the stats and car packages"
+  , {
+    # From stats::glm() examples section:
+    ## Dobson (1990) Page 93: Randomized Controlled Trial :
+    counts <- c(18,17,15,20,10,20,25,13,12)
+    outcome <- gl(3,1,9)
+    treatment <- gl(3,3)
+    data.frame(treatment, outcome, counts) # showing data
+    glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
+
+    chisq_out <- apa_print(anova(glm.D93, test = "Chisq"))
+    cp_out    <- apa_print(anova(glm.D93, test = "Cp"))
+    lrt_out   <- apa_print(anova(glm.D93, test = "LRT"))
+    rao_out   <- apa_print(anova(glm.D93, test = "Rao"))
+
+    expect_identical(
+      chisq_out$full_result
+      , list(
+        outcome = "$\\chi^2(2) = 5.45$, $p = .065$"
+        , treatment = "$\\chi^2(2) = 0.00$, $p > .999$"
+      )
+    )
+    expect_identical(
+      cp_out$full_result
+      , list(
+        outcome = "$C_p = 11.13$"
+        , treatment = "$C_p = 15.13$"
+      )
+    )
+    expect_identical(
+      lrt_out
+      , chisq_out
+    )
+    expect_identical(
+      rao_out$full_result
+      , list(
+        outcome = "$\\mathit{RS}(2) = 5.56$, $p = .062$"
+        , treatment = "$\\mathit{RS}(2) = 0.00$, $p > .999$"
+      )
+    )
+
+    car_lr_out   <- apa_print(car::Anova(glm.D93, type = 3, test.statistic = "LR"))
+    car_wald_out <- apa_print(car::Anova(glm.D93, type = 3, test.statistic = "Wald"))
+    car_f_out    <- apa_print(car::Anova(glm.D93, type = 3, test.statistic = "F"))
+
+    expect_identical(
+      car_lr_out$full_result
+      , list(
+        outcome = "$\\chi^2(2) = 5.45$, $p = .065$"
+        , treatment = "$\\chi^2(2) = 0.00$, $p > .999$"
+      )
+    )
+    expect_identical(
+      car_wald_out$full_result
+      , list(
+        Intercept = "$\\chi^2(1) = 317.37$, $p < .001$"
+        , outcome = "$\\chi^2(2) = 5.49$, $p = .064$"
+        , treatment = "$\\chi^2(2) = 0.00$, $p > .999$"
+      )
+    )
+    expect_identical(
+      car_f_out$full_result
+      , list(
+        outcome     = "$F(2, 4) = 2.11$, $p = .237$"
+        , treatment = "$F(2, 4) = 0.00$, $p > .999$"
+      )
+    )
+  }
+)
+
+
 # restore previous options
  options(op)