从 checkresiduals 函数中提取 p 值

Question

我正在使用forecast包进行预测。下面是我的预测结果。

#CODE 
library(forecast)
      DATA_SET<-data.frame(TEST=c(200,220,200,260,300,290,320,340,360,500,200,300,400,250,350,390,400,450,470,350,300,220,580,450,120,250,360,470)
                           )
      View(DATA_SET)

      # Making TS object
      TS_DATA_SET<-ts(DATA_SET,start=c(2010,1),frequency = 12)

      # Forecasting
      TS_FORECAST<-auto.arima(TS_DATA_SET)

现在我想从checkresiduals函数中提取p值，并将其存储到数据框中，

   #Checking residuals
   checkresiduals(TS_FORECAST, plot = FALSE)

##  Ljung-Box test
##
##   data:  Residuals from ARIMA(0,0,0) with non-zero mean
## Q* = 4.5113, df = 4.6, p-value = 0.4237
##
## Model df: 1.   Total lags used: 5.6

我正在尝试使用下面的代码，但是遇到了问题。

p-value<-data.frame(checkresiduals(TS_FORECAST, plot = FALSE))
p_value
#data frame with 0 columns and 0 rows

请问有人能帮我从checkresiduals函数中提取p值（p值=0.4237）到数据框中吗？

- silent_hunter

Answer 1

编辑:

我的第一种方法已经被实现在'checkresiduals()'函数中。现在该函数默认返回输出值。

---

旧答案:

不幸的是，函数checkresiduals()没有返回值，只是将它们用prints()打印出来。您可以通过写入不带括号的checkresiduals来查看该函数。或者您可以查看开发人员的github。

您可以通过在函数中放置一个return()来重写该函数。我只是复制并粘贴该函数，并将其插入到末尾：

 checkresiduals <- function(object, lag, df=NULL, test, plot=TRUE, ...) {
  showtest <- TRUE
  if (missing(test)) {
    if (is.element("lm", class(object))) {
      test <- "BG"
    } else {
      test <- "LB"
    }
    showtest <- TRUE
  }
  else if (test != FALSE) {
    test <- match.arg(test, c("LB", "BG"))
    showtest <- TRUE
  }
  else {
    showtest <- FALSE
  }

  # Extract residuals
  if (is.element("ts", class(object)) | is.element("numeric", class(object))) {
    residuals <- object
    object <- list(method = "Missing")
  }
  else {
    residuals <- residuals(object)
  }

  if (length(residuals) == 0L) {
    stop("No residuals found")
  }

  if ("ar" %in% class(object)) {
    method <- paste("AR(", object$order, ")", sep = "")
  } else if (!is.null(object$method)) {
    method <- object$method
  } else if ("HoltWinters" %in% class(object)) {
    method <- "HoltWinters"
  } else if ("StructTS" %in% class(object)) {
    method <- "StructTS"
  } else {
    method <- try(as.character(object), silent = TRUE)
    if ("try-error" %in% class(method)) {
      method <- "Missing"
    } else if (length(method) > 1 | base::nchar(method[1]) > 50) {
      method <- "Missing"
    }
  }
  if (method == "Missing") {
    main <- "Residuals"
  } else {
    main <- paste("Residuals from", method)
  }

  if (plot) {
    suppressWarnings(ggtsdisplay(residuals, plot.type = "histogram", main = main, ...))
  }

  # Check if we have the model
  if (is.element("forecast", class(object))) {
    object <- object$model
  }

  if (is.null(object) | !showtest) {
    return(invisible())
  }

  # Seasonality of data
  freq <- frequency(residuals)

  # Find model df
  if(grepl("STL \\+ ", method)){
    warning("The fitted degrees of freedom is based on the model used for the seasonally adjusted data.")
  }
  df <- modeldf(object)

  if (missing(lag)) {
    lag <- ifelse(freq > 1, 2 * freq, 10)
    lag <- min(lag, round(length(residuals)/5))
    lag <- max(df+3, lag)
  }

  if (!is.null(df)) {
    if (test == "BG") {
      # Do Breusch-Godfrey test
      BGtest <- lmtest::bgtest(object, order = lag)
      BGtest$data.name <- main
      print(BGtest)
      return(BGtest)
    }
    else {
      # Do Ljung-Box test
      LBtest <- Box.test(zoo::na.approx(residuals), fitdf = df, lag = lag, type = "Ljung")
      LBtest$method <- "Ljung-Box test"
      LBtest$data.name <- main
      names(LBtest$statistic) <- "Q*"
      print(LBtest)
      cat(paste("Model df: ", df, ".   Total lags used: ", lag, "\n\n", sep = ""))
      return(LBtest)
    }
  }

}

