在 textreg ARIMA 输出中包含观察数
Include number of observations in textreg ARIMA output
我运行 几个时间序列回归。我的示例使用不同的数据，并且我没有注意使这些模型正确 - 它们只是对我的问题的说明。抱歉，我使用的是下载的数据集，而不是创建的数据集。我希望没关系。
我的任务是打印各种模型的 texreg 输出，包括但不限于 ARIMA 预测。我需要包括所有模型的观察次数。 texreg 提取函数不包括 ARIMA 对象的观察次数 - 或者至少观察次数永远不会打印在 gof 统计信息中。
我想到了两个选项，我希望得到帮助来实现其中一个或两个：
- 向 gof 统计添加一个自定义行，其中包含观察次数。 应该 使用 custom.gof.rows 是可能的，正如所讨论的 here and here。但是，从 1.36.23 版开始，custom.gof.rows 并未出现在 texreg 包的一部分中。也许它是出于某种原因被取出来的？不确定。
如果有包含 custom.gof.rows 的 texreg 的秘密版本，任何人都可以 link 吗？
library('ggplot2')
library('forecast')
library('tseries')
library('texreg')
library('lmtest')
#data can be downloaded from: https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset
#Change to data path
path<-c("")
#Import data
daily_data = read.csv(paste(path,'day.csv',sep=""), header=TRUE, stringsAsFactors=FALSE)
#Mark a shift date
daily_data$shift<- ifelse(daily_data$instant>12,1,0)
#Create time series data
dt <- ts(daily_data$temp, frequency=12, start= c(2011, 1))
#Define input vars for ARIMA
(xvars <- as.matrix(daily_data[, c("instant", "shift")]))
#Basic ARIMA
a<- arima(dt, xreg=xvars)
#Auto-ARIMA
b<- auto.arima(dt, xreg=xvars)
##Where I want to include number of observations either automatically or using custom.gof.rows
screenreg(list(coeftest(b),a))
This is the output, which I had to link (sorry)
请注意，模型对象本身并不是原因（我认为）。我可以自己提取观察的数量并单独打印。
#Both models have number of observations in their objects that I can extract
nobs_a<-nobs(a)
nobs_b<-nobs(b)
cat("Number of obs ARIMA: ",nobs_a,"\n")
cat("Number of obs Auto-ARIMA: ",nobs_b,"\n")
##Custom.gof.rows doesn't appear to be in texreg version 1.36.23
screenreg((list(coeftest(b),a)), custom.gof.rows = list(nobs_a,nobs_b))
- 修改提取函数以包括观察的数量，以便自动包括它们，如扩展文档 here 所讨论的那样。这是我完全陌生的，但我在下面尝试修改 textreg extract.Arima 函数以使用 nobs 检索观察次数。
function (model, include.pvalues = FALSE, include.aic = TRUE,
include.loglik = TRUE, ...) { mask <- model$mask nam <- names(model$coef)
co <- model$coef sdev <-
sqrt(diag(model$var.coef)) if (include.pvalues == TRUE) {
t.rat <- rep(NA, length(mask))
t.rat[mask] <- co[mask]/sdev
pt <- 2 * pnorm(-abs(t.rat))
setmp <- rep(NA, length(mask))
setmp[mask] <- sdev } else {
pt <- numeric()
setmp <- sdev } gof <- numeric() gof.names <- character() gof.decimal <- logical() if
(include.aic == TRUE) {
aic <- AIC(model)
gof <- c(gof, aic)
gof.names <- c(gof.names, "AIC")
gof.decimal <- c(gof.decimal, TRUE) }
##This is the section I added - ripped more or less intact from the extract.lm function
if (include.nobs == TRUE) {
gof <- c(gof, nobs)
gof.names <- c(gof.names, "Num. obs.")
gof.decimal <- c(gof.decimal, FALSE) }
if (include.loglik == TRUE) {
lik <- model$loglik
gof <- c(gof, lik)
gof.names <- c(gof.names, "Log Likelihood")
gof.decimal <- c(gof.decimal, TRUE) }
tr <- createTexreg(coef.names = nam, coef = co, se = setmp,
pvalues = pt, gof.names = gof.names, gof = gof, gof.decimal = gof.decimal) >
return(tr) } <bytecode:
这会运行但会破坏 textreg 函数。如果我走错了路，我也很乐意 link 阅读有关一般更改提取函数的教程。您可能还会注意到我在上面的代码中使用 coeftest() 来提取 Auto-ARIMA 系数。我对为 Auto-ARIMA 预测对象编写一个 exract 函数非常感兴趣。我认为这无论如何都是必要的，因为 nobs 仅适用于 b 本身，而不适用于 coeftest(b)。不过那是题外话——我到了那儿就能弄明白。
任何人都可以通过其中一种或两种途径提供帮助，或者可以提出一种不同的方法来在 texreg 输出中包含观察值的数量吗？
非常感谢您的帮助。
custom.gof.rows 参数大约在一年前添加到 texreg、screenreg 等函数中，尚未找到进入 CRAN 的途径，但最终会。目前，如果您想使用此参数，可以从 GitHub 安装最新版本。您可以按如下方式进行（并且需要先安装 remotes 包）：
library(remotes)
install_github("leifeld/texreg")
custom.gof.rows 参数采用命名的向量列表。您只是简单地为 custom.gof.rows 参数提供了两个值作为单独的参数作为观察值。相反，您需要将两者都包装在一个向量中并为其附加一个名称，例如 "Num. obs."。您可以按如下方式执行此操作：
screenreg(list(coeftest(b), a),
custom.gof.rows = list("Num. obs." = c(nobs(b), nobs(a))))
就是说，无论如何，我现在已经将观察次数添加到 ARIMA 模型的 extract 方法中，因此您将不再需要此参数。
您尝试自行修改 extract 方法几乎成功了。您还需要在函数头中包含 include.nobs = TRUE 参数，并且您需要将该函数注册为通用 extract 函数的方法，如 article in the Journal of Statistical Software 中所述。
为了让您的生活更轻松，我为您做了这件事。因此，如果您按照上述说明从 GitHub 安装最新版本的 texreg，它应该会自动包含观察次数。
您正在使用的 arima 函数定义在 stats 包中。它 return 是一个 Arima 对象。我相应地更新了它的 extract 方法。您在上面使用的 auto.arima 函数是同一函数的包装器，并且 returns 是一个具有 class forecast_ARIMA 的对象（之前和其次，ARIMA ).这有点令人困惑，但了解这些区别很重要。我已经用观察次数为你更新了这两种方法，所以你的代码
screenreg(list(coeftest(b), a))
现在应该 return 以下输出：
======================================
Model 1 Model 2
--------------------------------------
ar1 0.96 ***
(0.04)
ar2 -0.30 ***
(0.05)
ar3 0.13 *
(0.05)
ar4 0.02
(0.05)
ar5 0.17 ***
(0.04)
intercept 0.43 ** 0.21 ***
(0.13) (0.05)
instant 0.00 0.00 *
(0.00) (0.00)
shift 0.02 0.25 ***
(0.05) (0.05)
--------------------------------------
AIC -440.26
BIC -421.88
Log Likelihood 224.13
Num. obs. 731
======================================
*** p < 0.001; ** p < 0.01; * p < 0.05
但是请注意，您正在使用 lmtest 包中的 coeftest 函数用于模型 b，它有一个单独的 class，因此 [=24] =] 方法，而不是 return 观察次数。要获得模型 b 的观察次数，您可以直接包含它，而不需要 coeftest:
screenreg(list(b, a), single.row = TRUE)
这 return 的输出如下：
=======================================================
Model 1 Model 2
-------------------------------------------------------
ar1 0.96 (0.04) ***
ar2 -0.30 (0.05) ***
ar3 0.13 (0.05) *
ar4 0.02 (0.05)
ar5 0.17 (0.04) ***
intercept 0.43 (0.13) ** 0.21 (0.05) ***
instant 0.00 (0.00) 0.00 (0.00) *
shift 0.02 (0.05) 0.25 (0.05) ***
-------------------------------------------------------
AIC -2183.71 -440.26
AICc -2183.46
BIC -2142.36 -421.88
Log Likelihood 1100.86 224.13
Num. obs. 731 731
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05
如果你坚持使用coeftest，你可以将extract函数的输出保存在一个中间对象中，然后操作这个对象。在下面的代码中，我从没有 coeftest 的版本中取出 GOF 块并将其注入到带有 coeftest 的版本中，然后将操作对象而不是原始对象交给 screenreg功能：
tr_b_coeftest <- extract(coeftest(b))
tr_b_plain <- extract(b)
tr_b_coeftest@gof.names <- tr_b_plain@gof.names
tr_b_coeftest@gof <- tr_b_plain@gof
tr_b_coeftest@gof.decimal <- tr_b_plain@gof.decimal
screenreg(list(tr_b_coeftest, a), single.row = TRUE)
这 return 的输出如下：
=======================================================
Model 1 Model 2
-------------------------------------------------------
ar1 0.96 (0.04) ***
ar2 -0.30 (0.05) ***
ar3 0.13 (0.05) *
ar4 0.02 (0.05)
ar5 0.17 (0.04) ***
intercept 0.43 (0.13) ** 0.21 (0.05) ***
instant 0.00 (0.00) 0.00 (0.00) *
shift 0.02 (0.05) 0.25 (0.05) ***
-------------------------------------------------------
AIC -2183.71 -440.26
AICc -2183.46
BIC -2142.36 -421.88
Log Likelihood 1100.86 224.13
Num. obs. 731 731
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05
我运行 几个时间序列回归。我的示例使用不同的数据，并且我没有注意使这些模型正确 - 它们只是对我的问题的说明。抱歉，我使用的是下载的数据集，而不是创建的数据集。我希望没关系。
我的任务是打印各种模型的 texreg 输出，包括但不限于 ARIMA 预测。我需要包括所有模型的观察次数。 texreg 提取函数不包括 ARIMA 对象的观察次数 - 或者至少观察次数永远不会打印在 gof 统计信息中。
我想到了两个选项，我希望得到帮助来实现其中一个或两个：
- 向 gof 统计添加一个自定义行，其中包含观察次数。 应该 使用 custom.gof.rows 是可能的，正如所讨论的 here and here。但是，从 1.36.23 版开始，custom.gof.rows 并未出现在 texreg 包的一部分中。也许它是出于某种原因被取出来的？不确定。
如果有包含 custom.gof.rows 的 texreg 的秘密版本，任何人都可以 link 吗？
library('ggplot2') library('forecast') library('tseries') library('texreg') library('lmtest') #data can be downloaded from: https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset #Change to data path path<-c("") #Import data daily_data = read.csv(paste(path,'day.csv',sep=""), header=TRUE, stringsAsFactors=FALSE) #Mark a shift date daily_data$shift<- ifelse(daily_data$instant>12,1,0) #Create time series data dt <- ts(daily_data$temp, frequency=12, start= c(2011, 1)) #Define input vars for ARIMA (xvars <- as.matrix(daily_data[, c("instant", "shift")])) #Basic ARIMA a<- arima(dt, xreg=xvars) #Auto-ARIMA b<- auto.arima(dt, xreg=xvars) ##Where I want to include number of observations either automatically or using custom.gof.rows screenreg(list(coeftest(b),a))
This is the output, which I had to link (sorry)
请注意，模型对象本身并不是原因（我认为）。我可以自己提取观察的数量并单独打印。
#Both models have number of observations in their objects that I can extract nobs_a<-nobs(a) nobs_b<-nobs(b) cat("Number of obs ARIMA: ",nobs_a,"\n") cat("Number of obs Auto-ARIMA: ",nobs_b,"\n") ##Custom.gof.rows doesn't appear to be in texreg version 1.36.23 screenreg((list(coeftest(b),a)), custom.gof.rows = list(nobs_a,nobs_b))
- 修改提取函数以包括观察的数量，以便自动包括它们，如扩展文档 here 所讨论的那样。这是我完全陌生的，但我在下面尝试修改 textreg extract.Arima 函数以使用 nobs 检索观察次数。
function (model, include.pvalues = FALSE, include.aic = TRUE, include.loglik = TRUE, ...) { mask <- model$mask nam <- names(model$coef) co <- model$coef sdev <- sqrt(diag(model$var.coef)) if (include.pvalues == TRUE) { t.rat <- rep(NA, length(mask)) t.rat[mask] <- co[mask]/sdev pt <- 2 * pnorm(-abs(t.rat)) setmp <- rep(NA, length(mask)) setmp[mask] <- sdev } else { pt <- numeric() setmp <- sdev } gof <- numeric() gof.names <- character() gof.decimal <- logical() if (include.aic == TRUE) { aic <- AIC(model) gof <- c(gof, aic) gof.names <- c(gof.names, "AIC") gof.decimal <- c(gof.decimal, TRUE) } ##This is the section I added - ripped more or less intact from the extract.lm function if (include.nobs == TRUE) { gof <- c(gof, nobs) gof.names <- c(gof.names, "Num. obs.") gof.decimal <- c(gof.decimal, FALSE) } if (include.loglik == TRUE) { lik <- model$loglik gof <- c(gof, lik) gof.names <- c(gof.names, "Log Likelihood") gof.decimal <- c(gof.decimal, TRUE) } tr <- createTexreg(coef.names = nam, coef = co, se = setmp, pvalues = pt, gof.names = gof.names, gof = gof, gof.decimal = gof.decimal) > return(tr) } <bytecode:
这会运行但会破坏 textreg 函数。如果我走错了路，我也很乐意 link 阅读有关一般更改提取函数的教程。您可能还会注意到我在上面的代码中使用 coeftest() 来提取 Auto-ARIMA 系数。我对为 Auto-ARIMA 预测对象编写一个 exract 函数非常感兴趣。我认为这无论如何都是必要的，因为 nobs 仅适用于 b 本身，而不适用于 coeftest(b)。不过那是题外话——我到了那儿就能弄明白。
任何人都可以通过其中一种或两种途径提供帮助，或者可以提出一种不同的方法来在 texreg 输出中包含观察值的数量吗？
非常感谢您的帮助。
custom.gof.rows 参数大约在一年前添加到 texreg、screenreg 等函数中，尚未找到进入 CRAN 的途径，但最终会。目前，如果您想使用此参数，可以从 GitHub 安装最新版本。您可以按如下方式进行（并且需要先安装 remotes 包）：
library(remotes)
install_github("leifeld/texreg")
custom.gof.rows 参数采用命名的向量列表。您只是简单地为 custom.gof.rows 参数提供了两个值作为单独的参数作为观察值。相反，您需要将两者都包装在一个向量中并为其附加一个名称，例如 "Num. obs."。您可以按如下方式执行此操作：
screenreg(list(coeftest(b), a),
custom.gof.rows = list("Num. obs." = c(nobs(b), nobs(a))))
就是说，无论如何，我现在已经将观察次数添加到 ARIMA 模型的 extract 方法中，因此您将不再需要此参数。
您尝试自行修改 extract 方法几乎成功了。您还需要在函数头中包含 include.nobs = TRUE 参数，并且您需要将该函数注册为通用 extract 函数的方法，如 article in the Journal of Statistical Software 中所述。
为了让您的生活更轻松，我为您做了这件事。因此，如果您按照上述说明从 GitHub 安装最新版本的 texreg，它应该会自动包含观察次数。
您正在使用的 arima 函数定义在 stats 包中。它 return 是一个 Arima 对象。我相应地更新了它的 extract 方法。您在上面使用的 auto.arima 函数是同一函数的包装器，并且 returns 是一个具有 class forecast_ARIMA 的对象（之前和其次，ARIMA ).这有点令人困惑，但了解这些区别很重要。我已经用观察次数为你更新了这两种方法，所以你的代码
screenreg(list(coeftest(b), a))
现在应该 return 以下输出：
======================================
Model 1 Model 2
--------------------------------------
ar1 0.96 ***
(0.04)
ar2 -0.30 ***
(0.05)
ar3 0.13 *
(0.05)
ar4 0.02
(0.05)
ar5 0.17 ***
(0.04)
intercept 0.43 ** 0.21 ***
(0.13) (0.05)
instant 0.00 0.00 *
(0.00) (0.00)
shift 0.02 0.25 ***
(0.05) (0.05)
--------------------------------------
AIC -440.26
BIC -421.88
Log Likelihood 224.13
Num. obs. 731
======================================
*** p < 0.001; ** p < 0.01; * p < 0.05
但是请注意，您正在使用 lmtest 包中的 coeftest 函数用于模型 b，它有一个单独的 class，因此 [=24] =] 方法，而不是 return 观察次数。要获得模型 b 的观察次数，您可以直接包含它，而不需要 coeftest:
screenreg(list(b, a), single.row = TRUE)
这 return 的输出如下：
=======================================================
Model 1 Model 2
-------------------------------------------------------
ar1 0.96 (0.04) ***
ar2 -0.30 (0.05) ***
ar3 0.13 (0.05) *
ar4 0.02 (0.05)
ar5 0.17 (0.04) ***
intercept 0.43 (0.13) ** 0.21 (0.05) ***
instant 0.00 (0.00) 0.00 (0.00) *
shift 0.02 (0.05) 0.25 (0.05) ***
-------------------------------------------------------
AIC -2183.71 -440.26
AICc -2183.46
BIC -2142.36 -421.88
Log Likelihood 1100.86 224.13
Num. obs. 731 731
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05
如果你坚持使用coeftest，你可以将extract函数的输出保存在一个中间对象中，然后操作这个对象。在下面的代码中，我从没有 coeftest 的版本中取出 GOF 块并将其注入到带有 coeftest 的版本中，然后将操作对象而不是原始对象交给 screenreg功能：
tr_b_coeftest <- extract(coeftest(b))
tr_b_plain <- extract(b)
tr_b_coeftest@gof.names <- tr_b_plain@gof.names
tr_b_coeftest@gof <- tr_b_plain@gof
tr_b_coeftest@gof.decimal <- tr_b_plain@gof.decimal
screenreg(list(tr_b_coeftest, a), single.row = TRUE)
这 return 的输出如下：
=======================================================
Model 1 Model 2
-------------------------------------------------------
ar1 0.96 (0.04) ***
ar2 -0.30 (0.05) ***
ar3 0.13 (0.05) *
ar4 0.02 (0.05)
ar5 0.17 (0.04) ***
intercept 0.43 (0.13) ** 0.21 (0.05) ***
instant 0.00 (0.00) 0.00 (0.00) *
shift 0.02 (0.05) 0.25 (0.05) ***
-------------------------------------------------------
AIC -2183.71 -440.26
AICc -2183.46
BIC -2142.36 -421.88
Log Likelihood 1100.86 224.13
Num. obs. 731 731
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05