Calculating P values for correlation map

Question

I have mapped the correlation of rainfall in Australia to different driver variables. stk_b7 - is monthly rainfall values from 1925-1946 in 5km grids. This is a bricked stack with the following dimensions -

stk_b7 
class: RasterStack 
dimensions : 681, 841, 572721, 264  (nrow, ncol, ncell, nlayers)
resolution : 0.05, 0.05  (x, y)

then I have a csv file mv_b7 with all of my driver variables eg. SOI, NINO3, NINO4 etc columns are the variable and rows are the months for the appropriate time period. Jan 1925, Feb 1925, ...

I have used the following code to calculate the correlation and map it.

myfun <- function(x) cor(as.vector(x),mv_b7$SOI)

#use the function with calc

r <- calc(stk_b7, myfun)

#map results of function applied to raster brick
rasterVis::gplot(r) + 
  geom_tile(aes(fill = value)) + 
   scale_fill_gradient2(low="blue",high="red",mid="white") +
  
  coord_equal() +
  theme_bw() +
  labs(x = "Longitude",
       y = "Latitude",
       fill = "Correlation(pearson)",
       title = "1900-2013 IPO-Precipitation Correlation",
       subtitle = "COBE SST data provided by NOAA, SILO gridded rainfall") +
  #guides(fill = guide_legend(override.aes = list(size = 2),ncol = 2)) +
  theme(plot.title = element_text(size = 12,hjust = 0.5,face = "bold"),
        plot.subtitle = element_text(size = 8,hjust = 0.5),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7),
        axis.title = element_text(size = 10),
        axis.text = element_text(size = 7))

This produces the correlation map as required. However now I need to calculate the corresponding P values. I am not sure if this is something that I would map or if I would just get one P value for the entire correlation map. I have tried using rcorr and cor.test with the same format but get the following error.

    myfun <- function(x) cor.test(as.vector(x),mv_b7$SOI)
    r <- calc(stk_b7, myfun)
    Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
    cannot use this function

Not sure if this has to do with NA's however am confused as to why the normal correlation will work but not the others. Any help on this is greatly appreciated

Answer 1

This is because the output of cor.test has multiple values, and you need to be explicit in using cor.test in order to get the p-value of that function call.

a <- rnorm(4) #Generating fake data
b <- rnorm(4)

cor.test(a,b)

Pearson's product-moment correlation

data:  a and b
t = -0.47545, df = 2, p-value = 0.6813
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9797036  0.9260328
sample estimates:
       cor 
-0.3186697 

To actually access the p-value, you can see that in the function documentation for cor.test, you will need to access the third value in the output.

cor.test(a,b)[3]
$p.value
[1] 0.6813303

So I would change your function to make sure you access the third value in your cor.test() call and it should then work.