你还需要从github文件中获得modeldf()函数。

modeldf <- function(object, ...){
  UseMethod("modeldf")
}

modeldf.Arima <- function(object, ...){
  length(object$coef)
}

使用这个解决方案，您可以使用原始的checkresiduals函数。现在，您可以通过以下方式调用p.value：

res_values <- checkresiduals(TS_FORECAST, plot = TRUE)
res_values$p.value

您也可以只使用Ljung-Box 和 Breusch-Godfrey test，而不需要使用checkresiduals()函数，因为这就是checkresiduals()函数的功能。

我认为编辑checkresiduals()函数是一种更方便的方法，这样您就可以像以前一样使用它。您可以将其粘贴到代码中，它应该会起作用。只需确保在调用函数之前声明modeldf()和modeldf().Arima。此外，它可以测试函数所做的任何测试。

---

选项2 因为它是可能的:

您可以使用capture.output()捕获输出。

capture.output(checkresiduals(TS_FORECAST, plot = FALSE))[5]

"Q* = 4.8322，df = 5，p-value = 0.4367"

使用grep命令可以在不更改函数的情况下提取p值。由于我不熟悉grep，因此无法对此任务提供适当的答案。

- mischva11

---

为了完整起见：我将我的回答更改作为Github提交。它应该在下一个更新中更改。 - mischva11

这个更新（来自@mischva11）现在是“checkresiduals”的一部分，所以它确实会返回值。如果需要的话，请参考我的有关转换为“data.frame”的答案。 - Brian Stamper

@Brian Stamper谢谢您，我更新了我的回答，我也喜欢您整洁的方法。 - mischva11

Answer 2

这里可以看到checkresiduals()的内部实现。

不幸的是，根据文档，它不返回任何值，因此无法轻松提取所需内容。

但你可以进行相同的计算（查看链接存储库中的第125行）：

Box.test(zoo::na.approx(TS_FORECAST$residuals), type = "Ljung")

为了访问p值，只需使用$p.value，在将输出分配给变量后使用。

请注意，在我的快速示例中，有些不同，因为我使用了默认值：

# r$p.value
# [1] 0.3678976

- RLave

Answer 3

如果您需要一个数据框，现在checkresiduals返回一个对象（感谢@mischva11），您可以使用broom包中的tidy函数将其转换为data.frame（实际上是一个tibble，broom是tidyverse的一部分，但这足以满足要求）。

> library(broom)
> p_value <- tidy(checkresiduals(TS_FORECAST, plot = FALSE))

        Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 0.87319, df = 3, p-value = 0.8319

Model df: 1.   Total lags used: 4

> p_value
# A tibble: 1 x 4
  statistic p.value parameter method        
      <dbl>   <dbl>     <dbl> <chr>         
1     0.873   0.832         3 Ljung-Box test

不幸的是，输出仍然被打印出来，如上所示，这可能不是您想要的结果。我发现避免这种情况的唯一方法是使用invisible(capture.output())。

invisible(capture.output(p_value <- tidy(checkresiduals(TS_FORECAST, plot = FALSE))))

- Brian Stamper