Retornando múltiplas saídas de um cálculo raster usando raster ou terra

Estou executando modelos de ajuste linear em uma série de rasters muito grandes. Frequentemente preciso retornar dois ou mais valores do modelo - por exemplo, o coeficiente e seu valor p associado. Em um caso não raster, eu faria uma função que retornasse todos os valores que eu queria como uma lista, assim:

#quick lm function to pull out coefficient and p.value
lm_fun <- function(x, y){
  lm <- lm(x ~ y)
  coef1 <- coef(summary(lm))[2,1]
  pval <- coef(summary(lm))[2,4]
  lm_results <- list(coef1 = coef1, pval = pval)
}

#sample data
dat <- c(runif(100, 2, 5), runif(100, 4, 10), runif(100, 6, 14))
timestep <- 1:length(dat)

lm_results <- lm_fun(dat, timestep)

No entanto, ao trabalhar com dados raster, uma função personalizada é aplicada ao raster por meio de calc (com o pacote raster) ou app (com o pacote terra) e retorna uma camada raster ou SpatRaster. Existe alguma maneira de modificar a função para produzir uma pilha raster, com saídas diferentes (por exemplo, coeficiente, valor p) em camadas diferentes? No momento, estou usando a abordagem abaixo, mas isso significa que estou executando o mesmo cálculo várias vezes para extrair partes diferentes dos resultados, o que é extremamente ineficiente computacionalmente. Estou aberto a usar raster ou terra.

# create three identical RasterLayer objects
r1 <- r2 <- r3 <- raster(nrow=100, ncol=100)
# Assign random cell values
values(r1) <- runif(ncell(r1), min=2, max=5)
values(r2) <- runif(ncell(r2), min=4, max=10) 
values(r3) <- runif(ncell(r3), min=6, max=14)
# combine three RasterLayer objects into a RasterStack
s <- stack(r1, r2, r3)
plot(s)

#calculate number of timesteps
nsteps <- dim(s)[3]
timestep <- 1:nsteps

#functions to calculate linear trend and associated p-value
trendfun <- function(x) { if (is.na(x[1])){ NA } else { lm(x ~ timestep)$coefficients[2]}}
pfun <- function(x) { if (is.na(x[1])){ NA } else { summary(lm(x ~ timestep))$coefficients[2,4]}}

# calculate trend and p-value using raster
library(raster)
s_trend <- calc(s, trendfun)
s_trend_pv <- calc(s, pfun)

#alternatively, calculate using terra package
library(terra)
s_rast <- rast(s)
s_trend <- app(s_rast, fun=trendfun)
s_pvalue <- app(s_rast, fun=pfun)