fix issue #567 (ergm), fix issue #572 (svyglm)

NAMESPACE

@@ -89,6 +89,7 @@ S3method(glance,survdiff)
 S3method(glance,survexp)
 S3method(glance,survfit)
 S3method(glance,survreg)
+S3method(glance,svyglm)
 S3method(glance,svyolr)
 S3method(glance,tbl_df)
 S3method(tidy,"NULL")
@@ -190,6 +191,7 @@ S3method(tidy,survdiff)
 S3method(tidy,survexp)
 S3method(tidy,survfit)
 S3method(tidy,survreg)
+S3method(tidy,svyglm)
 S3method(tidy,svyolr)
 S3method(tidy,table)
 S3method(tidy,ts)

R/ergm-tidiers.R

@@ -9,7 +9,7 @@
 #' @template param_confint
 #' @template param_exponentiate
 #' @template param_quick
-#' @param ... Additional arguments to pass to [ergm::summary.ergm()].
+#' @param ... Additional arguments to pass to [ergm::summary()].
 #'   **Cautionary note**: Mispecified arguments may be silently ignored.
 #' 
 #' @evalRd return_tidy(
@@ -53,7 +53,7 @@
 #' @export 
 #' @aliases ergm_tidiers
 #' @seealso [tidy()], [ergm::ergm()], [ergm::control.ergm()], 
