R/anova.R

In jackmwolf/pcsstools: Tools for Regression Using Pre-Computed Summary Statistics

# anova methods for pcsslm and pcsslmlist

# The original source for much of this came from R Core Team and contributors
# worldwide's `stats`` package.
#
# Below is the original license statement for R, which liscenses the stats
# package.
#
# This software is distributed under the terms of the GNU General
# Public License, either Version 2, June 1991 or Version 3, June 2007.
# The terms of version 2 of the license are in a file called COPYING
# which you should have received with
# this software and which can be displayed by RShowDoc("COPYING").
# Version 3 of the license can be displayed by RShowDoc("GPL-3").
#
# Copies of both versions 2 and 3 of the license can be found
# at https://www.R-project.org/Licenses/.
#
# A small number of files (the API header files listed in
#                          R_DOC_DIR/COPYRIGHTS) are distributed under the
# LESSER GNU GENERAL PUBLIC LICENSE, version 2.1 or later.
# This can be displayed by RShowDoc("LGPL-2.1"),
# or obtained at the URI given.
# Version 3 of the license can be displayed by RShowDoc("LGPL-3").
#
# 'Share and Enjoy.'
#
# GNU LESSER GENERAL PUBLIC LICENSE
# Version 3, 29 June 2007
#
# Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
#   Everyone is permitted to copy and distribute verbatim copies
# of this license document, but changing it is not allowed.
#
#
# This version of the GNU Lesser General Public License incorporates
# the terms and conditions of version 3 of the GNU General Public
# License, supplemented by the additional permissions listed below.
#
# 0. Additional Definitions.
#
# As used herein, "this License" refers to version 3 of the GNU Lesser
# General Public License, and the "GNU GPL" refers to version 3 of the GNU
# General Public License.
#
# "The Library" refers to a covered work governed by this License,
# other than an Application or a Combined Work as defined below.
#
# An "Application" is any work that makes use of an interface provided
# by the Library, but which is not otherwise based on the Library.
# Defining a subclass of a class defined by the Library is deemed a mode
# of using an interface provided by the Library.
#
# A "Combined Work" is a work produced by combining or linking an
# Application with the Library.  The particular version of the Library
# with which the Combined Work was made is also called the "Linked
# Version".
#
# The "Minimal Corresponding Source" for a Combined Work means the
# Corresponding Source for the Combined Work, excluding any source code
# for portions of the Combined Work that, considered in isolation, are
# based on the Application, and not on the Linked Version.
#
# The "Corresponding Application Code" for a Combined Work means the
# object code and/or source code for the Application, including any data
# and utility programs needed for reproducing the Combined Work from the
# Application, but excluding the System Libraries of the Combined Work.
#
# 1. Exception to Section 3 of the GNU GPL.
#
# You may convey a covered work under sections 3 and 4 of this License
# without being bound by section 3 of the GNU GPL.
#
# 2. Conveying Modified Versions.
#
# If you modify a copy of the Library, and, in your modifications, a
# facility refers to a function or data to be supplied by an Application
# that uses the facility (other than as an argument passed when the
#                         facility is invoked), then you may convey a copy of the modified
# version:
#
#   a) under this License, provided that you make a good faith effort to
# ensure that, in the event an Application does not supply the
# function or data, the facility still operates, and performs
# whatever part of its purpose remains meaningful, or
#
# b) under the GNU GPL, with none of the additional permissions of
# this License applicable to that copy.
#
# 3. Object Code Incorporating Material from Library Header Files.
#
# The object code form of an Application may incorporate material from
# a header file that is part of the Library.  You may convey such object
# code under terms of your choice, provided that, if the incorporated
# material is not limited to numerical parameters, data structure
# layouts and accessors, or small macros, inline functions and templates
# (ten or fewer lines in length), you do both of the following:
#
#   a) Give prominent notice with each copy of the object code that the
# Library is used in it and that the Library and its use are
# covered by this License.
#
# b) Accompany the object code with a copy of the GNU GPL and this license
# document.
#
# 4. Combined Works.
#
# You may convey a Combined Work under terms of your choice that,
# taken together, effectively do not restrict modification of the
# portions of the Library contained in the Combined Work and reverse
# engineering for debugging such modifications, if you also do each of
# the following:
#
#   a) Give prominent notice with each copy of the Combined Work that
# the Library is used in it and that the Library and its use are
# covered by this License.
#
# b) Accompany the Combined Work with a copy of the GNU GPL and this license
# document.
#
# c) For a Combined Work that displays copyright notices during
# execution, include the copyright notice for the Library among
# these notices, as well as a reference directing the user to the
# copies of the GNU GPL and this license document.
#
# d) Do one of the following:
#
#   0) Convey the Minimal Corresponding Source under the terms of this
# License, and the Corresponding Application Code in a form
# suitable for, and under terms that permit, the user to
# recombine or relink the Application with a modified version of
# the Linked Version to produce a modified Combined Work, in the
# manner specified by section 6 of the GNU GPL for conveying
# Corresponding Source.
#
# 1) Use a suitable shared library mechanism for linking with the
# Library.  A suitable mechanism is one that (a) uses at run time
# a copy of the Library already present on the user's computer
#        system, and (b) will operate properly with a modified version
#        of the Library that is interface-compatible with the Linked
#        Version.
#
#    e) Provide Installation Information, but only if you would otherwise
#    be required to provide such information under section 6 of the
#    GNU GPL, and only to the extent that such information is
#    necessary to install and execute a modified version of the
#    Combined Work produced by recombining or relinking the
#    Application with a modified version of the Linked Version. (If
#    you use option 4d0, the Installation Information must accompany
#    the Minimal Corresponding Source and Corresponding Application
#    Code. If you use option 4d1, you must provide the Installation
#    Information in the manner specified by section 6 of the GNU GPL
#    for conveying Corresponding Source.)
#
#   5. Combined Libraries.
#
#   You may place library facilities that are a work based on the
# Library side by side in a single library together with other library
# facilities that are not Applications and are not covered by this
# License, and convey such a combined library under terms of your
# choice, if you do both of the following:
#
#    a) Accompany the combined library with a copy of the same work based
#    on the Library, uncombined with any other library facilities,
#    conveyed under the terms of this License.
#
#    b) Give prominent notice with the combined library that part of it
#    is a work based on the Library, and explaining where to find the
#    accompanying uncombined form of the same work.
#
#   6. Revised Versions of the GNU Lesser General Public License.
#
#   The Free Software Foundation may publish revised and/or new versions
# of the GNU Lesser General Public License from time to time. Such new
# versions will be similar in spirit to the present version, but may
# differ in detail to address new problems or concerns.
#
#   Each version is given a distinguishing version number. If the
# Library as you received it specifies that a certain numbered version
# of the GNU Lesser General Public License "or any later version"
# applies to it, you have the option of following the terms and
# conditions either of that published version or of any later version
# published by the Free Software Foundation. If the Library as you
# received it does not specify a version number of the GNU Lesser
# General Public License, you may choose any version of the GNU Lesser
# General Public License ever published by the Free Software Foundation.
#
#   If the Library as you received it specifies that a proxy can decide
# whether future versions of the GNU Lesser General Public License shall
# apply, that proxy's public statement of acceptance of any version is
# permanent authorization for you to choose that version for the
# Library.


