如何从 R 中的 ca.po 函数中提取 p 值？

Posted 2023-02-14

【中文标题】如何从 R 中的 ca.po 函数中提取 p 值？

【英文标题】：How to extract p value from ca.po function in R?

【发布时间】：2022-01-17 14:54:42

【问题描述】：

我想获得两个 ca.po 模型的 p 值。谁能告诉我怎么做？

在？

library(data.table)
library(urca)
dt_xy = as.data.table(timeSeries::LPP2005REC[, 2:3])
res = urca::ca.po(dt_xy, type = "Pz", demean = demean, lag = "short")
summary(res)

还有结果。我在结果中标记了我需要的 p 值。

模型 1 p 值 = 0.9841

模型 2 p 值 = 0.1363

######################################## 
# Phillips and Ouliaris Unit Root Test # 
######################################## 

Test of type Pz 
detrending of series with constant and linear trend 

Response SPI :

Call:
lm(formula = SPI ~ zr + trd)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.036601 -0.003494  0.000243  0.004139  0.024975 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  9.702e-04  7.954e-04   1.220    0.223
zrSPI       -1.185e-02  5.227e-02  -0.227    0.821
zrSII       -3.037e-02  1.374e-01  -0.221    0.825
trd         -6.961e-07  3.657e-06  -0.190    0.849

Residual standard error: 0.007675 on 372 degrees of freedom
Multiple R-squared:  0.0004236, Adjusted R-squared:  -0.007637 
F-statistic: 0.05255 on 3 and 372 DF,  p-value: 0.9841 **<--- I need this p.value**


Response SII :

Call:
lm(formula = SII ~ zr + trd)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0096931 -0.0018105 -0.0002734  0.0017166  0.0115427 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -7.598e-05  3.012e-04  -0.252   0.8010  
zrSPI       -1.068e-02  1.979e-02  -0.540   0.5897  
zrSII       -9.574e-02  5.201e-02  -1.841   0.0664 .
trd          1.891e-06  1.385e-06   1.365   0.1730  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.002906 on 372 degrees of freedom
Multiple R-squared:  0.01476,   Adjusted R-squared:  0.006813 
F-statistic: 1.858 on 3 and 372 DF,  p-value: 0.1363 **<--- I need this p.value**



Value of test-statistic is: 857.4274 

Critical values of Pz are:
                  10pct    5pct     1pct
critical values 71.9586 81.3812 102.0167

【问题讨论】：

您尝试过tidy() 包中的tidy() 函数吗？

tidy() 将查看 res 对象中的内容。我在res 对象中找不到p.value。我不确定我的方向是否正确。

【参考方案1】：

您必须深入了解res 对象才能查看其属性以及那里可用的内容。

attributes(reg)
...
#> 
#> $testreg
#> Response SPI :
#> 
#> Call:
#> lm(formula = SPI ~ zr + trd)
#> 
...

返回一长串对象，但我们可以看到lm 的摘要在testreg 下被调用，我们可以看到这是res 的属性之一。我们也可以使用attr(res, "name")来访问res的属性，下面我们来看看testreg的组成部分。

names(attributes(res))
#>  [1] "z"         "type"      "model"     "lag"       "cval"      "res"      
#>  [7] "teststat"  "testreg"   "test.name" "class"
names(attr(res, "testreg"))
#> [1] "Response SPI" "Response SII"

如上所述，您正在寻找 2 个单独的 p 值，因此我们有两个单独的模型。让我们检索这些并看看它们是什么。

spi <- attr(res, "testreg")[["Response SPI"]]
sii <- attr(res, "testreg")[["Response SII"]]

class(spi)
#> [1] "summary.lm"

所以，它们每个都是summary.lm 对象。有很多关于如何从lm 或summary.lm 对象中提取p 值的文档，所以让我们使用here 描述的方法。

get_pval <- function(summary_lm) 
  pf(
    summary_lm$fstatistic[1L],
    summary_lm$fstatistic[2L],
    summary_lm$fstatistic[3L],
    lower.tail = FALSE
  )


get_pval(spi)
#>     value 
#> 0.9840898
get_pval(sii)
#>     value 
#> 0.1363474

你去吧，这就是你感兴趣的两个 p 值！

【讨论】：

感谢@caldwellst 从这个黑盒中挖掘出 p.values。