R2 incorrect for methods that use update(model, ~1) and have missing data

Fixes #803

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.13.0.1
+Version: 0.13.0.2
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

NEWS.md

@@ -8,6 +8,14 @@
   the full model), or can return singularity checks for each random effects term
   separately.
 
+## Bug fixes
+
+* Fixed issue with wrong computation of pseudo-R2 for some models where the
+  base-model (null model) was updated using the original data, which could
+  include missing values. Now the model frame is used, ensuring the correct
+  number of observations in the returned base-model, thus calculating the correct
+  log-likelihood and returning the correct pseudo-R2.
+
 # performance 0.13.0
 
 ## Breaking changes

R/r2_coxsnell.R

@@ -218,25 +218,45 @@ r2_coxsnell.svycoxph <- function(model, ...) {
 
 #' @export
 r2_coxsnell.multinom <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, ~1, trace = FALSE))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    trace = FALSE,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_coxsnell(model, l_base)
 }
 
 #' @export
 r2_coxsnell.clm2 <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, location = ~1, scale = ~1))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    location = ~1,
+    scale = ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_coxsnell(model, l_base)
 }
 
 #' @export
 r2_coxsnell.bayesx <- function(model, ...) {
-  junk <- utils::capture.output(l_base <- insight::get_loglikelihood(stats::update(model, ~1))) # nolint
+  junk <- utils::capture.output(
+    l_base <- insight::get_loglikelihood(stats::update(
+      model,
+      ~1,
+      data = insight::get_data(model, source = "mf")
+    ))
+  )
   .r2_coxsnell(model, l_base)
 }
 
 #' @export
 r2_coxsnell.clm <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, ~1))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   # if no loglik, return NA
   if (length(as.numeric(l_base)) == 0) {
     return(NULL)

R/r2_mcfadden.R

@@ -73,7 +73,11 @@ r2_mcfadden.glm <- function(model, verbose = TRUE, ...) {
     return(NULL)
   }
 
-  l_null <- insight::get_loglikelihood(stats::update(model, ~1))
+  l_null <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_mcfadden(model, l_null)
 }
 
@@ -162,21 +166,35 @@ r2_mcfadden.vglm <- function(model, ...) {
     insight::format_error("Can't get log-likelihood when `summ` is not zero.")
   }
 
-  l_null <- insight::get_loglikelihood(stats::update(model, ~1))
+  l_null <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_mcfadden(model, l_null)
 }
 
 
 #' @export
 r2_mcfadden.clm2 <- function(model, ...) {
-  l_null <- insight::get_loglikelihood(stats::update(model, location = ~1, scale = ~1))
+  l_null <- insight::get_loglikelihood(stats::update(
+    model,
+    location = ~1,
+    scale = ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_mcfadden(model, l_null)
 }
 
 
 #' @export
 r2_mcfadden.multinom <- function(model, ...) {
-  l_null <- insight::get_loglikelihood(stats::update(model, ~1, trace = FALSE))
+  l_null <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    trace = FALSE,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_mcfadden(model, l_null)
 }
 

R/r2_nagelkerke.R

@@ -170,19 +170,33 @@ r2_nagelkerke.negbinmfx <- r2_nagelkerke.logitmfx
 
 #' @export
 r2_nagelkerke.multinom <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, ~1, trace = FALSE))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    trace = FALSE,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_nagelkerke(model, l_base)
 }
 
 #' @export
 r2_nagelkerke.clm2 <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, location = ~1, scale = ~1))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    location = ~1,
+    scale = ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   .r2_nagelkerke(model, l_base)
 }
 
 #' @export
 r2_nagelkerke.clm <- function(model, ...) {
-  l_base <- insight::get_loglikelihood(stats::update(model, ~1))
+  l_base <- insight::get_loglikelihood(stats::update(
+    model,
+    ~1,
+    data = insight::get_data(model, source = "mf")
+  ))
   # if no loglik, return NA
   if (length(as.numeric(l_base)) == 0) {
     return(NULL)

tests/testthat/test-r2_nagelkerke.R

@@ -17,3 +17,36 @@ test_that("r2_nagelkerke", {
     }
   )
 })
+
+test_that("r2_nagelkerke, multinom, correct base-model with NA", {
+  skip_on_cran()
+  skip_if_not_installed("nnet")
+
+  n_obs <- 1000
+  softmax <- function(x){
+    exp(x - max(x)) / sum(exp(x - max(x)))
+  }
+  sample_y <- function(x){
+    sample(1:3, size = 1, prob = softmax(c(0.25 * x, -0.1 * x, 0 * x)))
+  }
+  set.seed(123)
+  sim_df <- data.frame(x = rnorm(n_obs, 0, 1), y = NA)
+
+  for(i in 1:nrow(sim_df)){
+    sim_df$y[i] <- sample_y(sim_df$x[i])
+  }
+
+  sim_df$x[1:500] <- NA
+  sim_df2 <- sim_df[!is.na(sim_df$x),]
+
+  m1 <- nnet::multinom(y ~ x, data = sim_df, trace = FALSE)
+  m2 <- nnet::multinom(y ~ x, data = sim_df2, trace = FALSE)
+
+  out1 <- r2_nagelkerke(m1)
+  out2 <- r2_nagelkerke(m2)
+  expect_equal(out1, out2, tolerance = 1e-4, ignore_attr = TRUE)
+
+  out1 <- r2_mcfadden(m1)
+  out2 <- r2_mcfadden(m2)
+  expect_equal(out1$R2, out2$R2, tolerance = 1e-4, ignore_attr = TRUE)
+})