-#'   [ergm::summary.ergm()]
+#'   [ergm::summary()]
 #' @family ergm tidiers
 tidy.ergm <- function(x, conf.int = FALSE, conf.level = .95,
                       exponentiate = FALSE, quick = FALSE, ...) {
@@ -62,7 +62,7 @@ tidy.ergm <- function(x, conf.int = FALSE, conf.level = .95,
     ret <- tibble(term = names(co), estimate = unname(co))
     return(process_ergm(ret, conf.int = FALSE, exponentiate = exponentiate))
   }
-  co <- ergm::summary.ergm(x, ...)$coefs
+  co <- ergm:::summary.ergm(x, ...)$coefs
 
   nn <- c("estimate", "std.error", "mcmc.error", "p.value")
   ret <- fix_data_frame(co, nn[1:ncol(co)])
@@ -103,26 +103,28 @@ tidy.ergm <- function(x, conf.int = FALSE, conf.level = .95,
 #' @family ergm tidiers
 glance.ergm <- function(x, deviance = FALSE, mcmc = FALSE, ...) {
   # will show appropriate warnings about standard errors, pseudolikelihood etc.
-  s <- ergm::summary.ergm(x, ...)
+  s <- ergm:::summary.ergm(x, ...)
   # dyadic (in)dependence and number of MCMLE iterations
   ret <- tibble(independence = s$independence, iterations = x$iterations)
   # log-likelihood
-  ret$logLik <- tryCatch(as.numeric(ergm::logLik.ergm(x)), error = function(e) NULL)
+  ret$logLik <- tryCatch(as.numeric(ergm:::logLik.ergm(x)), error = function(e) NULL)
   # null and residual deviance
   if (deviance & !is.null(ret$logLik)) {
-    dyads <- ergm::get.miss.dyads(x$constrained, x$constrained.obs)
-    dyads <- statnet.common::NVL(dyads, network::network.initialize(1))
-    dyads <- network::network.edgecount(dyads)
-    dyads <- network::network.dyadcount(x$network, FALSE) - dyads
+    # thanks to Pavel N. Krivitsky (@krivit)
+    # https://github.com/tidymodels/broom/issues/567
+    if (packageVersion("ergm") < "3.10") {
+      dyads <- sum(as.rlebdm(x$constrained, x$constrained.obs, which = "informative"))
+    } else {
+      dyads <- stats::nobs(x)
+    }
 
-    ret$null.deviance <- ergm::logLikNull(x)
+    ret$null.deviance <- ergm:::logLikNull(x)
     ret$null.deviance <- ifelse(is.na(ret$null.deviance), 0, -2 * ret$null.deviance)
     ret$df.null <- dyads
 
     ret$residual.deviance <- -2 * ret$logLik
     ret$df.residual <- dyads - length(x$coef)
   }
-
   ret$AIC <- tryCatch(stats::AIC(x), error = function(e) NULL)
   ret$BIC <- tryCatch(stats::BIC(x), error = function(e) NULL)
 

R/survey-tidiers.R

@@ -1,3 +1,116 @@
+#' @templateVar class svyglm
+#' @template title_desc_tidy_lm_wrapper
+#'
+#' @param x A `svyglm` object returned from [survey::svyglm()].
+#'
+#' @export
+#' @family lm tidiers
+#' @seealso [survey::svyglm()], [stats::glm()]
+# [NOTES]
+# - partly copies tidy.glm, which itself copies tidy.lm
+# - depends on process_lm, but that might change in the future
+tidy.svyglm <- function(x, conf.int = FALSE, conf.level = .95,
+                        exponentiate = FALSE, ...) {
+
+  s <- survey:::summary.svyglm(x)
+  ret <- tidy.summary.lm(s)
+
+  process_lm(ret, x,
+             conf.int = conf.int, conf.level = conf.level,
+             exponentiate = exponentiate
+  )
+
+}
+
+#' @templateVar class svyglm
+#' @template title_desc_glance
+#'
+#' @param x A `svyglm` object returned from [survey::svyglm()].
+#' @param maximal A `svyglm` object corresponding to the maximal model against
+#' which to compute the BIC: see Lumley and Scott (2015) for details. Defaults
+#' to \code{x}, which is equivalent to not using a maximal model.
+#' @template param_unused_dots
+#'
+#' @evalRd return_glance(
+#'   "null.deviance",
+#'   "df.null",
+#'   "AIC",
+#'   "BIC,
+#'   "deviance",
+#'   "df.residual"
+#' )
+#'
+#' @examples
+#'
+#' library(survey)
+#' 
+#' set.seed(123)
+#' data(api)
+#' 
+#' # survey design
+#' dstrat <-
+#'   svydesign(
+#'     id =  ~ 1,
+#'     strata =  ~ stype,
+#'     weights =  ~ pw,
+#'     data = apistrat,
+#'     fpc =  ~ fpc
+#'   )
+#' 
+#' # model
+#' m <- survey::svyglm(
+#'   formula = sch.wide ~ ell + meals + mobility,
+#'   design = dstrat,
+#'   family = quasibinomial()
+#' )
+#' 
+#' glance(m)
+#'
+#' @references Lumley T, Scott A (2015). AIC and BIC for modelling with complex 
+#' survey data. *Journal of Survey Statistics and Methodology*, 3(1).
+#' <https://doi.org/10.1093/jssam/smu021>.
+#' 
+#' @export
+#' @family lm tidiers
+#' @seealso [survey::svyglm()], [stats::glm()], [survey::anova.svyglm]
+glance.svyglm <- function(x, maximal = x, ...) {
+
+  ret <- with(x, tibble::tibble(null.deviance, df.null))
+
+  # log-likelihood does not apply (returns deviance instead)
+  #
+  # > ret$logLik <- tryCatch(as.numeric(stats::logLik(x)), error = function(e) NULL)
+  # Warning in logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
+
+  # AIC
+  #
+  #   equivalent to stats::AIC(x)
+  #   not always directly computed by svyglm, e.g. if family = quasibinomial()
+  #
+  ret$AIC <- survey:::AIC.svyglm(x)["AIC"]
+
+  # BIC
+  #
+  #   equivalent to stats::BIC(x, maximal)
+  #
+  ret$BIC <- survey:::dBIC(x, maximal)["BIC"]
+
+  # deviance
+  #
+  #   equivalent to stats::deviance(x)
+  #
+  ret$deviance <- x$deviance
+
+  # df.residual
+  #
+  #   equivalent to stats::df.residual(x)
+  #
+  ret$df.residual <- x$df.residual
+
+  ret
+
+}
+
 #' @templateVar class svyolr
 #' @template title_desc_tidy
 #' 
@@ -44,4 +157,3 @@ glance.svyolr <- function(x, ...) {
     nobs = stats::nobs(x)
   )
 }
-