Извлечь p-значение из функции проверки остатков - PullRequest
Новая
       4

Извлечь p-значение из функции проверки остатков

4 голосов
/ 18 апреля 2019

Я прогнозирую с пакетом прогноза. Ниже приведены результаты моего прогноза. 

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

Итак, теперь я хочу извлечь значение 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

Так кто-нибудь может мне помочь, как извлечь p-значение (p-значение = 0,4237) из функции checkresiduals в data.frame?

Ответы [ 2 ]

3 голосов
/ 18 апреля 2019

К сожалению, функция 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)
    }
  }

}

вам также нужна функция modeldf() из файла github 

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

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

С этим решением вы используете свою оригинальную функцию проверки остатков. Теперь вы можете вызвать 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 = 0,4367"

С помощью команды grep должна быть возможность извлечь значение p без изменения функции. Так как я не знаком с grep, я не могу дать правильный ответ на эту задачу.

3 голосов
/ 18 апреля 2019

Здесь вы можете увидеть внутреннюю часть checkresiduals().

К сожалению, согласно документам, он не возвращает значение, поэтому вы не можете легко извлечь то, что вынужно.

Но вы можете сделать тот же расчет (посмотрите на строку 125 в связанном репо): 

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

Чтобы получить доступ к p-значению, просто используйте $p.value после назначениявывод в переменную.

Обратите внимание, что в моем быстром примере это немного отличается, потому что я использовал значения по умолчанию: 

# r$p.value
# [1] 0.3678976
...