En écologie on se retrouve, très souvent, face à des processus qui sont liés à l’espace mais aussi qui évoluent dans le temps. Nous sommes maintenant habitués à l’analyse spatiale qui est devenu à la portée du plus grand nombre grâce aux SIG, tels QGis. Nous sommes beaucoup moins habitués à l’analyse temporelle qui est moins prise en compte par les SIG.

Analyse spatio-temporelle

Que dire alors de l’analyse spatio-temporelle, c’est à dire de l’analyse simultanée d’un processus variant dans l’espace et dans le temps? Comme pour l’analyse spatiale, un premier niveau pour aborder le processus qui nous intéresse est l’analyse exploratoire des données. Un deuxième niveau nous amène à essayer de comprendre les paramètres qui contrôlent le processus pour remplir les trous dans les données et réduire le « bruit » de nos observations. Enfin, au plus haut du traitement, nous pouvons modéliser statistiquement le processus et obtenir un outil qui nous permette une vision complétée de nos données (boucher les trous dans l’espace et le temps de nos observations), de trier les paramètres du processus et de quantifier leur importance, et finalement, de prévoir l’état futur de notre processus.

Si nous prenons seulement l’aspect spatial des données on comprendra mieux le sujet. Nous disposons, par exemple, de données sur une population d’oiseaux dont le territoire progresse d’année en année. D’un point de vue spatial nous aurons une distribution du nombre d’individus. Le premier pas pour traiter ces données sera l’analyse exploratoire des données (histogrammes, qqplots, polygones de Voronoï,…). On aura alors une idée de l’adéquation de notre échantillonnage et le processus sous-jacent qui nous intéresse.

Le deuxième pas sera de passer de nos échantillons ponctuels à des données surfaciques, soit par des interpolations simples (IDW, Spline,…) soit par des méthodes statistiques (Krigeage,..). A ce point on disposera d’une image de notre processus, nettoyée des erreurs de mesure et du « bruit » des observations.

Modélisation des données

On pourra alors utiliser des méthodes de modélisation telles que les modèles linéaires, les modèles linéaires généralisés (GLM), les modèles linéaires mixtes généralisés (GLMM),les modèles additifs généralisés (GAM), les modèles mixtes additifs généralisés (GAMM)… qui nous permettent d’inclure nos connaissances scientifiques sur le processus et d’obtenir des prédictions moins biaisées.

Les outils que nous aborderons dans cet article se situent un cran plus loin que ces modèles (GLM, GAM,…) pour traiter des données spatio-temporelles. En effet, les différents modèles linéaires correspondent à ce que l’on appelle l’approche descriptive de la modélisation spatio-temporelle. Pour faire simple, on définit la corrélation entre deux valeurs dans l’espace et le temps.

La modélisation IDE correspond à ce que l’on appelle l’approche dynamique. Le modèle considère que la valeur actuelle à un endroit donné est égale à un « facteur de propagation » de la valeur à ce point au moment précédent à laquelle on ajoute une « erreur d’innovation ».

Bref, on est en plein dans un domaine qu’on ne peut aborder à la légère. Il faut bien sacrifier du temps à comprendre ce que l’on fait, pourquoi et comment on le fait. Dans cet article j’éviterai les explications mathématiques. Vous trouverez les liens qui renvoient vers ces explications. Pour rester dans le domaine habituel des SIG je ferai des simplifications pas toujours exactes à 100% du point de vue mathématique, je m’en excuse d’avance au près des mathématiciens.

Les modèles d’équation intégro-différentielle (IDE) décrivent la dépendance conditionnelle entre le processus spatial à un moment futur et le processus au moment présent par l’intermédiaire d’un opérateur intégral.

La non linéarité ou la dépendance temporelle dans la dynamique du modèle est souvent intégrée en permettant aux paramètres de varier temporellement, ou en ajustant un modèle avec un opérateur linéaire invariant temporellement dans une fenêtre coulissante. Les deux procédures ont tendance à être excellentes à des fins de prévision sur de courts horizons temporels, mais elles prennent généralement beaucoup de temps de calcul.

En termes moins techniques, les modèles IDE sont des outils très puissant pour prévoir l’état suivant d’un système spatio-temporel à partir des données des états précédents et de l’état actuel.

Pour savoir ce qu’est IDE

Pour comprendre les modèles spatio-temporels dynamiques, reportez-vous au chapitre 5 de l’ouvrage Spatio-Temporal Statistics with R de Christopher K. Wikle, Andrew Zammit-Mangion, et Noel Cressie, piublié par Chapman & Hall/CRC et disponible en téléchargement ici.

La documentation de la bibliothèque R d’IDE est disponible ici.

