# Mapping multiple trends with confidence

A tutorial to compute trends by groups and plot/map the results
R
datavisualization
OSM
spatial
modelling
Author

Michaël

Published

2019-02-06

Modified

2023-11-25

We will use `dplyr::nest` to create a list-column and will apply a model (with `purrr::map`) to each row, then we will extract each slope and its p value with `map` and `broom::tidy`.

## Setup

``````library(tidyverse)
library(httr)
library(sf)
library(janitor)
library(fs)
library(broom)
library(scales)
library(rnaturalearth)
library(glue)``````

## Data

Map data. Départements polygons from modified Adminexpress

``````dep <- st_read("~/data/adminexpress/adminexpress_cog_simpl_000_2022.gpkg",
"departement")``````
``````Reading layer `departement' from data source
using driver `GPKG'
Simple feature collection with 103 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -63.15272 ymin: -21.38883 xmax: 55.83661 ymax: 51.08876
Geodetic CRS:  WGS 84``````

Population data by département from INSEE.

``````if (!file_exists("pop.xlsx")) {
GET("https://www.insee.fr/fr/statistiques/fichier/2012713/TCRD_004.xlsx",
write_disk("pop.xlsx"))
}

pop <- read_xlsx("pop.xlsx", skip = 3, na = "nd") |>
clean_names() |>
rename(insee_dep = x1,
dep = x2,
x2023 = x2023_p) |>
select(-part_dans_la_population_francaise_p_en_percent) |>
pivot_longer(starts_with("x"), names_to = "annee", values_to = "pop") |>
mutate(annee = parse_number(annee))``````

## Plot

``````trends |>
ggplot(aes(fct_reorder(insee_dep, estimate), estimate,
color = if_else(p.value <= .05,
if_else(estimate >= 0, "positive", "négative"),
"n.s."))) +
geom_point() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .5) +
scale_color_manual(name = "tendance",
values = c("positive" = "red",
"n.s." = "lightgray",
"négative" = "blue")) +
scale_y_continuous(labels = number_format(big.mark = " ", style_positive = "plus")) +
coord_flip() +
labs(title = "Évolution des populations départementales françaises",
subtitle = glue("{year_min}-{year_max}"),
x = "dép.",
y = bquote(Delta[population] ~ (habitant %.% an^{-1})),
caption = glue("r.iresmi.net - {Sys.Date()}
données INSEE")) +
guides(color = guide_legend(reverse = TRUE)) +
theme(plot.caption = element_text(size = 7))``````

## Map

``````pop_dep <- dep |>
inner_join(trends, join_by(insee_dep)) |>
inner_join(filter(pop,
annee == year_max,
!insee_dep %in% c("F", "M")),
join_by(insee_dep)) |>
st_transform("EPSG:2154")

mean_fr <- trends |>
summarise(mean(estimate, na.rm = TRUE)) |>
pull()

world <- ne_countries(scale = "medium", returnclass = "sf") |>
filter(continent == "Europe") |>
st_transform("EPSG:2154")

pop_dep |>
ggplot() +
geom_sf(data = world, fill = "grey97", color = 0) +
geom_sf(color = "lightgrey", fill = "floralwhite", size = .2) +
stat_sf_coordinates(data = filter(pop_dep, p.value > .05),
aes(size = pop),
fill = "lightgrey", color = "lightgrey", shape = 21, alpha = 0.8) +
stat_sf_coordinates(data = filter(pop_dep, p.value <= .05),
aes(size = pop, fill = estimate),
color = "lightgrey", shape = 21, alpha = 0.8) +
coord_sf(xlim = c(160000, 1190000),
ylim = c(6085000, 7070000)) +
Delta ~ population[.(year_min)-.(year_max)],
(habitant %.% an^{-1}))),
moy. == .(sprintf("%+i", round(mean_fr))) ~ ", en gris : n.s.")),
labels = number_format(big.mark = " ", style_positive = "plus"),
low = "darkblue", mid = "white", high = "darkred", midpoint = mean_fr) +
scale_size_area(name = glue("habitants ({year_max})"),
labels = number_format(big.mark = " "),
max_size = 10) +
guides(size = guide_legend(reverse = TRUE)) +
labs(title = "Évolution des populations départementales françaises",
subtitle = glue("Métropole, {year_min}-{year_max}"),
x = "",
y = "",
caption = glue("r.iresmi.net - {Sys.Date()}
données INSEE {year_max}
fond cartographique : contributeurs Openstreetmap 2014, Naturalearth")) +
theme_void() +
theme(plot.caption = element_text(size = 7),
legend.text.align = 1)``````