【问题标题】:Extract p-value from checkresiduals function从检查残差函数中提取 p 值
【发布时间】:2019-04-18 08:27:34
【问题描述】:

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

#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-value (p-value = 0.4237) 到 data.frame 中吗?

【问题讨论】:

    标签: r forecasting


    【解决方案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 filemodeldf() 函数

    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-BoxBreusch-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 值 = 0.4367”

    使用 grep 命令应该可以在不更改函数的情况下提取 p 值。由于我不熟悉 grep,因此我无法就此任务提供正确的答案。

    【讨论】:

    • 为了完整性:我将我的答案更改为 Github 提交。应该会在下一次更新中更改。
    • 此更新(来自@mischva11)现在是checkresiduals 的一部分,因此它现在确实返回值。如果需要,请参阅我关于转换为 data.frame 的回答。
    • @Brian Stamper 谢谢布赖恩,我更新了我的答案,我也喜欢你的整洁方法
    【解决方案2】:

    Here你可以看到checkresiduals()的内部。

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

    但是您可以进行相同的计算(查看链接回购中的第 125 行):

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

    要访问 p 值,只需在将输出分配给变量后使用 $p.value

    请注意,在我的快速示例中,它有点不同,因为我使用了默认值:

    # r$p.value
    # [1] 0.3678976
    

    【讨论】:

      【解决方案3】:

      如果您特别需要一个数据框,现在checkresiduals 返回一个对象(感谢@mischva11),您可以使用broom 包中的tidy 函数将其转换为data.frame(实际上是一个tibblebroomtidyverse 的一部分,但这应该足够了)。

      > 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))))
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2012-09-27
        • 1970-01-01
        • 2020-05-12
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多