How to extract p value from ca.po function in R?-CodePudding

I want to get the p-value of both ca.po models. Can someone show me how?

at?

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)

And the results. I marked the p-values I need in the result.

Model 1 p-value = 0.9841

Model 2 p-value = 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

CodePudding user response：

You have to dig into the res object to see its attributes and what's available there.

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

A long list of objects is returned, but we can see what is looking like the summary of lm being called under testreg, which we can see is one of the attributes of res. We can also access attributes of res using attr(res, "name"), so let's look at the components of 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"

As you noted above, you're looking for 2 separate p-values, so makes since we have two separate models. Let's retrieve these and look at what they are.

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

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

So, each of them is a summary.lm object. There's lots of documentation on how to extract p-values from lm or summary.lm objects, so let's use the method described 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

And there you go, those are the two p-values you were interested in!