Tracer les pentes estimées, comme dans la question, est une bonne chose à faire. Plutôt que de filtrer par importance, bien que - ou conjointement avec elle -, pourquoi ne pas tracer une mesure de la mesure dans laquelle chaque régression correspond aux données? Pour cela, l'erreur quadratique moyenne de la régression est facilement interprétée et significative.
À titre d'exemple, le Rcode ci-dessous génère une série temporelle de 11 rasters, effectue les régressions et affiche les résultats de trois manières: sur la ligne du bas, sous forme de grilles distinctes des pentes estimées et des erreurs quadratiques moyennes; sur la rangée supérieure, comme la superposition de ces grilles avec les véritables pentes sous-jacentes (que vous n’aurez en pratique jamais, mais que la simulation informatique permet de comparer). La superposition, parce qu’elle utilise la couleur pour une variable (pente estimée) et la luminosité pour une autre (MSE), n’est pas facile à interpréter dans cet exemple particulier, mais avec les cartes séparées sur la ligne du bas, elle peut être utile et intéressante.
(Ignorez les légendes superposées de la superposition. Notez également que le jeu de couleurs de la carte "Pentes vraies" n’est pas tout à fait le même que celui des cartes des pentes estimées: une erreur aléatoire fait en sorte que certaines des pentes gamme plus extrême que les pentes réelles. C’est un phénomène général lié à la régression vers la moyenne .)
En passant, ce n’est pas le moyen le plus efficace d’effectuer un grand nombre de régressions pour la même série de fois: la matrice de projection peut plutôt être précalculée et appliquée à chaque "pile" de pixels plus rapidement que de la recalculer pour chaque régression. Mais cela n'a pas d'importance pour cette petite illustration.
# 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))
