Graficar las pendientes estimadas, como en la pregunta, es una gran cosa para hacer. Sin embargo, en lugar de filtrar por significado, o en conjunto con él, ¿por qué no trazar una medida de qué tan bien se ajusta cada regresión a los datos? Para esto, el error cuadrático medio de la regresión se interpreta fácilmente y es significativo.
Como ejemplo, el código R
a continuación genera una serie de tiempo de 11 rásteres, realiza las regresiones y muestra los resultados de tres formas: en la fila inferior, como cuadrículas separadas de pendientes estimadas y errores de la media al cuadrado; en la fila superior, como la superposición de esas cuadrículas junto con las verdaderas pendientes subyacentes (que en la práctica nunca tendrá, pero la simulación por computadora permite la comparación). La superposición, porque utiliza el color para una variable (pendiente estimada) y la luminosidad para otra (MSE), no es fácil de interpretar en este ejemplo en particular, pero junto con los mapas separados en la fila inferior puede ser útil e interesante.

(Ignorelasleyendassuperpuestasenlasuperposición.Tengaencuenta,también,queelesquemadecolorparaelmapade"Pendientes verdaderas" no es exactamente el mismo que para los mapas de pendientes estimadas: el error aleatorio causa algunas de las pendientes estimadas para abarcar un rango más extremo que las pendientes reales. Este es un fenómeno general relacionado con la regresión hacia la media .)
Por cierto, esta no es la forma más eficiente de hacer un gran número de regresiones para el mismo conjunto de veces: en cambio, la proyección La matriz se puede calcular previamente y aplicar a cada "pila" de píxeles más rápidamente que para calcularla en cada regresión. Pero eso no importa para esta pequeña ilustración.
# Specify the extent in space and time.
#
n.row <- 60; n.col <- 100; n.time <- 11
#
# Generate data.
#
set.seed(17)
sd.err <- outer(1:n.row, 1:n.col, function(x,y) 5 * ((1/2 - y/n.col)^2 + (1/2 - x/n.row)^2))
e <- array(rnorm(n.row * n.col * n.time, sd=sd.err), dim=c(n.row, n.col, n.time))
beta.1 <- outer(1:n.row, 1:n.col, function(x,y) sin((x/n.row)^2 - (y/n.col)^3)*5) / n.time
beta.0 <- outer(1:n.row, 1:n.col, function(x,y) atan2(y, n.col-x))
times <- 1:n.time
y <- array(outer(as.vector(beta.1), times) + as.vector(beta.0),
dim=c(n.row, n.col, n.time)) + e
#
# Perform the regressions.
#
regress <- function(y) {
fit <- lm(y ~ times)
return(c(fit$coeff[2], summary(fit)$sigma))
}
system.time(b <- apply(y, c(1,2), regress))
#
# Plot the results.
#
library(raster)
plot.raster <- function(x, ...) plot(raster(x, xmx=n.col, ymx=n.row), ...)
par(mfrow=c(2,2))
plot.raster(b[1,,], main="Slopes with errors")
plot.raster(b[2,,], add=TRUE, alpha=.5, col=gray(255:0/256))
plot.raster(beta.1, main="True slopes")
plot.raster(b[1,,], main="Estimated slopes")
plot.raster(b[2,,], main="Mean squared errors", col=gray(255:0/256))