terra::regress takes a SpatRaster with multiple layers and returns a raster with the itercept and slopes for each cell. I use this method to estimate the trend in climatic time series surfaces.
It is possible to also get the P-value of the slope for each cell? And can this be extended to other non-parametric metrics such as Theil-sen and Mann-kendall?
Here is a minimal, self-contained, reproducible example (please always include one).
library(terra)
r <- rast(system.file("ex/logo.tif", package="terra"))
x <- regress(r, 1:nlyr(r))
As you point out, that does not return p-values. To get the p-value you would need to write your own function. Something like this:
f <- function(Y, X) {
apply(Y, 1, function(y) {
m <- lm(y ~ X)
s <- summary(m)
as.vector(coefficients(s)[2, c(1,4)])
})
}
out <- app(r, fun=f, X=1:3, wopt=list(names=c("slope", "pvalue")))
This is much slower though. I could probably implement that in regress but I am unsure how useful that would be. After all, this is not hypothesis testing. The slope might be informative enough?