I have single band rasters of data values averaged per year from 1986-2023. I would like to run a trend anlysis that gives me the linear regressions of each pixel across the years. I am not sure where to start but I can forsee the problem that if I run a function on the raster stack the year data is lost as it is only saved in the image file name and doesn’t appear in the stack dataframe.

I have so far attempted to run

fun <- function (x) { if (is.na(x[ 1 ])) else {cochrane.orcutt(lm(x[ 1 :a] ~ x[( 1 +a):(a+b)]))$coefficients[ 2 ] }} slope <- calc(stack, fun)