Dans l’article précédent, nous avons considéré le cas d’un noyau spatialement invariant, c’est-à-dire le cas où le noyau m(s,r;θp) n’est qu’une fonction de r s . Dans cet article, nous considérons le cas où un ou plusieurs des paramètres θ sont variables spatialement. De tels modèles sont nécessaires lorsque le processus spatio-temporel présente, par exemple, une dérive qui varie selon l’espace.
Mais il faut des données…
Pour modéliser de telles données, nous avons besoin d’un grand nombre de points temporels,disons au moins 15. Cela est important, car il est difficile d’obtenir des estimations raisonnables des paramètres spatialement répartis, à moins que les données couvrent une grande partie du domaine spatial pendant au moins quelques points de temps consécutifs.
Les noyaux variables spatialement peuvent être introduits en spécifiant l’argument kernel_basis dans l’appel à IDE. Les fonctions de base utilisées par IDE sont de la même classe que celles utilisées par le package FRK (fixed rank krigging).
Rappel :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.
C’est quoi des fonctions de base?
Une image valant plus qu’un long discours, voyons un exemple:
On a un phénomène représenté par la courbe rouge de l’image de droite. On peut essayer de trouver une seule fonction qui nous donne le bon y pour chaque x de la courbe, mais vous présentez que ce n’est pas simple.
Une autre manière d’aborder et de résoudre le problème, c’est de prendre une série de fonctions, toutes égales, mais qui sont décalées sur l’axe des x.
Si vous observez la figure de gauche, chaque courbe, à sa base couvre une vingtaine de valeurs de x.
On va avoir donc la première courbe, représentée par une fonction (dans cette image ce sont des fonctions bicarrées) f1(x), puis une deuxième f2(x), etc.
Il nous suffit de trouver un facteur ∝1 , ∝2, … à multiplier pour chaque fonction et à les additionner pour construire notre courbe rouge :
Y(x)= ∝1.f1(x) + ∝2.f2(x) + . . . + ∝n.fn(x)
On doit donc faire trois choix pour ajuster nos fonctions de base à notre processus:
- le type de fonction
- la taille de la base à utiliser que nous appellerons l’ouverture (aperture)
- le nombre de résolutions. On verra plus loin à quoi cela correspond.
Les types de fonction
Parmi les nombreux types de fonctions existantes, le package FRK de R, qui est celui qui est utilisé pour la modélisation IDE, implémente quatre types de fonctions:
Les fonction bicarrées
Une petite formule pour définir ce type de fonction:
Et une image pour mieux comprendre:
Dans la formule, R représente l’ouverture de la fonction.
Les fonctions Gaussiennes
Ces fonctions correspondent à la formule:
Et son image correspondante:
Dans la formule, σ, la déviation standard, représente l’ouverture de la fonction.
Les fonctions exp (exponentielles)
Ces fonctions correspondent à la formule:
Et son image correspondante:
Dans la formule, τ, la valeur d’inflexion, représente l’ouverture de la fonction.
Les fonctions Matern32
Ces fonctions correspondent à la formule:
Et son image correspondante:
Dans la formule, κ, la valeur d’échelle, représente l’ouverture de la fonction.
Ouverture des fonctions
Les images des fonctions précédentes permettent de voir le profil de celles-ci. Au niveau spatial (x/y) elles correspondent à des cercles ou des ellipses.
On doit s’assurer que les bases des fonctions couvrent complètement la zone d’étude et se superposent suffisamment pour que au moins cette surface corresponde à la moitié de la fonction adjacente. Si non, il y a des valeurs du processus qui ne pourraient pas être représentées, même en affectant un facteur multiplicatif plus fort.
Une fonction de R permet de tracer les emprises des fonctions de base retenues:
Cette image montre les emprises de 12 fonctions de base, réparties régulièrement. Notez que la zone d’étude est plus petite que le cadre de l’image.
Par défaut, l’aperture des fonctions est de 1. Selon le type de fonction, on peut ne pas retrouver une répartition telle que sur cette figure.
Si vous choisissez les fonctions bicarrées, vous n’aurez rien à faire. Par contre, pour les autres familles, il faudra ajuster l’aperture.
En changeant seulement le type de famille de bisquare à Gaussian, on obtient la distribution suivante:
Il faut augmenter l’aperture (scale_aperture) à 3 pour avoir la distribution suivante:
Le nombre de résolutions
Dans les images précédentes, on suppose qu’il y a un processus à ajuster sur l’ensemble de la zone d’étude. Ce processus est à l’origine des variations des paramètres du noyau.
Mais il se peut qu’on ait des raisons de penser qu’il n’y a pas qu’un seul processus mais, par exemple deux, mais d’échelle différente. On aurait un processus fin, et un processus plus global. Dans ce cas nous pouvons définir deux niveaux de résolution.
Par exemple:
L’ajustement des paramètres sera beaucoup plus fin, mais n’oubliez pas que tout ceci se paye au niveau du temps de traitement. Évitez la surenchère qu niveau des paramètres de modélisation, car si une demie douzaine de fonctions suffisent pour estimer les bons paramètres, en ajouter une centaine ne donnera pas un meilleur résultat, mais vous aurez le temps de partir en weekend avant d’obtenir le résultat du calcul.
Le script ide_auto_basis
Le script ide_auto_basis vous permet de tester les différentes valeurs, en produisant les images ci-dessus.
Si vous ne l’avez pas téléchargé à partir de l’article précédent, voici le lien pour le faire:
(Voir les différences des versions)
Si non, vous avez ici le code du script
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
basis=selection bisquare;Gaussian;exp;Matern32
nresolution= number 1
scale_aperture=number 1
amplitude=selection constant_basis();auto_basis
aperture=selection constant_basis();auto_basis
xshift=selection constant_basis();auto_basis
yshift=selection constant_basis();auto_basis
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 »)
library(« FRK »)
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 (basis==0) {base= »bisquare »}
if (basis==1) {base= »Gaussian »}
if (basis==2) {base= »exp »}
if (basis==3) {base= »Matern32″}
mbasis_1<-auto_basis(manifold=plane(),
data= mSTIDF,
nres=nresolution,
scale_aperture=scale_aperture,
type=base)
« mbasis »
mbasis_1
« kernel basis »
if (amplitude==0) {
kernel_basis<-list(thetam1=constant_basis())
}else {
kernel_basis<-list(thetam1=mbasis_1)}
if (aperture==0 ) {
kernel_basis[2]<-(thetam2=constant_basis())
}else{
kernel_basis[2]<-(thetam2=mbasis_1)}
if (xshift==0 ) {
kernel_basis[3]<-(thetam3= constant_basis())
}else {
kernel_basis[3]<-(thetam3= mbasis_1)}
if (yshift==0) {
kernel_basis[4]<-(thetam4= constant_basis())
}else {
kernel_basis[4]<-(thetam4= mbasis_1)}
kernel_basis
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,
kernel_basis= kernel_basis)
« 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 »)
p1<-show_kernel(fit_results_obs$IDEmodel,scale=0.2)
plot(p1)
c1<-cor(block_df$cnt,block_df$IDE_pred, method= »pearson »)
titre<-paste(« Bloc-Predictions vs observations. Correlation= »,c1 )
p2<-ggplot(block_df,aes(cnt, IDE_pred)) +
ggtitle(titre) + geom_point() +
stat_smooth(method= »lm », se=FALSE)
plot(p2)
c2<-cor(alea_df$cnt,alea_df$IDE_pred, method= »pearson »)
titre<-paste(« Aleatoire-Predictions vs observations. Correlation= »,c2 )
p3<-ggplot(alea_df,aes(cnt, IDE_pred)) +
ggtitle(titre) + geom_point() +
stat_smooth(method= »lm », se=FALSE)
plot(p3)
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<- -2loglikIDE + 2pIDE
bic<- -2loglikIDE + pIDElog(m)
Criteria[nom_du_fichier, « AIC »] <- aic
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()
}
Quand vous l’exécutez, vous obtenez une sortie avec le dessin des emprises de la fonction de base, que vous pouvez charger dans votre navigateur :
Une fois trouvés les bons paramètres, vous pourrez utiliser le script ide_modele_variant.
Le script ide_modele_variant
Le lien de téléchargement est le même que pour le script précédent.
Le code du script est le suivant:
##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
##basis=selection bisquare;Gaussian;exp;Matern32
##nresolution= number 1
##scale_aperture=number 1
##amplitude=selection constant_basis();auto_basis
##aperture=selection constant_basis();auto_basis
##xshift=selection constant_basis();auto_basis
##yshift=selection constant_basis();auto_basis
##output= output vector
library(IDE)
library(dplyr)
library(« sp »)
library(« spacetime »)
library(« sf »)
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)
data$t<-as.Date(data$t,origin= »0000-01-01″)
mSTIDF <- stConstruct(x = data, # data set
space = c(« X », « Y »), # spatial fields
time = t) # time field
if (basis==0) {base= »bisquare »}
if (basis==1) {base= »Gaussian »}
if (basis==2) {base= »exp »}
if (basis==3) {base= »Matern32″}
mbasis_1<-auto_basis(manifold=plane(),
data= mSTIDF,
nres=nresolution,
scale_aperture=scale_aperture,
type=base)
> »kernel basis »
if (amplitude==0) {
kernel_basis<-list(thetam1=constant_basis())
}else {
kernel_basis<-list(thetam1=mbasis_1)}
if (aperture==0 ) {
kernel_basis[2]<-(thetam2=constant_basis())
}else{
kernel_basis[2]<-(thetam2=mbasis_1)}
if (xshift==0 ) {
kernel_basis[3]<-(thetam3= constant_basis())
}else {
kernel_basis[3]<-(thetam3= mbasis_1)}
if (yshift==0) {
kernel_basis[4]<-(thetam4= constant_basis())
}else {
kernel_basis[4]<-(thetam4= mbasis_1)}
>kernel_basis
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 + X + Y ,
data = mSTIDF,
dt = as.difftime(1, units = unit),
grid_size = grid_size,
forecast=forecast,kernel_basis=kernel_basis)
> »Ajustement du modele… »
fit_results <- fit.IDE(IDEmodel,method= »DEoptim »,itermax=itermax)
show_kernel(fit_results$IDEmodel)
kp<-fit_results$IDEmodel$get(« k ») %>% unlist()
> »kernel parameters: «
>kp
co<-coef(fit_results$IDEmodel)
> »regression coefficients: «
>coef(fit_results$IDEmodel)
sg<-fit_results$IDEmodel$get(« sigma2_eps »)
> »sigma2_eps : »
>sg
abs_ev <- eigen(fit_results$IDEmodel$get(« M »))$values %>%abs()
> »complex eigenvalues of the evolution matrix M »
> summary(abs_ev)
> »Realisation des predictions… »
ST_grid <- predict(fit_results$IDEmodel)
ST_grid_df <- ST_grid %>%
data.frame() %>%
mutate(IDE_predZse = sqrt(IDE_predse^2 +
fit_results$IDEmodel$get(« 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)
spdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
output<-st_as_sf(spdfN)
La première partie des paramètres est la même que celle du modèle invariant spatialement qui ont été décrits dans l’article précédent.
Les paramètres basis, nresolution et scale_aperture sont ceux obtenus avec le script ide_auto_basis.
Voyons maintenant les paramètres de définition de la variation spatiale du noyau.
Paramètres theta
Maintenant, rappelons-nous que θp1 corresponds à l’amplitude du noyau, θp2 à l’échelle ou ouverture, θp3 à la « dérive » de décalage directionnel X et θp4 au décalage de direction Y (dérive).
L’amplitude du noyau variant spatialement est donnée par θp1(s) et contrôle la stationnarité temporelle, le paramètre échelle de longueur (variance) variant spatialement θp2(s) correspond à un paramètre échelle du noyau (ouverture) (c.-à-d. que la largeur du noyau augmente quand θp2 augmente) et la moyenne des paramètres de décalage θp3(s) et θp4(s) correspondent à un décalage spatialement variable du noyau par rapport aux emplacements.
Pour chacun de ces quatre paramètres on peut choisir individuellement de les considérer comme invariant spatialement ou bien variant spatialement.
Dans le premier cas, il faut choisir l’option constant_basis. Si on souhaite le considérer comme variant spatialement et l’ajuster avec les froncions de base définies plus haut, il faut sélectionner auto_basis.
Il n’est pas facile d’expliquer le cas d’un noyau variant empaillement, tellement ceci est lié aux types de données et au type de processus sous-jacent.
Il faut bien comprendre qu’il ne s’agit pas d’ajuster la distribution spatiale comme on le ferait pour un krigeage. On ne considère pas la distribution des valeurs des observations à l’intérieur d’un pas de temps.
On considère ce qui arrive entre le moment t-1 et le moment t. Supposons qu’il s’agisse d’une population de poissons dans une rivière. Après la reproduction les individus auront une tendance à descendre le courant plutôt qu’à le remonter. Il y aura alors une dérive spatiale du paramètre θp3(s), le valeur du noyau correspondant à un grand nombre d’individus au moment t-1 sera légèrement en amont du même noyau au moment t.
Nous verrons dans l’exemple qui suit les sorties graphiques associées à ceux paramètres.
Données exemple
Chargez la couche MOcarolinawren dans QGis.
Dans la boîte de traitements->R->modele IDE double cliquez sur ide modele variant
Dans notre exemple nous demandons deux résolutions et nous considérons qu’il y a une dérive du noyau sur les deux axes, x et y.
La différence entre une et deux résolutions est que, dans le premier cas on cherchera à ajuster 12 fonctions par paramètre theta variable et que dans le deuxième cas on aura 120 fonctions par paramètre… avec ce que cela implique en temps de calcul.
Valeur de theta 3 avec un nombre de résolutions de 1
[[3]]
Number of basis functions: 12
Number of resolutions: 1
Type of manifold: plane
Dimension of manifold: 2
First basis function:
function (s) { stopifnot(ncol(s) == dimensions(manifold)) y <- distance(manifold, s, c) (1 - (y/R)^2)^2 * (y < R) }
Valeur de theta 3 avec un nombre de résolutions de 2
[[3]]
Number of basis functions: 120
Number of resolutions: 2
Type of manifold: plane
Dimension of manifold: 2
First basis function:
function (s) { stopifnot(ncol(s) == dimensions(manifold)) y <- distance(manifold, s, c) (1 - (y/R)^2)^2 * (y < R) }Sortie dans le fichier R Plots
Que vous choisissiez temporaire ou permanent, vous verrez un tracé dans le fichier R Plots.
Vecteurs de variation du noyau
Dans le tracé on peut visualiser l’advection qui a été générée par les fonctions de base choisies. La longueur des vecteurs est proportionnelle à la variation induite.
Rappelez-vous que vous pouvez utiliser le script ide_auto_basis pour avoir les emprises des fonctions de base.
Couche résultante
Dans la couche en sortie on retrouve deux attributs, Ypred et Ypredse.
Le premier correspond aux prédictions de valeurs du processus Y, le deuxième à la déviation standard (l’erreur).
Pour pouvoir travailler sur cette sortie, il est nécessaire de filtrer les couches pour ne retenir qu’un pas de temps. Dans nos données exemple nous avons 21 pas de temps, et dans notre sortie nous en avons 22 (les 21 des données d’origine plus un pas de temps futur).
L’image suivante montre le résultat de Ypred pour le pas de temps 21, le dernier de nos données d’origine.
Et l’image qui suit la carte des erreurs (Ypredse).
A titre de comparaison, le traitement avec le script invariant, nous a donné pour le même pas de temps, comme valeurs de Ypred
Vous remarquerez que, nos seulement la distribution n’est pas la même, mais que les valeurs maximales ici sont autour de 28, tandis qu’avec le modèle variant ont observe des valeurs de 37.
Au niveau de la carte d’erreur, le résultat du modèle invariant est:
Il est plus difficile de comparer les deux sorties d’erreur, le modèle invariant ayant une amplitude beaucoup plus grande.
En tout état de cause, en ce qui concerne le sujet de cet article, nous sommes arrivés à la fin. Vous avez tous les éléments pour utiliser les scripts et obtenir les résultats du modèle. Mais en réalité, c’est maintenant que le gros du travail commence.
Lequel de ces deux modèles faut-il retenir? Et si vous essayer un modèle variant avec d’autres paramètres theta, ou si vous décidez d’utiliser une autre famille de fonctions de base, vous serez confrontés à d’autres résultats de modèle, différents mais à première vue plausibles.
Le travail qui commence ici est tout simplement la validation des différents modèles et la comparaison de ces modèles.
Ce n’est pas le sujet de cet article, mais celui d’un prochain.
