{"id":9240,"date":"2021-02-16T10:04:41","date_gmt":"2021-02-16T09:04:41","guid":{"rendered":"https:\/\/www.sigterritoires.fr\/?p=9240"},"modified":"2021-06-07T18:33:52","modified_gmt":"2021-06-07T16:33:52","slug":"modelisation-spatio-temporelle-avec-une-equation-integro-differentielle-ide-avec-qgis-et-r","status":"publish","type":"post","link":"https:\/\/www.sigterritoires.fr\/index.php\/modelisation-spatio-temporelle-avec-une-equation-integro-differentielle-ide-avec-qgis-et-r\/","title":{"rendered":"Mod\u00e9lisation spatio-temporelle avec une \u00e9quation int\u00e9gro-diff\u00e9rentielle (IDE) avec QGis et R"},"content":{"rendered":"\n

En \u00e9cologie on se retrouve, tr\u00e8s souvent, face<\/a> \u00e0 des processus qui sont li\u00e9s \u00e0 l’espace mais aussi qui \u00e9voluent dans le temps. Nous sommes maintenant habitu\u00e9s \u00e0 l’analyse spatiale qui est devenu \u00e0 la port\u00e9e du plus grand nombre gr\u00e2ce aux SIG<\/a>, tels QGis. Nous sommes beaucoup moins habitu\u00e9s \u00e0 l’analyse temporelle qui est moins prise en compte par les SIG.<\/p>\n\n\n\n

Analyse spatio-temporelle<\/h2>\n\n\n\n
Que dire alors de l’analyse spatio-temporelle, c’est \u00e0 dire de l’analyse simultan\u00e9e 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\u00e9resse est l’analyse exploratoire des donn\u00e9es. Un deuxi\u00e8me niveau nous am\u00e8ne \u00e0 essayer de comprendre les param\u00e8tres qui contr\u00f4lent le processus pour remplir les trous dans les donn\u00e9es et r\u00e9duire le \u00ab\u00a0bruit\u00a0\u00bb de nos observations. Enfin, au plus haut du traitement, nous pouvons mod\u00e9liser statistiquement le processus et obtenir un outil qui nous permette une vision compl\u00e9t\u00e9e de nos donn\u00e9es (boucher les trous dans l’espace et le temps de nos observations), de trier les param\u00e8tres du processus et de quantifier leur importance, et finalement, de pr\u00e9voir l’\u00e9tat futur de notre processus.<\/p>\n\n\n\n\n\n\n\n
Si nous prenons seulement l’aspect spatial des donn\u00e9es on comprendra mieux le sujet. Nous disposons, par exemple, de donn\u00e9es sur une population d’oiseaux dont le territoire progresse d’ann\u00e9e en ann\u00e9e. D’un point de vue spatial nous aurons une distribution du nombre d’individus. Le premier pas pour traiter ces donn\u00e9es sera l’analyse exploratoire des donn\u00e9es (histogrammes, qqplots, polygones de Vorono\u00ef,…). On aura alors une id\u00e9e de l’ad\u00e9quation de notre \u00e9chantillonnage et le processus sous-jacent qui nous int\u00e9resse. <\/p>\n\n\n\n
Le deuxi\u00e8me pas sera de passer de nos \u00e9chantillons ponctuels \u00e0 des donn\u00e9es surfaciques, soit par des interpolations simples (IDW, Spline,…) soit par des m\u00e9thodes statistiques (Krigeage,..). A ce point on disposera d’une image de notre processus, nettoy\u00e9e des erreurs de mesure et du \u00ab\u00a0bruit\u00a0\u00bb des observations.<\/p>\n\n\n\n
Mod\u00e9lisation des donn\u00e9es<\/h2>\n\n\n\n
On pourra alors utiliser des m\u00e9thodes de mod\u00e9lisation telles que les mod\u00e8les lin\u00e9aires, les mod\u00e8les lin\u00e9aires g\u00e9n\u00e9ralis\u00e9s (GLM), les mod\u00e8les lin\u00e9aires mixtes g\u00e9n\u00e9ralis\u00e9s (GLMM),les mod\u00e8les additifs g\u00e9n\u00e9ralis\u00e9s (GAM), les mod\u00e8les mixtes additifs g\u00e9n\u00e9ralis\u00e9s (GAMM)… qui nous permettent d’inclure nos connaissances scientifiques sur le processus et d’obtenir des pr\u00e9dictions moins biais\u00e9es. <\/p>\n\n\n\n
Les outils que nous aborderons dans cet article se situent un cran plus loin que ces mod\u00e8les (GLM, GAM,…) pour traiter des donn\u00e9es spatio-temporelles. En effet, les diff\u00e9rents mod\u00e8les lin\u00e9aires correspondent \u00e0 ce que l’on appelle l’approche descriptive de la mod\u00e9lisation spatio-temporelle. Pour faire simple, on d\u00e9finit la corr\u00e9lation entre deux valeurs dans l’espace et le temps. <\/p>\n\n\n\n
La mod\u00e9lisation IDE correspond \u00e0 ce que l’on appelle l’approche dynamique. Le mod\u00e8le consid\u00e8re que la valeur actuelle \u00e0 un endroit donn\u00e9 est \u00e9gale \u00e0 un \u00ab\u00a0facteur de propagation\u00a0\u00bb de la valeur \u00e0 ce point au moment pr\u00e9c\u00e9dent \u00e0 laquelle on ajoute une \u00ab\u00a0erreur d’innovation\u00a0\u00bb.<\/p>\n\n\n\n
Bref, on est en plein dans un domaine qu’on ne peut aborder \u00e0 la l\u00e9g\u00e8re. Il faut bien sacrifier du temps \u00e0 comprendre ce que l’on fait, pourquoi et comment on le fait. Dans cet article j’\u00e9viterai les explications math\u00e9matiques. Vous trouverez les liens qui renvoient vers ces explications. Pour rester dans le domaine habituel des SIG je ferai des simplifications pas toujours exactes \u00e0 100% du point de vue math\u00e9matique, je m’en excuse d’avance au pr\u00e8s des math\u00e9maticiens. <\/p>\n\n\n\n
Les mod\u00e8les d\u2019\u00e9quation int\u00e9gro-diff\u00e9rentielle (IDE) d\u00e9crivent la d\u00e9pendance conditionnelle entre le processus spatial \u00e0 un moment futur et le processus au moment pr\u00e9sent par l\u2019interm\u00e9diaire d\u2019un op\u00e9rateur int\u00e9gral. <\/p>\n\n\n\n
La non lin\u00e9arit\u00e9 ou la d\u00e9pendance temporelle dans la dynamique du mod\u00e8le est souvent int\u00e9gr\u00e9e en permettant aux param\u00e8tres de varier temporellement, ou en ajustant un mod\u00e8le avec un op\u00e9rateur lin\u00e9aire invariant temporellement dans une fen\u00eatre coulissante. Les deux proc\u00e9dures ont tendance \u00e0 \u00eatre excellentes \u00e0 des fins de pr\u00e9vision sur de courts horizons temporels, mais elles prennent g\u00e9n\u00e9ralement beaucoup de temps de calcul.<\/p>\n\n\n\n
En termes moins techniques, les mod\u00e8les IDE sont des outils tr\u00e8s puissant pour pr\u00e9voir l’\u00e9tat suivant d’un syst\u00e8me spatio-temporel \u00e0 partir des donn\u00e9es des \u00e9tats pr\u00e9c\u00e9dents et de l’\u00e9tat actuel.<\/p>\n\n\n\n
Pour savoir ce qu’est IDE<\/h2>\n\n\n\n
Pour comprendre les mod\u00e8les 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\u00e9 par Chapman & Hall\/CRC et disponible en t\u00e9l\u00e9chargement ici<\/a>.<\/p>\n\n\n\n
La documentation de la biblioth\u00e8que R d’IDE est disponible ici.<\/a><\/p>\n\n\n\n
Les traitements QGis et R<\/h2>\n\n\n\n
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\u00e9mentent IDE.<\/p>\n\n\n\n
Par contre, un package R est d\u00e9di\u00e9 \u00e0 la mise en place de mod\u00e8les IDE. Si vous \u00eates un habitu\u00e9 de R, vous pouvez r\u00e9cup\u00e9rer les donn\u00e9es spatio-temporelles de QGis et les traiter dans RStudio.<\/p>\n\n\n\n
Nous avons vu comment migrer les donn\u00e9es de QGis vers RStudio dans cet article.<\/a><\/p>\n\n\n\n
Pour faciliter l’acc\u00e8s aux mod\u00e8les IDE directement dans QGis, nous avons d\u00e9velopp\u00e9 deux scripts R directement utilisables \u00e0 partir des traitements QGis.<\/p>\n\n\n\n
Pr\u00e9alablement, vous devez vous assurer que les traitements QGis 3.x sont bien configur\u00e9s pour l’utilisation de scripts R.<\/p>\n\n\n\n
Vous trouverez comment faire dans cet article.<\/a><\/p>\n\n\n\n
T\u00e9l\u00e9chargements des scripts R et des donn\u00e9es d’exemple<\/h1>\n\n\n\n
Cliquez sur le bouton \u00ab\u00a0T\u00e9l\u00e9charger\u00a0\u00bb pour r\u00e9cup\u00e9rer un fichier zipp\u00e9, contenant deux scripts R et un fichier de donn\u00e9es qui vous permettra de suivre la d\u00e9marche de cet article.<\/p>\n\n\n\n
Si non, vous avez ici le code du script<\/p>\n\n\n
<\/div>
Code du script ide_modele_invariant<\/div>
<\/div><\/div>
<\/p>\n
modele IDE=group<\/h2>\n
Layer=vector<\/h2>\n
t=Field Layer<\/h2>\n
val=Field Layer<\/h2>\n
grid_size=number 50<\/h2>\n
units=selection secs;mins;hours;days;weeks<\/h2>\n
forecast=number 1<\/h2>\n
itermax=number 100<\/h2>\n
effectuer_la_validation_du_modele=boolean TRUE<\/h2>\n
pas_de_temps_de_validation=number<\/h2>\n
pourcentage_points_aleatoires_de_validation=number 10<\/h2>\n
repertoire_resultats=optional folder<\/h2>\n
nom_du_fichier= optional string<\/h2>\n
Resultat_grille= output vector<\/h2>\n
Resultat_observations= output vector<\/h2>\n
library(IDE)
\nlibrary(dplyr)
\nlibrary(\u00ab\u00a0sp\u00a0\u00bb)
\nlibrary(\u00ab\u00a0spacetime\u00a0\u00bb)
\nlibrary(\u00ab\u00a0sf\u00a0\u00bb)
\nlibrary(\u00ab\u00a0verification\u00a0\u00bb)
\nlibrary(\u00ab\u00a0scoringRules\u00a0\u00bb)
\nlibrary(\u00ab\u00a0tidyr\u00a0\u00bb)
\nlibrary(\u00ab\u00a0ggplot2\u00a0\u00bb)
\nlibrary(\u00ab\u00a0gridExtra\u00a0\u00bb)
\ndf<-as.data.frame(Layer)
\ncoords<-st_coordinates(Layer)
\ncrs<-st_crs(Layer)$proj4string
\ndfc<-cbind(df,coords)
\ndata<-dplyr::select(dfc,t,all_of(val),X,Y)<\/p>\n
names(data)<-c(\u00ab\u00a0t\u00a0\u00bb,\u00a0\u00bbcnt\u00a0\u00bb,\u00a0\u00bbX\u00a0\u00bb,\u00a0\u00bbY\u00a0\u00bb)
\ndata$cnt<-as.double(data$cnt)
\ndata$X<-as.double(data$X)
\ndata$Y<-as.double(data$Y)
\nis.date <- function(x) inherits(x, ‘Date’)
\nis.convertible.to.date <- function(x) !is.na(as.Date(as.character(x),origin=\u00a0\u00bb0000-01-01″))
\nif (!sapply(data$t[1], is.date)&& !is.character(data$t[1])) {
\n if (near(data$t[1], as.integer(data$t[1]))){
\n data$t<-as.Date(data$t,origin=\u00a0\u00bb0000-01-01″)
\n trans=1
\n }
\n}
\nif (!sapply(data$t[1], is.date)&& is.character(data$t[1])) {
\n if (sapply(data$t[1], is.convertible.to.date)){
\n data$t<-as.Date(data$t)
\n trans=2
\n }
\n}
\nif (effectuer_la_validation_du_modele) {
\ntv<-pas_de_temps_de_validation
\nt0<-min(data$t)
\ntmax<-max(data$t)
\nif (tv<=t0) {
\n tv<-t0+1
\n }
\nif (tv>tmax) {
\ntv<-tmax-1
\n}
\nif (tv==tmax && forecast==0){
\nforecast=1
\n}<\/p>\n
paste(\u00ab\u00a0pas de temps retenu pour validation=\u00a0\u00bb,tv)
\n}
\nmSTIDF <- stConstruct(x = data, # data set
\nspace = c(\u00ab\u00a0X\u00a0\u00bb, \u00ab\u00a0Y\u00a0\u00bb), # spatial fields
\ntime = t) # time field
\nif (effectuer_la_validation_du_modele) {<\/p>\n
\u00e9chantillonnage pour validation<\/h1>\n
block pas de temps<\/h1>\n
mtot<-length(mSTIDF)
\nmSTIDF$timeY <- format(time(mSTIDF), \u00ab\u00a0%Y-%m-%d\u00a0\u00bb)
\nif (trans==1) {
\nvt<-as.Date(tv,origin=\u00a0\u00bb0000-01-01″)
\n}else{
\nvt<-as.Date(tv)
\n}
\nvalblock_idx<-which(mSTIDF$timeY%in%c(as.character(vt)))
\npaste(\u00ab\u00a0Le bloc de validation du pas de temps comporte\u00a0\u00bb,length(valblock_idx),\u00a0\u00bbobservations\u00a0\u00bb)
\nobs_idx <- setdiff(1:mtot, valblock_idx)<\/p>\n
pourcentage d’observations<\/h1>\n
pc<-pourcentage_points_aleatoires_de_validation\/100
\nset.seed(1)
\nvalrandom_idx <- sample(obs_idx,
\n pc * length(obs_idx),
\n replace = FALSE) %>% sort()
\npaste(\u00ab\u00a0L’echantillon aleatoire de validation comporte\u00a0\u00bb,length(valrandom_idx),\u00a0\u00bbobservations\u00a0\u00bb)
\nobs_idx <- setdiff(obs_idx, valrandom_idx)
\npaste(\u00ab\u00a0Le mod\u00e8le sera calcul\u00e9 \u00e0 partir de\u00a0\u00bb,length(obs_idx),\u00a0\u00bbobservations\u00a0\u00bb)
\nobs<-mSTIDF[obs_idx,]
\nblock<-mSTIDF[valblock_idx,]
\nrandom<-mSTIDF[valrandom_idx,]
\n}else{
\nobs<-mSTIDF
\n}
\nif (units==0) {unit=\u00a0\u00bbsecs\u00a0\u00bb}
\nif (units==1) {unit=\u00a0\u00bbmins\u00a0\u00bb}
\nif (units==2) {unit=\u00a0\u00bbhours\u00a0\u00bb}
\nif (units==3) {unit=\u00a0\u00bbdays\u00a0\u00bb}
\nif (units==4) {unit=\u00a0\u00bbweeks\u00a0\u00bb}
\n\u00ab\u00a0calcul du modele IDE\u2026\u00a0\u00bb
\nIDEmodel <- IDE(f = cnt ~ 1 ,
\n data = obs,
\n dt = as.difftime(1, units = unit),
\n grid_size = grid_size,
\n forecast=forecast)
\n\u00ab\u00a0Ajustement du modele\u2026\u00a0\u00bb
\nfit_results_obs <- fit.IDE(IDEmodel,method=\u00a0\u00bbDEoptim\u00a0\u00bb,itermax=itermax)
\nkp<-fit_results_obs$IDEmodel$get(\u00ab\u00a0k\u00a0\u00bb) %>% unlist()
\n\u00ab\u00a0kernel parameters: \u00ab\u00a0
\nkp
\nco<-coef(fit_results_obs$IDEmodel)
\n\u00ab\u00a0regression coefficients: \u00ab\u00a0
\ncoef(fit_results_obs$IDEmodel)
\nsigma2_eps<-fit_results_obs$IDEmodel$get(\u00ab\u00a0sigma2_eps\u00a0\u00bb)
\n\u00ab\u00a0sigma2_eps :\u00a0\u00bb
\nsigma2_eps
\nabs_ev <- eigen(fit_results_obs$IDEmodel$get(\u00ab\u00a0M\u00a0\u00bb))$values %>%abs()
\n\u00ab\u00a0complex eigenvalues of the evolution matrix M\u00a0\u00bb
\nsummary(abs_ev)
\n\u00ab\u00a0Realisation des predictions\u2026\u00a0\u00bb
\n\u00ab\u00a0Grille complete\u00a0\u00bb
\nST_grid <- predict(fit_results_obs$IDEmodel)
\nST_grid_df <- ST_grid %>%
\ndata.frame() %>%
\nmutate(IDE_pred = ST_grid$Ypred,
\nIDE_predse = ST_grid$Ypredse,
\nIDE_predZse = sqrt(ST_grid$Ypredse^2 +
\nsigma2_eps))
\nST_grid_df$t<-as.integer(format(ST_grid_df$t,\u00a0\u00bb%d\u00a0\u00bb))-1
\ncoords<-dplyr::select(ST_grid_df,X,Y)
\nattr<-dplyr::select(ST_grid_df,-X,-Y,-IDE_pred,-IDE_predse)
\nnames(attr)<-c(\u00ab\u00a0t\u00a0\u00bb,\u00a0\u00bbVal\u00a0\u00bb,\u00a0\u00bbPredict\u00a0\u00bb,\u00a0\u00bbErreurY\u00a0\u00bb,\u00a0\u00bbErreurZ\u00a0\u00bb)
\nspdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
\nResultat_grille<-st_as_sf(spdfN)<\/p>\n
\u00ab\u00a0Previsions observations\u00a0\u00bb
\npred_IDE_obs <- predict(fit_results_obs$IDEmodel,
\nnewdata = obs)<\/p><\/blockquote>\n
obs_df <- obs %>%
\ndata.frame() %>%
\nmutate(IDE_pred = pred_IDE_obs$Ypred,
\nIDE_predse = pred_IDE_obs$Ypredse,
\nIDE_predZse = sqrt(IDE_predse^2 +
\nsigma2_eps))<\/p>\n
obs_df$t<-as.integer(format(obs_df$time,\u00a0\u00bb%d\u00a0\u00bb))-1
\ncoords<-dplyr::select(obs_df,X,Y)
\nattr<-dplyr::select(obs_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
\nnames(attr)<-c(\u00ab\u00a0ID\u00a0\u00bb,\u00a0\u00bbt\u00a0\u00bb,\u00a0\u00bbVal\u00a0\u00bb,\u00a0\u00bbPredict\u00a0\u00bb,\u00a0\u00bbErreurY\u00a0\u00bb,\u00a0\u00bbErreurZ\u00a0\u00bb)
\nattr$t<-as.integer(attr$t)
\nspdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
\nResultat_observations<-st_as_sf(spdfN)
\nif (effectuer_la_validation_du_modele) {
\nif (repertoire_resultats==\u00a0\u00bb) {
\nrepertoire_resultats=(paste(path.expand(‘~’),\u00a0\u00bb\\Documents\u00a0\u00bb,sep=\u00a0\u00bb\u00a0\u00bb))
\n}<\/p>\n
\u00ab\u00a0Predictions des donnees de validation\u00a0\u00bb<\/p>\n
pr\u00e9dictions pour block pas de temps<\/h1>\n
pred_IDE_block <- predict(fit_results_obs$IDEmodel,
\n newdata = block)
\nblock_df <- block %>%
\ndata.frame() %>%
\nmutate(IDE_pred = pred_IDE_block$Ypred,
\nIDE_predse = pred_IDE_block$Ypredse,
\nIDE_predZse = sqrt(IDE_predse^2 +
\nsigma2_eps))
\nblock_df$t<-as.integer(format(block_df$time,\u00a0\u00bb%d\u00a0\u00bb))-1
\ncoords<-dplyr::select(block_df,X,Y)
\nattr<-dplyr::select(block_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
\nnames(attr)<-c(\u00ab\u00a0ID\u00a0\u00bb,\u00a0\u00bbt\u00a0\u00bb,\u00a0\u00bbVal\u00a0\u00bb,\u00a0\u00bbPredict\u00a0\u00bb,\u00a0\u00bbErreurY\u00a0\u00bb,\u00a0\u00bbErreurZ\u00a0\u00bb)
\nattr$t<-as.integer(attr$t)
\nspdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
\nsfn<-st_as_sf(spdfN)
\nnomfic<-paste(repertoire_resultats,\u00a0\u00bb\/\u00a0\u00bb,nom_du_fichier,\u00a0\u00bb_block.shp\u00a0\u00bb,sep=\u00a0\u00bb\u00a0\u00bb)
\nst_write(sfn, nomfic, driver=\u00a0\u00bbESRI Shapefile<\/a>\u00a0\u00bb,append=FALSE)<\/p>\n
pr\u00e9dictions pour pourcentage \u00e9chantillons<\/h1>\n
pred_IDE_random <- predict(fit_results_obs$IDEmodel,
\n newdata = random)
\nalea_df <- random %>%
\ndata.frame() %>%
\nmutate(IDE_pred = pred_IDE_random$Ypred,
\nIDE_predse = pred_IDE_random$Ypredse,
\nIDE_predZse = sqrt(IDE_predse^2 +
\nsigma2_eps))
\nalea_df$t<-as.integer(format(alea_df$time,\u00a0\u00bb%d\u00a0\u00bb))-1
\ncoords<-dplyr::select(alea_df,X,Y)
\nattr<-dplyr::select(alea_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
\nnames(attr)<-c(\u00ab\u00a0ID\u00a0\u00bb,\u00a0\u00bbt\u00a0\u00bb,\u00a0\u00bbVal\u00a0\u00bb,\u00a0\u00bbPredict\u00a0\u00bb,\u00a0\u00bbErreurY\u00a0\u00bb,\u00a0\u00bbErreurZ\u00a0\u00bb)
\nattr$t<-as.integer(attr$t)
\nspdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
\nsfn<-st_as_sf(spdfN)
\nnomfic<-paste(repertoire_resultats,\u00a0\u00bb\/\u00a0\u00bb,nom_du_fichier,\u00a0\u00bb_alea.shp\u00a0\u00bb,sep=\u00a0\u00bb\u00a0\u00bb)
\nst_write(sfn, nomfic, driver=\u00a0\u00bbESRI Shapefile\u00a0\u00bb,append=FALSE)
\npaste(\u00ab\u00a0Les graphiques sont dans le repertoire\u00a0\u00bb,repertoire_resultats)<\/p><\/blockquote>\n
setwd(repertoire_resultats)
\npdf(\u00ab\u00a0Rplots.pdf\u00a0\u00bb)
\nc1<-cor(block_df$cnt,block_df$IDE_pred, method=\u00a0\u00bbpearson\u00a0\u00bb)
\ntitre<-paste(\u00ab\u00a0Bloc-Predictions vs observations. Correlation=\u00a0\u00bb,c1 )
\np1<-ggplot(block_df,aes(cnt, IDE_pred)) +
\n ggtitle(titre) + geom_point() +
\n stat_smooth(method=\u00a0\u00bblm\u00a0\u00bb, se=FALSE)
\nplot(p1)
\nc2<-cor(alea_df$cnt,alea_df$IDE_pred, method=\u00a0\u00bbpearson\u00a0\u00bb)
\ntitre<-paste(\u00ab\u00a0Aleatoire-Predictions vs observations. Correlation=\u00a0\u00bb,c2 )
\np2<-ggplot(alea_df,aes(cnt, IDE_pred)) +
\n ggtitle(titre) + geom_point() +
\n stat_smooth(method=\u00a0\u00bblm\u00a0\u00bb, se=FALSE)
\nplot(p2)
\nBias <- function(x,y) mean(x – y) # x: val. obs.
\nPCV <- function(x,y) mean((x – y)^2) # y: predictions
\nSCV <- function(x,y,v) mean((x – y)^2 \/ v) # v: pred. variances
\nCRPS <- function(x, y, s) verification::crps(x, cbind(y, s))$CRPS
\nDiagblock <- block_df %>% summarise(
\nBias_bloc = Bias(IDE_pred, cnt),
\nPCV_bloc = PCV(IDE_pred, cnt),
\nSCV_bloc = SCV(IDE_pred, cnt, IDE_predZse^2),
\nCRPS_bloc = CRPS(cnt, IDE_pred, IDE_predZse)
\n)
\nDiagrandom <- alea_df %>% summarise(
\nBias_alea = Bias(IDE_pred, cnt),
\nPCV_alea = PCV(IDE_pred, cnt),
\nSCV_alea = SCV(IDE_pred, cnt, IDE_predZse^2),
\nCRPS_alea = CRPS(cnt, IDE_pred, IDE_predZse)
\n)<\/p>\n
\u00ab\u00a0Diagnostics de la validation par bloc temporel\u00a0\u00bb
\nDiagblock
\n\u00ab\u00a0Diagnostics de la validation par echantillonnage aleatoire\u00a0\u00bb
\nDiagrandom
\nnval<-length(block)
\npred_IDE_block <- predict(fit_results_obs$IDEmodel,
\n newdata = block,
\n covariances = TRUE)
\nVeps <- diag(rep(sigma2_eps, nval))
\nL_IDE <- t(chol(pred_IDE_block$Cov + Veps))
\nIntIDE <- coef(fit_results_obs$IDEmodel)
\nnsim <- 100
\nE <- matrix(rnorm(nvalnsim), nval, nsim)
\nSims_IDE <- IntIDE + pred_IDE_block$newdata$Ypred + L_IDE %% E
\nes1<-es_sample(block$cnt, dat = as.matrix(Sims_IDE))
\npaste(\u00ab\u00a0Valeur de es_sample bloc=\u00a0\u00bb,es1)
\ndistances <- block %>%
\ncoordinates() %>%
\ndist() %>%
\nas.matrix()
\nweights <- 0.5^distances
\nvs1<-vs_sample(block$cnt, dat = as.matrix(Sims_IDE),
\nw = weights, p = 1)
\npaste(\u00ab\u00a0Valeur de vs_sample bloc=\u00a0\u00bb,vs1)<\/p><\/blockquote>\n
nval<-length(random)
\npred_IDE_random <- predict(fit_results_obs$IDEmodel,
\nnewdata = random,
\ncovariances = TRUE)
\nVeps <- diag(rep(sigma2_eps, nval))
\nL_IDE <- t(chol(pred_IDE_random$Cov + Veps))
\nIntIDE <- coef(fit_results_obs$IDEmodel)
\nE <- matrix(rnorm(nvalnsim), nval, nsim)
\nSims_IDE <- IntIDE + pred_IDE_random$newdata$Ypred + L_IDE %<\/em>% E
\nes2<-es_sample(random$cnt, dat = as.matrix(Sims_IDE))<\/p>\n
paste(\u00ab\u00a0Valeur de es_sample al\u00e9a=\u00a0\u00bb,es2)
\ndistances <- random %>%
\ncoordinates() %>%
\ndist() %>%
\nas.matrix()
\nweights <- 0.5^distances
\nvs2<-vs_sample(random$cnt, dat = as.matrix(Sims_IDE),
\n w = weights, p = 1)
\npaste(\u00ab\u00a0Valeur de vs_sample al\u00e9a=\u00a0\u00bb,vs2)
\nloglikIDE <- -fit_results_obs$IDEmodel$negloglik()
\npIDE <- 7
\nm<-length(obs)
\nCriteria <- data.frame(AIC = c(0), BIC = c(0),
\n row.names = c( nom_du_fichier))
\naic<- -2loglikIDE + 2pIDE
\nCriteria[nom_du_fichier, \u00ab\u00a0AIC\u00a0\u00bb] <- aic
\nbic<- -2loglikIDE + pIDElog(m)
\nCriteria[nom_du_fichier, \u00ab\u00a0BIC\u00a0\u00bb] <- bic
\nCriteria
\nplot.new()
\ngrid.table(Criteria)
\nBias_bloc = Bias(block_df$IDE_pred, block_df$cnt)
\nPCV_bloc = PCV(block_df$IDE_pred, block_df$cnt)
\nSCV_bloc = SCV(block_df$IDE_pred, block_df$cnt, block_df$IDE_predZse^2)
\nCRPS_bloc = CRPS(block_df$cnt, block_df$IDE_pred, block_df$IDE_predZse)
\nBias_alea = Bias(alea_df$IDE_pred, alea_df$cnt)
\nPCV_alea = PCV(alea_df$IDE_pred, alea_df$cnt)
\nSCV_alea = SCV(alea_df$IDE_pred, alea_df$cnt, alea_df$IDE_predZse^2)
\nCRPS_alea = CRPS(alea_df$cnt, alea_df$IDE_pred, alea_df$IDE_predZse)
\nresume <- matrix(c(\u00ab\u00a0Correlation\u00a0\u00bb,c1,c2,\u00a0\u00bbBiais\u00a0\u00bb,Bias_bloc,Bias_alea,\u00a0\u00bbPCV\u00a0\u00bb,PCV_bloc,PCV_alea,\u00a0\u00bbSCV\u00a0\u00bb,SCV_bloc,SCV_alea,\u00a0\u00bbCRPS\u00a0\u00bb,CRPS_bloc,CRPS_alea,\u00a0\u00bbES\u00a0\u00bb,es1,es2,\u00a0\u00bbVS\u00a0\u00bb,vs1,vs2),ncol=3,byrow=TRUE)
\ncolnames(resume) <- c(\u00ab\u00a0Score\u00a0\u00bb,\u00a0\u00bbBloc t\u00a0\u00bb,\u00a0\u00bbAleatoire\u00a0\u00bb)
\nresume <- as.table(resume)
\nresume
\nplot.new()
\ngrid.table(resume,rows=NULL)
\ndev.off()
\n}else{
\npred_IDE_obs <- predict(fit_results_obs$IDEmodel,
\n newdata = obs)
\nobs_df <- obs %>%
\ndata.frame() %>%
\nmutate(IDE_pred = pred_IDE_obs$Ypred,
\nIDE_predse = pred_IDE_obs$Ypredse,
\nIDE_predZse = sqrt(IDE_predse^2 +
\nsigma2_eps))
\nobs_df$t<-as.integer(format(obs_df$time,\u00a0\u00bb%d\u00a0\u00bb))-1
\ncoords<-dplyr::select(obs_df,X,Y)
\nattr<-dplyr::select(obs_df,sp.ID,t,cnt,IDE_pred,IDE_predse,IDE_predZse)
\nnames(attr)<-c(\u00ab\u00a0ID\u00a0\u00bb,\u00a0\u00bbt\u00a0\u00bb,\u00a0\u00bbVal\u00a0\u00bb,\u00a0\u00bbPredict\u00a0\u00bb,\u00a0\u00bbErreurY\u00a0\u00bb,\u00a0\u00bbErreurZ\u00a0\u00bb)
\nattr$t<-as.integer(attr$t)
\nspdfN <- SpatialPointsDataFrame(coords = coords, data = attr,proj4string=CRS(crs))
\nsfn<-st_as_sf(spdfN)
\n} <\/div><\/div><\/blockquote>\n\n\n\n
ideQgis_v1-4<\/a>T\u00e9l\u00e9charger<\/a><\/div>\n\n\n\n
(Voir les diff\u00e9rentes versions<\/a> )<\/p>\n\n\n\n
Les donn\u00e9es<\/h3>\n\n\n\n
Les donn\u00e9es sont extraites des D\u00e9nombrements du Relev\u00e9 des oiseaux nicheurs (BBS). Ces donn\u00e9es proviennent du Relev\u00e9 des oiseaux nicheurs de l\u2019Am\u00e9rique du Nord.
Nous consid\u00e9rons en particulier les d\u00e9nombrements annuels du troglodyte de Caroline (Thryothorus ludovicianus) pour la p\u00e9riode de 1967 \u00e0 2014. L\u2019unit\u00e9 d\u2019\u00e9chantillonnage du BBS est un parcours routier d\u2019environ 39,2 km. Dans chaque unit\u00e9 d\u2019\u00e9chantillonnage, les observateurs b\u00e9n\u00e9voles font 50 arr\u00eats et comptent les oiseaux pendant une p\u00e9riode de 3 minutes lorsqu\u2019ils ex\u00e9cutent leurs itin\u00e9raires (habituellement en juin). Il y a plus de 4000 routes dans l\u2019enqu\u00eate nord-am\u00e9ricaine, mais toutes ne sont pas disponibles chaque ann\u00e9e. Aux fins des analyses pr\u00e9sent\u00e9es ici, nous consid\u00e9rons que les d\u00e9nombrements des parcours totaux se produisent chaque ann\u00e9e (pendant la saison de reproduction) et d\u00e9finissons l\u2019emplacement spatial de chaque parcours comme \u00e9tant le centro\u00efde de la route. Ainsi, nous consid\u00e9rons que les donn\u00e9es sont discr\u00e8tes dans le temps, g\u00e9ostatistiques et irr\u00e9guli\u00e8res dans l\u2019espace, et non gaussiennes dans le sens qu\u2019elles sont des comptages.<\/p>\n\n\n\n
Copiez les fichiers shapefile dans un r\u00e9pertoire de donn\u00e9es et chargez la couche dans QGis.<\/p>\n\n\n\n
$\"\"$ <\/a><\/figure>\n\n\n\n
Les scripts<\/h3>\n\n\n\n
Comme leur nom l’indique, il y a deux types de mod\u00e8les IDE:<\/p>\n\n\n\n
ceux que l’on consid\u00e8re avec un noyau invariant dans l’espace : le script ide_modele_invariant.rsx<\/li>
ceux que l’on consid\u00e8re avec un noyau qui varie dans l’espace: le script ide_modele_variant.rsx<\/li><\/ul>\n\n\n\n
Nous verrons \u00e0 quoi cela correspond un peu plus bas.<\/p>\n\n\n\n
Ces scripts sont \u00e0 copier, avec leur fichier d’aide, dans le r\u00e9pertoire C:\\Users\\xxxxx\\AppData\\Roaming\\QGIS\\QGIS3\\profiles\\default\\processing\\rscripts<\/p>\n\n\n\n
Si QGis est configur\u00e9 pour les scripts R vous trouverez ces deux scripts dans la fen\u00eatre Traitements->R->modele IDE<\/p>\n\n\n\n
$\"\"$ <\/a><\/figure>\n\n\n\n
Les \u00e9l\u00e9ments du mod\u00e8le IDE<\/h1>\n\n\n\n
Le mod\u00e8le IDE est d\u00e9fini par <\/p>\n\n\n\n
$\"\"$ <\/a><\/figure>\n\n\n\n
On va essayer de traduire ceci en langage courant:<\/p>\n\n\n\n
La valeur du processus qui nous int\u00e9resse (Y<\/em><\/strong>) \u00e0 un endroit donn\u00e9 s<\/em><\/strong> et \u00e0 un temps donn\u00e9 t<\/em><\/strong> r\u00e9sulte de trois \u00e9l\u00e9ments:<\/p>\n\n\n\n
la valeur prise par ce qu’on appellera le noyau du mod\u00e8le (m(s,x;\u03b8p))<\/li>
la valeur du processus au moment pr\u00e9c\u00e9dent t-1<\/em><\/strong><\/li>
une \u00ab\u00a0erreur\u00a0\u00bb ou \u00ab\u00a0bruit\u00a0\u00bb ind\u00e9pendant du processus \u03b7t<\/li><\/ul>\n\n\n\n $\"\"$ <\/a><\/figure>\n\n\n\n
Observations ou valeurs du processus?<\/h2>\n\n\n\n
Tout d’abord, sachez que certaines notations ont leur importance. Ici on travaille pour obtenir Yt. Dans la litt\u00e9rature vous verrez deux types de notation: Yt ou Zt. Ceci a son importance. Quand on parle de Z on parle de vos donn\u00e9es observ\u00e9es. Quand on parle de Y on parle du processus sous-jacent. La diff\u00e9rence vient du fait que ce que l’on observe (Z) est toujours une image biais\u00e9e du processus sous-jacent (Y), ne serait-ce que par les erreurs de mesure.<\/p>\n\n\n\n
Le noyau du mod\u00e8le<\/h2>\n\n\n\n
Dans le noyau du mod\u00e8le m(s, x; \u03b8p) la valeur m que l’on va appliquer d\u00e9pend de trois \u00e9l\u00e9ments:<\/p>\n\n\n\n
le plus simple est s<\/em><\/strong>, c’est \u00e0 dire la localisation dans l’espace du point qui nous int\u00e9resse<\/li>
puis, nous avons \u03b8p qui correspond \u00e0 quatre param\u00e8tres qui d\u00e9terminent la redistribution du poids qui sera affect\u00e9 \u00e0 la valeur de Yt-1. Pour simplifier, si on a un processus avec une croissance de 10% entre t-1 et t, \u03b8p aura une valeur de 1.1 (ceci est une grosssse simplification)<\/li>
le troisi\u00e8me \u00e9l\u00e9ment est x, la variation spatiale, et il va falloir s’arr\u00eater un peu pour bien le comprendre.<\/li><\/ul>\n\n\n\n
La variation spatiale<\/h2>\n\n\n\n
Reprenons les poids \u03b8p 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 \u00eatre le cas. Le pourcentage de diff\u00e9rence entre les deux pas de temps peut r\u00e9sulter d’une fonction d\u00e9pendante de l’endroit o\u00f9 l’on se trouve (s). Dans ce cas, la valeur du noyau d\u00e9pendra toujours des param\u00e8tres \u03b8p, mais aussi de la valeur de la fonction x \u00e0 l’endroit donn\u00e9 (s). On dira alors que notre noyau est variant.<\/p>\n\n\n\n
Les param\u00e8tres \u03b8p<\/h2>\n\n\n\n
Le noyau du mod\u00e8le IDE peut \u00eatre associ\u00e9e \u00e0 une courbe gaussienne. Le premier des quatre param\u00e8tres theta_1 correspond \u00e0 la hauteur (amplitude) de la courbe en cloche. Le deuxi\u00e8me, theta 2, correspond \u00e0 l’aperture de la cloche la rendant plus pointue ou plus \u00e9vas\u00e9e. Les deux autres param\u00e8tres, theta 3 et theta 4 correspondent \u00e0 un d\u00e9calage sur la dimension des x et des y, respectivement, en fonction de l’endroit o\u00f9 l’on se trouve.<\/p>\n\n\n\n
Le moment est venu de se jeter \u00e0 l’eau. On va voir deux types de simulation. La premi\u00e8re en consid\u00e9rant que le noyau du mod\u00e8le ne varie pas spatialement, la deuxi\u00e8me en consid\u00e9rant que le noyau d\u00e9pend de l’endroit o\u00f9 l’on se trouve et donc qu’il varie spatialement.<\/p>\n\n\n\n
Exemple de simulation avec un noyau spatialement invariant<\/h1>\n\n\n\n
$\"\"$ <\/a><\/figure><\/div>\n\n\n\n
Le mod\u00e8le d\u2019\u00e9quation int\u00e9gro-diff\u00e9rentielle (IDE) est : (1)construit \u00e0 l\u2019aide de la fonction \u00ab IDE \u00bb, (2) ajust\u00e9 \u00e0 l\u2019aide de la fonction \u00ab IDE.fit \u00bb et (3)utilis\u00e9 pour la pr\u00e9vision \u00e0 l\u2019aide de la fonction \u00ab predict \u00bb.<\/p><\/blockquote><\/figure>\n\n\n\n
Si vous ouvrez dans un \u00e9diteur de texte le script ide_modele_invariant.rsx vous verrez ces trois \u00e9tapes.<\/p>\n\n\n\n
Construction du mod\u00e8le IDE<\/h2>\n\n\n\n
Tout d\u2019abord on d\u00e9finit les caract\u00e9ristiques du mod\u00e8le. Dans notre cas:<\/p>\n\n\n\n
IDEmodel <- IDE(f = cnt ~ 1 ,\n data = mSTIDF,\n dt = as.difftime(1, units = unit),\n grid_size = grid_size,\n forecast=forecast)<\/code><\/pre>\n\n\n\nEn plus de la \u00ab\u00a0formule\u00a0\u00bb (f=) qui d\u00e9termine quels sont les colonnes de donn\u00e9es \u00e0 mod\u00e9liser (un attribut de la table, plus les coordonn\u00e9es de chaque point), la couche de donn\u00e9es (data=) mise \u00e0 un format R de donn\u00e9es spatio-temporelles (STIDF), nous avons \u00e0 rentrer trois param\u00e8tres.<\/p>\n\n\n\n Le pas de temps dt et les unit\u00e9s<\/h2>\n\n\n\nContrairement \u00e0 ce que l’on peut penser, ce n’est pas si simple d’indiquer le pas de temps de nos donn\u00e9es!!!<\/p>\n\n\n\n Les outils de mod\u00e9lisation n\u00e9cessitent des pas de temps de m\u00eame longueur. Jusque-l\u00e0, tout para\u00eet simple. Sauf que si nos donn\u00e9es correspondent \u00e0 des observations mensuelles, le probl\u00e8me est que les mois n’ont pas la m\u00eame longueur… de 28 \u00e0 31 jours. Et si nos donn\u00e9es sont annuelles, les ann\u00e9es n’ont pas la m\u00eame longueur (bissextiles ou pas). Ceci fait que les unit\u00e9s de pas de temps acceptables pour les mod\u00e8les spatio-temporels vont de la seconde \u00e0 la semaine, pas plus.<\/p>\n\n\n\n On est alors oblig\u00e9s de biaiser un peu pour traiter des donn\u00e9es mensuelles ou annuelles.<\/p>\n\n\n\n Le moyen le plus simple pour ne pas avoir \u00e0 jongler avec nos donn\u00e9es est de calculer un num\u00e9ro d’ordre de nos observations est d’utiliser cet ordre comme param\u00e8tre temporel. <\/span><\/strong><\/p>\n\n\n\n Supposons que nous avons de donn\u00e9es 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\u00eame, au lieu d’utiliser la vraie date, on \u00e9tablit une \u00e9chelle de 1 \u00e0 X correspondante.<\/p>\n\n\n\n Dans le cas o\u00f9 vous utiliserez le num\u00e9ro d’ordre, il faut d\u00e9finir les unit\u00e9s avec la valeur \u00ab\u00a0days\u00a0\u00bb<\/span><\/strong>.<\/p>\n\n\n\n Il n’est pas possible de pr\u00e9voir tous les formats possibles de date et d’heure. Si vous rencontrez des probl\u00e8mes, calculez un num\u00e9ro d’ordre et d\u00e9finissez les unit\u00e9s comme \u00e9tant des jours. Ceci ne modifie en rien les r\u00e9sultats.<\/p>\n\n\n\n
La grille de mod\u00e9lisation<\/h2>\n\n\n\n
Nos donn\u00e9es sont sous forme de points d\u2019\u00e9chantillonnage r\u00e9partis sur une zone d’\u00e9tude. Pour la mise en place du mod\u00e8le nous d\u00e9finissons un maillage carr\u00e9 sur lequel on \u00e9valuera les param\u00e8tres du mod\u00e8le et ensuite on calculera les valeurs de Y. Le param\u00e8tre grid_size d\u00e9termine le nombre de mailles de la grille. N’oubliez pas que le nombre de cellules d\u00e9terminera le temps de traitement… et il est d\u00e9j\u00e0 tr\u00e8s long sans faire appel \u00e0 des grilles trop grandes.<\/p>\n\n\n\n
Les pr\u00e9visions<\/h2>\n\n\n\n
Le mod\u00e8le que nous mettons en place va nous permettre de calculer les pr\u00e9dictions pour chaque cellule de la grille et chaque pas de temps pr\u00e9sent dans nos donn\u00e9es. Mais il peut nous permettre de pr\u00e9voir les valeurs pour les prochains pas de temps dans le futur. La valeur du param\u00e8tre forecast <\/em><\/strong>indique le nombre de pas de temps, \u00e0 partir du dernier de nos observations, \u00e0 calculer. Si nos donn\u00e9es sont annuelles et vont de 1994 \u00e0 2010, forecast= 2 fera que le mod\u00e8le calcule les donn\u00e9es des ann\u00e9es 2011 et 2012.<\/p>\n\n\n\n
Dans ces scripts nous n’avons pas impl\u00e9ment\u00e9, pour ne pas surcharger inutilement l’interface, un autre param\u00e8tre possible : hindcast. Celui-ci ag\u00eet comme forecast mais pour la p\u00e9riode ant\u00e9rieure \u00e0 nos donn\u00e9es. Si on ajoute dans le code, apr\u00e8s forecast \u00ab\u00a0hindcast=2\u00a0\u00bb, le mod\u00e8le calculera les donn\u00e9es pour les deux ann\u00e9es ant\u00e9rieures \u00e0 nos donn\u00e9es, 1992 et 1993.<\/p>\n\n\n\n
Si vous souhaitez utiliser l’option hindcast, il suffit de modifier le code du script avec la valeur que vous souhaitez.<\/p>\n\n\n\n
`IDEmodel <- IDE(f = cnt ~ 1 ,\n data = mSTIDF,\n dt = as.difftime(1, units = unit),\n grid_size = grid_size,\n forecast=forecast,hindcast=1)<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n`Ajustement du mod\u00e8le IDE<\/h2>\n\n\n\nMaintenant que nous avons d\u00e9fini quels param\u00e8tres on veut utiliser dans notre mod\u00e8le, on \u00ab\u00a0ajuste\u00a0\u00bb le mod\u00e8le \u00e0 partir des observations. Ceci est fait avec la commande fit.IDE<\/p>\n\n\n\n fit_results <- fit.IDE(IDEmodel,method=\"DEoptim\",itermax=itermax)<\/code><\/pre>\n\n\n\nLe mod\u00e8le IDE calcul\u00e9 de l\u2019objet, ici l’objet fit_results, contient les estimations initiales des param\u00e8tres, ainsi que les pr\u00e9dictions pour tous les pas de temps des estimations des param\u00e8tres. Les param\u00e8tres dans ce cas sont la variance de l’erreur de mesure, la variance de la perturbation al\u00e9atoire \u03b7t(s)(dont la structure de covariance est fixe), les param\u00e8tres du noyau, et les coefficients de r\u00e9gression \u03b2.<\/p>\n\n\n\n Nous parlons bien des valeurs des diff\u00e9rents param\u00e8tres du mod\u00e8le pour chaque pas de temps et non des valeurs du processus \u00e0 chaque maille de la grille et chaque pas de temps.<\/p>\n\n\n\n Ce calcul-l\u00e0 sera fait dans l’\u00e9tape suivante.<\/p>\n\n\n\n Dans cette \u00e9tape, le seul param\u00e8tre sur lequel vous pouvez intervenir est le nombre d’it\u00e9rations pour l’ajustement du mod\u00e8le. Il n’y a pas de recette, tout d\u00e9pend de la complexit\u00e9 du mod\u00e8le et des observations disponibles.<\/p>\n\n\n\n Par contre, il faut savoir que le temps de calcul de cette \u00e9tape est tr\u00e8s,tr\u00e8s long. Avec un ordinateur plut\u00f4t puissant, pour les donn\u00e9es d’exemple fournies,le temps de calcul pour une centaine d’it\u00e9rations tourne autour d’une heure ou deux…<\/p>\n\n\n\n Les r\u00e9sultats de chaque it\u00e9ration sont affich\u00e9s dans l’onglet \u00ab\u00a0Journal\u00a0\u00bb de l’interface et dans la sortie console du script. Vous pourrez voir l\u00e0 \u00e0 partir de quelle it\u00e9ration on a atteint les valeurs finales, et si vous devez r\u00e9-ex\u00e9cuter le calcul, vous pouvez limiter le nombre d’it\u00e9rations \u00e0 ce chiffre-l\u00e0.<\/p>\n\n\n\n Dans le journal vous aurez aussi le d\u00e9tail des diff\u00e9rents calculs des param\u00e8tres.<\/p>\n\n\n\n Calcul des pr\u00e9dictions sur la maille d\u00e9finie<\/h2>\n\n\n\nLa derni\u00e8re \u00e9tape consiste \u00e0 utiliser l’ajustement r\u00e9alis\u00e9 et pr\u00e9dire pour chaque n\u0153ud de la maille spatiale d\u00e9finie dans la premi\u00e8re \u00e9tape et pour chaque pas de temps souhait\u00e9s la valeur de ,notre processus Y.<\/p>\n\n\n\n Ceci est fait avec la fonction \u00ab\u00a0predict\u00a0\u00bb<\/p>\n\n\n\n ST_grid <- predict(fit_results$IDEmodel)<\/code><\/pre>\n\n\n\nLes lignes suivantes de code, apr\u00e8s celle-ci, transforment le format grid de R en une couche de points au format d\u00e9fini par l’utilisateur dans le champ \u00ab\u00a0output\u00a0\u00bb.<\/p>\n\n\n\n Donn\u00e9es exemple<\/h1>\n\n\n\nPassons \u00e0 la partie pratique, m\u00eame s’il n’y a pas beaucoup de manipulations et beaucoup de patience \u00e0 avoir.<\/p>\n\n\n\n Chargez la couche MOcarolinawren dans QGis.<\/p>\n\n\n\n Dans la bo\u00eete de traitements->R->modele IDE double cliquez sur ide modele invariant<\/p>\n\n\n\n <\/a><\/figure>\n\n\n\nDans le champ t<\/em><\/strong> s\u00e9lectionnez l’attribut t<\/em><\/strong> qui correspond \u00e0 un num\u00e9ro d’ordre de l’ann\u00e9e. Ses valeurs sont comprises entre 1 et 21.<\/p>\n\n\n\n Dans le champ val <\/em><\/strong>s\u00e9lectionnez l’attribut cnt <\/em><\/strong>qui correspond aux comptages d’oiseaux.<\/p>\n\n\n\n Dans le champ grid_size, laissez la valeur par d\u00e9faut (50). On aura alors une grille en sortie de 50 x 50 points.<\/p>\n\n\n\n Dans le champ units<\/em><\/strong>, comme nous l’avons indiqu\u00e9 plus haut, on ne peut pas d\u00e9finir le pas de temps en ann\u00e9es. On dit alors qu’il s’ag\u00eet de jours (days<\/em><\/strong>) car le fait que tous les quatre ans il y ait un jour de diff\u00e9rence (1\/365\u00e8me de plus dans la longueur du pas) est s\u00fbrement n\u00e9gligeable.<\/p>\n\n\n\n Dans le champ forecast<\/em><\/strong>, laissez la valeur par d\u00e9faut (1). Dans la grille de sortie on aura alors un pas de temps suppl\u00e9mentaire, avec la valeur t<\/em><\/strong> \u00e0 22 qui correspondra aux pr\u00e9dictions pour l’ann\u00e9e suivante \u00e0 la derni\u00e8re dans les observations.<\/p>\n\n\n\n Si vous \u00eates impatient, r\u00e9duisez le nombre d’it\u00e9rations par d\u00e9faut. Si non, armez vous de patience et laissez ce chiffre.<\/p>\n\n\n\n D\u00e9finissez la couche en sortie, au format que vous pr\u00e9f\u00e9rez (shp, gpkg,…) et un fichier pour les sorties de la console. Ce dernier est optionnel mais il contient tous les param\u00e8tres de mod\u00e9lisation.<\/p>\n\n\n\n Cliquez sur Ex\u00e9cuter<\/em><\/strong>… bon courage.<\/p>\n\n\n\nR\u00e9sultats de l’exemple<\/h1>\n\n\n\nQuatre heures apr\u00e8s, nous obtenons une couche qui s’affiche dans QGis:<\/p>\n\n\n\n<\/a><\/figure>\n\n\n\n<\/p>\n\n\n\n Elle contient notre grille de pr\u00e9dictions de taille 30×50.<\/p>\n\n\n\n Avant de voir plus en d\u00e9tail les r\u00e9sultats, voyons quelles informations sont fournies par la sortie console du script. Vous le trouvez dans l’onglet \u00ab\u00a0Journal<\/strong>\u00a0\u00bb de la fen\u00eatre du script, en tant que lien dans le panneau \u00ab\u00a0Visualiseur de r\u00e9sultats<\/strong>\u00a0\u00bb et, si vous avez d\u00e9fini un fichier permanent dans le champ \u00ab\u00a0R Console Output [optionnel]<\/strong>\u00ab\u00a0, dans ce fichier.<\/p>\n\n\n\n\n[1] \"calcul du modele IDE...\"\n[1] \"Ajustement du modele...\"\nIteration: 1 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397\nIteration: 2 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397\nIteration: 3 bestvalit: 1855.495069 bestmemit: 0.011295 -0.470852 1.598254 1.059066 3.813749 -0.653397\n\t\t.\n\t\t.\n\t\t.\nIteration: 99 bestvalit: 1800.983447 bestmemit: 0.052054 -0.350126 -0.459083 2.066409 3.177588 -0.728638\nIteration: 100 bestvalit: 1800.983447 bestmemit: 0.052054 -0.350126 -0.459083 2.066409 3.177588 -0.728638\n[1] \"kernel parameters: \"\npar1 par2 par3 par4\n52.05424217 0.03015948 -0.45908287 2.06640855\n[1] \"regression coefficients: \"\nIntercept X Y\n70.3594335 -0.3226177 -2.4435913\n[1] \"sigma2_eps :\"\npar5\n23.98881\n[1] \"complex eigenvalues of the evolution matrix M\"\nMin. 1st Qu. Median Mean 3rd Qu. Max.\n0.00000 0.00000 0.00301 0.18171 0.34210 0.98394\n[1] \"Realisation des predictions...\"\nWarning message:\nIn abbreviate_shapefile_names(obj) :\nField names abbreviated for ESRI Shapefile driver\n<\/code><\/pre>\n\n\n\nParam\u00e8tres descriptifs<\/h2>\n\n\n\nLes sorties kernel parameters, regression coefficients et sigma2_eps indiquent les valeurs obtenues pour l’ajustement du mod\u00e8le. Ils n’ont pas de signification par rapport \u00e0 la \u00ab\u00a0justesse\u00a0\u00bb de l’ajustement, sauf \u00e0 les conna\u00eetre dans la r\u00e9alit\u00e9. Ceci n’est pas g\u00e9n\u00e9ralement le cas, mais on peut utiliser IDE pour tester des simulations de mod\u00e8les pour lesquelles nous fixons ces param\u00e8tres. On peut alors comparer l’ajustement des param\u00e8tres par rapport aux valeurs fournies aux simulations.<\/p>\n\n\n\nValeurs propres de la matrice M<\/h2>\n\n\n\nPar contre, il y a une sortie qui nous permet de v\u00e9rifier que le r\u00e9sultat peut \u00eatre pris en compte et passer \u00e0 l’\u00e9tape suivante de toute mod\u00e9lisation : la validation des r\u00e9sultats.<\/p>\n\n\n\n Comme nous l’avons vu plus haut, la formule de l’\u00e9quation int\u00e9gro-diff\u00e9rentielle est<\/p>\n\n\n\n<\/a><\/figure>\n\n\n\nMais, une fois le mod\u00e8le calcul\u00e9 et ajust\u00e9, on peut l’\u00e9crire sous la forme:<\/p>\n\n\n\n<\/a><\/figure><\/div>\n\n\n\no\u00f9 M est une matrice qu’on peut appeler \u00ab\u00a0matrice de transition\u00a0\u00bb, 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<\/sub> \u00e0 Y_{t<\/sub>.<\/p>\n\n\n\n}} On peut donc calculer les valeurs propres (eigenvalues) de cette matrice. Ceci permet de savoir si le mod\u00e8le calcul\u00e9 est stable ou pas.<\/p>\n\n\n\n Dans notre sortie, la synth\u00e8se de ces valeurs propres appara\u00eet sous le titre \u00ab\u00a0complex eigenvalues of the evolution matrix M<\/strong>\u00ab\u00a0.<\/p>\n\n\n\n Nous notons que si la valeur maximale des valeurs propres est >=1 , alors le mod\u00e8le, est instable, et Yt grandira sans limite lorsque t augmente. Inversement, si le module maximal de toutes les valeurs propres est inf\u00e9rieur \u00e0 1, alors le mod\u00e8le est stable.<\/p>\n\n\n\n Dans notre exemple, la valeur maximale des valeurs propres est de 0.98394. Le mod\u00e8le peut donc \u00eatre consid\u00e9r\u00e9 comme stable.<\/p>\n\n\n\nLa couche r\u00e9sultat<\/h3>\n\n\n\nA la fin du traitement vous aurez une couche de points contenant la grille de pr\u00e9dictions:<\/p>\n\n\n\n<\/a><\/figure>\n\n\n\nL’attribut t correspond au pas de temps. Vous trouverez les 21 pas originaux plus le pas 22 qui correspond aux pr\u00e9dictions du prochain pas (forecast).<\/p>\n\n\n\n Ypred <\/em><\/strong>correspond \u00e0 la valeur calcul\u00e9e par le mod\u00e8le, Ypredse<\/strong> \u00e0 la valeur de l’erreur standard dans le domaine du processus (Y) pour le point et le pas de temps. YpredZ<\/strong> correspond \u00e0 l’erreur standard dans le domaine des observations (Z) pour le point et le pas de temps.<\/p>\n\n\n\n En appliquant le m\u00eame filtre sur les donn\u00e9es en entr\u00e9e et la grille en sortie sur l’attribut t<\/em><\/strong>, vous pourrez voir la correspondance:<\/p>\n\n\n\n<\/a><\/figure>\n\n\n\nDans le prochain article nous verrons l’utilisation du deuxi\u00e8me script, avec un noyau variant spatialement.<\/p>\n","protected":false},"excerpt":{"rendered":" En \u00e9cologie on se retrouve, tr\u00e8s souvent, face \u00e0 des processus qui sont li\u00e9s \u00e0 l’espace mais aussi qui \u00e9voluent dans le temps. Nous sommes maintenant habitu\u00e9s \u00e0 l’analyse spatiale qui est devenu \u00e0 la port\u00e9e…<\/p>\n","protected":false},"author":1,"featured_media":9297,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"give_campaign_id":0,"_bbp_topic_count":0,"_bbp_reply_count":0,"_bbp_total_topic_count":0,"_bbp_total_reply_count":0,"_bbp_voice_count":0,"_bbp_anonymous_reply_count":0,"_bbp_topic_count_hidden":0,"_bbp_reply_count_hidden":0,"_bbp_forum_subforum_count":0,"sfsi_plus_gutenberg_text_before_share":"","sfsi_plus_gutenberg_show_text_before_share":"","sfsi_plus_gutenberg_icon_type":"","sfsi_plus_gutenberg_icon_alignemt":"","sfsi_plus_gutenburg_max_per_row":"","_monsterinsights_skip_tracking":false,"_monsterinsights_sitenote_active":false,"_monsterinsights_sitenote_note":"","_monsterinsights_sitenote_category":0,"jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[62,1641],"tags":[340,1659,1658,51,1657,1627,170],"class_list":["post-9240","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-qgis-2","category-scripts-r","tag-ide","tag-integro-differentielle","tag-invariant","tag-modele","tag-noyau","tag-script-r","tag-spatio-temporel"],"aioseo_notices":[],"jetpack_featured_media_url":"https:\/\/i0.wp.com\/www.sigterritoires.fr\/wp-content\/uploads\/2021\/02\/ide5b.png?fit=587%2C503&ssl=1","jetpack_shortlink":"https:\/\/wp.me\/p6XU0A-2p2","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/posts\/9240","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/comments?post=9240"}],"version-history":[{"count":0,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/posts\/9240\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/media\/9297"}],"wp:attachment":[{"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/media?parent=9240"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/categories?post=9240"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.sigterritoires.fr\/index.php\/wp-json\/wp\/v2\/tags?post=9240"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}