Les traitements QGis et R

Au jour d’aujourd’hui, il n’y a pas de traitement disponible directement avec QGis. Ni Grass, ni SAGA ni les autres fournisseurs de traitements habituels de QGis implémentent IDE.

Par contre, un package R est dédié à la mise en place de modèles IDE. Si vous êtes un habitué de R, vous pouvez récupérer les données spatio-temporelles de QGis et les traiter dans RStudio.

Nous avons vu comment migrer les données de QGis vers RStudio dans cet article.

Pour faciliter l’accès aux modèles IDE directement dans QGis, nous avons développé deux scripts R directement utilisables à partir des traitements QGis.

Préalablement, vous devez vous assurer que les traitements QGis 3.x sont bien configurés pour l’utilisation de scripts R.

Vous trouverez comment faire dans cet article.

Téléchargements des scripts R et des données d’exemple

Cliquez sur le bouton « Télécharger » pour récupérer un fichier zippé, contenant deux scripts R et un fichier de données qui vous permettra de suivre la démarche de cet article.

Si non, vous avez ici le code du script

Code du script ide_modele_invariant

modele IDE=group

Layer=vector

t=Field Layer

val=Field Layer

grid_size=number 50

units=selection secs;mins;hours;days;weeks

forecast=number 1

itermax=number 100

effectuer_la_validation_du_modele=boolean TRUE

pas_de_temps_de_validation=number

pourcentage_points_aleatoires_de_validation=number 10

repertoire_resultats=optional folder

nom_du_fichier= optional string

Resultat_grille= output vector

Resultat_observations= output vector

library(IDE)
library(dplyr)
library(« sp »)
library(« spacetime »)
library(« sf »)
library(« verification »)
library(« scoringRules »)
library(« tidyr »)
library(« ggplot2 »)
library(« gridExtra »)
df<-as.data.frame(Layer)
coords<-st_coordinates(Layer)
crs<-st_crs(Layer)$proj4string
dfc<-cbind(df,coords)
data<-dplyr::select(dfc,t,all_of(val),X,Y)

names(data)<-c(« t », »cnt », »X », »Y »)
data$cnt<-as.double(data$cnt)
data$X<-as.double(data$X)
data$Y<-as.double(data$Y)
is.date <- function(x) inherits(x, ‘Date’)
is.convertible.to.date <- function(x) !is.na(as.Date(as.character(x),origin= »0000-01-01″))
if (!sapply(data$t[1], is.date)&& !is.character(data$t[1])) {
if (near(data$t[1], as.integer(data$t[1]))){
data$t<-as.Date(data$t,origin= »0000-01-01″)
trans=1
}
}
if (!sapply(data$t[1], is.date)&& is.character(data$t[1])) {
if (sapply(data$t[1], is.convertible.to.date)){
data$t<-as.Date(data$t)
trans=2
}
}
if (effectuer_la_validation_du_modele) {
tv<-pas_de_temps_de_validation
t0<-min(data$t)
tmax<-max(data$t)
if (tv<=t0) {
tv<-t0+1
}
if (tv>tmax) {
tv<-tmax-1
}
if (tv==tmax && forecast==0){
forecast=1
}

paste(« pas de temps retenu pour validation= »,tv)
}
mSTIDF <- stConstruct(x = data, # data set
space = c(« X », « Y »), # spatial fields
time = t) # time field
if (effectuer_la_validation_du_modele) {

échantillonnage pour validation

block pas de temps

mtot<-length(mSTIDF)
mSTIDF$timeY <- format(time(mSTIDF), « %Y-%m-%d »)
if (trans==1) {
vt<-as.Date(tv,origin= »0000-01-01″)
}else{
vt<-as.Date(tv)
}
valblock_idx<-which(mSTIDF$timeY%in%c(as.character(vt)))
paste(« Le bloc de validation du pas de temps comporte »,length(valblock_idx), »observations »)
obs_idx <- setdiff(1:mtot, valblock_idx)

pourcentage d’observations

pc<-pourcentage_points_aleatoires_de_validation/100
set.seed(1)
valrandom_idx <- sample(obs_idx,
pc * length(obs_idx),
replace = FALSE) %>% sort()
paste(« L’echantillon aleatoire de validation comporte »,length(valrandom_idx), »observations »)
obs_idx <- setdiff(obs_idx, valrandom_idx)
paste(« Le modèle sera calculé à partir de »,length(obs_idx), »observations »)
obs<-mSTIDF[obs_idx,] block<-mSTIDF[valblock_idx,] random<-mSTIDF[valrandom_idx,] }else{
obs<-mSTIDF
}
if (units==0) {unit= »secs »}
if (units==1) {unit= »mins »}
if (units==2) {unit= »hours »}
if (units==3) {unit= »days »}
if (units==4) {unit= »weeks »}
« calcul du modele IDE… »
IDEmodel <- IDE(f = cnt ~ 1 ,
data = obs,
dt = as.difftime(1, units = unit),
grid_size = grid_size,
forecast=forecast)
« Ajustement du modele… »
fit_results_obs <- fit.IDE(IDEmodel,method= »DEoptim »,itermax=itermax)
kp<-fit_results_obs$IDEmodel$get(« k ») %>% unlist()
« kernel parameters: « 
kp
co<-coef(fit_results_obs$IDEmodel)
« regression coefficients: « 
coef(fit_results_obs$IDEmodel)
sigma2_eps<-fit_results_obs$IDEmodel$get(« sigma2_eps »)
« sigma2_eps : »
sigma2_eps
abs_ev <- eigen(fit_results_obs$IDEmodel$get(« M »))$values %>%abs()
« complex eigenvalues of the evolution matrix M »
summary(abs_ev)
« Realisation des predictions… »
« Grille complete »
ST_grid <- predict(fit_results_obs$IDEmodel)
ST_grid_df <- ST_grid %>%
data.frame() %>%
mutate(IDE_pred = ST_grid$Ypred,
IDE_predse = ST_grid$Ypredse,
IDE_predZse = sqrt(ST_grid$Ypredse^2 +
sigma2_eps))
ST_grid_df$t<-as.integer(format(ST_grid_df$t, »%d »))-1
coords<-dplyr::select(ST_grid_df,X,Y)
attr<-dplyr::select(ST_grid_df,-X,-Y,-IDE_pred,-IDE_predse)
names(attr)<-c(« t », »Val », »Predict », »ErreurY », »ErreurZ »)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
Resultat_grille<-st_as_sf(spdfN)

