Ce document présente la méthodologie détaillée d’ajustement du modèle de projection des emplois de l’ARTM. Ce modèle s’appuie sur un travail préliminaire des chercheurs Cédric Brunelle (INRS) et Jean Dubé (Université Laval).
Dernière compilation
2026-07-22
1 Introduction
Le modèle a pour objectif de prédire le nombre de personnes actives de 15 ans ou plus travaillant dans les différentes ZAT du territoire de L’ARTM. Dans le reste du document, cette variable sera aussi appelée nombre d’emplois.
Pour chaque ZAT, nous disposons des données historiques de 1996, 2001, 2006, 2011, 2016 et 2021. Il s’agit donc d’une structure en panel, car nous avons de multiples valeurs qui se suivent dans le temps pour une même observation. Il s’agit aussi d’un panel large, car nous avons plus de ZAT que d’années.
1.1 Formulation originale du modèle
Le modèle original a été réalisé dans le logiciel Stata avec la fonction xtabond. Elle implémente l’estimateur d’Arellano-Bond basé sur la méthode des moments.
Ce modèle introduit la dimension temporelle dans l’équation de régression en ajoutant des lags de la variable dépendante. Cependant, puisque nous sommes dans un contexte de données par panel, cela soulève un problème d’endogénéité. En effet, nous risquons d’avoir des corrélations importantes entre les effets propres à chaque ZAT et la variable dépendante spatialement décalée.
\(\epsilon_{it}\) un terme d’erreur distribué normalement et idiosyncratique
Pour un tel modèle, l’estimation par OLS pose un problème majeur, car les deux termes d’erreurs \(\alpha_i\) et \(\epsilon_{it}\) sont confondus. La corrélation entre \(y_{it-1}\) et \(\alpha_i\) pose alors un problème d’endogénéité.
Pour contourner ce problème, l’estimateur d’Arellano-Bond modélise plutôt la différence des variables décalées temporellement.
\[
E\left[y_{i, t-s} \cdot \Delta \varepsilon_{i t}\right]=0 \quad \text { for } s \geq 2
\]
En d’autres termes, l’endogénéité est contrôlée par des versions décalées plus profondes de la variable dépendante. Une telle approche n’est pas idéale dans notre contexte, car nous avons assez peu de réplications temporelles dans nos données. Concrètement, utiliser un lag minimal de 2 impose de nous passer de deux périodes dans nos données, soit près du tiers des données.
1.2 Reformulation proposée
Nous proposons de reformuler ce modèle pour gagner un peu de flexibilité et limiter certains effets problématiques observés dans les résultats des projections du modèle original.
Notre formulation suppose que le nombre d’emplois par ZAT est la somme de cinq éléments:
L’effet propre de chaque ZAT. En effet, chaque ZAT a son volume d’emploi original, une moyenne qui lui est propre. Nous ne disposons que de 5 observations par ZAT, nous ne pouvons donc pas estimer cet effet comme un effet fixe. Nous allons donc l’estimer comme un effet aléatoire, soit une déviation de l’intercept général propre à chaque ZAT et provenant d’une distribution normale. Cette estimation repose sur la mutualisation de l’information entre ZAT.
Les effets fixes, soit la contribution d’un ensemble de termes tels que la présence d’une station de métro, de train, d’un hôpital, etc. Ces effets peuvent changer dans le temps et sont identiques pour chaque ZAT. Cependant, puisque chaque ZAT à son effet propre et que la fonction box-cox (utilisée pour transformer la variable Y, voir plus loin) produit des effets quasi multiplicatifs, ces effets fixes ont des impacts plus grands pour des ZAT déjà plus grandes.
Les effets non linéaires. Ils ont un rôle similaire aux effets fixes, mais un lien non linéaire est estimé plutôt qu’un simple coefficient, notamment pour capturer des effets de seuils.
La tendance temporelle. Il s’agit d’une tendance régionale indépendante de la ZAT. Nous avons notamment constaté que l’emploi a tendance à croitre dans la région à cause de la croissance économique et de la population. Il s’agit donc d’une tendance d’arrière-plan. Nous faisons varier cette tendance par sous-région pour capter des dynamismes différenciés. Nous nous attendons à avoir un effet positif (croissance), que nous allons modéliser avec une transformation logarithmique pour représenter une croissance qui ralentit.
Notre formulation suppose que chaque ZAT est caractérisée par une structure temporelle résiduelle suivant un processus de type ARMA(1). Ce processus tient compte des dynamiques de long terme et de court terme au sein d’une ZAT. Il s’agit d’un processus stationnaire qui n’est donc pas responsable de l’éventuelle croissance au sein de la ZAT. Notons aussi que la variance n’est pas supposée homogène, mais est autorisée à varier par grande région.
Chaque ZAT combine donc deux termes temporels: une tendance régionale atténuée localement et une tendance locale stationnaire.
Reading layer `empm23_region_analyse' from data source
`C:\Users\jgelb\OneDrive - ARTM\Donnees ouvertes\Données et analytique - DO-017_Projections_emplois\Ressources\ajustement du modèle\data\enquete.gpkg'
using driver `GPKG'
Simple feature collection with 8 features and 13 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 222446.5 ymin: 4994150 xmax: 342561 ymax: 5111729
Projected CRS: NAD83(CSRS) / MTM zone 8
Code
# Création d'un dictionnaire associant chaque ZAT à sa régiondico_reg_ori <- zat %>%st_centroid() %>%st_intersection(regions['nomreg2'])# Idem, mais pour les super-zatdico_super <- zat %>%st_centroid() %>%st_intersection(super_zat) %>%st_intersection(regions['nomreg']) %>%st_drop_geometry() %>%fselect(c('id_zat', 'id', 'nomreg')) %>%rename(c('idzat'='id_zat', 'idzone'='id'))# lecture du fichier de données structurées partagé par les chercheurs des # université Laval et INRSbase_data <-read_dta('data/ScenarioProjection_Etages.dta')# suppression d'une colonne inutilebase_data$`_Y1`<-NULL# création de variables complémentaires et nettoyage de certaines colonnesbase_data <- base_data %>%mutate(ln_Area_Commercial =log(Total_Area_Commercial), ln_Area_Bureau =log(Total_Area_Bureau) ) %>%mutate_if(is.numeric, ~replace(., is.na(.), 0)) %>%mutate_if(is.numeric, ~replace(., !is.finite(.), 0)) %>%arrange(id_zat, AN) %>%group_by(id_zat) %>%mutate(ylnlag =lag(ln_Y1),ylag =lag(Y1),Y_1 =NULL,Z_RESID =NULL ) %>%ungroup()# correction manuelle d'une ZAT avec un nombre d'emplois invraissemblable en 1996test1 <- base_data$id_zat ==860test2 <- base_data$AN ==1996base_data[ (test1 & test2), 'Y1'] <-100# suppression des "_" des noms de colonnes, ce caractère peut poser des # problèmes avec l'outil de modélisation STANold_names <-names(base_data)new_names <-gsub('_', '', old_names, fixed =TRUE)names(base_data) <- new_names# jointure avec les données sur les superZATbase_data <-left_join( base_data, dico_super, by =c('idzat'))# nous remettons certaines données dans leur échelle linéaire# pour faciliter l'estimation du modèle, nous allons convertir nos m2 de # surface de plancher en hectaresbase_data$SUMPLANCHERTRANSPAREA <- base_data$SUMPLANCHERTRANSPAREA /10000base_data$SUMPLANCHERMANUFAREA <- base_data$SUMPLANCHERMANUFAREA /10000base_data$SUMPLANCHERCOMMERCIALAREA <- base_data$SUMPLANCHERCOMMERCIALAREA /10000# on définit des limites supérieures pour ces variables avec des observations# extrèmes qui ont tendance à casser les modèlesbase_data$SUMPLANCHERCOMMERCIALAREA <-pmin(base_data$SUMPLANCHERCOMMERCIALAREA,150)base_data$SUMPLANCHERMANUFAREA <-pmin(base_data$SUMPLANCHERMANUFAREA,75)base_data$SUMPLANCHERTRANSPAREA <-pmin(base_data$SUMPLANCHERTRANSPAREA, 20)base_data$lnAreaCommercial <-ifelse( base_data$SUMPLANCHERCOMMERCIALAREA >0, log(base_data$SUMPLANCHERCOMMERCIALAREA), 0)base_data$lnAreaManuf <-ifelse( base_data$SUMPLANCHERMANUFAREA >0, log(base_data$SUMPLANCHERMANUFAREA), 0)base_data$lnAreaTransp <-ifelse( base_data$SUMPLANCHERTRANSPAREA, log(base_data$SUMPLANCHERTRANSPAREA), 0)base_data$lnTotal10Etages <-ifelse( base_data$Total10Etages, log(base_data$Total10Etages), 0)# nous convertissons la variable AN en interval régulier de type integerbase_data$AN_int <-as.integer(as.factor(base_data$AN))# nous recodons la variable nombre de bâtiments plus de 10 étages pour # pouvoir lui attribuer un effet non linéaire monotonique dans le modèlebase_data$Total10EtagesMO <-case_when( base_data$Total10Etages <3~ base_data$Total10Etages, base_data$Total10Etages >=3& base_data$Total10Etages <6~3, base_data$Total10Etages >=6& base_data$Total10Etages <11~4, base_data$Total10Etages >=11& base_data$Total10Etages <21~5, base_data$Total10Etages >=21~6)# definition de la transformation boxcoxboxcox <-function(y, lambda) {if (lambda ==0) {return(log(y)) } else {return((y^lambda -1) / lambda) }}# et sa fonction de transformation inverseboxcox.inv <-function(yt, lambda) {if (lambda ==0) {return(exp(yt)) } else {return((lambda * yt +1)^(1/ lambda)) }}
2.1 Structuration des données de Métro et train
Pour chaque ZAT, nous allons calculer pour chaque année le nombre de stations de métro en service dans un rayon de 500m autours de la ZAT. Cette variable permettra donc de discriminer les ZAT proches du métro, mais aussi de couvrir les cas d’ouverture de stations de métro (trois stations ouvertes après 1996).
Nous allons aussi ajouter les stations de la ligne bleue projetées pour l’horizon 2031. Elles seront donc incluses dans les projections.
Code
# chargement des données des stations de métrostation_metro <-st_read("data/reseau_transport.gpkg",layer ="reseau_metro_station", quiet =TRUE)# intersection avec les ZATdico_metro <- zat %>%st_buffer(250) %>%st_intersection(station_metro['nom']) %>%st_drop_geometry() %>%group_by(id_zat,nom) %>%slice_head(n =1) %>%ungroup() %>%group_by(id_zat) %>%summarise(nbmetro =n(),metro_names =paste0(nom, collapse =';') )base_data <-left_join(base_data, dico_metro, by =c('idzat'='id_zat'))base_data$nbmetro <-ifelse(is.na(base_data$nbmetro), 0, base_data$nbmetro)# on marque comme fermées ces trois stations ouvertes en 1996names_2007 <-c("Cartier", "De la Concorde", "Montmorency")for(n in names_2007){ base_data$nbmetro <-ifelse( (grepl(n,base_data$metro_names,fixed =TRUE)) & (base_data$AN <2007), 0, base_data$nbmetro )}# on ajoute les stations futures de la ligne bleuestation_bleu <-st_read('data/transport_projete.gpkg',layer ='metro_projete_l5_station_20221114',quiet =TRUE ) %>%st_transform(st_crs(zat)) %>%fselect(c('a_station')) %>%st_buffer(250) %>%st_intersection(zat['id_zat']) %>%st_drop_geometry() %>%group_by(id_zat,a_station) %>%slice_head(n =1) %>%ungroup() %>%group_by(id_zat) %>%summarise(total =n())bonus <- station_bleu$total[match(base_data$idzat, station_bleu$id_zat)]bonus <-ifelse(is.na(bonus), 0, bonus)base_data$nbmetro <-ifelse( (base_data$AN %in%c(2031, 2036, 2041)), base_data$nbmetro + bonus, base_data$nbmetro)
Nous effectuons la même opération pour les stations de train.
Code
# lecture des donnéestrain_stations <-st_read('data/reseau_transport.gpkg',layer ="stations_train",quiet = T) %>%st_transform(st_crs(zat)) %>%fselect(c('valide_de', 'valide_a', 'NO_STATION'))# Pour chaque station, nous devons déterminer les années pendant lesquelles# elle était validetrain_stations$valide_de <-ifelse( train_stations$valide_de ==0, 1980, train_stations$valide_de)train_stations$station_id <-1:nrow(train_stations)dico_train <- train_stations %>%st_buffer(250) %>%st_intersection(zat[c('id_zat')]) %>%st_drop_geometry() %>%group_by(id_zat, station_id) %>%slice_head(n =1) %>%ungroup()dico_train$start <-as.numeric(substr(as.character(dico_train$valide_de),1,4))dico_train$end <-as.numeric(substr(as.character(dico_train$valide_a),1,4))train_stations$start <-as.numeric(substr(as.character(train_stations$valide_de),1,4))train_stations$start <-as.numeric(substr(as.character(train_stations$valide_de),1,4))base_data$nbtrain <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_trains <- dico_train %>%filter( (start <= year) & (end >= year)) %>%group_by(id_zat) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_trains$nb[match(base_data$idzat, ok_trains$id_zat)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$nbtrain <-ifelse( base_data$AN == year, vals, base_data$nbtrain )}# validation cartographiquemaps <-lapply(years, function(year){ sub_data <-subset(base_data, (base_data$AN == year) & (base_data$nbtrain >0)) sub_train <-subset(train_stations, train_stations$start < year) sub_data$geometry <- zat$geometry[match(sub_data$idzat, zat$id_zat)] sub_data <-st_as_sf(sub_data, crs =st_crs(zat)) map <-tm_shape(sub_data) +tm_borders('black') +tm_shape(sub_train) +tm_dots('red') +tm_title(year)return(map)})tmap_arrange(maps, ncol =3)
Nous allons également ajouter les stations du REM. Il faut ajouter l’existant et le projeté.
Code
rem_stations_existing <-st_read('data/reseau_transport.gpkg',layer ='reseau_rem_station', quiet =TRUE)%>%st_transform(st_crs(zat)) %>%fselect(c('id')) %>%rename(any_of(c('NO_STATION'='id'))) %>%mutate(valide_de =2023,valide_a =99999999 )%>%st_set_geometry('geometry')rem_stations_projete <-st_read('data/transport_projete.gpkg',layer ='reseau_rem_projete_station_20231011',quiet =TRUE)%>%st_drop_geometry() %>%st_as_sf(coords =c('point_x', 'point_y'), crs =2950) %>%st_transform(st_crs(zat)) %>%fselect(c('id', 'antenne')) %>%rename(any_of(c('NO_STATION'='id'))) %>%mutate(NO_STATION = NO_STATION +200,valide_de =ifelse(grepl('A2', antenne, fixed =TRUE),2027,2025),valide_a =99999999 ) %>%mutate(antenne =NULL ) %>%st_set_geometry('geometry')rem_stations <-rbind(rem_stations_existing, rem_stations_projete)# Pour chaque station, nous devons déterminer les années pendant lesquelles# elle était validerem_stations$valide_de <-ifelse( rem_stations$valide_de ==0, 1980, rem_stations$valide_de)rem_stations$station_id <-1:nrow(rem_stations)dico_rem <- rem_stations %>%st_buffer(250) %>%st_intersection(zat[c('id_zat')]) %>%st_drop_geometry() %>%group_by(id_zat, station_id) %>%slice_head(n =1) %>%ungroup()dico_rem$start <-as.numeric(substr(as.character(dico_rem$valide_de),1,4))dico_rem$end <-as.numeric(substr(as.character(dico_rem$valide_a),1,4))rem_stations$start <-as.numeric(substr(as.character(rem_stations$valide_de),1,4))rem_stations$start <-as.numeric(substr(as.character(rem_stations$valide_de),1,4))base_data$nbrem <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_rem <- dico_rem %>%filter( (start <= year) & (end >= year)) %>%group_by(id_zat) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_rem$nb[match(base_data$idzat, ok_rem$id_zat)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$nbrem <-ifelse( base_data$AN == year, vals, base_data$nbrem )}# validation cartographiquemaps <-lapply(years, function(year){ sub_data <-subset(base_data, (base_data$AN == year) & (base_data$nbrem >0)) sub_rem <-subset(rem_stations, rem_stations$start < year)if(nrow(sub_data) >0){ sub_data$geometry <- zat$geometry[match(sub_data$idzat, zat$id_zat)] sub_data <-st_as_sf(sub_data, crs =st_crs(zat)) map <-tm_shape(sub_data) +tm_borders('black') +tm_shape(sub_rem) +tm_dots('red') +tm_title(year)return(map) }else{return(NULL) }})maps[sapply(maps, is.null)] <-NULLtmap_arrange(maps, ncol =3)
Nous allons calculer pour chaque super zone le nombre total de stations de métro de train et de REM. Ceci permettra notamment d’ajouter des effets régionaux pour régulariser le modèle et représenter un potentiel effet diffusé des stations.
Code
dico_rem_zone <- rem_stations %>%st_buffer(250) %>%st_intersection(super_zat[c('id')]) %>%rename(c('idzone'='id')) %>%st_drop_geometry() %>%group_by(idzone, station_id) %>%slice_head(n =1) %>%ungroup()dico_rem_zone$start <-as.numeric(substr(as.character(dico_rem_zone$valide_de),1,4))dico_rem_zone$end <-as.numeric(substr(as.character(dico_rem_zone$valide_a),1,4))base_data$totalremzone <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_rem <- dico_rem_zone %>%filter( (start <= year) & (end >= year)) %>%group_by(idzone) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_rem$nb[match(base_data$idzone, ok_rem$idzone)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$totalremzone <-ifelse( base_data$AN == year, vals, base_data$totalremzone )}dico_train_zone <- train_stations %>%st_buffer(250) %>%st_intersection(super_zat[c('id')]) %>%rename(c('idzone'='id')) %>%st_drop_geometry() %>%group_by(idzone, station_id) %>%slice_head(n =1) %>%ungroup()dico_train_zone$start <-as.numeric(substr(as.character(dico_train_zone$valide_de),1,4))dico_train_zone$end <-as.numeric(substr(as.character(dico_train_zone$valide_a),1,4))base_data$totaltrainzone <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_trains <- dico_train_zone %>%filter( (start <= year) & (end >= year)) %>%group_by(idzone) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_trains$nb[match(base_data$idzone, ok_trains$idzone)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$totaltrainzone <-ifelse( base_data$AN == year, vals, base_data$totaltrainzone )}dico_metro_zone <- super_zat %>%st_buffer(250) %>%st_intersection(station_metro['nom']) %>%st_drop_geometry() %>%group_by(id,nom) %>%slice_head(n =1) %>%ungroup() %>%group_by(id) %>%summarise(totalmetrozone =n(),metro_names_zone =paste0(nom, collapse =';') )base_data <-left_join(base_data, dico_metro_zone, by =c('idzone'='id'))base_data$totalmetrozone <-ifelse(is.na(base_data$totalmetrozone), 0, base_data$totalmetrozone)names_2007 <-c("Cartier", "De la Concorde", "Montmorency")for(n in names_2007){ base_data$totalmetrozone <-ifelse( (grepl(n,base_data$metro_names_zone,fixed =TRUE)) & (base_data$AN <2007), base_data$totalmetrozone -1, base_data$totalmetrozone )}station_bleu <-st_read('data/transport_projete.gpkg',layer ='metro_projete_l5_station_20221114',quiet =TRUE) %>%st_transform(st_crs(super_zat)) %>%fselect(c('a_station')) %>%st_buffer(250) %>%st_intersection(super_zat['id']) %>%st_drop_geometry() %>%group_by(id,a_station) %>%slice_head(n =1) %>%ungroup() %>%group_by(id) %>%summarise(total =n())bonus <- station_bleu$total[match(base_data$idzone, station_bleu$id)]bonus <-ifelse(is.na(bonus), 0, bonus)base_data$totalmetrozone <-ifelse( (base_data$AN %in%c(2031, 2036, 2041)), base_data$totalmetrozone + bonus, base_data$totalmetrozone)
2.2 Structuration des données d’autoroutes
Sur le même principe que la section prédécente, nous allons ajouter des variables sur l’accumulation de km linéaires d’autoroute par ZAT. Nous avons pour cela codé manuellement les années d’ouverture des différents tronçons autoroutiers en nous basant sur la source suivante : Répertoire des autoroutes du Québec.
Code
autoroutes <-st_read('data/autoroutes.gpkg', quiet = T) %>%st_transform(st_crs(zat))base_data$autoroutesKM <-0for (year inseq(1996, 2041, by =5)){ inter_zat <- autoroutes %>%filter(valide_de < year) %>%st_intersection(zat['id_zat']) %>%mutate(length_km =as.numeric(st_length(geom)) /1000 ) %>%st_drop_geometry() %>%group_by(id_zat) %>%summarise(autoroutesKM =sum(length_km) ) values <- inter_zat$autoroutesKM[match(base_data$idzat, inter_zat$id_zat)] values <-ifelse(is.na(values), 0, values) base_data$autoroutesKM <-ifelse( base_data$AN == year, values, base_data$autoroutesKM )}# nous créons une variable ordinale pour le kilométrage autoroutier. Elle # pourra être utilisée comme un effet monotone non linéaire dans le modèlebase_data$autoroutesMO <-case_when( base_data$autoroutesKM ==0~0, base_data$autoroutesKM >0& base_data$autoroutesKM <=2~1, base_data$autoroutesKM >2& base_data$autoroutesKM <=5~2, base_data$autoroutesKM >5& base_data$autoroutesKM <=10~3, base_data$autoroutesKM >10& base_data$autoroutesKM <=15~4, base_data$autoroutesKM >5~5)
On finalise la structuration des données ci-dessous et on exporte le tableau final pour ne pas avoir à relancer cette étape de structuration à chaque fois.
Code
# on termine la structuration des données et on les exportesbase_data$ANint <- (base_data$AN -1996) /5+1base_data$lnANint <-log(base_data$ANint)base_data$nomreg <-as.factor(base_data$nomreg )base_data$idzat <-as.factor(base_data$idzat )# nous retirons les données de 2021 du fait de leur mauvaise qualitébase_data2 <-subset(base_data, base_data$AN <2021)saveRDS(base_data, "artefacts/base_data_for_projection.rds")
2.3 Détermination de λ
Considérant la nature strictement positive et fortement asymétrique de la variable à modéliser, les chercheurs des Universités INRS et Laval avaient proposé d’utiliser la fonction de transformation log. Bien que ce choix soit justifié dans ce contexte, nous avions observé dans l’analyse des prédictions du modèle original que certaines ZAT voyaient leur nombre d’emplois exploser, probablement à cause de l’effet de la fonction log utilisée pour transformer le nombre d’emplois.
Nous proposons donc d’utiliser la transformation BoxCox qui est une généralisation de la fonction log. Elle nécessite de déterminer un paramètre \(\lambda\) contrôlant le caractère exponentiel de la transformation.
Puisque nous souhaitons obtenir un effet moins explosif que celui de la fonction log, nous devons trouver la bonne valeur de \(\lambda\) avec $ 0 $.
Pour cela, nous allons ajuster le modèle complet avec nos différentes transformations. Pour comparer les modèles, nous allons nous baser sur trois mesures: le rmse, le R2 bayésien et le graphique des posterior-predictive checks.
Le code ci-dessous prend beaucoup de temps, car nous devons ajuster 30 modèles très complexes.
Code
# déterminer la bonne valeur de lambdavals <-seq(0,0.3, 0.01)for (v in vals){ base_data2[[paste0('bcY',v)]] <-boxcox(base_data2$Y1, v)}
Une fois que les modèles sont ajustés, nous pouvons déterminer le meilleur d’entre eux en nous basant sur le R2 bayésien ainsi que le RMSE:
Code
vals <-seq(0,0.3, 0.01)models <-readRDS('artefacts/lambda_selection2.rds')comp_df <-data.frame(lambda = vals,r2 =sapply(models, function(x){median(x[[1]]$criteria$loo_R2)}),rmse =sapply(models, function(x){x[[4]]}))# on libère un peu de mémoire viverm(models)gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 5925362 316.5 16905717 902.9 11909575 636.1
Vcells 20389695 155.6 401317675 3061.9 487190703 3717.0
Il apparait assez clairement que le RMSE et le R2 s’améliorent significativement avec l’augmentation de \(\lambda\). Le gain semble maximal lorsque \(\lambda > 0.2\). Au delà de cette limite, les gains deviennent anecdotiques. Nous allons donc retenir le modèle avec \(\lambda = 0.2\)
2.4 Ajustement des modèles
Dans cette section, nous réalisons l’ajustement des modèles.
Nous ajustons trois modèles différents, chacun introduisant une petite variation sur la prise en compte de la tendance temporelle:
lnANint + (1 + lnANint | nomreg) : cette formulation implique une tendance qui ralentit (log sur le temps) et dont l’effet est différent d’une région à l’autre. Chaque région à son propre effet moyen et sa propre vitesse de croissance. Il existe aussi une corrélation par région entre l’effet moyen et la vitesse de croissance.
lnANint + (1 | nomreg) + (0 + lnANint | idzone/nomreg) : cette formulation implique une tendance qui ralentit (log sur le temps) et dont l’effet est différent d’une super-zat à l’autre et d’une région à l’autre, mais avec un principe hiérarchique entre les super-zat et les régions. Chaque région a aussi son propre effet moyen.
ANint + (1 | nomreg) + (0 + ANint | idzone/nomreg) : idem, mais sans la variation qui ralentit (pas de log sur le temps).
Les prédictions finales seront obtenues en croisant les prédictions de ces différents modèles par moyennisation bayésienne.
Code
# nous allons intégrer des priors de régularisation, ainsi que des priors # sur les effets positifs des effets fixespriors1 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 1)", class ="sd", coef ='lnANint', group ='nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)# on intègre la croissance au niveau des regionsfit3 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint ++ (1+ lnANint | nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors1,#control = list(adapt_delta = 0.95),refresh =10 )saveRDS(fit3, 'model_lnANINT_RE.rds')# on intègre la croissance au niveau des superzat et des zonespriors3 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 0.25)", class ="sd", coef ='lnANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='lnANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit5 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint + (1| nomreg) + (0+ lnANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,prior = priors3,iter =5000,control =list(adapt_delta =0.9),refresh =10 )saveRDS(fit5, 'model_lnANINT_REZONES.rds')priors4 <-prior("normal(0.25, 0.1)", class ="b", coef ='ANint') +prior("normal(0, 0.25)", class ="sd", coef ='ANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='ANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit6 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + ANint + (1| nomreg) + (0+ ANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,prior = priors4,iter =5000,control =list(adapt_delta =0.9),refresh =10 )saveRDS(fit6, 'model_ANINT_REZONES.rds')fit3 <-add_criterion(fit3, criterion =c('loo', 'waic', 'loo_R2'))fit6 <-add_criterion(fit6, criterion =c('loo', 'waic', 'loo_R2'))fit5 <-add_criterion(fit5, criterion =c('loo', 'waic', 'loo_R2'))save(fit3, fit5, fit6, file ='artefacts/fitted_gen2.rda')
2.4.1 Regard sur la qualité d’ajustement
Avant d’analyser le contenu des modèles, nous allons vérifier que les conditions d’application sont respectées et qu’ils parviennent à capter une partie significative du phénomène étudié.
Pour faciliter l’analyse, nous allons baser ce diagnostic sur un seul des trois modèles.
Code
load('artefacts/fitted_gen2.rda')# extraction d'échantillon des distributions postérieurs prédictives# application de la méthode de combinaison des prédictions par Bayesian # averaging.posterior_preds1 <-predict(fit3, summary =FALSE, cores =10, ndraws =2500)posterior_preds2 <-predict(fit5, summary =FALSE, cores =10, ndraws =2500)posterior_preds3 <-predict(fit6, summary =FALSE, cores =10, ndraws =2500)posterior_preds <-rbind(posterior_preds1, posterior_preds2, posterior_preds3)posterior_preds_trans <-boxcox.inv(posterior_preds, 0.2)pred_df <-data.frame(y = fit5$data$bcY0.2,y_real =boxcox.inv(fit5$data$bcY0.2, 0.2),ypred =colQuantiles(posterior_preds, probs =0.5),ypred_low =colQuantiles(posterior_preds, probs =0.2),ypred_high =colQuantiles(posterior_preds, probs =0.8),ypred_real =colQuantiles(posterior_preds_trans, probs =0.5),ypred_real_low =colQuantiles(posterior_preds_trans, probs =0.2),ypred_real_high =colQuantiles(posterior_preds_trans, probs =0.8))plot1 <-ggplot(pred_df) +geom_point(aes(x = y, y = ypred), size =0.1) +labs(x ='bc(emploi)', y ='bc(emploi prédit)')+geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed')+theme_bw() +coord_fixed()plot2 <-ggplot(pred_df) +geom_point(aes(x = y_real, y = ypred_real), size =0.1) +geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed') +labs(x ='emploi', y ='emploi prédit')+theme_bw() +coord_fixed()plot1plot2
La figure comparant les prédictions et les valeurs réelles montre un très bon ajustement, que ce soit dans l’échelle transformée ou dans l’échelle originale des données.
Les Root Median Squared Error calculés à l’échelle des régions varient sensiblement autour du RMSE moyen régional. Nous notons tout de même que nous avons des meilleures prédictions pour la Couronne Sud et des prédictions légèrement moins bonnes pour Montréal Centre et Ouest.
Les prédictions ont des intervalles de confiance qui grandissent avec le nombre d’emploi. Cependant, en dehors de quelques valeurs extrêmes, nous avons dans l’ensemble un ajustement très raisonable.
Code
pp_check(fit3, ndraws =25) +xlim(0,40)
Code
pp_check(fit5, ndraws =25) +xlim(0,40)
Code
pp_check(fit6, ndraws =25) +xlim(0,40)
La figure ci-dessus indique clairement que le modèle est capable de générer des données reproduisant la distribution des données originales.
Code
summary(fit3$criteria$loo_R2)
R2
Min. :0.8386
1st Qu.:0.8535
Median :0.8565
Mean :0.8564
3rd Qu.:0.8594
Max. :0.8709
Code
summary(fit5$criteria$loo_R2)
R2
Min. :0.8369
1st Qu.:0.8549
Median :0.8580
Mean :0.8579
3rd Qu.:0.8609
Max. :0.8745
Code
summary(fit6$criteria$loo_R2)
R2
Min. :0.8366
1st Qu.:0.8557
Median :0.8589
Mean :0.8587
3rd Qu.:0.8620
Max. :0.8751
Le R2 conditionnel du modèle est très élevé (0.86).
Il est intéressant de constater qu’une partie non négligeable de la variance provient tout de même des effets fixes (0.65). Ceci signifie que nos variables utilisées pour prédire le nombre d’emplois sont responsables d’une partie non négligeable de sa variation.
Parmis les effets aléatoires, il apparait clairement que l’effet propre des ZAT joue un rôle important dans la qualité de prédiction du modèle. Les effets régionaux et les ajustements locaux des tendances temporelles jouent un rôle plus mineur.
Code
library(gamlss, quietly =TRUE)resids <-residuals(fit5)resid_mod <-gamlss(resids[,1] ~1, family = TF)qts_obs <-quantile(resids[,1], probs =seq(0.01,0.99,0.01))qts_theo <-qTF(seq(0.01,0.99,0.01),mu = resid_mod$mu.coefficients[[1]],sigma =exp(resid_mod$sigma.coefficients[[1]]),nu =exp(resid_mod$nu.coefficients[[1]]) )df_plot <-data.frame(qts_obs = qts_obs,qts_theo = qts_theo)plot1 <-ggplot(df_plot) +geom_abline(slope =1, intercept =0, color ='red') +geom_point(aes(x = qts_obs, y = qts_theo)) +coord_fixed() +theme_bw() +labs(x ='quantiles observés', y ='quantiles attendus', title ='QQplot des résidus')resid_df <-data.frame(resids)plot2 <-ggplot(resid_df) +geom_histogram(aes(x = Estimate, y = ..density..),bins =50, color ='white') +stat_function(fun = dTF, args =list(mu = resid_mod$mu.coefficients[[1]],sigma =exp(resid_mod$sigma.coefficients[[1]]),nu =exp(resid_mod$nu.coefficients[[1]]) ),color ='red', size =1 ) +theme_bw() +labs(title ='Distribution des résidus')plot1plot2
GAMLSS-RS iteration 1: Global Deviance = 27869.27
GAMLSS-RS iteration 2: Global Deviance = 27543.16
GAMLSS-RS iteration 3: Global Deviance = 27477.71
GAMLSS-RS iteration 4: Global Deviance = 27465.26
GAMLSS-RS iteration 5: Global Deviance = 27462.93
GAMLSS-RS iteration 6: Global Deviance = 27462.5
GAMLSS-RS iteration 7: Global Deviance = 27462.42
GAMLSS-RS iteration 8: Global Deviance = 27462.41
GAMLSS-RS iteration 9: Global Deviance = 27462.4
GAMLSS-RS iteration 10: Global Deviance = 27462.4
Les résidus naturels du modèles suivent bien une distribution de Student tel que spécifié.
Partial autocorrelations of series 'x', by lag
1 2 3 4
-0.070 -0.297 -0.535 -1.180
Notre modèle parvient bien à retirer l’autocorrélation temporelle d’ordre 1, mais il reste de l’autocorrélation négative d’ordre 2 et 3. Ceci est problématique car certains paramètres du modèle pourraient ne pas être bien estimés. D’un point de vue pratique, cela signifie que pour une période temporelle t, les résidus du modèle tendent à être négativement associés avec les résidus à t-1. Cependant, considérant que nous n’avons que 6 périodes de temps dans les données, cet enjeu est beaucoup moins problématique.
Nous pouvons aussi regarder l’autocorrélation spatiale de nos résidus à chaque année.
Code
dfs <-split(df, df$AN)for(x in dfs){ year <- (unique(x$AN) -1 )*5+1996 x$idzat <-as.numeric(as.character(x$idzat)) x_space <-left_join(zat, x, by =c('id_zat'='idzat')) x_space <-subset(x_space, !is.na(x_space$resid)) queen <-poly2nb(x_space) WQueen <-nb2listw(queen, style ='W', zero.policy =TRUE) M <-moran.mc(x_space$resid, listw = WQueen, nsim =999) text <-paste0(year, ', moran I = ',round(M$statistic[[1]],2)) map <-tm_shape(x_space) +tm_polygons('resid', fill.scale =tm_scale_intervals(style ='jenks', n =7, values ='RdBu',midpoint =0 ), lwd =0.1) +tm_title(text)print(map)}
Les niveaux d’autocorrélation spatiale des résidus sont négligeables, ce qui est important pour ce type de variable fortement spatialement autocorrélée.
2.4.2 Validation croisée
Puisque nos modèles semblent être correctement ajustés, nous allons également tenté de déterminer dans quelle mesure ils sont pertinents pour effectuer des prédictions.
Pour cela, nous allons simplement réajuster des modèles identiques, mais en retirant les données de 2016. L’idée est alors de voir si nous aurions eu des prédictions réalistes comparativement à ce qui a été effectivement observé en 2016.
Code
base_data3 <-subset(base_data2, base_data2$AN <2016)priors1 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 1)", class ="sd", coef ='lnANint', group ='nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)priors3 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 0.25)", class ="sd", coef ='lnANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='lnANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)priors4 <-prior("normal(0.25, 0.1)", class ="b", coef ='ANint') +prior("normal(0, 0.25)", class ="sd", coef ='ANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='ANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit3_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint ++ (1+ lnANint | nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors1,#control = list(adapt_delta = 0.95),refresh =10 )fit4_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint + (1| nomreg) + (0+ lnANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors3,#control = list(adapt_delta = 0.95),refresh =10 )fit5_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + ANint + (1| nomreg) + (0+ ANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg ),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,prior = priors4,iter =3000,#control = list(adapt_delta = 0.95),refresh =10 )save(fit3_cv, fit4_cv, fit5_cv, file ='artefacts/fitted_gen1_cv.rda')
Maintenant que nous avons nos modèles, regardons leur prédiction en 2016
Code
load('artefacts/fitted_gen1_cv.rda')data_2016 <- base_data2 %>%fsubset(Y1 >20) %>%fsubset(idzat %in% fit5_cv$data$idzat)data_2016$idzat <-as.factor(as.integer(as.character(data_2016$idzat)))# on effectue les prédictions en utilisant les trois modèles. Nous allons # ainsi récupérer 1200 posterior samples que nous pourrons combiner pour # obtenir la prédiction finale par moyennisation bayésiennepred_2016_fit3 <-posterior_predict(fit3_cv, newdata = data_2016, ndraws =1200)pred_2016_fit4 <-posterior_predict(fit4_cv, newdata = data_2016, ndraws =1200)pred_2016_fit5 <-posterior_predict(fit5_cv, newdata = data_2016, ndraws =1200)pred_2016 <-rbind(pred_2016_fit3, pred_2016_fit4, pred_2016_fit5)pred_2016_trans <-boxcox.inv(pred_2016_fit5, 0.2)data_2016$pred <-colQuantiles(pred_2016_trans, probs =0.5)data_2016$pred_low <-colQuantiles(pred_2016_trans, probs =0.2)data_2016$pred_high <-colQuantiles(pred_2016_trans, probs =0.8)plot1 <-ggplot(data_2016) +geom_errorbar(aes(x = Y1,ymin = pred_low, ymax = pred_high),linewidth =0.1,alpha =0.3 ) +geom_point(aes(x = Y1, y = pred), size =0.4) +geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed') +labs(x ='emploi', y ='emploi prédit')+theme_bw() +coord_fixed()ggsave(plot = plot1, device ='png', filename ='cv_2016.png',path ='images/', width =20, height =20, units ='cm', dpi =200)knitr::include_graphics('images/cv_2016.png')
La figure ci-dessus illustre que le modèle parvient à faire une prédiction crédible pour 2016, surtout lorsque l’on considère l’incertitude avec notre intervalle de confiance.
Ci-dessous, nous calculons le Root Median Squared Error du modèle complet et du modèle de cross-validation. Si les deux RMdSE sont proches, cela signifie que le modèle complet n’a pas fait d’overfitting pour 2016.
Il est intéressant de constater que les RMdSE obtenus avec le modèle complet et avec le modèle sans les données de 2016 sont très similaires. Ceci signifique le modèle complet ne fait pas d’over-fitting, et que l’on peut raisonnablement espérer que sa qualité d’ajustement actuelle se répercutera sur les projections.
2.4.3 Analyse des tendances captées par les modèles
2.4.3.1 Tendances temporelles
Nous avons ajusté ci-dessus trois versions différentes du modèle incluant des variations dans la façon d’inclure la tendance temporelle. Dans le premier modèle, nous représentons le fait que chaque ZAT de la région suit une tendance centrale non linéaire (lnANint), et que cette tendance varie par région d’analyse (1 + lnANint | nomreg). Les figures ci-dessous illustrent les tendances modélisées.
Plus spécifiquement, elles illustrent le parcours d’une ZAT identique dans différentes régions. Elles doivent donc être interprétées de façon relative plutôt qu’absolue.
Code
load('artefacts/fitted_gen2.rda')# tendance apprises par la première version du modèle# pour cela, nous allons fixer toutes les valeurs d'une observation# et la duppliquer pour chaque zone, région et annéedf1 <- fit3$data[1,]pred_data <-expand.grid(nomreg =unique(fit3$data$nomreg),AN =seq(1996,2041,5))pred_data$AN_int <- (pred_data$AN -1996) /5+1for(var innames(df1)){if( !(var %in%c('nomreg', 'AN_int'))){ pred_data[[var]] <- df1[[var]] }}pred_data$lnANint <-log(pred_data$AN_int)pred_data$AN <- (pred_data$AN_int-1) *5+1996preds <-posterior_epred(fit3, pred_data, incl_autocor =FALSE)pred_data$pred <-colQuantiles(preds, probs =0.5)pred_data$pred_low <-colQuantiles(preds, probs =0.3)pred_data$pred_high <-colQuantiles(preds, probs =0.7)# représentation dans l'échelle transforméeplot1 <-ggplot(pred_data) +geom_line(aes(x = AN,y = pred_low,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN,y = pred_high,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN, y = pred, color = nomreg, group = nomreg), linewidth =1) +labs(x ='année',y ='bc(emplois)' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5)))# représentation dans l'échelle originalepreds_trans <-boxcox.inv(preds, 0.2)pred_data$pred_trans <-colQuantiles(preds_trans, probs =0.5)pred_data$pred_trans_low <-colQuantiles(preds_trans, probs =0.3)pred_data$pred_trans_high <-colQuantiles(preds_trans, probs =0.7)plot2 <-ggplot(pred_data) +geom_line(aes(x = AN,y = pred_trans_low,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN,y = pred_trans_high,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN, y = pred_trans, color = nomreg, group = nomreg), linewidth =1) +labs(x ='année',y ='emplois' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5)))# représentation dans une échelle relativepred_data <- pred_data %>%ungroup() %>%fsubset(AN ==1996) %>%fselect(c('nomreg', 'pred_trans')) %>%rename(c('pred_trans_1996'='pred_trans')) %>%right_join(pred_data) %>%mutate(pred_prt = pred_trans / pred_trans_1996 *100 )plot3 <-ggplot(pred_data) +geom_line(aes(x = AN, y = pred_prt, color = nomreg, group = nomreg))+labs(x ='année',y ='emplois (base 100 1996)' ) +theme_bw() +scale_x_continuous(breaks =c(seq(1996,2041,5)))plot1plot2plot3
Les trois figures ci-dessus permettent de constater que le modèle prévoit une croissance pour chacune des régions. La croissance la plus faible serait prévue pour Montréal Centre, suivi ensuite par Montréal Est et Montréal Ouest. Les couronnes connaitraient quant à elle les croissances les plus élevées.
Le second modèle va un peu plus loin en faisant varier le taux de croissance par région et par Super-Zat.
Cette seconde version du modèle permet de capturer l’hétérogénéité au niveau des Super-Zat. Nous pouvons notamment constater que certaines d’entre elles sont marquées par des croissances plus faibles, voire même des décroissances.
Le troisième modèle inclut la tendance temporelle avec un coefficient linéaire. Il produit donc des croissances plus importantes, surtout dans le long terme.
Il serait probablement plus pertinent de combiner les prédictions des trois modèles pour limiter l’effet problématique des décroissances.
Cette approche se nomme bayesian model averaging. Elle permet notamment de combiner l’incertitude des modèles ainsi que les éléments qu’ils ont appris respectivement.
2.5 Regard sur les projections
Nous pouvons à présent porter un regard sur les prédictions réalisées par les modèles.
Code
# réalisation des projections base_data$lnANint <-log(base_data$AN_int)base_data$ANint <-log(base_data$AN_int)base_data_pred <-subset(base_data, base_data$idzat %in% fit5$data$idzat)base_data_pred$idzat <- base_data_pred$idzat %>%as.character() %>%as.integer() %>%as.factor()preds <-pp_average(fit3, fit5, fit6, weights ='loo',summary =FALSE,ndraws =2500,cores =10,newdata = base_data_pred )# les draws issus des distributions postérieures doivent être retransformés # avec l'inverse de notre fonction BoxCox pour ensuite pouvoir calculer nos # intervalles de confiancepreds_trans <-boxcox.inv(preds, 0.2)# il se peut que nous obtenions des NA, ce qui est causé par des valeurs très# petites que nous remplaçons donc par la valeur 0preds_trans[is.na(preds_trans)] <-0# nous pouvons ensuite calculer pour chaque simulation (draw) les totaux# d'emploispreds_trans <-t(preds_trans)preds_trans <-data.frame(preds_trans)names(preds_trans) <-paste0('pred',1:ncol(preds_trans))pred_names <-names(preds_trans)preds2 <-cbind(preds_trans, base_data_pred[c('AN', 'idzat', 'idzone', 'nomreg')])# calcul des totaux pour l'ensemble du territoirepred_artm <- preds2 %>%fgroup_by(AN) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%fselect(pred_names) %>%as.matrix() %>%rowQuantiles(probs =c(0.1,0.5,0.9)) %>%as.data.frame() %>%cbind(unique(preds2$AN)) %>%frename(c('est_low', 'est_med', 'est_high', 'AN'))pred_artm$est_mean <- preds2 %>%fgroup_by(AN) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%fselect(pred_names) %>%as.matrix() %>%rowMeans2()# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by ='AN')pred_artm$total_reel <-ifelse(pred_artm$AN >2021, NA, pred_artm$total_reel)ggplot(pred_artm) +geom_line(aes(x = AN, y = total_reel)) +geom_point(aes(x = AN, y = total_reel)) +geom_ribbon(aes(x = AN, ymin = est_low, ymax = est_high), alpha =0.5) +geom_line(aes(x = AN, y = est_med), color ='red') +geom_point(aes(x = AN, y = est_med)) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2021, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)") +theme_bw() +scale_y_continuous(labels =function(x){x/1000})
Il est intéressant de constater qu’une fois que les prédictions sont agrégées, le modèle tend à sur-prédire le nombre d’emploi dans la région. Il faudra tenir compte de ce phénomème dans nos choix finaux de scénarios de projection.
Code
# calcul des totaux par régionpred_artm <- preds2 %>%fgroup_by(AN, nomreg) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%pivot(how ='longer', values = pred_names) %>%fgroup_by(AN, nomreg) %>%fsummarise(est_low =fquantile(value, probs =0.1),est_med =fquantile(value, probs =0.5),est_high =fquantile(value, probs =0.9) )# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN, nomreg) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by =c('AN', 'nomreg'))pred_artm$total_reel <-ifelse(pred_artm$AN >2016, NA, pred_artm$total_reel)pred_artm$nomreg <-as.factor(pred_artm$nomreg)ggplot(pred_artm) +geom_line(aes(x = AN, y = total_reel)) +geom_point(aes(x = AN, y = total_reel)) +geom_ribbon(aes(x = AN, ymin = est_low, ymax = est_high), alpha =0.5) +geom_line(aes(x = AN, y = est_med), color ='red') +geom_point(aes(x = AN, y = est_med)) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2016, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)") +theme_bw() +scale_y_continuous(labels =function(x){x/1000}) +facet_wrap(vars(nomreg), scales ='free')
Une fois nos projections éclatées par régions, nous pouvons constater que l’écart dans les projections au niveau régional s’explique vraisemblablement par notre sur-prédiction dans les régions Montréal Est et Ouest. En effet, ces deux régions ont connu une baisse significative de leur nombre d’emplois entre 2006 et 2011.
Code
# calcul des totaux par superzatpred_artm <- preds2 %>%fgroup_by(AN, nomreg, idzone) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%pivot(how ='longer', values = pred_names) %>%fgroup_by(AN, nomreg, idzone) %>%fsummarise(est_low =fquantile(value, probs =0.1),est_med =fquantile(value, probs =0.5),est_high =fquantile(value, probs =0.9) )# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN, nomreg, idzone) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by =c('AN', 'nomreg', 'idzone'))pred_artm$total_reel <-ifelse(pred_artm$AN >2016, NA, pred_artm$total_reel)pred_artm$nomreg <-as.factor(pred_artm$nomreg)pred_artm$idzone <-as.character(pred_artm$idzone)for(reg inunique(pred_artm$nomreg)){ sub_data <-subset(pred_artm, pred_artm$nomreg == reg) this_plot <-ggplot(sub_data) +geom_ribbon(aes(x = AN,ymin = est_low,ymax = est_high, group = idzone), alpha =0.25) +geom_line(aes(x = AN, y = total_reel,group = idzone, color = idzone),linewidth =1) +geom_point(aes(x = AN, y = total_reel, color = idzone)) +geom_line(aes(x = AN, y = est_med, group = idzone, color = idzone),linetype ='dashed',linewidth =1) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2016, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)",title = reg ) +theme_bw() +scale_y_continuous(labels =function(x){x/1000})print(this_plot)}
Les figures ci-dessus illustrent comment le modèle lisse les projections pour conserver les tendances observées à long terme. Pour chacune des régions, nous observons que les tendances par super-zat projetées semblent poursuivre les tendances observées dans les données historiques. Nous remarquons tout de même quelques super-zat dans lesquelles les prédictions tendent à être assez éloignées des valeurs effectivement observées dans les données historiques.
2.6 Projections
L’objectif principal du modèle développé ici est d’estimer le nombre d’emplois par ZAT pour les horizons 2026, 2031, 2036 et 2041. L’idée est notamment de se doter d’un scénario de référence, ainsi que d’un scénario optimiste et d’un scénario pessimiste.
2.6.1 Déterminer les paramètres de projection
Le modèle de régression bayésien que nous avons estimé permet de déterminer la distribution des prédictions que nous effectuons. Cette prédiction comporte donc l’ensemble de l’incertitude sur nos projections. L’idéal serait donc de conserver cette incertitude complète et de la transférer aux modèles de déplacements. Cependant, puisque c’est peu réaliste compte tenu des temps de calcul de ces modèles, nous allons essayer d’extraire des point wise estimates de nos projections.
Pour sélectionner les meilleurs point wise estimates, nous allons regarder pour le modèle sans les données de 2016, lesquels auraient donné les meilleurs résultats.
Code
library(bayestestR)library(ggpubr)load('artefacts/fitted_gen1_cv.rda')# extraction des prédictions de notre modèledata_2016 <- base_data2 %>%fsubset(Y1 >20) %>%fsubset(idzat %in% fit5_cv$data$idzat)data_2016$idzat <-as.factor(as.integer(as.character(data_2016$idzat)))preds_comp1 <-posterior_predict(fit3_cv, newdata = data_2016,summary =FALSE, cores =10)preds_comp2 <-posterior_predict(fit4_cv, newdata = data_2016,summary =FALSE, cores =10)preds_comp3 <-posterior_predict(fit5_cv, newdata = data_2016,summary =FALSE, cores =10)pred_2016 <-rbind(preds_comp1, preds_comp2, preds_comp3)pred_2016_trans <-boxcox.inv(pred_2016, 0.2)# sélection d'un ensemble de prédictions pour visualiser leur distributionobs <-sample(1:ncol(pred_2016_trans), size =25, replace =FALSE)plots <-lapply(obs, function(i){ x <- pred_2016_trans[,i] target <- data_2016$Y1[[i]] x <- x[x <quantile(x, probs =0.99)] x <- x[x >0] plot <-ggplot() +geom_histogram(aes(x = x), color ='white', bins =50) +geom_vline(xintercept = target, color ='red',linetype ='dashed', size =1) +labs(x ='', y ='')return(plot)})ggarrange(plotlist = plots, nrow =5, ncol =5)
La figure ci-dessus indique clairement que les distributions des prédictions d’emplois tendent à être asymétriques. Il est donc important de choisir quel percentile de la distribution nous allons utiliser pour nos prédiction. Nous allons notamment regarder les valeurs aux alentours de la médiane.
Code
mae <-function(x,y){mean(abs(x-y))}mean_error <-function(x,y){mean(x-y)}rmse <-function(x,y){sqrt(mean((x-y)**2))}smape <-function(x,y){mean( abs(abs(x-y) / ((x+y)/2)) )}pcts <-seq(0.25,0.75, by =0.025)preds <-colQuantiles(pred_2016_trans, probs = pcts)df_scores2 <-sapply(1:ncol(preds), function(i){ x <- preds[,i] y <- data_2016$Y1 reg <-lm(x ~0+y) vals <-c(mean_error(x,y),mae(x,y),rmse(x,y),smape(x,y), reg$coefficients[[1]] )return(vals)})df_scores2 <-data.frame(df_scores2)df_scores2$metric <-c('erreur moyenne','erreur absolue moyenne','rmse','smape', 'coeff' )df_scores2 <-pivot(df_scores2, how ='longer', ids ='metric')df_scores2$pct <-gsub('X','',df_scores2$variable,fixed =TRUE)df_scores2$pct <- pcts[as.integer(df_scores2$pct)]df_scores2$metric <-as.factor(df_scores2$metric)ggplot(df_scores2) +geom_line(aes(x = pct,y = value )) +facet_wrap(vars(metric), scales ='free')
Il semblerait qu’une valeur comprise entre les 50e et 60e percentile permettrait d’obtenir les prédictions donnant un bon compromis entre nos différentes métriques d’ajustement.
Code
df_scores2_wide <- df_scores2 %>%pivot(how ='wider', values ='value', names ='metric')# le meilleur coefficient est celui le plus proche de 1test <-abs(df_scores2_wide$coeff-1) ==min(abs(df_scores2_wide$coeff-1))best_coeff <- df_scores2_wide$pct[test] # la meilleure erreur absolue moyenne et proche de 0test <- df_scores2_wide$`erreur absolue moyenne`==min(df_scores2_wide$`erreur absolue moyenne`)best_mae <- df_scores2_wide$pct[test] # la meilleur erreur moyenne est la plus proche de 0test <-abs(df_scores2_wide$`erreur moyenne`) ==min(abs(df_scores2_wide$`erreur moyenne`))best_me <- df_scores2_wide$pct[test]# le meilleur RMSE est le plus petittest <- df_scores2_wide$rmse ==min(df_scores2_wide$rmse)best_rmse <- df_scores2_wide$pct[test]# le meilleur smape est le plus petittest <- df_scores2_wide$smape ==min(df_scores2_wide$smape)best_smape <- df_scores2_wide$pct[test]df_best_pcts <-data.frame(score =c('coeff', 'mae', 'me', 'rmse', 'smape'),best_ctl =c(best_coeff, best_mae, best_me, best_rmse, best_smape))knitr::kable(df_best_pcts)
score
best_ctl
coeff
0.625
mae
0.550
me
0.575
rmse
0.575
smape
0.500
Nous proposons de retenir la valeur de 0.58 pour le quantile de prédiction. Elle offre un bon compromis entre nos différentes mesures d’erreur.
Nous devons encore choisir nos bornes supérieures et inférieures qui délimiteront nos scénarios pessimiste et optimiste.
Pour notre scénario pessimiste, nous proposons de retenir le percentile garantissant que la croissance générale de la région ne passe pas en dessous des totaux de 2016.
il semblerait que le 54e percentile permette d’éviter une projection de décroissance. Nous pourrions donc l’utiliser comme notre valeur basse. Comme valeur haute, nous proposons de prendre un écart de 10 points de percentile avec la valeur basse.
Code
qtls2 <-colQuantiles(preds_total, probs =c(0.54,0.58,0.64))base_data_pred$pred <- qtls2[,2]base_data_pred$low <- qtls2[,1]base_data_pred$high <- qtls2[,3]this_preds <- base_data_pred %>%group_by(AN) %>%summarise(total_y =sum(Y1),total_med =sum(pred),total_low =sum(low),total_high =sum(high) )real <-subset(this_preds, this_preds$AN <2021)real_2016 <-max(real$total_y)proj <- this_preds %>%filter(AN >2016) %>%add_row(AN =2016,total_med = real_2016,total_low = real_2016,total_high = real_2016 )this_plot <-ggplot() +geom_line(data = real, mapping =aes(x = AN, y = total_y)) +geom_point(data = real, mapping =aes(x = AN, y = total_y)) +geom_ribbon(data = proj, mapping =aes(x = AN, ymin = total_low,ymax = total_high ),alpha =0.4) +geom_line(data = proj, mapping =aes(x = AN, y = total_med)) +geom_point(data = proj, mapping =aes(x = AN, y = total_med)) +geom_vline(xintercept =2016, linetype ='dashed') +theme_bw() +scale_x_continuous(breaks =seq(1996, 2041, by =5)) +labs(x ='année', y ="volume d'emplois")this_plot
Code
# par regionthis_preds <- base_data_pred %>%group_by(AN, nomreg) %>%summarise(total_y =sum(Y1),total_med =sum(pred),total_low =sum(low),total_high =sum(high) )real <-subset(this_preds, this_preds$AN <2021)real_2016 <-subset(real, real$AN ==2016)real_2016$total_med <- real_2016$total_yreal_2016$total_low <- real_2016$total_yreal_2016$total_high <- real_2016$total_yproj <- this_preds %>%filter(AN >2016) %>%rbind(real_2016)this_plot <-ggplot() +geom_line(data = real, mapping =aes(x = AN, y = total_y, group = nomreg, color = nomreg)) +geom_point(data = real, mapping =aes(x = AN, y = total_y, color = nomreg)) +geom_ribbon(data = proj, mapping =aes(x = AN, ymin = total_low,ymax = total_high, group = nomreg),alpha =0.2) +geom_line(data = proj, mapping =aes(x = AN, y = total_med ,group = nomreg, color = nomreg)) +geom_point(data = proj, mapping =aes(x = AN, y = total_med, color = nomreg))this_plot
2.6.2 Appliquer les projections
Nous pouvons maintenant appliquer les projections en appliquant nos différentes observations.
Une question reste à aborder : quel sera l’impact de l’ouverture des stations de REM sur l’emploi. Pour le moment, elles ont été codées comme des stations de train. Or, nous avons pu observer dans les effets fixes que le train avait un impact proche de zéro.
Nous pourrions formuler l’hypothèse qu’une station de REM aura un impact plus proche de l’ouverture d’une station de métro que de train. Dans ce cas, nous pourrions compter une station de REM comme une fraction de métro. Nous allons coder l’impact d’une station de REM comme la moitiée de celui d’une station de métro.
Code
load('artefacts/fitted_gen2.rda')# on s'assure la reproductibilité des samplesset.seed(666)# réalisation des projections base_data$lnANint <-log(base_data$AN_int)base_data$ANint <- base_data$AN_intbase_data_pred <-subset(base_data, base_data$idzat %in% fit5$data$idzat)base_data_pred$idzat <- base_data_pred$idzat %>%as.character() %>%as.integer() %>%as.factor()#_______________# les stations de REM sont traitées comme la moitiée des stations de métrobase_data_pred2 <- base_data_predbase_data_pred2$nbmetro <- base_data_pred$nbmetro + base_data_pred$nbrem *0.5base_data_pred2$totalmetrozone <- base_data_pred$totalmetrozone + base_data_pred$totalremzone *0.5#_______________# Combinaison des projections des trois modèlespreds1 <-posterior_predict(fit3, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds2 <-posterior_predict(fit5, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds3 <-posterior_predict(fit6, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds_total <-rbind(preds1, preds2, preds3)# les draws issus des distributions postérieures doivent être retransformés # avec l'inverse de notre fonction BoxCox pour ensuite pouvoir calculer nos # intervalles de confiancepreds_total_trans <-boxcox.inv(preds_total, 0.2)# il se peut que nous obtenions des NA, ce qui est causé par des valeurs très# petites que nous remplaçons donc par la valeur 0preds_total_trans[is.na(preds_total_trans)] <-0qtls <-colQuantiles(preds_total_trans, probs =c(0.54,0.58,0.64))gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 23808325 1271.6 46312850 2473.4 46312850 2473.4
Vcells 1820340595 13888.1 3045465411 23235.1 2537821176 19362.1
Code
qtls <-data.frame(qtls)names(qtls) <-c('s1_low', 's1_med', 's1_high')projections <-cbind( base_data_pred[c('AN', 'idzat', 'Y1', 'idzone', 'nomreg')], qtls)# on marque les données projetée avec NA pour bien distinguer les données # originales des projections projections$Y1 <-ifelse( projections$AN >2016, NA, projections$Y1)saveRDS(projections, 'artefacts/projections.rds')# nous exportons les projections par ZAT en CSVwrite.csv(projections, 'projection_emplois_ZAT_ARTM.csv', row.names = F)# nous exportons aussi des échantillons des valeurs postérieures pour permettre # des agrégations probabilistes# nous allons retenir uniquement 1500 échantillons postérieursids <-sample(1:nrow(preds_total_trans), size =1500, , replace = F)posteriors_export <- preds_total_trans[ids,] %>%t() %>%data.frame() names(posteriors_export) <-paste0("sample_", 1:ncol(posteriors_export))posteriors_export <-cbind(base_data_pred2[c('AN', 'idzat', 'Y1', 'idzone', 'nomreg')], posteriors_export)write.csv(posteriors_export, 'projection_emplois_ZAT_ARTM_samples.csv', row.names = F)
Terminons en montrant comment utiliser les échantillons postérieurs pour calculer des intervalles de confiances pour des projections agrégées.
Code
post_samples <-read.csv('projection_emplois_ZAT_ARTM_samples.csv')#_______________________________# example avec agrégation régionale# on extrait les données historiques connuespart1 <- post_samples %>%fsubset(AN <2021) %>%fgroup_by(AN) %>%fsummarise(emploi =fsum(Y1) )# on extrait un intervalle de confiance pour les données # projetée (15% - 85%). Pour cela, il faut calculer la somme régionale # d'emploi pour chacun des échantillon bayésien, puis effectuer la # l'agrégation pour obtenir les intervales de confiances# on identifie les colonnes contenant les échantillonssamples_cols <-names(post_samples)samples_cols <- samples_cols[grepl('sample_', samples_cols)]part2 <- post_samples %>%fsubset(AN >=2021) %>%# on applique uniquement aux données projetéesfselect(c("AN", samples_cols)) %>%# on sélectionner les colonnes nécessairespivot(how ='longer', values = samples_cols) %>%# on pivote au format longfgroup_by(AN,variable) %>%# pour chaque année, et chaque échantillonfsummarise( # on calcule la somme des emplois dan la région emploi =fsum(value) ) %>%ungroup() %>%fgroup_by(AN) %>%# les sommes (1500) obtenues sont organisées par fsummarise( # années afin de pouvoir calculer les intervalles low_est =fquantile(emploi, probs =0.15), # de confiancemed_est =fquantile(emploi, probs =0.5),high_est =fquantile(emploi, probs =0.85) )# pour assurer la continuité de la figure finale, nous devons # ajouter une ligne dans le tableau part2part2 <-rbind(part2, data.frame(AN =2016, low_est =last(part1$emploi),med_est =last(part1$emploi),high_est =last(part1$emploi) ))ggplot() +geom_line(data = part1, mapping =aes(x = AN, y = emploi)) +geom_point(data = part1, mapping =aes(x = AN, y = emploi)) +geom_ribbon(data = part2, mapping =aes(x = AN, ymin = low_est,ymax = high_est ),alpha =0.4) +geom_line(data = part2, mapping =aes(x = AN, y = med_est)) +geom_point(data = part2, mapping =aes(x = AN, y = med_est)) +geom_vline(xintercept =2016, linetype ='dashed') +theme_bw() +scale_x_continuous(breaks =seq(1996, 2041, by =5)) +labs(x ='année', y ="volume d'emplois")
Code source
---title: "Modèle de projection des emplois"subtitle: "Format HTML"author: "J. GELB -- RVD -- ARTM"date: '2025-06-03'format: html: code-fold: show code-tools: true css: ["css/theme_ARTM.css"] embed-resources: true lang: fr number-sections: true toc: true toc-collapsed: true toc-depth: 3 toc-float: true toc-location: 'left'---Objectif: Ce document présente la méthodologie détaillée d'ajustement du modèle de projection des emplois de l'ARTM. Ce modèle s'appuie sur un travail préliminaire des chercheurs Cédric Brunelle (INRS) et Jean Dubé (Université Laval).Dernière compilation: `r Sys.Date()`# IntroductionLe modèle a pour objectif de prédire le nombre de personnes actives de 15 ans ou plus travaillant dans les différentes ZAT du territoire de L'ARTM. Dans le reste du document, cette variable sera aussi appelée *nombre d'emplois*.Pour chaque ZAT, nous disposons des données historiques de 1996, 2001, 2006, 2011, 2016 et 2021. Il s'agit donc d'une structure en panel, car nous avons de multiples valeurs qui se suivent dans le temps pour une même observation. Il s'agit aussi d'un panel large, car nous avons plus de ZAT que d'années.## Formulation originale du modèleLe modèle original a été réalisé dans le logiciel Stata avec la fonction `xtabond`. Elle implémente l'estimateur d'Arellano-Bond basé sur la méthode des moments.Ce modèle introduit la dimension temporelle dans l'équation de régression en ajoutant des lags de la variable dépendante. Cependant, puisque nous sommes dans un contexte de données par panel, cela soulève un problème d'endogénéité. En effet, nous risquons d'avoir des corrélations importantes entre les effets propres à chaque ZAT et la variable dépendante spatialement décalée.$$y_{i t}=\beta_1 y_{i, t-1}+\beta_2 x_{i t}+\alpha_i+\varepsilon_{i t}$$avec : * $\alpha_i$ une erreur propre à la ZAT *i** $\epsilon_{it}$ un terme d'erreur distribué normalement et idiosyncratiquePour un tel modèle, l'estimation par OLS pose un problème majeur, car les deux termes d'erreurs $\alpha_i$ et $\epsilon_{it}$ sont confondus. La corrélation entre $y_{it-1}$ et $\alpha_i$ pose alors un problème d'endogénéité.Pour contourner ce problème, l'estimateur d'Arellano-Bond modélise plutôt la différence des variables décalées temporellement.$$\Delta y_{i t}=\beta_1 \Delta y_{i, t-1}+\beta_2 \Delta x_{i t}+\Delta \varepsilon_{i t}$$Le tout, combiné avec une condition d'estimation:$$E\left[y_{i, t-s} \cdot \Delta \varepsilon_{i t}\right]=0 \quad \text { for } s \geq 2$$En d'autres termes, l'endogénéité est contrôlée par des versions décalées plus profondes de la variable dépendante. Une telle approche n'est pas idéale dans notre contexte, car nous avons assez peu de réplications temporelles dans nos données. Concrètement, utiliser un lag minimal de 2 impose de nous passer de deux périodes dans nos données, soit près du tiers des données.## Reformulation proposéeNous proposons de reformuler ce modèle pour gagner un peu de flexibilité et limiter certains effets problématiques observés dans les résultats des projections du modèle original.Notre formulation suppose que le nombre d'emplois par ZAT est la somme de cinq éléments: 1. L'effet propre de chaque ZAT. En effet, chaque ZAT a son volume d'emploi original, une moyenne qui lui est propre. Nous ne disposons que de 5 observations par ZAT, nous ne pouvons donc pas estimer cet effet comme un effet fixe. Nous allons donc l'estimer comme un effet aléatoire, soit une déviation de l'intercept général propre à chaque ZAT et provenant d'une distribution normale. Cette estimation repose sur la mutualisation de l'information entre ZAT.2. Les effets fixes, soit la contribution d'un ensemble de termes tels que la présence d'une station de métro, de train, d'un hôpital, etc. Ces effets peuvent changer dans le temps et sont identiques pour chaque ZAT. Cependant, puisque chaque ZAT à son effet propre et que la fonction box-cox (utilisée pour transformer la variable Y, voir plus loin) produit des effets quasi multiplicatifs, ces effets fixes ont des impacts plus grands pour des ZAT déjà plus grandes.3. Les effets non linéaires. Ils ont un rôle similaire aux effets fixes, mais un lien non linéaire est estimé plutôt qu'un simple coefficient, notamment pour capturer des effets de seuils.4. La tendance temporelle. Il s'agit d'une tendance régionale indépendante de la ZAT. Nous avons notamment constaté que l'emploi a tendance à croitre dans la région à cause de la croissance économique et de la population. Il s'agit donc d'une tendance d'arrière-plan. Nous faisons varier cette tendance par sous-région pour capter des dynamismes différenciés. Nous nous attendons à avoir un effet positif (croissance), que nous allons modéliser avec une transformation logarithmique pour représenter une croissance qui ralentit.5. Notre formulation suppose que chaque ZAT est caractérisée par une structure temporelle résiduelle suivant un processus de type ARMA(1). Ce processus tient compte des dynamiques de long terme et de court terme au sein d'une ZAT. Il s'agit d'un processus stationnaire qui n'est donc pas responsable de l'éventuelle croissance au sein de la ZAT. Notons aussi que la variance n'est pas supposée homogène, mais est autorisée à varier par grande région.Chaque ZAT combine donc deux termes temporels: une tendance régionale atténuée localement et une tendance locale stationnaire. $$\begin{split} bc(Y_{zt}, \lambda) &\sim MVStudent(\mu_{zt}, \Sigma, \nu) \\ \mu_{zt} &= \beta_0 +\beta X_{zt} + \sum^k_{j = 1} f_j(\zeta_{ztj}) + \eta_{z} + (\beta_{an} + \gamma_{reg})\times AN \\ \eta &\sim Normal(0,\sigma_z) \\ \gamma &\sim Normal(0,\sigma_{reg}) \\ \Sigma &= diag(\sigma) \cdot R_{ARMA(1,1)} \cdot diag(\sigma)\\ \sigma &= exp(\alpha_0 + \alpha X)\end{split}$$# Chargement des donnéesDans cette première section, nous chargeons simplement les données pour notre modèle et nous structurons certaines variables complémentaires.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falselibrary(sf, quietly =TRUE)library(spdep, quietly =TRUE)library(dplyr, quietly =TRUE)library(collapse, quietly =TRUE)library(tmap, quietly =TRUE)library(DescTools, quietly =TRUE)library(ggplot2, quietly =TRUE)library(tidyr, quietly =TRUE)library(haven, quietly =TRUE)library(brms, quietly =TRUE)library(tidybayes, quietly =TRUE)library(DHARMa, quietly =TRUE)library(matrixStats, quietly =TRUE)# chargement des ZAT et super ZATzat <-st_read('data/limite_analyse.gpkg',layer ="zat_zones_2018", quiet =TRUE) %>%mutate(id_zat = id ) %>%fselect(c("id_zat")) %>%st_set_geometry('geometry')super_zat <-st_read('data/limite_analyse.gpkg',layer ="zat_super_zones_2018", quiet =TRUE) %>%st_transform(st_crs(zat)) %>%fselect(c("id")) %>%st_set_geometry('geometry')# Chargement des régions d'analyseregions <-st_read("data/enquete.gpkg", layer ="empm23_region_analyse") %>%mutate(nomreg2 = nom,nomreg =ifelse(nom =="MONTREAL CENTRE-VILLE", "MONTREAL CENTRE", nom) ) %>%st_transform(st_crs(zat)) %>%st_set_geometry('geometry')# Création d'un dictionnaire associant chaque ZAT à sa régiondico_reg_ori <- zat %>%st_centroid() %>%st_intersection(regions['nomreg2'])# Idem, mais pour les super-zatdico_super <- zat %>%st_centroid() %>%st_intersection(super_zat) %>%st_intersection(regions['nomreg']) %>%st_drop_geometry() %>%fselect(c('id_zat', 'id', 'nomreg')) %>%rename(c('idzat'='id_zat', 'idzone'='id'))# lecture du fichier de données structurées partagé par les chercheurs des # université Laval et INRSbase_data <-read_dta('data/ScenarioProjection_Etages.dta')# suppression d'une colonne inutilebase_data$`_Y1`<-NULL# création de variables complémentaires et nettoyage de certaines colonnesbase_data <- base_data %>%mutate(ln_Area_Commercial =log(Total_Area_Commercial), ln_Area_Bureau =log(Total_Area_Bureau) ) %>%mutate_if(is.numeric, ~replace(., is.na(.), 0)) %>%mutate_if(is.numeric, ~replace(., !is.finite(.), 0)) %>%arrange(id_zat, AN) %>%group_by(id_zat) %>%mutate(ylnlag =lag(ln_Y1),ylag =lag(Y1),Y_1 =NULL,Z_RESID =NULL ) %>%ungroup()# correction manuelle d'une ZAT avec un nombre d'emplois invraissemblable en 1996test1 <- base_data$id_zat ==860test2 <- base_data$AN ==1996base_data[ (test1 & test2), 'Y1'] <-100# suppression des "_" des noms de colonnes, ce caractère peut poser des # problèmes avec l'outil de modélisation STANold_names <-names(base_data)new_names <-gsub('_', '', old_names, fixed =TRUE)names(base_data) <- new_names# jointure avec les données sur les superZATbase_data <-left_join( base_data, dico_super, by =c('idzat'))# nous remettons certaines données dans leur échelle linéaire# pour faciliter l'estimation du modèle, nous allons convertir nos m2 de # surface de plancher en hectaresbase_data$SUMPLANCHERTRANSPAREA <- base_data$SUMPLANCHERTRANSPAREA /10000base_data$SUMPLANCHERMANUFAREA <- base_data$SUMPLANCHERMANUFAREA /10000base_data$SUMPLANCHERCOMMERCIALAREA <- base_data$SUMPLANCHERCOMMERCIALAREA /10000# on définit des limites supérieures pour ces variables avec des observations# extrèmes qui ont tendance à casser les modèlesbase_data$SUMPLANCHERCOMMERCIALAREA <-pmin(base_data$SUMPLANCHERCOMMERCIALAREA,150)base_data$SUMPLANCHERMANUFAREA <-pmin(base_data$SUMPLANCHERMANUFAREA,75)base_data$SUMPLANCHERTRANSPAREA <-pmin(base_data$SUMPLANCHERTRANSPAREA, 20)base_data$lnAreaCommercial <-ifelse( base_data$SUMPLANCHERCOMMERCIALAREA >0, log(base_data$SUMPLANCHERCOMMERCIALAREA), 0)base_data$lnAreaManuf <-ifelse( base_data$SUMPLANCHERMANUFAREA >0, log(base_data$SUMPLANCHERMANUFAREA), 0)base_data$lnAreaTransp <-ifelse( base_data$SUMPLANCHERTRANSPAREA, log(base_data$SUMPLANCHERTRANSPAREA), 0)base_data$lnTotal10Etages <-ifelse( base_data$Total10Etages, log(base_data$Total10Etages), 0)# nous convertissons la variable AN en interval régulier de type integerbase_data$AN_int <-as.integer(as.factor(base_data$AN))# nous recodons la variable nombre de bâtiments plus de 10 étages pour # pouvoir lui attribuer un effet non linéaire monotonique dans le modèlebase_data$Total10EtagesMO <-case_when( base_data$Total10Etages <3~ base_data$Total10Etages, base_data$Total10Etages >=3& base_data$Total10Etages <6~3, base_data$Total10Etages >=6& base_data$Total10Etages <11~4, base_data$Total10Etages >=11& base_data$Total10Etages <21~5, base_data$Total10Etages >=21~6)# definition de la transformation boxcoxboxcox <-function(y, lambda) {if (lambda ==0) {return(log(y)) } else {return((y^lambda -1) / lambda) }}# et sa fonction de transformation inverseboxcox.inv <-function(yt, lambda) {if (lambda ==0) {return(exp(yt)) } else {return((lambda * yt +1)^(1/ lambda)) }}```## Structuration des données de Métro et trainPour chaque ZAT, nous allons calculer pour chaque année le nombre de stations de métro en service dans un rayon de 500m autours de la ZAT. Cette variable permettra donc de discriminer les ZAT proches du métro, mais aussi de couvrir les cas d'ouverture de stations de métro (trois stations ouvertes après 1996).Nous allons aussi ajouter les stations de la ligne bleue projetées pour l'horizon 2031. Elles seront donc incluses dans les projections.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# chargement des données des stations de métrostation_metro <-st_read("data/reseau_transport.gpkg",layer ="reseau_metro_station", quiet =TRUE)# intersection avec les ZATdico_metro <- zat %>%st_buffer(250) %>%st_intersection(station_metro['nom']) %>%st_drop_geometry() %>%group_by(id_zat,nom) %>%slice_head(n =1) %>%ungroup() %>%group_by(id_zat) %>%summarise(nbmetro =n(),metro_names =paste0(nom, collapse =';') )base_data <-left_join(base_data, dico_metro, by =c('idzat'='id_zat'))base_data$nbmetro <-ifelse(is.na(base_data$nbmetro), 0, base_data$nbmetro)# on marque comme fermées ces trois stations ouvertes en 1996names_2007 <-c("Cartier", "De la Concorde", "Montmorency")for(n in names_2007){ base_data$nbmetro <-ifelse( (grepl(n,base_data$metro_names,fixed =TRUE)) & (base_data$AN <2007), 0, base_data$nbmetro )}# on ajoute les stations futures de la ligne bleuestation_bleu <-st_read('data/transport_projete.gpkg',layer ='metro_projete_l5_station_20221114',quiet =TRUE ) %>%st_transform(st_crs(zat)) %>%fselect(c('a_station')) %>%st_buffer(250) %>%st_intersection(zat['id_zat']) %>%st_drop_geometry() %>%group_by(id_zat,a_station) %>%slice_head(n =1) %>%ungroup() %>%group_by(id_zat) %>%summarise(total =n())bonus <- station_bleu$total[match(base_data$idzat, station_bleu$id_zat)]bonus <-ifelse(is.na(bonus), 0, bonus)base_data$nbmetro <-ifelse( (base_data$AN %in%c(2031, 2036, 2041)), base_data$nbmetro + bonus, base_data$nbmetro)```Nous effectuons la même opération pour les stations de train.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# lecture des donnéestrain_stations <-st_read('data/reseau_transport.gpkg',layer ="stations_train",quiet = T) %>%st_transform(st_crs(zat)) %>%fselect(c('valide_de', 'valide_a', 'NO_STATION'))# Pour chaque station, nous devons déterminer les années pendant lesquelles# elle était validetrain_stations$valide_de <-ifelse( train_stations$valide_de ==0, 1980, train_stations$valide_de)train_stations$station_id <-1:nrow(train_stations)dico_train <- train_stations %>%st_buffer(250) %>%st_intersection(zat[c('id_zat')]) %>%st_drop_geometry() %>%group_by(id_zat, station_id) %>%slice_head(n =1) %>%ungroup()dico_train$start <-as.numeric(substr(as.character(dico_train$valide_de),1,4))dico_train$end <-as.numeric(substr(as.character(dico_train$valide_a),1,4))train_stations$start <-as.numeric(substr(as.character(train_stations$valide_de),1,4))train_stations$start <-as.numeric(substr(as.character(train_stations$valide_de),1,4))base_data$nbtrain <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_trains <- dico_train %>%filter( (start <= year) & (end >= year)) %>%group_by(id_zat) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_trains$nb[match(base_data$idzat, ok_trains$id_zat)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$nbtrain <-ifelse( base_data$AN == year, vals, base_data$nbtrain )}# validation cartographiquemaps <-lapply(years, function(year){ sub_data <-subset(base_data, (base_data$AN == year) & (base_data$nbtrain >0)) sub_train <-subset(train_stations, train_stations$start < year) sub_data$geometry <- zat$geometry[match(sub_data$idzat, zat$id_zat)] sub_data <-st_as_sf(sub_data, crs =st_crs(zat)) map <-tm_shape(sub_data) +tm_borders('black') +tm_shape(sub_train) +tm_dots('red') +tm_title(year)return(map)})tmap_arrange(maps, ncol =3)```Nous allons également ajouter les stations du REM. Il faut ajouter l'existant et le projeté.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falserem_stations_existing <-st_read('data/reseau_transport.gpkg',layer ='reseau_rem_station', quiet =TRUE)%>%st_transform(st_crs(zat)) %>%fselect(c('id')) %>%rename(any_of(c('NO_STATION'='id'))) %>%mutate(valide_de =2023,valide_a =99999999 )%>%st_set_geometry('geometry')rem_stations_projete <-st_read('data/transport_projete.gpkg',layer ='reseau_rem_projete_station_20231011',quiet =TRUE)%>%st_drop_geometry() %>%st_as_sf(coords =c('point_x', 'point_y'), crs =2950) %>%st_transform(st_crs(zat)) %>%fselect(c('id', 'antenne')) %>%rename(any_of(c('NO_STATION'='id'))) %>%mutate(NO_STATION = NO_STATION +200,valide_de =ifelse(grepl('A2', antenne, fixed =TRUE),2027,2025),valide_a =99999999 ) %>%mutate(antenne =NULL ) %>%st_set_geometry('geometry')rem_stations <-rbind(rem_stations_existing, rem_stations_projete)# Pour chaque station, nous devons déterminer les années pendant lesquelles# elle était validerem_stations$valide_de <-ifelse( rem_stations$valide_de ==0, 1980, rem_stations$valide_de)rem_stations$station_id <-1:nrow(rem_stations)dico_rem <- rem_stations %>%st_buffer(250) %>%st_intersection(zat[c('id_zat')]) %>%st_drop_geometry() %>%group_by(id_zat, station_id) %>%slice_head(n =1) %>%ungroup()dico_rem$start <-as.numeric(substr(as.character(dico_rem$valide_de),1,4))dico_rem$end <-as.numeric(substr(as.character(dico_rem$valide_a),1,4))rem_stations$start <-as.numeric(substr(as.character(rem_stations$valide_de),1,4))rem_stations$start <-as.numeric(substr(as.character(rem_stations$valide_de),1,4))base_data$nbrem <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_rem <- dico_rem %>%filter( (start <= year) & (end >= year)) %>%group_by(id_zat) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_rem$nb[match(base_data$idzat, ok_rem$id_zat)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$nbrem <-ifelse( base_data$AN == year, vals, base_data$nbrem )}# validation cartographiquemaps <-lapply(years, function(year){ sub_data <-subset(base_data, (base_data$AN == year) & (base_data$nbrem >0)) sub_rem <-subset(rem_stations, rem_stations$start < year)if(nrow(sub_data) >0){ sub_data$geometry <- zat$geometry[match(sub_data$idzat, zat$id_zat)] sub_data <-st_as_sf(sub_data, crs =st_crs(zat)) map <-tm_shape(sub_data) +tm_borders('black') +tm_shape(sub_rem) +tm_dots('red') +tm_title(year)return(map) }else{return(NULL) }})maps[sapply(maps, is.null)] <-NULLtmap_arrange(maps, ncol =3)```Nous allons calculer pour chaque super zone le nombre total de stations de métro de train et de REM. Ceci permettra notamment d'ajouter des effets régionaux pour régulariser le modèle et représenter un potentiel effet diffusé des stations.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedico_rem_zone <- rem_stations %>%st_buffer(250) %>%st_intersection(super_zat[c('id')]) %>%rename(c('idzone'='id')) %>%st_drop_geometry() %>%group_by(idzone, station_id) %>%slice_head(n =1) %>%ungroup()dico_rem_zone$start <-as.numeric(substr(as.character(dico_rem_zone$valide_de),1,4))dico_rem_zone$end <-as.numeric(substr(as.character(dico_rem_zone$valide_a),1,4))base_data$totalremzone <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_rem <- dico_rem_zone %>%filter( (start <= year) & (end >= year)) %>%group_by(idzone) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_rem$nb[match(base_data$idzone, ok_rem$idzone)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$totalremzone <-ifelse( base_data$AN == year, vals, base_data$totalremzone )}dico_train_zone <- train_stations %>%st_buffer(250) %>%st_intersection(super_zat[c('id')]) %>%rename(c('idzone'='id')) %>%st_drop_geometry() %>%group_by(idzone, station_id) %>%slice_head(n =1) %>%ungroup()dico_train_zone$start <-as.numeric(substr(as.character(dico_train_zone$valide_de),1,4))dico_train_zone$end <-as.numeric(substr(as.character(dico_train_zone$valide_a),1,4))base_data$totaltrainzone <-0years <-seq(1996, 2041, by =5)# pour chaque année de notre analysefor (year in years){# nous gardons uniquement les stations actives et nous comptons # combien il y en a par ZAT ok_trains <- dico_train_zone %>%filter( (start <= year) & (end >= year)) %>%group_by(idzone) %>%summarise(nb =n())# on extrait ensuite ce comptage vals <- ok_trains$nb[match(base_data$idzone, ok_trains$idzone)] vals <-ifelse(is.na(vals), 0, vals)# on attribue ce comptage à l'année en question dans les données panel base_data$totaltrainzone <-ifelse( base_data$AN == year, vals, base_data$totaltrainzone )}dico_metro_zone <- super_zat %>%st_buffer(250) %>%st_intersection(station_metro['nom']) %>%st_drop_geometry() %>%group_by(id,nom) %>%slice_head(n =1) %>%ungroup() %>%group_by(id) %>%summarise(totalmetrozone =n(),metro_names_zone =paste0(nom, collapse =';') )base_data <-left_join(base_data, dico_metro_zone, by =c('idzone'='id'))base_data$totalmetrozone <-ifelse(is.na(base_data$totalmetrozone), 0, base_data$totalmetrozone)names_2007 <-c("Cartier", "De la Concorde", "Montmorency")for(n in names_2007){ base_data$totalmetrozone <-ifelse( (grepl(n,base_data$metro_names_zone,fixed =TRUE)) & (base_data$AN <2007), base_data$totalmetrozone -1, base_data$totalmetrozone )}station_bleu <-st_read('data/transport_projete.gpkg',layer ='metro_projete_l5_station_20221114',quiet =TRUE) %>%st_transform(st_crs(super_zat)) %>%fselect(c('a_station')) %>%st_buffer(250) %>%st_intersection(super_zat['id']) %>%st_drop_geometry() %>%group_by(id,a_station) %>%slice_head(n =1) %>%ungroup() %>%group_by(id) %>%summarise(total =n())bonus <- station_bleu$total[match(base_data$idzone, station_bleu$id)]bonus <-ifelse(is.na(bonus), 0, bonus)base_data$totalmetrozone <-ifelse( (base_data$AN %in%c(2031, 2036, 2041)), base_data$totalmetrozone + bonus, base_data$totalmetrozone)```## Structuration des données d'autoroutesSur le même principe que la section prédécente, nous allons ajouter des variables sur l'accumulation de km linéaires d'autoroute par ZAT. Nous avons pour cela codé manuellement les années d'ouverture des différents tronçons autoroutiers en nous basant sur la source suivante : [Répertoire des autoroutes du Québec](https://www.transports.gouv.qc.ca/fr/projets-infrastructures/info-reseau-routier/repertoire-autoroutes/Pages/repertoire-des-autoroutes-du-Quebec.aspx).```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseautoroutes <-st_read('data/autoroutes.gpkg', quiet = T) %>%st_transform(st_crs(zat))base_data$autoroutesKM <-0for (year inseq(1996, 2041, by =5)){ inter_zat <- autoroutes %>%filter(valide_de < year) %>%st_intersection(zat['id_zat']) %>%mutate(length_km =as.numeric(st_length(geom)) /1000 ) %>%st_drop_geometry() %>%group_by(id_zat) %>%summarise(autoroutesKM =sum(length_km) ) values <- inter_zat$autoroutesKM[match(base_data$idzat, inter_zat$id_zat)] values <-ifelse(is.na(values), 0, values) base_data$autoroutesKM <-ifelse( base_data$AN == year, values, base_data$autoroutesKM )}# nous créons une variable ordinale pour le kilométrage autoroutier. Elle # pourra être utilisée comme un effet monotone non linéaire dans le modèlebase_data$autoroutesMO <-case_when( base_data$autoroutesKM ==0~0, base_data$autoroutesKM >0& base_data$autoroutesKM <=2~1, base_data$autoroutesKM >2& base_data$autoroutesKM <=5~2, base_data$autoroutesKM >5& base_data$autoroutesKM <=10~3, base_data$autoroutesKM >10& base_data$autoroutesKM <=15~4, base_data$autoroutesKM >5~5)```On finalise la structuration des données ci-dessous et on exporte le tableau final pour ne pas avoir à relancer cette étape de structuration à chaque fois.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# on termine la structuration des données et on les exportesbase_data$ANint <- (base_data$AN -1996) /5+1base_data$lnANint <-log(base_data$ANint)base_data$nomreg <-as.factor(base_data$nomreg )base_data$idzat <-as.factor(base_data$idzat )# nous retirons les données de 2021 du fait de leur mauvaise qualitébase_data2 <-subset(base_data, base_data$AN <2021)saveRDS(base_data, "artefacts/base_data_for_projection.rds")```## Détermination de λConsidérant la nature strictement positive et fortement asymétrique de la variable à modéliser, les chercheurs des Universités INRS et Laval avaient proposé d'utiliser la fonction de transformation *log*. Bien que ce choix soit justifié dans ce contexte, nous avions observé dans l'analyse des prédictions du modèle original que certaines ZAT voyaient leur nombre d'emplois exploser, probablement à cause de l'effet de la fonction *log* utilisée pour transformer le nombre d'emplois.Nous proposons donc d'utiliser la transformation BoxCox qui est une généralisation de la fonction *log*. Elle nécessite de déterminer un paramètre $\lambda$ contrôlant le caractère exponentiel de la transformation. Puisque nous souhaitons obtenir un effet moins explosif que celui de la fonction *log*, nous devons trouver la bonne valeur de $\lambda$ avec $ 0 \leq \lambda \leq 1 $.Pour cela, nous allons ajuster le modèle complet avec nos différentes transformations. Pour comparer les modèles, nous allons nous baser sur trois mesures: le rmse, le R2 bayésien et le graphique des posterior-predictive checks.Le code ci-dessous prend beaucoup de temps, car nous devons ajuster 30 modèles très complexes. ```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| eval: true# déterminer la bonne valeur de lambdavals <-seq(0,0.3, 0.01)for (v in vals){ base_data2[[paste0('bcY',v)]] <-boxcox(base_data2$Y1, v)}``````{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| eval: falselibrary(future)library(future.apply)plan(multisession(workers =15))models <-future_lapply(vals, function(v){ formula <-bf(as.formula(paste0('bcY',v, ' ~ s(lnAreaCommercial, k = 3) + s(SUMPLANCHERMANUFAREA, k = 3) + mo(Total10EtagesMO) + s(SUMPLANCHERTRANSPAREA, k = 3) + mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint + (1|idzat) + (1 + lnANint | nomreg) + arma(time = AN_int, p = 1, q = 1, gr = idzat, cov = TRUE) '))) fit1 <-brm( formula,data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =1,cores =1,iter =3000,#control = list(adapt_delta = 0.95),refresh =10 ) fit1 <-add_criterion(fit1, c("loo", "loo_R2")) this_plot <-pp_check(fit1) preds <-predict(fit1, robust =TRUE) rmse <-sqrt(mean((boxcox.inv(fit1$data[,1],v) -boxcox.inv(preds[,1],v))**2))return(list(fit1,preds,this_plot,rmse) )})saveRDS(models, 'artefacts/lambda_selection2.rds')```Une fois que les modèles sont ajustés, nous pouvons déterminer le meilleur d'entre eux en nous basant sur le R2 bayésien ainsi que le RMSE: ```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsevals <-seq(0,0.3, 0.01)models <-readRDS('artefacts/lambda_selection2.rds')comp_df <-data.frame(lambda = vals,r2 =sapply(models, function(x){median(x[[1]]$criteria$loo_R2)}),rmse =sapply(models, function(x){x[[4]]}))# on libère un peu de mémoire viverm(models)gc()ggplot(comp_df) +geom_point(aes(x = lambda, y = r2, size = rmse)) +geom_text(aes(x = lambda, y = r2-0.02, label =round(rmse)), size =3) +ylim(0.7,0.9) +theme_bw()```Il apparait assez clairement que le RMSE et le R2 s'améliorent significativement avec l'augmentation de $\lambda$. Le gain semble maximal lorsque $\lambda > 0.2$. Au delà de cette limite, les gains deviennent anecdotiques. Nous allons donc retenir le modèle avec $\lambda = 0.2$## Ajustement des modèlesDans cette section, nous réalisons l'ajustement des modèles. Nous ajustons trois modèles différents, chacun introduisant une petite variation sur la prise en compte de la tendance temporelle: 1. `lnANint + (1 + lnANint | nomreg)` : cette formulation implique une tendance qui ralentit (log sur le temps) et dont l'effet est différent d'une région à l'autre. Chaque région à son propre effet moyen et sa propre vitesse de croissance. Il existe aussi une corrélation par région entre l'effet moyen et la vitesse de croissance.2. `lnANint + (1 | nomreg) + (0 + lnANint | idzone/nomreg)` : cette formulation implique une tendance qui ralentit (log sur le temps) et dont l'effet est différent d'une super-zat à l'autre et d'une région à l'autre, mais avec un principe hiérarchique entre les super-zat et les régions. Chaque région a aussi son propre effet moyen.3. `ANint + (1 | nomreg) + (0 + ANint | idzone/nomreg)` : idem, mais sans la variation qui ralentit (pas de log sur le temps).Les prédictions finales seront obtenues en croisant les prédictions de ces différents modèles par moyennisation bayésienne.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| eval: false# nous allons intégrer des priors de régularisation, ainsi que des priors # sur les effets positifs des effets fixespriors1 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 1)", class ="sd", coef ='lnANint', group ='nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)# on intègre la croissance au niveau des regionsfit3 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint ++ (1+ lnANint | nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors1,#control = list(adapt_delta = 0.95),refresh =10 )saveRDS(fit3, 'model_lnANINT_RE.rds')# on intègre la croissance au niveau des superzat et des zonespriors3 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 0.25)", class ="sd", coef ='lnANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='lnANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit5 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint + (1| nomreg) + (0+ lnANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,prior = priors3,iter =5000,control =list(adapt_delta =0.9),refresh =10 )saveRDS(fit5, 'model_lnANINT_REZONES.rds')priors4 <-prior("normal(0.25, 0.1)", class ="b", coef ='ANint') +prior("normal(0, 0.25)", class ="sd", coef ='ANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='ANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit6 <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10EtagesMO) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MCGILLHOSP + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + ANint + (1| nomreg) + (0+ ANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE), sigma ~ nomreg),data =subset(base_data2, base_data2$Y1 >20),family =student(),chains =4,cores =4,prior = priors4,iter =5000,control =list(adapt_delta =0.9),refresh =10 )saveRDS(fit6, 'model_ANINT_REZONES.rds')fit3 <-add_criterion(fit3, criterion =c('loo', 'waic', 'loo_R2'))fit6 <-add_criterion(fit6, criterion =c('loo', 'waic', 'loo_R2'))fit5 <-add_criterion(fit5, criterion =c('loo', 'waic', 'loo_R2'))save(fit3, fit5, fit6, file ='artefacts/fitted_gen2.rda')```### Regard sur la qualité d'ajustementAvant d'analyser le contenu des modèles, nous allons vérifier que les conditions d'application sont respectées et qu'ils parviennent à capter une partie significative du phénomène étudié.Pour faciliter l'analyse, nous allons baser ce diagnostic sur un seul des trois modèles. ```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| layout-ncol: 2load('artefacts/fitted_gen2.rda')# extraction d'échantillon des distributions postérieurs prédictives# application de la méthode de combinaison des prédictions par Bayesian # averaging.posterior_preds1 <-predict(fit3, summary =FALSE, cores =10, ndraws =2500)posterior_preds2 <-predict(fit5, summary =FALSE, cores =10, ndraws =2500)posterior_preds3 <-predict(fit6, summary =FALSE, cores =10, ndraws =2500)posterior_preds <-rbind(posterior_preds1, posterior_preds2, posterior_preds3)posterior_preds_trans <-boxcox.inv(posterior_preds, 0.2)pred_df <-data.frame(y = fit5$data$bcY0.2,y_real =boxcox.inv(fit5$data$bcY0.2, 0.2),ypred =colQuantiles(posterior_preds, probs =0.5),ypred_low =colQuantiles(posterior_preds, probs =0.2),ypred_high =colQuantiles(posterior_preds, probs =0.8),ypred_real =colQuantiles(posterior_preds_trans, probs =0.5),ypred_real_low =colQuantiles(posterior_preds_trans, probs =0.2),ypred_real_high =colQuantiles(posterior_preds_trans, probs =0.8))plot1 <-ggplot(pred_df) +geom_point(aes(x = y, y = ypred), size =0.1) +labs(x ='bc(emploi)', y ='bc(emploi prédit)')+geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed')+theme_bw() +coord_fixed()plot2 <-ggplot(pred_df) +geom_point(aes(x = y_real, y = ypred_real), size =0.1) +geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed') +labs(x ='emploi', y ='emploi prédit')+theme_bw() +coord_fixed()plot1plot2```La figure comparant les prédictions et les valeurs réelles montre un très bon ajustement, que ce soit dans l'échelle transformée ou dans l'échelle originale des données.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsepred_df$nomreg <- fit5$data$nomregdf1 <- pred_df %>%mutate(error = (y_real - ypred_real)**2 ) %>%group_by(nomreg) %>%summarise(rmse =sqrt(median(error)) )df2 <-data.frame(nomreg ='TOTAL',rmse =sqrt(median((pred_df$y_real - pred_df$ypred_real)**2)))knitr::kable(rbind(df1, df2), digits =0)```Les Root Median Squared Error calculés à l'échelle des régions varient sensiblement autour du RMSE moyen régional. Nous notons tout de même que nous avons des meilleures prédictions pour la Couronne Sud et des prédictions légèrement moins bonnes pour Montréal Centre et Ouest.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedfs <-lapply(1:12, function(i){ df <- pred_df[sample(1:nrow(pred_df), size =250, replace = F),] df$i <- ireturn(df)})dfs <-do.call(rbind, dfs)dfs$i <-as.factor(dfs$i)ggplot(dfs) +geom_errorbar(aes(x = y_real,ymin = ypred_real_low, ymax = ypred_real_high),linewidth =0.1 ) +geom_point(aes(x = y_real, y = ypred_real), size =1) +geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed') +labs(x ='emploi', y ='emploi prédit')+theme_bw() +coord_fixed() +facet_wrap(vars(i), nrow =4, ncol =3)```Les prédictions ont des intervalles de confiance qui grandissent avec le nombre d'emploi. Cependant, en dehors de quelques valeurs extrêmes, nous avons dans l'ensemble un ajustement très raisonable.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsepp_check(fit3, ndraws =25) +xlim(0,40)pp_check(fit5, ndraws =25) +xlim(0,40)pp_check(fit6, ndraws =25) +xlim(0,40)```La figure ci-dessus indique clairement que le modèle est capable de générer des données reproduisant la distribution des données originales.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsesummary(fit3$criteria$loo_R2)summary(fit5$criteria$loo_R2)summary(fit6$criteria$loo_R2)```Le R2 conditionnel du modèle est très élevé (0.86).```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseR2_fixe <-bayes_R2(fit5,summary =TRUE,robust =TRUE,re_formula =NA )R2_fixe```Il est intéressant de constater qu'une partie non négligeable de la variance provient tout de même des effets fixes (0.65). Ceci signifie que nos variables utilisées pour prédire le nombre d'emplois sont responsables d'une partie non négligeable de sa variation. ```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# effet spécifique des ZATR2_part_zat <-bayes_R2(fit5,summary =TRUE,robust =TRUE,re_formula =~ (1| idzat) )R2_part_zat# effet spécifique des régionsR2_part_reg <-bayes_R2(fit5,summary =TRUE,robust =TRUE,re_formula =~ (1| nomreg) )R2_part_reg# effet spécifique des trends temporelles par zoneR2_part_temp <-bayes_R2(fit5,summary =TRUE,robust =TRUE,re_formula =~ (0+ lnANint | idzone/nomreg) )R2_part_temp```Parmis les effets aléatoires, il apparait clairement que l'effet propre des ZAT joue un rôle important dans la qualité de prédiction du modèle. Les effets régionaux et les ajustements locaux des tendances temporelles jouent un rôle plus mineur.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| layout-ncol: 2library(gamlss, quietly =TRUE)resids <-residuals(fit5)resid_mod <-gamlss(resids[,1] ~1, family = TF)qts_obs <-quantile(resids[,1], probs =seq(0.01,0.99,0.01))qts_theo <-qTF(seq(0.01,0.99,0.01),mu = resid_mod$mu.coefficients[[1]],sigma =exp(resid_mod$sigma.coefficients[[1]]),nu =exp(resid_mod$nu.coefficients[[1]]) )df_plot <-data.frame(qts_obs = qts_obs,qts_theo = qts_theo)plot1 <-ggplot(df_plot) +geom_abline(slope =1, intercept =0, color ='red') +geom_point(aes(x = qts_obs, y = qts_theo)) +coord_fixed() +theme_bw() +labs(x ='quantiles observés', y ='quantiles attendus', title ='QQplot des résidus')resid_df <-data.frame(resids)plot2 <-ggplot(resid_df) +geom_histogram(aes(x = Estimate, y = ..density..),bins =50, color ='white') +stat_function(fun = dTF, args =list(mu = resid_mod$mu.coefficients[[1]],sigma =exp(resid_mod$sigma.coefficients[[1]]),nu =exp(resid_mod$nu.coefficients[[1]]) ),color ='red', size =1 ) +theme_bw() +labs(title ='Distribution des résidus')plot1plot2```Les résidus naturels du modèles suivent bien une distribution de Student tel que spécifié.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedf <- fit5$datadf$resid <- resids[,1]collapse::psacf(df$resid, df$idzat, df$AN_int,plot =FALSE )collapse::pspacf(df$resid, df$idzat, df$AN_int,plot =FALSE )```Notre modèle parvient bien à retirer l'autocorrélation temporelle d'ordre 1, mais il reste de l'autocorrélation négative d'ordre 2 et 3. Ceci est problématique car certains paramètres du modèle pourraient ne pas être bien estimés. D'un point de vue pratique, cela signifie que pour une période temporelle *t*, les résidus du modèle tendent à être négativement associés avec les résidus à *t-1*. Cependant, considérant que nous n'avons que 6 périodes de temps dans les données, cet enjeu est beaucoup moins problématique.Nous pouvons aussi regarder l'autocorrélation spatiale de nos résidus à chaque année.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| layout-ncol: 2dfs <-split(df, df$AN)for(x in dfs){ year <- (unique(x$AN) -1 )*5+1996 x$idzat <-as.numeric(as.character(x$idzat)) x_space <-left_join(zat, x, by =c('id_zat'='idzat')) x_space <-subset(x_space, !is.na(x_space$resid)) queen <-poly2nb(x_space) WQueen <-nb2listw(queen, style ='W', zero.policy =TRUE) M <-moran.mc(x_space$resid, listw = WQueen, nsim =999) text <-paste0(year, ', moran I = ',round(M$statistic[[1]],2)) map <-tm_shape(x_space) +tm_polygons('resid', fill.scale =tm_scale_intervals(style ='jenks', n =7, values ='RdBu',midpoint =0 ), lwd =0.1) +tm_title(text)print(map)}```Les niveaux d'autocorrélation spatiale des résidus sont négligeables, ce qui est important pour ce type de variable fortement spatialement autocorrélée.### Validation croiséePuisque nos modèles semblent être correctement ajustés, nous allons également tenté de déterminer dans quelle mesure ils sont pertinents pour effectuer des prédictions.Pour cela, nous allons simplement réajuster des modèles identiques, mais en retirant les données de 2016. L'idée est alors de voir si nous aurions eu des prédictions réalistes comparativement à ce qui a été effectivement observé en 2016.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| eval: falsebase_data3 <-subset(base_data2, base_data2$AN <2016)priors1 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 1)", class ="sd", coef ='lnANint', group ='nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)priors3 <-prior("normal(0.25, 0.1)", class ="b", coef ='lnANint') +prior("normal(0, 0.25)", class ="sd", coef ='lnANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='lnANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)priors4 <-prior("normal(0.25, 0.1)", class ="b", coef ='ANint') +prior("normal(0, 0.25)", class ="sd", coef ='ANint', group ='idzone') +prior("normal(0, 0.5)", class ="sd", coef ='ANint', group ='idzone:nomreg') +prior("normal(0.75, 1)", class ="b", coef ='MAJORHOSPITAL') +prior("normal(0.75, 1)", class ="b", coef ='POWERCENTRE') +prior("normal(0.75, 1)", class ="b", coef ='WMANUF15K') +prior("normal(0.75, 1)", class ="b", coef ='nbmetro') +prior("normal(0.75, 1)", class ="b", coef ='nbtrain') +prior("normal(0.15, 0.5)", class ="b", coef ='totalmetrozone') +prior("normal(0.15, 0.5)", class ="b", coef ='totaltrainzone') +prior("normal(0, 0.15)", dpar ='sigma', class ='b') +prior("gamma(3,0.75)", class ='nu') +prior('normal(0.7,0.25)', class ='ar', lb =-1, ub =1) +prior('normal(0,0.35)', class ='ma', lb =-1, ub =1)fit3_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint ++ (1+ lnANint | nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors1,#control = list(adapt_delta = 0.95),refresh =10 )fit4_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + lnANint + (1| nomreg) + (0+ lnANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,iter =3000,prior = priors3,#control = list(adapt_delta = 0.95),refresh =10 )fit5_cv <-brm(bf(bcY0.2~s(lnAreaCommercial, k =3) +s(SUMPLANCHERMANUFAREA, k =3) +mo(Total10Etages) +s(SUMPLANCHERTRANSPAREA, k =3) +mo(autoroutesMO) + PCTLOTSAGRIEXTRAREA + MAJORHOSPITAL + POWERCENTRE + WMANUF15K + nbmetro + nbtrain + totalmetrozone + totaltrainzone + ANint + (1| nomreg) + (0+ ANint | idzone/nomreg) + (1|idzat) +arma(time = AN_int, gr = idzat, cov =TRUE),sigma =~ nomreg ),data =subset(base_data3, base_data3$Y1 >20),family =student(),chains =4,cores =4,prior = priors4,iter =3000,#control = list(adapt_delta = 0.95),refresh =10 )save(fit3_cv, fit4_cv, fit5_cv, file ='artefacts/fitted_gen1_cv.rda')```Maintenant que nous avons nos modèles, regardons leur prédiction en 2016```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseload('artefacts/fitted_gen1_cv.rda')data_2016 <- base_data2 %>%fsubset(Y1 >20) %>%fsubset(idzat %in% fit5_cv$data$idzat)data_2016$idzat <-as.factor(as.integer(as.character(data_2016$idzat)))# on effectue les prédictions en utilisant les trois modèles. Nous allons # ainsi récupérer 1200 posterior samples que nous pourrons combiner pour # obtenir la prédiction finale par moyennisation bayésiennepred_2016_fit3 <-posterior_predict(fit3_cv, newdata = data_2016, ndraws =1200)pred_2016_fit4 <-posterior_predict(fit4_cv, newdata = data_2016, ndraws =1200)pred_2016_fit5 <-posterior_predict(fit5_cv, newdata = data_2016, ndraws =1200)pred_2016 <-rbind(pred_2016_fit3, pred_2016_fit4, pred_2016_fit5)pred_2016_trans <-boxcox.inv(pred_2016_fit5, 0.2)data_2016$pred <-colQuantiles(pred_2016_trans, probs =0.5)data_2016$pred_low <-colQuantiles(pred_2016_trans, probs =0.2)data_2016$pred_high <-colQuantiles(pred_2016_trans, probs =0.8)plot1 <-ggplot(data_2016) +geom_errorbar(aes(x = Y1,ymin = pred_low, ymax = pred_high),linewidth =0.1,alpha =0.3 ) +geom_point(aes(x = Y1, y = pred), size =0.4) +geom_abline(intercept =0, slope =1, color ='red', linetype ='dashed') +labs(x ='emploi', y ='emploi prédit')+theme_bw() +coord_fixed()ggsave(plot = plot1, device ='png', filename ='cv_2016.png',path ='images/', width =20, height =20, units ='cm', dpi =200)knitr::include_graphics('images/cv_2016.png')```La figure ci-dessus illustre que le modèle parvient à faire une prédiction crédible pour 2016, surtout lorsque l'on considère l'incertitude avec notre intervalle de confiance.Ci-dessous, nous calculons le Root Median Squared Error du modèle complet et du modèle de cross-validation. Si les deux RMdSE sont proches, cela signifie que le modèle complet n'a pas fait d'overfitting pour 2016.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# calcul du RMSE du modèle completpreds_comp1 <-posterior_predict(fit3, summary =FALSE, cores =10)preds_comp2 <-posterior_predict(fit5, summary =FALSE, cores =10)preds_comp3 <-posterior_predict(fit6, summary =FALSE, cores =10)preds_comp <-rbind(preds_comp1, preds_comp2, preds_comp3)preds_comp <-boxcox.inv(preds_comp, 0.2)Y <-boxcox.inv(fit5$data$bcY0.2, 0.2)rmse_total <-sapply(1:nrow(preds_comp), function(i){sqrt(median((preds_comp[i,] - Y)**2))})median(rmse_total)# calcul du RMSE du modèle de cross-validationpreds_cv1 <-posterior_predict(fit3_cv, summary =FALSE, cores =10)preds_cv2 <-posterior_predict(fit4_cv, summary =FALSE, cores =10)preds_cv3 <-posterior_predict(fit5_cv, summary =FALSE, cores =10)preds_cv <-rbind(preds_cv1,preds_cv2,preds_cv3)preds_cv <-boxcox.inv(preds_cv, 0.2)Y <-boxcox.inv(fit5_cv$data$bcY0.2, 0.2)rmse_cv <-sapply(1:nrow(preds_cv), function(i){sqrt(median((preds_cv[i,] - Y)**2))})median(rmse_cv)```Il est intéressant de constater que les RMdSE obtenus avec le modèle complet et avec le modèle sans les données de 2016 sont très similaires. Ceci signifique le modèle complet ne fait pas d'over-fitting, et que l'on peut raisonnablement espérer que sa qualité d'ajustement actuelle se répercutera sur les projections.### Analyse des tendances captées par les modèles#### Tendances temporellesNous avons ajusté ci-dessus trois versions différentes du modèle incluant des variations dans la façon d'inclure la tendance temporelle. Dans le premier modèle, nous représentons le fait que chaque ZAT de la région suit une tendance centrale non linéaire (`lnANint`), et que cette tendance varie par région d'analyse (`1 + lnANint | nomreg`). Les figures ci-dessous illustrent les tendances modélisées.Plus spécifiquement, elles illustrent le parcours d'une ZAT identique dans différentes régions. Elles doivent donc être interprétées de façon relative plutôt qu'absolue.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false#| layout-ncol: 3load('artefacts/fitted_gen2.rda')# tendance apprises par la première version du modèle# pour cela, nous allons fixer toutes les valeurs d'une observation# et la duppliquer pour chaque zone, région et annéedf1 <- fit3$data[1,]pred_data <-expand.grid(nomreg =unique(fit3$data$nomreg),AN =seq(1996,2041,5))pred_data$AN_int <- (pred_data$AN -1996) /5+1for(var innames(df1)){if( !(var %in%c('nomreg', 'AN_int'))){ pred_data[[var]] <- df1[[var]] }}pred_data$lnANint <-log(pred_data$AN_int)pred_data$AN <- (pred_data$AN_int-1) *5+1996preds <-posterior_epred(fit3, pred_data, incl_autocor =FALSE)pred_data$pred <-colQuantiles(preds, probs =0.5)pred_data$pred_low <-colQuantiles(preds, probs =0.3)pred_data$pred_high <-colQuantiles(preds, probs =0.7)# représentation dans l'échelle transforméeplot1 <-ggplot(pred_data) +geom_line(aes(x = AN,y = pred_low,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN,y = pred_high,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN, y = pred, color = nomreg, group = nomreg), linewidth =1) +labs(x ='année',y ='bc(emplois)' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5)))# représentation dans l'échelle originalepreds_trans <-boxcox.inv(preds, 0.2)pred_data$pred_trans <-colQuantiles(preds_trans, probs =0.5)pred_data$pred_trans_low <-colQuantiles(preds_trans, probs =0.3)pred_data$pred_trans_high <-colQuantiles(preds_trans, probs =0.7)plot2 <-ggplot(pred_data) +geom_line(aes(x = AN,y = pred_trans_low,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN,y = pred_trans_high,group = nomreg, color = nomreg), linetype ='dashed') +geom_line(aes(x = AN, y = pred_trans, color = nomreg, group = nomreg), linewidth =1) +labs(x ='année',y ='emplois' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5)))# représentation dans une échelle relativepred_data <- pred_data %>%ungroup() %>%fsubset(AN ==1996) %>%fselect(c('nomreg', 'pred_trans')) %>%rename(c('pred_trans_1996'='pred_trans')) %>%right_join(pred_data) %>%mutate(pred_prt = pred_trans / pred_trans_1996 *100 )plot3 <-ggplot(pred_data) +geom_line(aes(x = AN, y = pred_prt, color = nomreg, group = nomreg))+labs(x ='année',y ='emplois (base 100 1996)' ) +theme_bw() +scale_x_continuous(breaks =c(seq(1996,2041,5)))plot1plot2plot3```Les trois figures ci-dessus permettent de constater que le modèle prévoit une croissance pour chacune des régions. La croissance la plus faible serait prévue pour Montréal Centre, suivi ensuite par Montréal Est et Montréal Ouest. Les couronnes connaitraient quant à elle les croissances les plus élevées.Le second modèle va un peu plus loin en faisant varier le taux de croissance par région et par Super-Zat.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedf1 <- fit5$data[1,]pred_data <- fit5$data %>%group_by(idzone, nomreg) %>%slice_head(n =1) %>%fselect(c('idzone', 'nomreg'))dfs <-lapply(seq(1996, 2041, by =5), function(i){ df2 <- pred_data df2$AN <- ireturn(df2)})pred_data <-do.call(rbind,dfs)pred_data$AN_int <- (pred_data$AN -1996) /5+1for(var innames(df1)){if( !(var %in%c('idzone', 'nomreg', 'AN_int'))){ pred_data[[var]] <- df1[[var]] }}pred_data$lnANint <-log(pred_data$AN_int)preds <-posterior_epred(fit5, pred_data, incl_autocor =FALSE)preds <-boxcox.inv(preds, 0.2)pred_data$pred <-colQuantiles(preds, probs =0.5)pred_data$id <-paste0(pred_data$idzone, pred_data$nomreg)ggplot(pred_data) +geom_line(aes(x = AN, y = pred, color = nomreg, group = id)) +labs(x ='année',y ='emplois' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5)))+ylim(5500,9000)```Cette seconde version du modèle permet de capturer l'hétérogénéité au niveau des Super-Zat. Nous pouvons notamment constater que certaines d'entre elles sont marquées par des croissances plus faibles, voire même des décroissances.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedf1 <- fit6$data[1,]pred_data <- fit6$data %>%group_by(idzone, nomreg) %>%slice_head(n =1) %>%fselect(c('idzone', 'nomreg'))dfs <-lapply(seq(1996, 2041, by =5), function(i){ df2 <- pred_data df2$AN <- ireturn(df2)})pred_data <-do.call(rbind,dfs)pred_data$AN_int <- (pred_data$AN -1996) /5+1for(var innames(df1)){if( !(var %in%c('idzone', 'nomreg', 'AN_int'))){ pred_data[[var]] <- df1[[var]] }}pred_data$ANint <- pred_data$AN_intpreds <-posterior_epred(fit6, pred_data, incl_autocor =FALSE)preds <-boxcox.inv(preds, 0.2)pred_data$pred <-colQuantiles(preds, probs =0.5)pred_data$id <-paste0(pred_data$idzone, pred_data$nomreg)ggplot(pred_data) +geom_line(aes(x = AN, y = pred, color = nomreg, group = id)) +labs(x ='année',y ='emplois' ) +theme_bw()+scale_x_continuous(breaks =c(seq(1996,2041,5))) +ylim(5500,10000)```Le troisième modèle inclut la tendance temporelle avec un coefficient linéaire. Il produit donc des croissances plus importantes, surtout dans le long terme. Il serait probablement plus pertinent de combiner les prédictions des trois modèles pour limiter l'effet problématique des décroissances.Cette approche se nomme *bayesian model averaging*. Elle permet notamment de combiner l'incertitude des modèles ainsi que les éléments qu'ils ont appris respectivement.## Regard sur les projectionsNous pouvons à présent porter un regard sur les prédictions réalisées par les modèles.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# réalisation des projections base_data$lnANint <-log(base_data$AN_int)base_data$ANint <-log(base_data$AN_int)base_data_pred <-subset(base_data, base_data$idzat %in% fit5$data$idzat)base_data_pred$idzat <- base_data_pred$idzat %>%as.character() %>%as.integer() %>%as.factor()preds <-pp_average(fit3, fit5, fit6, weights ='loo',summary =FALSE,ndraws =2500,cores =10,newdata = base_data_pred )# les draws issus des distributions postérieures doivent être retransformés # avec l'inverse de notre fonction BoxCox pour ensuite pouvoir calculer nos # intervalles de confiancepreds_trans <-boxcox.inv(preds, 0.2)# il se peut que nous obtenions des NA, ce qui est causé par des valeurs très# petites que nous remplaçons donc par la valeur 0preds_trans[is.na(preds_trans)] <-0# nous pouvons ensuite calculer pour chaque simulation (draw) les totaux# d'emploispreds_trans <-t(preds_trans)preds_trans <-data.frame(preds_trans)names(preds_trans) <-paste0('pred',1:ncol(preds_trans))pred_names <-names(preds_trans)preds2 <-cbind(preds_trans, base_data_pred[c('AN', 'idzat', 'idzone', 'nomreg')])# calcul des totaux pour l'ensemble du territoirepred_artm <- preds2 %>%fgroup_by(AN) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%fselect(pred_names) %>%as.matrix() %>%rowQuantiles(probs =c(0.1,0.5,0.9)) %>%as.data.frame() %>%cbind(unique(preds2$AN)) %>%frename(c('est_low', 'est_med', 'est_high', 'AN'))pred_artm$est_mean <- preds2 %>%fgroup_by(AN) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%fselect(pred_names) %>%as.matrix() %>%rowMeans2()# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by ='AN')pred_artm$total_reel <-ifelse(pred_artm$AN >2021, NA, pred_artm$total_reel)ggplot(pred_artm) +geom_line(aes(x = AN, y = total_reel)) +geom_point(aes(x = AN, y = total_reel)) +geom_ribbon(aes(x = AN, ymin = est_low, ymax = est_high), alpha =0.5) +geom_line(aes(x = AN, y = est_med), color ='red') +geom_point(aes(x = AN, y = est_med)) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2021, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)") +theme_bw() +scale_y_continuous(labels =function(x){x/1000})```Il est intéressant de constater qu'une fois que les prédictions sont agrégées, le modèle tend à sur-prédire le nombre d'emploi dans la région. Il faudra tenir compte de ce phénomème dans nos choix finaux de scénarios de projection.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# calcul des totaux par régionpred_artm <- preds2 %>%fgroup_by(AN, nomreg) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%pivot(how ='longer', values = pred_names) %>%fgroup_by(AN, nomreg) %>%fsummarise(est_low =fquantile(value, probs =0.1),est_med =fquantile(value, probs =0.5),est_high =fquantile(value, probs =0.9) )# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN, nomreg) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by =c('AN', 'nomreg'))pred_artm$total_reel <-ifelse(pred_artm$AN >2016, NA, pred_artm$total_reel)pred_artm$nomreg <-as.factor(pred_artm$nomreg)ggplot(pred_artm) +geom_line(aes(x = AN, y = total_reel)) +geom_point(aes(x = AN, y = total_reel)) +geom_ribbon(aes(x = AN, ymin = est_low, ymax = est_high), alpha =0.5) +geom_line(aes(x = AN, y = est_med), color ='red') +geom_point(aes(x = AN, y = est_med)) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2016, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)") +theme_bw() +scale_y_continuous(labels =function(x){x/1000}) +facet_wrap(vars(nomreg), scales ='free')```Une fois nos projections éclatées par régions, nous pouvons constater que l'écart dans les projections au niveau régional s'explique vraisemblablement par notre sur-prédiction dans les régions Montréal Est et Ouest. En effet, ces deux régions ont connu une baisse significative de leur nombre d'emplois entre 2006 et 2011.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: false# calcul des totaux par superzatpred_artm <- preds2 %>%fgroup_by(AN, nomreg, idzone) %>%fsummarise(across(pred_names,list(tot = fsum))) %>%pivot(how ='longer', values = pred_names) %>%fgroup_by(AN, nomreg, idzone) %>%fsummarise(est_low =fquantile(value, probs =0.1),est_med =fquantile(value, probs =0.5),est_high =fquantile(value, probs =0.9) )# ajout des vrais totaux pour comparaisonreal_tots <- base_data_pred %>%group_by(AN, nomreg, idzone) %>%summarise(total_reel =sum(Y1))pred_artm <-left_join(pred_artm, real_tots, by =c('AN', 'nomreg', 'idzone'))pred_artm$total_reel <-ifelse(pred_artm$AN >2016, NA, pred_artm$total_reel)pred_artm$nomreg <-as.factor(pred_artm$nomreg)pred_artm$idzone <-as.character(pred_artm$idzone)for(reg inunique(pred_artm$nomreg)){ sub_data <-subset(pred_artm, pred_artm$nomreg == reg) this_plot <-ggplot(sub_data) +geom_ribbon(aes(x = AN,ymin = est_low,ymax = est_high, group = idzone), alpha =0.25) +geom_line(aes(x = AN, y = total_reel,group = idzone, color = idzone),linewidth =1) +geom_point(aes(x = AN, y = total_reel, color = idzone)) +geom_line(aes(x = AN, y = est_med, group = idzone, color = idzone),linetype ='dashed',linewidth =1) +scale_x_continuous(breaks =seq(1996, 2041, by =5) ) +geom_vline(xintercept =2016, linetype ='dashed') +labs(x ='année',y ="volume d'emploi (milliers)",title = reg ) +theme_bw() +scale_y_continuous(labels =function(x){x/1000})print(this_plot)}```Les figures ci-dessus illustrent comment le modèle lisse les projections pour conserver les tendances observées à long terme. Pour chacune des régions, nous observons que les tendances par super-zat projetées semblent poursuivre les tendances observées dans les données historiques. Nous remarquons tout de même quelques super-zat dans lesquelles les prédictions tendent à être assez éloignées des valeurs effectivement observées dans les données historiques.## ProjectionsL'objectif principal du modèle développé ici est d'estimer le nombre d'emplois par ZAT pour les horizons 2026, 2031, 2036 et 2041. L'idée est notamment de se doter d'un scénario de référence, ainsi que d'un scénario optimiste et d'un scénario pessimiste.### Déterminer les paramètres de projectionLe modèle de régression bayésien que nous avons estimé permet de déterminer la distribution des prédictions que nous effectuons. Cette prédiction comporte donc l'ensemble de l'incertitude sur nos projections. L'idéal serait donc de conserver cette incertitude complète et de la transférer aux modèles de déplacements. Cependant, puisque c'est peu réaliste compte tenu des temps de calcul de ces modèles, nous allons essayer d'extraire des *point wise estimates* de nos projections.Pour sélectionner les meilleurs *point wise estimates*, nous allons regarder pour le modèle sans les données de 2016, lesquels auraient donné les meilleurs résultats.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falselibrary(bayestestR)library(ggpubr)load('artefacts/fitted_gen1_cv.rda')# extraction des prédictions de notre modèledata_2016 <- base_data2 %>%fsubset(Y1 >20) %>%fsubset(idzat %in% fit5_cv$data$idzat)data_2016$idzat <-as.factor(as.integer(as.character(data_2016$idzat)))preds_comp1 <-posterior_predict(fit3_cv, newdata = data_2016,summary =FALSE, cores =10)preds_comp2 <-posterior_predict(fit4_cv, newdata = data_2016,summary =FALSE, cores =10)preds_comp3 <-posterior_predict(fit5_cv, newdata = data_2016,summary =FALSE, cores =10)pred_2016 <-rbind(preds_comp1, preds_comp2, preds_comp3)pred_2016_trans <-boxcox.inv(pred_2016, 0.2)# sélection d'un ensemble de prédictions pour visualiser leur distributionobs <-sample(1:ncol(pred_2016_trans), size =25, replace =FALSE)plots <-lapply(obs, function(i){ x <- pred_2016_trans[,i] target <- data_2016$Y1[[i]] x <- x[x <quantile(x, probs =0.99)] x <- x[x >0] plot <-ggplot() +geom_histogram(aes(x = x), color ='white', bins =50) +geom_vline(xintercept = target, color ='red',linetype ='dashed', size =1) +labs(x ='', y ='')return(plot)})ggarrange(plotlist = plots, nrow =5, ncol =5)```La figure ci-dessus indique clairement que les distributions des prédictions d'emplois tendent à être asymétriques. Il est donc important de choisir quel percentile de la distribution nous allons utiliser pour nos prédiction. Nous allons notamment regarder les valeurs aux alentours de la médiane.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsemae <-function(x,y){mean(abs(x-y))}mean_error <-function(x,y){mean(x-y)}rmse <-function(x,y){sqrt(mean((x-y)**2))}smape <-function(x,y){mean( abs(abs(x-y) / ((x+y)/2)) )}pcts <-seq(0.25,0.75, by =0.025)preds <-colQuantiles(pred_2016_trans, probs = pcts)df_scores2 <-sapply(1:ncol(preds), function(i){ x <- preds[,i] y <- data_2016$Y1 reg <-lm(x ~0+y) vals <-c(mean_error(x,y),mae(x,y),rmse(x,y),smape(x,y), reg$coefficients[[1]] )return(vals)})df_scores2 <-data.frame(df_scores2)df_scores2$metric <-c('erreur moyenne','erreur absolue moyenne','rmse','smape', 'coeff' )df_scores2 <-pivot(df_scores2, how ='longer', ids ='metric')df_scores2$pct <-gsub('X','',df_scores2$variable,fixed =TRUE)df_scores2$pct <- pcts[as.integer(df_scores2$pct)]df_scores2$metric <-as.factor(df_scores2$metric)ggplot(df_scores2) +geom_line(aes(x = pct,y = value )) +facet_wrap(vars(metric), scales ='free')```Il semblerait qu'une valeur comprise entre les 50e et 60e percentile permettrait d'obtenir les prédictions donnant un bon compromis entre nos différentes métriques d'ajustement.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsedf_scores2_wide <- df_scores2 %>%pivot(how ='wider', values ='value', names ='metric')# le meilleur coefficient est celui le plus proche de 1test <-abs(df_scores2_wide$coeff-1) ==min(abs(df_scores2_wide$coeff-1))best_coeff <- df_scores2_wide$pct[test] # la meilleure erreur absolue moyenne et proche de 0test <- df_scores2_wide$`erreur absolue moyenne`==min(df_scores2_wide$`erreur absolue moyenne`)best_mae <- df_scores2_wide$pct[test] # la meilleur erreur moyenne est la plus proche de 0test <-abs(df_scores2_wide$`erreur moyenne`) ==min(abs(df_scores2_wide$`erreur moyenne`))best_me <- df_scores2_wide$pct[test]# le meilleur RMSE est le plus petittest <- df_scores2_wide$rmse ==min(df_scores2_wide$rmse)best_rmse <- df_scores2_wide$pct[test]# le meilleur smape est le plus petittest <- df_scores2_wide$smape ==min(df_scores2_wide$smape)best_smape <- df_scores2_wide$pct[test]df_best_pcts <-data.frame(score =c('coeff', 'mae', 'me', 'rmse', 'smape'),best_ctl =c(best_coeff, best_mae, best_me, best_rmse, best_smape))knitr::kable(df_best_pcts)```Nous proposons de retenir la valeur de 0.58 pour le quantile de prédiction. Elle offre un bon compromis entre nos différentes mesures d'erreur. Nous devons encore choisir nos bornes supérieures et inférieures qui délimiteront nos scénarios pessimiste et optimiste.Pour notre scénario pessimiste, nous proposons de retenir le percentile garantissant que la croissance générale de la région ne passe pas en dessous des totaux de 2016.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseload('artefacts/fitted_gen2.rda')# réalisation des projections base_data$lnANint <-log(base_data$AN_int)base_data$ANint <- base_data$AN_intbase_data_pred <-subset(base_data, base_data$idzat %in% fit5$data$idzat)base_data_pred$idzat <- base_data_pred$idzat %>%as.character() %>%as.integer() %>%as.factor()#_______________# obtention des échantillons pour nos distributions à posterioripreds1 <-posterior_predict(fit3, newdata = base_data_pred,summary =FALSE, cores =10, ndraws =1200)preds2 <-posterior_predict(fit5, newdata = base_data_pred,summary =FALSE, cores =10, ndraws =1200)preds3 <-posterior_predict(fit6, newdata = base_data_pred,summary =FALSE, cores =10, ndraws =1200)preds_total <-rbind(preds1, preds2, preds3)preds_total <-boxcox.inv(preds_total, 0.2)qtls <-colQuantiles(preds_total, probs =0.58-seq(0,0.15, by =0.01))total_2011 <-sum(base_data_pred$Y1[base_data_pred$AN ==2011])print(total_2011)total_2016 <-sum(base_data_pred$Y1[base_data_pred$AN ==2016])print(total_2016)print(colSums(qtls[base_data_pred$AN ==2021,]))```il semblerait que le 54e percentile permette d'éviter une projection de décroissance. Nous pourrions donc l'utiliser comme notre valeur basse. Comme valeur haute, nous proposons de prendre un écart de 10 points de percentile avec la valeur basse.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseqtls2 <-colQuantiles(preds_total, probs =c(0.54,0.58,0.64))base_data_pred$pred <- qtls2[,2]base_data_pred$low <- qtls2[,1]base_data_pred$high <- qtls2[,3]this_preds <- base_data_pred %>%group_by(AN) %>%summarise(total_y =sum(Y1),total_med =sum(pred),total_low =sum(low),total_high =sum(high) )real <-subset(this_preds, this_preds$AN <2021)real_2016 <-max(real$total_y)proj <- this_preds %>%filter(AN >2016) %>%add_row(AN =2016,total_med = real_2016,total_low = real_2016,total_high = real_2016 )this_plot <-ggplot() +geom_line(data = real, mapping =aes(x = AN, y = total_y)) +geom_point(data = real, mapping =aes(x = AN, y = total_y)) +geom_ribbon(data = proj, mapping =aes(x = AN, ymin = total_low,ymax = total_high ),alpha =0.4) +geom_line(data = proj, mapping =aes(x = AN, y = total_med)) +geom_point(data = proj, mapping =aes(x = AN, y = total_med)) +geom_vline(xintercept =2016, linetype ='dashed') +theme_bw() +scale_x_continuous(breaks =seq(1996, 2041, by =5)) +labs(x ='année', y ="volume d'emplois")this_plot# par regionthis_preds <- base_data_pred %>%group_by(AN, nomreg) %>%summarise(total_y =sum(Y1),total_med =sum(pred),total_low =sum(low),total_high =sum(high) )real <-subset(this_preds, this_preds$AN <2021)real_2016 <-subset(real, real$AN ==2016)real_2016$total_med <- real_2016$total_yreal_2016$total_low <- real_2016$total_yreal_2016$total_high <- real_2016$total_yproj <- this_preds %>%filter(AN >2016) %>%rbind(real_2016)this_plot <-ggplot() +geom_line(data = real, mapping =aes(x = AN, y = total_y, group = nomreg, color = nomreg)) +geom_point(data = real, mapping =aes(x = AN, y = total_y, color = nomreg)) +geom_ribbon(data = proj, mapping =aes(x = AN, ymin = total_low,ymax = total_high, group = nomreg),alpha =0.2) +geom_line(data = proj, mapping =aes(x = AN, y = total_med ,group = nomreg, color = nomreg)) +geom_point(data = proj, mapping =aes(x = AN, y = total_med, color = nomreg))this_plot```### Appliquer les projectionsNous pouvons maintenant appliquer les projections en appliquant nos différentes observations.Une question reste à aborder : quel sera l'impact de l'ouverture des stations de REM sur l'emploi. Pour le moment, elles ont été codées comme des stations de train. Or, nous avons pu observer dans les effets fixes que le train avait un impact proche de zéro. Nous pourrions formuler l'hypothèse qu'une station de REM aura un impact plus proche de l'ouverture d'une station de métro que de train. Dans ce cas, nous pourrions compter une station de REM comme une fraction de métro. Nous allons coder l'impact d'une station de REM comme la moitiée de celui d'une station de métro.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falseload('artefacts/fitted_gen2.rda')# on s'assure la reproductibilité des samplesset.seed(666)# réalisation des projections base_data$lnANint <-log(base_data$AN_int)base_data$ANint <- base_data$AN_intbase_data_pred <-subset(base_data, base_data$idzat %in% fit5$data$idzat)base_data_pred$idzat <- base_data_pred$idzat %>%as.character() %>%as.integer() %>%as.factor()#_______________# les stations de REM sont traitées comme la moitiée des stations de métrobase_data_pred2 <- base_data_predbase_data_pred2$nbmetro <- base_data_pred$nbmetro + base_data_pred$nbrem *0.5base_data_pred2$totalmetrozone <- base_data_pred$totalmetrozone + base_data_pred$totalremzone *0.5#_______________# Combinaison des projections des trois modèlespreds1 <-posterior_predict(fit3, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds2 <-posterior_predict(fit5, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds3 <-posterior_predict(fit6, newdata = base_data_pred2,summary =FALSE, cores =10, ndraws =2500)preds_total <-rbind(preds1, preds2, preds3)# les draws issus des distributions postérieures doivent être retransformés # avec l'inverse de notre fonction BoxCox pour ensuite pouvoir calculer nos # intervalles de confiancepreds_total_trans <-boxcox.inv(preds_total, 0.2)# il se peut que nous obtenions des NA, ce qui est causé par des valeurs très# petites que nous remplaçons donc par la valeur 0preds_total_trans[is.na(preds_total_trans)] <-0qtls <-colQuantiles(preds_total_trans, probs =c(0.54,0.58,0.64))gc()qtls <-data.frame(qtls)names(qtls) <-c('s1_low', 's1_med', 's1_high')projections <-cbind( base_data_pred[c('AN', 'idzat', 'Y1', 'idzone', 'nomreg')], qtls)# on marque les données projetée avec NA pour bien distinguer les données # originales des projections projections$Y1 <-ifelse( projections$AN >2016, NA, projections$Y1)saveRDS(projections, 'artefacts/projections.rds')# nous exportons les projections par ZAT en CSVwrite.csv(projections, 'projection_emplois_ZAT_ARTM.csv', row.names = F)# nous exportons aussi des échantillons des valeurs postérieures pour permettre # des agrégations probabilistes# nous allons retenir uniquement 1500 échantillons postérieursids <-sample(1:nrow(preds_total_trans), size =1500, , replace = F)posteriors_export <- preds_total_trans[ids,] %>%t() %>%data.frame() names(posteriors_export) <-paste0("sample_", 1:ncol(posteriors_export))posteriors_export <-cbind(base_data_pred2[c('AN', 'idzat', 'Y1', 'idzone', 'nomreg')], posteriors_export)write.csv(posteriors_export, 'projection_emplois_ZAT_ARTM_samples.csv', row.names = F)```Terminons en montrant comment utiliser les échantillons postérieurs pour calculer des intervalles de confiances pour des projections agrégées.```{r}#| code-fold: true#| out-width: '100%'#| fig.align: 'center'#| warning: false#| message: falsepost_samples <-read.csv('projection_emplois_ZAT_ARTM_samples.csv')#_______________________________# example avec agrégation régionale# on extrait les données historiques connuespart1 <- post_samples %>%fsubset(AN <2021) %>%fgroup_by(AN) %>%fsummarise(emploi =fsum(Y1) )# on extrait un intervalle de confiance pour les données # projetée (15% - 85%). Pour cela, il faut calculer la somme régionale # d'emploi pour chacun des échantillon bayésien, puis effectuer la # l'agrégation pour obtenir les intervales de confiances# on identifie les colonnes contenant les échantillonssamples_cols <-names(post_samples)samples_cols <- samples_cols[grepl('sample_', samples_cols)]part2 <- post_samples %>%fsubset(AN >=2021) %>%# on applique uniquement aux données projetéesfselect(c("AN", samples_cols)) %>%# on sélectionner les colonnes nécessairespivot(how ='longer', values = samples_cols) %>%# on pivote au format longfgroup_by(AN,variable) %>%# pour chaque année, et chaque échantillonfsummarise( # on calcule la somme des emplois dan la région emploi =fsum(value) ) %>%ungroup() %>%fgroup_by(AN) %>%# les sommes (1500) obtenues sont organisées par fsummarise( # années afin de pouvoir calculer les intervalles low_est =fquantile(emploi, probs =0.15), # de confiancemed_est =fquantile(emploi, probs =0.5),high_est =fquantile(emploi, probs =0.85) )# pour assurer la continuité de la figure finale, nous devons # ajouter une ligne dans le tableau part2part2 <-rbind(part2, data.frame(AN =2016, low_est =last(part1$emploi),med_est =last(part1$emploi),high_est =last(part1$emploi) ))ggplot() +geom_line(data = part1, mapping =aes(x = AN, y = emploi)) +geom_point(data = part1, mapping =aes(x = AN, y = emploi)) +geom_ribbon(data = part2, mapping =aes(x = AN, ymin = low_est,ymax = high_est ),alpha =0.4) +geom_line(data = part2, mapping =aes(x = AN, y = med_est)) +geom_point(data = part2, mapping =aes(x = AN, y = med_est)) +geom_vline(xintercept =2016, linetype ='dashed') +theme_bw() +scale_x_continuous(breaks =seq(1996, 2041, by =5)) +labs(x ='année', y ="volume d'emplois")```