#' ANOVA for linear models fit using PCSS
#' 
#' Compute an analysis of variance table for one or more linear model fitted
#' using PCSS.
#' 
#' @param object,... objects of class \code{pcsslm}.
#' @param test a character string specifying the test statistic to be used. Can 
#'   be one of \code{"F"}, \code{"Chisq"} or \code{"Cp"}, with partial matching 
#'   allowed, or \code{NULL} for no test.
#' @param scale numeric. An estimate of the noise variance \eqn{\sigma^2}. If 
#'   zero this will be estimated from the largest model considered.
#' 
#' @author R Core Team and contributors worldwide. Modified by Jack Wolf
#' @method anova pcsslm
#' 
#' @importFrom stats formula stat.anova
#' 
#' @return An object of class \code{"anova"} inheriting from class 
#'   \code{"data.frame"}.
#' 
#' @export
anova.pcsslm <- function(object, ...) {
  if (length(list(object, ...)) > 1L) {
    return(anova.pcsslmlist(object, ...))
  }

  ssr <- object$`Sum Sq`["SSE"]
  mss <- NA

  # if (ssr < 1e-10 * mss) {
  #   warning("ANOVA F-tests on an essentially perfect fit are unreliable")
  # }

  dfr <- object$df[2]
  p <- object$df[1]
  if (p > 0L) {
    ss <- c(object$`Sum Sq`["SSR"], ssr)
    df <- c(p - 1, dfr)

    tlabels <- c(attr(object$terms, "term.labels"))
    tlabels <- tlabels[tlabels %in% names(which(!object$aliased))]
    tlabels <- paste(tlabels, collapse = ", ")
  }
  else {
    ss <- ssr
    df <- dfr
    tlabels <- character()
  }
  ms <- ss / df
  f <- ms / (ssr / dfr)
  P <- pf(f, df, dfr, lower.tail = FALSE)
  table <- data.frame(df, ss, ms, f, P)
  table[length(P), 4:5] <- NA
  dimnames(table) <- list(c(tlabels, "Residuals"), c(
    "Df",
    "Sum Sq", "Mean Sq", "F value", "Pr(>F)"
  ))

  structure(table,
    heading = c(
      "Analysis of Variance Table\n",
      paste("Response:", deparse(formula(object)[[2L]]))
    ),
    class = c("anova", "data.frame")
  )
}

#' @rdname anova.pcsslm 
anova.pcsslmlist <- function(object, ..., scale = 0, test = "F") {
  objects <- list(object, ...)
  responses <- as.character(lapply(objects, function(x) deparse(x$terms[[2L]])))
  sameresp <- responses == responses[1L]
  if (!all(sameresp)) {
    objects <- objects[sameresp]
    warning(gettextf(
      "models with response %s removed because response differs from model 1",
      sQuote(deparse(responses[!sameresp]))
    ), domain = NA)
  }

  nmodels <- length(objects)
  if (nmodels == 1) {
    return(anova.pcsslm(object))
  }
  resdf <- as.numeric(lapply(objects, function(x) x$df[2]))
  resdev <- as.numeric(lapply(objects, function(x) x$`Sum Sq`["SSE"]))
  table <- data.frame(resdf, resdev, c(NA, -diff(resdf)), c(
    NA,
    -diff(resdev)
  ))
  variables <- lapply(objects, function(x) {
    paste(deparse(formula(x)),
      collapse = "\n"
    )
  })
  dimnames(table) <- list(1L:nmodels, c(
    "Res.Df", "RSS",
    "Df", "Sum of Sq"
  ))
  title <- "Analysis of Variance Table\n"
  topnote <- paste0("Model ", format(1L:nmodels), ": ",
    variables,
    collapse = "\n"
  )
  if (!is.null(test)) {
    bigmodel <- order(resdf)[1L]
    scale <- if (scale > 0) {
      scale
    } else {
      resdev[bigmodel] / resdf[bigmodel]
    }
    table <- stat.anova(
      table = table, test = test, scale = scale,
      df.scale = resdf[bigmodel], n = sum(objects[[1]]$df[1:2])
    )
  }
  structure(table, heading = c(title, topnote), class = c(
    "anova",
    "data.frame"
  ))
}