« Previsions observations »
pred_IDE_obs <- predict(fit_results_obs$IDEmodel,
newdata = obs)

obs_df <- obs %>%
data.frame() %>%
mutate(IDE_pred = pred_IDE_obs$Ypred,
IDE_predse = pred_IDE_obs$Ypredse,
IDE_predZse = sqrt(IDE_predse^2 +
sigma2_eps))

obs_df$t<-as.integer(format(obs_df$time, »%d »))-1
coords<-dplyr::select(obs_df,X,Y)
attr<-dplyr::select(obs_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
names(attr)<-c(« ID », »t », »Val », »Predict », »ErreurY », »ErreurZ »)
attr$t<-as.integer(attr$t)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
Resultat_observations<-st_as_sf(spdfN)
if (effectuer_la_validation_du_modele) {
if (repertoire_resultats== ») {
repertoire_resultats=(paste(path.expand(‘~’), »\Documents »,sep= » »))
}

« Predictions des donnees de validation »

prédictions pour block pas de temps

pred_IDE_block <- predict(fit_results_obs$IDEmodel,
newdata = block)
block_df <- block %>%
data.frame() %>%
mutate(IDE_pred = pred_IDE_block$Ypred,
IDE_predse = pred_IDE_block$Ypredse,
IDE_predZse = sqrt(IDE_predse^2 +
sigma2_eps))
block_df$t<-as.integer(format(block_df$time, »%d »))-1
coords<-dplyr::select(block_df,X,Y)
attr<-dplyr::select(block_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
names(attr)<-c(« ID », »t », »Val », »Predict », »ErreurY », »ErreurZ »)
attr$t<-as.integer(attr$t)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
sfn<-st_as_sf(spdfN)
nomfic<-paste(repertoire_resultats, »/ »,nom_du_fichier, »_block.shp »,sep= » »)
st_write(sfn, nomfic, driver= »ESRI Shapefile »,append=FALSE)

prédictions pour pourcentage échantillons

pred_IDE_random <- predict(fit_results_obs$IDEmodel,
newdata = random)
alea_df <- random %>%
data.frame() %>%
mutate(IDE_pred = pred_IDE_random$Ypred,
IDE_predse = pred_IDE_random$Ypredse,
IDE_predZse = sqrt(IDE_predse^2 +
sigma2_eps))
alea_df$t<-as.integer(format(alea_df$time, »%d »))-1
coords<-dplyr::select(alea_df,X,Y)
attr<-dplyr::select(alea_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
names(attr)<-c(« ID », »t », »Val », »Predict », »ErreurY », »ErreurZ »)
attr$t<-as.integer(attr$t)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
sfn<-st_as_sf(spdfN)
nomfic<-paste(repertoire_resultats, »/ »,nom_du_fichier, »_alea.shp »,sep= » »)
st_write(sfn, nomfic, driver= »ESRI Shapefile »,append=FALSE)
paste(« Les graphiques sont dans le repertoire »,repertoire_resultats)

setwd(repertoire_resultats)
pdf(« Rplots.pdf »)
c1<-cor(block_df$cnt,block_df$IDE_pred, method= »pearson »)
titre<-paste(« Bloc-Predictions vs observations. Correlation= »,c1 )
p1<-ggplot(block_df,aes(cnt, IDE_pred)) +
ggtitle(titre) + geom_point() +
stat_smooth(method= »lm », se=FALSE)
plot(p1)
c2<-cor(alea_df$cnt,alea_df$IDE_pred, method= »pearson »)
titre<-paste(« Aleatoire-Predictions vs observations. Correlation= »,c2 )
p2<-ggplot(alea_df,aes(cnt, IDE_pred)) +
ggtitle(titre) + geom_point() +
stat_smooth(method= »lm », se=FALSE)
plot(p2)
Bias <- function(x,y) mean(x – y) # x: val. obs.
PCV <- function(x,y) mean((x – y)^2) # y: predictions
SCV <- function(x,y,v) mean((x – y)^2 / v) # v: pred. variances
CRPS <- function(x, y, s) verification::crps(x, cbind(y, s))$CRPS
Diagblock <- block_df %>% summarise(
Bias_bloc = Bias(IDE_pred, cnt),
PCV_bloc = PCV(IDE_pred, cnt),
SCV_bloc = SCV(IDE_pred, cnt, IDE_predZse^2),
CRPS_bloc = CRPS(cnt, IDE_pred, IDE_predZse)
)
Diagrandom <- alea_df %>% summarise(
Bias_alea = Bias(IDE_pred, cnt),
PCV_alea = PCV(IDE_pred, cnt),
SCV_alea = SCV(IDE_pred, cnt, IDE_predZse^2),
CRPS_alea = CRPS(cnt, IDE_pred, IDE_predZse)
)

« Diagnostics de la validation par bloc temporel »
Diagblock
« Diagnostics de la validation par echantillonnage aleatoire »
Diagrandom
nval<-length(block)
pred_IDE_block <- predict(fit_results_obs$IDEmodel,
newdata = block,
covariances = TRUE)
Veps <- diag(rep(sigma2_eps, nval))
L_IDE <- t(chol(pred_IDE_block$Cov + Veps))
IntIDE <- coef(fit_results_obs$IDEmodel)
nsim <- 100
E <- matrix(rnorm(nval*nsim), nval, nsim)
Sims_IDE <- IntIDE + pred_IDE_block$newdata$Ypred + L_IDE %*% E
es1<-es_sample(block$cnt, dat = as.matrix(Sims_IDE))
paste(« Valeur de es_sample bloc= »,es1)
distances <- block %>%
coordinates() %>%
dist() %>%
as.matrix()
weights <- 0.5^distances
vs1<-vs_sample(block$cnt, dat = as.matrix(Sims_IDE),
w = weights, p = 1)
paste(« Valeur de vs_sample bloc= »,vs1)

nval<-length(random)
pred_IDE_random <- predict(fit_results_obs$IDEmodel,
newdata = random,
covariances = TRUE)
Veps <- diag(rep(sigma2_eps, nval))
L_IDE <- t(chol(pred_IDE_random$Cov + Veps))
IntIDE <- coef(fit_results_obs$IDEmodel)
E <- matrix(rnorm(nvalnsim), nval, nsim)
Sims_IDE <- IntIDE + pred_IDE_random$newdata$Ypred + L_IDE %% E
es2<-es_sample(random$cnt, dat = as.matrix(Sims_IDE))

paste(« Valeur de es_sample aléa= »,es2)
distances <- random %>%
coordinates() %>%
dist() %>%
as.matrix()
weights <- 0.5^distances
vs2<-vs_sample(random$cnt, dat = as.matrix(Sims_IDE),
w = weights, p = 1)
paste(« Valeur de vs_sample aléa= »,vs2)
loglikIDE <- -fit_results_obs$IDEmodel$negloglik()
pIDE <- 7
m<-length(obs)
Criteria <- data.frame(AIC = c(0), BIC = c(0),
row.names = c( nom_du_fichier))
aic<- -2*loglikIDE + 2*pIDE
Criteria[nom_du_fichier, « AIC »] <- aic
bic<- -2*loglikIDE + pIDE*log(m)
Criteria[nom_du_fichier, « BIC »] <- bic
Criteria
plot.new()
grid.table(Criteria)
Bias_bloc = Bias(block_df$IDE_pred, block_df$cnt)
PCV_bloc = PCV(block_df$IDE_pred, block_df$cnt)
SCV_bloc = SCV(block_df$IDE_pred, block_df$cnt, block_df$IDE_predZse^2)
CRPS_bloc = CRPS(block_df$cnt, block_df$IDE_pred, block_df$IDE_predZse)
Bias_alea = Bias(alea_df$IDE_pred, alea_df$cnt)
PCV_alea = PCV(alea_df$IDE_pred, alea_df$cnt)
SCV_alea = SCV(alea_df$IDE_pred, alea_df$cnt, alea_df$IDE_predZse^2)
CRPS_alea = CRPS(alea_df$cnt, alea_df$IDE_pred, alea_df$IDE_predZse)
resume <- matrix(c(« Correlation »,c1,c2, »Biais »,Bias_bloc,Bias_alea, »PCV »,PCV_bloc,PCV_alea, »SCV »,SCV_bloc,SCV_alea, »CRPS »,CRPS_bloc,CRPS_alea, »ES »,es1,es2, »VS »,vs1,vs2),ncol=3,byrow=TRUE)
colnames(resume) <- c(« Score », »Bloc t », »Aleatoire »)
resume <- as.table(resume)
resume
plot.new()
grid.table(resume,rows=NULL)
dev.off()
}else{
pred_IDE_obs <- predict(fit_results_obs$IDEmodel,
newdata = obs)
obs_df <- obs %>%
data.frame() %>%
mutate(IDE_pred = pred_IDE_obs$Ypred,
IDE_predse = pred_IDE_obs$Ypredse,
IDE_predZse = sqrt(IDE_predse^2 +
sigma2_eps))
obs_df$t<-as.integer(format(obs_df$time, »%d »))-1
coords<-dplyr::select(obs_df,X,Y)
attr<-dplyr::select(obs_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
names(attr)<-c(« ID », »t », »Val », »Predict », »ErreurY », »ErreurZ »)
attr$t<-as.integer(attr$t)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
sfn<-st_as_sf(spdfN)
}

(Voir les différentes versions )

Les données

Les données sont extraites des Dénombrements du Relevé des oiseaux nicheurs (BBS). Ces données proviennent du Relevé des oiseaux nicheurs de l’Amérique du Nord.
Nous considérons en particulier les dénombrements annuels du troglodyte de Caroline (Thryothorus ludovicianus) pour la période de 1967 à 2014. L’unité d’échantillonnage du BBS est un parcours routier d’environ 39,2 km. Dans chaque unité d’échantillonnage, les observateurs bénévoles font 50 arrêts et comptent les oiseaux pendant une période de 3 minutes lorsqu’ils exécutent leurs itinéraires (habituellement en juin). Il y a plus de 4000 routes dans l’enquête nord-américaine, mais toutes ne sont pas disponibles chaque année. Aux fins des analyses présentées ici, nous considérons que les dénombrements des parcours totaux se produisent chaque année (pendant la saison de reproduction) et définissons l’emplacement spatial de chaque parcours comme étant le centroïde de la route. Ainsi, nous considérons que les données sont discrètes dans le temps, géostatistiques et irrégulières dans l’espace, et non gaussiennes dans le sens qu’elles sont des comptages.

Copiez les fichiers shapefile dans un répertoire de données et chargez la couche dans QGis.

Les scripts

Comme leur nom l’indique, il y a deux types de modèles IDE:

• ceux que l’on considère avec un noyau invariant dans l’espace : le script ide_modele_invariant.rsx
• ceux que l’on considère avec un noyau qui varie dans l’espace: le script ide_modele_variant.rsx

Nous verrons à quoi cela correspond un peu plus bas.

Ces scripts sont à copier, avec leur fichier d’aide, dans le répertoire C:\Users\xxxxx\AppData\Roaming\QGIS\QGIS3\profiles\default\processing\rscripts

Si QGis est configuré pour les scripts R vous trouverez ces deux scripts dans la fenêtre Traitements->R->modele IDE

Les éléments du modèle IDE

Le modèle IDE est défini par

On va essayer de traduire ceci en langage courant:

La valeur du processus qui nous intéresse (Y) à un endroit donné s et à un temps donné t résulte de trois éléments:

• la valeur prise par ce qu’on appellera le noyau du modèle (m(s,x;θp))
• la valeur du processus au moment précédent t-1
• une « erreur » ou « bruit » indépendant du processus ηt

Observations ou valeurs du processus?

Tout d’abord, sachez que certaines notations ont leur importance. Ici on travaille pour obtenir Yt. Dans la littérature vous verrez deux types de notation: Yt ou Zt. Ceci a son importance. Quand on parle de Z on parle de vos données observées. Quand on parle de Y on parle du processus sous-jacent. La différence vient du fait que ce que l’on observe (Z) est toujours une image biaisée du processus sous-jacent (Y), ne serait-ce que par les erreurs de mesure.

Le noyau du modèle

Dans le noyau du modèle m(s, x; θp) la valeur m que l’on va appliquer dépend de trois éléments:

• le plus simple est s, c’est à dire la localisation dans l’espace du point qui nous intéresse
• puis, nous avons θp qui correspond à quatre paramètres qui déterminent la redistribution du poids qui sera affecté à la valeur de Yt-1. Pour simplifier, si on a un processus avec une croissance de 10% entre t-1 et t, θp aura une valeur de 1.1 (ceci est une grosssse simplification)
• le troisième élément est x, la variation spatiale, et il va falloir s’arrêter un peu pour bien le comprendre.

La variation spatiale

Reprenons les poids θp de notre exemple. Si partout dans notre espace de travail la croissance du processus Y est de 10%, on peut dire que notre noyau est spatialement invariant. Mais ceci peut ne pas être le cas. Le pourcentage de différence entre les deux pas de temps peut résulter d’une fonction dépendante de l’endroit où l’on se trouve (s). Dans ce cas, la valeur du noyau dépendra toujours des paramètres θp, mais aussi de la valeur de la fonction x à l’endroit donné (s). On dira alors que notre noyau est variant.

Les paramètres θp

Le noyau du modèle IDE peut être associée à une courbe gaussienne. Le premier des quatre paramètres theta_1 correspond à la hauteur (amplitude) de la courbe en cloche. Le deuxième, theta 2, correspond à l’aperture de la cloche la rendant plus pointue ou plus évasée. Les deux autres paramètres, theta 3 et theta 4 correspondent à un décalage sur la dimension des x et des y, respectivement, en fonction de l’endroit où l’on se trouve.

Le moment est venu de se jeter à l’eau. On va voir deux types de simulation. La première en considérant que le noyau du modèle ne varie pas spatialement, la deuxième en considérant que le noyau dépend de l’endroit où l’on se trouve et donc qu’il varie spatialement.

Exemple de simulation avec un noyau spatialement invariant

Si vous ouvrez dans un éditeur de texte le script ide_modele_invariant.rsx vous verrez ces trois étapes.

Construction du modèle IDE

Tout d’abord on définit les caractéristiques du modèle. Dans notre cas:

IDEmodel <- IDE(f = cnt ~  1  ,
                data = mSTIDF,
                dt = as.difftime(1, units = unit),
                grid_size = grid_size,
                forecast=forecast)

En plus de la « formule » (f=) qui détermine quels sont les colonnes de données à modéliser (un attribut de la table, plus les coordonnées de chaque point), la couche de données (data=) mise à un format R de données spatio-temporelles (STIDF), nous avons à rentrer trois paramètres.

Le pas de temps dt et les unités

Contrairement à ce que l’on peut penser, ce n’est pas si simple d’indiquer le pas de temps de nos données!!!

Les outils de modélisation nécessitent des pas de temps de même longueur. Jusque-là, tout paraît simple. Sauf que si nos données correspondent à des observations mensuelles, le problème est que les mois n’ont pas la même longueur… de 28 à 31 jours. Et si nos données sont annuelles, les années n’ont pas la même longueur (bissextiles ou pas). Ceci fait que les unités de pas de temps acceptables pour les modèles spatio-temporels vont de la seconde à la semaine, pas plus.

On est alors obligés de biaiser un peu pour traiter des données mensuelles ou annuelles.

Le moyen le plus simple pour ne pas avoir à jongler avec nos données est de calculer un numéro d’ordre de nos observations est d’utiliser cet ordre comme paramètre temporel.

Supposons que nous avons de données annuelles entre 2001 et 2019. On calcule un champ avec une valeur de 1 pour 2001, 2 pour 2002, etc. Pour des mois, on fait de même, au lieu d’utiliser la vraie date, on établit une échelle de 1 à X correspondante.

Dans le cas où vous utiliserez le numéro d’ordre, il faut définir les unités avec la valeur « days ».

Il n’est pas possible de prévoir tous les formats possibles de date et d’heure. Si vous rencontrez des problèmes, calculez un numéro d’ordre et définissez les unités comme étant des jours. Ceci ne modifie en rien les résultats.

La grille de modélisation

Nos données sont sous forme de points d’échantillonnage répartis sur une zone d’étude. Pour la mise en place du modèle nous définissons un maillage carré sur lequel on évaluera les paramètres du modèle et ensuite on calculera les valeurs de Y. Le paramètre grid_size détermine le nombre de mailles de la grille. N’oubliez pas que le nombre de cellules déterminera le temps de traitement… et il est déjà très long sans faire appel à des grilles trop grandes.

Les prévisions

Le modèle que nous mettons en place va nous permettre de calculer les prédictions pour chaque cellule de la grille et chaque pas de temps présent dans nos données. Mais il peut nous permettre de prévoir les valeurs pour les prochains pas de temps dans le futur. La valeur du paramètre forecast indique le nombre de pas de temps, à partir du dernier de nos observations, à calculer. Si nos données sont annuelles et vont de 1994 à 2010, forecast= 2 fera que le modèle calcule les données des années 2011 et 2012.

Dans ces scripts nous n’avons pas implémenté, pour ne pas surcharger inutilement l’interface, un autre paramètre possible : hindcast. Celui-ci agît comme forecast mais pour la période antérieure à nos données. Si on ajoute dans le code, après forecast « hindcast=2 », le modèle calculera les données pour les deux années antérieures à nos données, 1992 et 1993.

Si vous souhaitez utiliser l’option hindcast, il suffit de modifier le code du script avec la valeur que vous souhaitez.

IDEmodel <- IDE(f = cnt ~  1  ,
                data = mSTIDF,
                dt = as.difftime(1, units = unit),
                grid_size = grid_size,
                forecast=forecast,hindcast=1)

Ajustement du modèle IDE

Maintenant que nous avons défini quels paramètres on veut utiliser dans notre modèle, on « ajuste » le modèle à partir des observations. Ceci est fait avec la commande fit.IDE

fit_results <- fit.IDE(IDEmodel,method="DEoptim",itermax=itermax)

Le modèle IDE calculé de l’objet, ici l’objet fit_results, contient les estimations initiales des paramètres, ainsi que les prédictions pour tous les pas de temps des estimations des paramètres. Les paramètres dans ce cas sont la variance de l’erreur de mesure, la variance de la perturbation aléatoire ηt(s)(dont la structure de covariance est fixe), les paramètres du noyau, et les coefficients de régression β.

Nous parlons bien des valeurs des différents paramètres du modèle pour chaque pas de temps et non des valeurs du processus à chaque maille de la grille et chaque pas de temps.

Ce calcul-là sera fait dans l’étape suivante.

Dans cette étape, le seul paramètre sur lequel vous pouvez intervenir est le nombre d’itérations pour l’ajustement du modèle. Il n’y a pas de recette, tout dépend de la complexité du modèle et des observations disponibles.

Par contre, il faut savoir que le temps de calcul de cette étape est très,très long. Avec un ordinateur plutôt puissant, pour les données d’exemple fournies,le temps de calcul pour une centaine d’itérations tourne autour d’une heure ou deux…

Les résultats de chaque itération sont affichés dans l’onglet « Journal » de l’interface et dans la sortie console du script. Vous pourrez voir là à partir de quelle itération on a atteint les valeurs finales, et si vous devez ré-exécuter le calcul, vous pouvez limiter le nombre d’itérations à ce chiffre-là.

Dans le journal vous aurez aussi le détail des différents calculs des paramètres.

Calcul des prédictions sur la maille définie

La dernière étape consiste à utiliser l’ajustement réalisé et prédire pour chaque nœud de la maille spatiale définie dans la première étape et pour chaque pas de temps souhaités la valeur de ,notre processus Y.

Ceci est fait avec la fonction « predict »

ST_grid <- predict(fit_results$IDEmodel)

Les lignes suivantes de code, après celle-ci, transforment le format grid de R en une couche de points au format défini par l’utilisateur dans le champ « output ».

Données exemple

Passons à la partie pratique, même s’il n’y a pas beaucoup de manipulations et beaucoup de patience à avoir.

Chargez la couche MOcarolinawren dans QGis.

Dans la boîte de traitements->R->modele IDE double cliquez sur ide modele invariant

Dans le champ t sélectionnez l’attribut t qui correspond à un numéro d’ordre de l’année. Ses valeurs sont comprises entre 1 et 21.

Dans le champ val sélectionnez l’attribut cnt qui correspond aux comptages d’oiseaux.

Dans le champ grid_size, laissez la valeur par défaut (50). On aura alors une grille en sortie de 50 x 50 points.

Dans le champ units, comme nous l’avons indiqué plus haut, on ne peut pas définir le pas de temps en années. On dit alors qu’il s’agît de jours (days) car le fait que tous les quatre ans il y ait un jour de différence (1/365ème de plus dans la longueur du pas) est sûrement négligeable.

Dans le champ forecast, laissez la valeur par défaut (1). Dans la grille de sortie on aura alors un pas de temps supplémentaire, avec la valeur t à 22 qui correspondra aux prédictions pour l’année suivante à la dernière dans les observations.

Si vous êtes impatient, réduisez le nombre d’itérations par défaut. Si non, armez vous de patience et laissez ce chiffre.

Définissez la couche en sortie, au format que vous préférez (shp, gpkg,…) et un fichier pour les sorties de la console. Ce dernier est optionnel mais il contient tous les paramètres de modélisation.

Cliquez sur Exécuter… bon courage.

Résultats de l’exemple

Quatre heures après, nous obtenons une couche qui s’affiche dans QGis:

Elle contient notre grille de prédictions de taille 30×50.

Avant de voir plus en détail les résultats, voyons quelles informations sont fournies par la sortie console du script. Vous le trouvez dans l’onglet « Journal » de la fenêtre du script, en tant que lien dans le panneau « Visualiseur de résultats » et, si vous avez défini un fichier permanent dans le champ « R Console Output [optionnel]« , dans ce fichier.

[1] "calcul du modele IDE..."
[1] "Ajustement du modele..."
Iteration: 1 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397
Iteration: 2 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397
Iteration: 3 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397
		.
		.
		.
Iteration: 99 bestvalit: 1800.983447 bestmemit: 0.052054 -0.350126 -0.459083 2.066409 3.177588 -0.728638
Iteration: 100 bestvalit: 1800.983447 bestmemit: 0.052054 -0.350126 -0.459083 2.066409 3.177588 -0.728638
[1] "kernel parameters: "
par1 par2 par3 par4
52.05424217 0.03015948 -0.45908287 2.06640855
[1] "regression coefficients: "
Intercept X Y
70.3594335 -0.3226177 -2.4435913
[1] "sigma2_eps :"
par5
23.98881
[1] "complex eigenvalues of the evolution matrix M"
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00000 0.00000 0.00301 0.18171 0.34210 0.98394
[1] "Realisation des predictions..."
Warning message:
In abbreviate_shapefile_names(obj) :
Field names abbreviated for ESRI Shapefile driver

Paramètres descriptifs

Les sorties kernel parameters, regression coefficients et sigma2_eps indiquent les valeurs obtenues pour l’ajustement du modèle. Ils n’ont pas de signification par rapport à la « justesse » de l’ajustement, sauf à les connaître dans la réalité. Ceci n’est pas généralement le cas, mais on peut utiliser IDE pour tester des simulations de modèles pour lesquelles nous fixons ces paramètres. On peut alors comparer l’ajustement des paramètres par rapport aux valeurs fournies aux simulations.

Valeurs propres de la matrice M

Par contre, il y a une sortie qui nous permet de vérifier que le résultat peut être pris en compte et passer à l’étape suivante de toute modélisation : la validation des résultats.

Comme nous l’avons vu plus haut, la formule de l’équation intégro-différentielle est

Mais, une fois le modèle calculé et ajusté, on peut l’écrire sous la forme:

où M est une matrice qu’on peut appeler « matrice de transition », de dimensions n x n (dans notre cas 50 x 50), avec les valeurs pour chaque cellule de la grille qui permettent de passer de la valeur de Y_t-1 à Y_t.

On peut donc calculer les valeurs propres (eigenvalues) de cette matrice. Ceci permet de savoir si le modèle calculé est stable ou pas.

Dans notre sortie, la synthèse de ces valeurs propres apparaît sous le titre « complex eigenvalues of the evolution matrix M« .

Nous notons que si la valeur maximale des valeurs propres est >=1 , alors le modèle, est instable, et Yt grandira sans limite lorsque t augmente. Inversement, si le module maximal de toutes les valeurs propres est inférieur à 1, alors le modèle est stable.

Dans notre exemple, la valeur maximale des valeurs propres est de 0.98394. Le modèle peut donc être considéré comme stable.

La couche résultat

A la fin du traitement vous aurez une couche de points contenant la grille de prédictions:

L’attribut t correspond au pas de temps. Vous trouverez les 21 pas originaux plus le pas 22 qui correspond aux prédictions du prochain pas (forecast).

Ypred correspond à la valeur calculée par le modèle, Ypredse à la valeur de l’erreur standard dans le domaine du processus (Y) pour le point et le pas de temps. YpredZ correspond à l’erreur standard dans le domaine des observations (Z) pour le point et le pas de temps.

En appliquant le même filtre sur les données en entrée et la grille en sortie sur l’attribut t, vous pourrez voir la correspondance:

Dans le prochain article nous verrons l’utilisation du deuxième script, avec un noyau variant spatialement.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois