Trazar las pendientes estimadas, como en la pregunta, es una gran cosa que hacer. Sin embargo, en lugar de filtrar por importancia, o junto con él, ¿por qué no trazar una medida de qué tan bien cada regresión se ajusta a los datos? Para esto, el error cuadrático medio de la regresión se interpreta fácilmente y es significativo.
Como ejemplo, el Rsiguiente código genera una serie temporal de 11 rásteres, realiza las regresiones y muestra los resultados de tres maneras: en la fila inferior, como cuadrículas separadas de pendientes estimadas y errores cuadrados medios; 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 le ofrece una comparación). La superposición, porque usa color para una variable (pendiente estimada) y claridad 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.
(Por favor, ignore las leyendas superpuestas en la superposición. Tenga en cuenta también que el esquema de color para el mapa "Verdaderas pendientes" no es exactamente el mismo que para los mapas de pendientes estimadas: un error aleatorio hace que algunas de las pendientes estimadas abarquen un rango más extremo que las pendientes verdaderas. 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 una gran cantidad de regresiones para el mismo conjunto de tiempos: en cambio, la matriz de proyección puede calcularse previamente y aplicarse a cada "pila" de píxeles más rápidamente que volver a calcularla para 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))
