changepoint icon indicating copy to clipboard operation
changepoint copied to clipboard

logLik() should return an object of class "logLik"

Open beanumber opened this issue 10 months ago • 5 comments

Thanks for the great work on this package!

The behavior of changepoint::logLik.cpt() is problematic for three reasons:

  1. It returns a numeric vector of length 2 (instead of the expected 1)
  2. It returns an object of class double instead of an object of class logLik, and thus the resulting object doesn't have the attributes that logLik objects should have.
  3. It doesn't return the actual log-likelihood (but rather -2 times the log-likelihood)!

This means that other generic functions already defined in stats like AIC() and BIC() don't work as expected.

library(changepoint)
#> Loading required package: zoo
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
#> Successfully loaded changepoint package version 2.2.4
#>  See NEWS for details of changes.
x <- cpt.meanvar(wave.c44137, penalty = "AIC")

# current behavior
logLik(x)
#>      -2*logLik -2*Loglike+pen 
#>       215528.5       215532.5
str(logLik(x))
#>  Named num [1:2] 215529 215533
#>  - attr(*, "names")= chr [1:2] "-2*logLik" "-2*Loglike+pen"
AIC(x)
#> Error in UseMethod("logLik"): no applicable method for 'logLik' applied to an object of class "cpt"
BIC(x)
#> Error in UseMethod("logLik"): no applicable method for 'logLik' applied to an object of class "cpt"

Created on 2024-04-03 with reprex v2.1.0

I think it would be better if changepoint::logLik.cpt() returned a logLik object with the appropriate attributes and values. Something like this should do the trick:

x <- changepoint::cpt.meanvar(changepoint::wave.c44137, penalty = "AIC")

logLik.cpt <- function(object, ...) {
  y <- changepoint::likelihood(object) |>
    suppressWarnings()
  ll <- -y[1] / 2
  attr(ll, "df") <- length(object@cpts)
  attr(ll, "nobs") <- length([email protected])
  class(ll) <- "logLik"
  return(ll)
}

# preferred behavior
logLik(x)
#> 'log Lik.' -107764.3 (df=2)
str(logLik(x))
#> Class 'logLik' : -107764 (df=2)
attributes(logLik(x))
#> $names
#> [1] "-2*logLik"
#> 
#> $df
#> [1] 2
#> 
#> $nobs
#> [1] 63651
#> 
#> $class
#> [1] "logLik"
AIC(x)
#> [1] 215532.5
BIC(x)
#> [1] 215550.7

Created on 2024-04-03 with reprex v2.1.0

beanumber avatar Apr 03 '24 21:04 beanumber