Companion Notebook for An agent-based model to investigate the effects of urban segregation around the clock on inequalities in health behaviour
Clémentine Cottineau-Mugadza, Julien Perret, Romain Reuillon, Sébastien Rey-Coyrehourcq, Julie Vallée
2025-10-07
Libraries
library(tidyverse)
library(reshape2)
library(sf)
library(spdep)
library(tmap)
library(cowplot)
library(here)
library(OasisR) ## package OK for version R.4.4.3 (windows)Synthetic population - H24 library (section 2.1)
Import data for synthetic population located in ‘night’ cell
pop0, in ‘day’ cell pop1 and in ‘evening’ cell
pop2
pop_0 <- st_read( "data/population.gpkg",
layer="Population_0", quiet = TRUE)
pop_1 <- st_read("data/population.gpkg",
layer="Population_1", quiet = TRUE)
pop_2 <- st_read("data/population.gpkg",
layer="Population_2", quiet = TRUE)Create new colums : gender (2 groups) ; Age (3 groups); Educ (3 Groups)
pop_0 <- pop_0 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
pop = men + women,
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3'))),
pct_educ3 = (educ3 / pop) * 100,
pct_age3 = (age3 / pop) * 100,
pct_women = (women / pop) * 100
)
pop_1 <- pop_1 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)
pop_2 <- pop_2 %>%
rowwise() %>%
mutate(
men = sum(across(starts_with('pop_1'))),
women = sum(across(starts_with('pop_2'))),
age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
educ1 = sum(across(ends_with('_1'))),
educ2 = sum(across(ends_with('_2'))),
educ3 = sum(across(ends_with('_3')))
)Distribution of the synthetic population in night cells in the Paris Region (Figure 1)
Total Population density
bks_pop <- c(1, 10, 100, 1000, 10000, Inf)
tm_shape(pop_0) + tm_fill(col="pop", style="fixed",
breaks = bks_pop, palette="Greys") +
tm_legend(outside=TRUE)
Proportion of people with 2+ years of higher education
bks_edu <- c(0, 3.5, 7.7, 11.4, 17.2, 100)
tm_shape(pop_0) + tm_fill(col="pct_educ3", style="fixed",
breaks = bks_edu, palette="Reds") +
tm_legend(outside=TRUE)
Proportion of people aged 60+
bks_age <- c(0, 12.5, 17.1, 20.6, 25.1, 100)
tm_shape(pop_0) + tm_fill(col="pct_age3", style="fixed",
breaks = bks_age, palette="Blues") +
tm_legend(outside=TRUE)
Proportion of women
bks_sex <- c(0, 45.8, 50, 51.6, 54.5, 100)
tm_shape(pop_0) + tm_fill(col="pct_women", style="fixed",
breaks = bks_sex, palette="Greens") +
tm_legend(outside=TRUE)
Duncan Dissimilarity Index (Appendix 3)
Duncan’s segregation index is one-group form of dissimilarity index DI and measures the unevenness of a group distribution compared to the rest of the population. It can be interpreted as the share of the group that would have to move to achieve an even distribution compared to the rest of the population. You should not include a column with total population, because this will be interpreted as a group.
In night cells
dfpop_0 <- pop_0 %>% st_drop_geometry()
pop_0_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_0[,2:19]))
colnames(pop_0_Duncan_18categ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_18categ) <- paste (colnames(dfpop_0[ , 2:19]))
pop_0_Duncan_sex<- as.data.frame(ISDuncan(dfpop_0[,20:21]))
colnames(pop_0_Duncan_sex)[1] <- "Duncan_night"
rownames(pop_0_Duncan_sex) <- paste (colnames(dfpop_0[ , 20:21]))
pop_0_Duncan_age<- as.data.frame(ISDuncan(dfpop_0[,22:24]))
colnames(pop_0_Duncan_age)[1] <- "Duncan_night"
rownames(pop_0_Duncan_age) <- paste (colnames(dfpop_0[ , 22:24]))
pop_0_Duncan_educ<- as.data.frame(ISDuncan(dfpop_0[,25:27]))
colnames(pop_0_Duncan_educ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_educ) <- paste (colnames(dfpop_0[ , 25:27]))
pop_0_all_Duncan <- bind_rows(pop_0_Duncan_18categ, pop_0_Duncan_sex, pop_0_Duncan_age, pop_0_Duncan_educ)
pop_0_all_Duncan_round <- pop_0_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))In day cells
dfpop_1 <- pop_1 %>% st_drop_geometry()
pop_1_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_1[,2:19]))
colnames(pop_1_Duncan_18categ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_18categ) <- paste (colnames(dfpop_1[ , 2:19]))
pop_1_Duncan_sex<- as.data.frame(ISDuncan(dfpop_1[,20:21]))
colnames(pop_1_Duncan_sex)[1] <- "Duncan_day"
rownames(pop_1_Duncan_sex) <- paste (colnames(dfpop_1[ , 20:21]))
pop_1_Duncan_age<- as.data.frame(ISDuncan(dfpop_1[,22:24]))
colnames(pop_1_Duncan_age)[1] <- "Duncan_day"
rownames(pop_1_Duncan_age) <- paste (colnames(dfpop_1[ , 22:24]))
pop_1_Duncan_educ<- as.data.frame(ISDuncan(dfpop_1[,25:27]))
colnames(pop_1_Duncan_educ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_educ) <- paste (colnames(dfpop_1[ , 25:27]))
pop_1_all_Duncan <- bind_rows(pop_1_Duncan_18categ, pop_1_Duncan_sex, pop_1_Duncan_age, pop_1_Duncan_educ)
pop_1_all_Duncan_round <- pop_1_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))In evening cells
dfpop_2 <- pop_2 %>% st_drop_geometry()
pop_2_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_2[,2:19]))
colnames(pop_2_Duncan_18categ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_18categ) <- paste (colnames(dfpop_2[ , 2:19]))
pop_2_Duncan_sex<- as.data.frame(ISDuncan(dfpop_2[,20:21]))
colnames(pop_2_Duncan_sex)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_sex) <- paste (colnames(dfpop_2[ , 20:21]))
pop_2_Duncan_age<- as.data.frame(ISDuncan(dfpop_2[,22:24]))
colnames(pop_2_Duncan_age)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_age) <- paste (colnames(dfpop_2[ , 22:24]))
pop_2_Duncan_educ<- as.data.frame(ISDuncan(dfpop_2[,25:27]))
colnames(pop_2_Duncan_educ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_educ) <- paste (colnames(dfpop_2[ , 25:27]))
pop_2_all_Duncan <- bind_rows(pop_2_Duncan_18categ, pop_2_Duncan_sex, pop_2_Duncan_age, pop_2_Duncan_educ)
pop_2_all_Duncan_round <- pop_2_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))Build table
pop_all_Duncan_round <- bind_cols(pop_0_all_Duncan_round, pop_1_all_Duncan_round, pop_2_all_Duncan_round)
pop_all_Duncan_round ## Duncan_night Duncan_day Duncan_evening
## pop_1_1_1 0.152 0.549 0.548
## pop_1_1_2 0.082 0.442 0.423
## pop_1_1_3 0.257 0.570 0.550
## pop_1_2_1 0.182 0.394 0.372
## pop_1_2_2 0.091 0.391 0.345
## pop_1_2_3 0.274 0.437 0.386
## pop_1_3_1 0.184 0.564 0.553
## pop_1_3_2 0.121 0.665 0.698
## pop_1_3_3 0.284 0.528 0.560
## pop_2_1_1 0.138 0.558 0.559
## pop_2_1_2 0.088 0.445 0.396
## pop_2_1_3 0.275 0.492 0.493
## pop_2_2_1 0.160 0.371 0.362
## pop_2_2_2 0.086 0.364 0.320
## pop_2_2_3 0.263 0.375 0.373
## pop_2_3_1 0.150 0.501 0.481
## pop_2_3_2 0.131 0.573 0.592
## pop_2_3_3 0.295 0.570 0.605
## men 0.033 0.230 0.159
## women 0.033 0.230 0.159
## age1 0.076 0.297 0.262
## age2 0.052 0.280 0.224
## age3 0.095 0.370 0.345
## educ1 0.241 0.349 0.361
## educ2 0.067 0.252 0.236
## educ3 0.322 0.358 0.404Moran’s index of spatial autocorrelation (Appendix 3)
In night cells
dfpop_0_prop <- dfpop_0/dfpop_0$pop
pop_0_geom <- pop_0[, -1:-28]
pop_0_prop <- bind_cols(pop_0_geom, dfpop_0_prop)
nb0 <- poly2nb(pop_0_prop, queen=TRUE)
lw0 <- nb2listw(nb0, style="W", zero.policy=TRUE)
list <- colnames(dfpop_0_prop[,2:27])
pop_0_moran <- data.frame(matrix(ncol = 2))
colnames(pop_0_moran) <- c("pop_night_IMoran" , "pop_night_Moranpvalue")
for (i in list)
{
pop_0_i <- pop_0_prop[,i]
colnames(pop_0_i)[1] ="ok"
pop_0_moran[i, 1] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [3]
pop_0_moran[i, 2] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [2]
}
pop_0_moran <- pop_0_moran[-1,] # Delete empty first line
pop_0_moran_round <- pop_0_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_0_moran_round## pop_night_IMoran pop_night_Moranpvalue
## pop_1_1_1 0.022 0.000
## pop_1_1_2 0.035 0.000
## pop_1_1_3 0.069 0.000
## pop_1_2_1 0.100 0.000
## pop_1_2_2 0.003 0.327
## pop_1_2_3 0.198 0.000
## pop_1_3_1 0.027 0.000
## pop_1_3_2 0.013 0.026
## pop_1_3_3 0.068 0.000
## pop_2_1_1 0.020 0.002
## pop_2_1_2 0.040 0.000
## pop_2_1_3 0.031 0.000
## pop_2_2_1 0.066 0.000
## pop_2_2_2 0.010 0.074
## pop_2_2_3 0.110 0.000
## pop_2_3_1 0.023 0.000
## pop_2_3_2 0.029 0.000
## pop_2_3_3 0.050 0.000
## men 0.018 0.004
## women 0.018 0.004
## age1 0.076 0.000
## age2 0.026 0.000
## age3 0.068 0.000
## educ1 0.303 0.000
## educ2 0.071 0.000
## educ3 0.367 0.000In day cells
dfpop_1_prop <- dfpop_1/dfpop_1$pop
pop_1_geom <- pop_1[, -1:-28]
pop_1_prop <- bind_cols(pop_1_geom, dfpop_1_prop)
nb1 <- poly2nb(pop_1_prop, queen=TRUE)
lw1 <- nb2listw(nb1, style="W", zero.policy=TRUE)
list <- colnames(dfpop_1_prop[,2:27])
pop_1_moran <- data.frame(matrix(ncol = 2))
colnames(pop_1_moran) <- c( "pop_day_IMoran" , "pop_day_Moranpvalue")
for (i in list)
{
pop_1_i <- pop_1_prop[,i]
colnames(pop_1_i)[1] ="ok"
pop_1_moran[i, 1] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [3]
pop_1_moran[i, 2] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [2]
}
pop_1_moran <- pop_1_moran[-1,] # Delete empty first line
pop_1_moran_round <- pop_1_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_1_moran_round## pop_day_IMoran pop_day_Moranpvalue
## pop_1_1_1 0.182 0
## pop_1_1_2 0.197 0
## pop_1_1_3 0.184 0
## pop_1_2_1 0.085 0
## pop_1_2_2 0.060 0
## pop_1_2_3 0.062 0
## pop_1_3_1 0.187 0
## pop_1_3_2 0.331 0
## pop_1_3_3 0.110 0
## pop_2_1_1 0.189 0
## pop_2_1_2 0.126 0
## pop_2_1_3 0.104 0
## pop_2_2_1 0.060 0
## pop_2_2_2 0.036 0
## pop_2_2_3 0.055 0
## pop_2_3_1 0.086 0
## pop_2_3_2 0.325 0
## pop_2_3_3 0.225 0
## men 0.150 0
## women 0.150 0
## age1 0.175 0
## age2 0.239 0
## age3 0.279 0
## educ1 0.226 0
## educ2 0.201 0
## educ3 0.172 0In evening cells
dfpop_2_prop <- dfpop_2/dfpop_2$pop
pop_2_geom <- pop_2[, -1:-28]
pop_2_prop <- bind_cols(pop_2_geom, dfpop_2_prop)
nb2 <- poly2nb(pop_2_prop, queen=TRUE)
lw2 <- nb2listw(nb2, style="W", zero.policy=TRUE)
list <- colnames(dfpop_2_prop[,2:27])
pop_2_moran <- data.frame(matrix(ncol = 2))
colnames(pop_2_moran) <- c( "pop_evening_IMoran" , "pop_evening_Moranpvalue")
for (i in list)
{
pop_2_i <- pop_2_prop[,i]
colnames(pop_2_i)[1] ="ok"
pop_2_moran[i, 1] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [3]
pop_2_moran[i, 2] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [2]
}
pop_2_moran <- pop_2_moran[-1,] # Delete empty first line
pop_2_moran_round <- pop_2_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_2_moran_round## pop_evening_IMoran pop_evening_Moranpvalue
## pop_1_1_1 0.125 0
## pop_1_1_2 0.193 0
## pop_1_1_3 0.113 0
## pop_1_2_1 0.088 0
## pop_1_2_2 0.078 0
## pop_1_2_3 0.074 0
## pop_1_3_1 0.196 0
## pop_1_3_2 0.351 0
## pop_1_3_3 0.107 0
## pop_2_1_1 0.192 0
## pop_2_1_2 0.112 0
## pop_2_1_3 0.077 0
## pop_2_2_1 0.072 0
## pop_2_2_2 0.057 0
## pop_2_2_3 0.054 0
## pop_2_3_1 0.095 0
## pop_2_3_2 0.352 0
## pop_2_3_3 0.232 0
## men 0.148 0
## women 0.148 0
## age1 0.145 0
## age2 0.262 0
## age3 0.306 0
## educ1 0.201 0
## educ2 0.193 0
## educ3 0.166 0Combine Moran tables pop_0, pop_1, pop_2
Analysis of empirical data about dietary habits (section 2.2)
Summary by sociodemographic groups (cf. table 1)
data <- read.csv("data/sociodemo_distribution_healthy.csv", stringsAsFactors = F)
head(data)## Sex Age Edu n_idf N_2002 N_2008 conso_5_2002 conso_5_2008
## 1 1 1 1 493871 52 148 0.01923077 0.02702703
## 2 1 1 2 496964 97 223 0.03092784 0.04932735
## 3 1 1 3 199476 84 137 0.03571429 0.07299270
## 4 1 2 1 1077670 241 355 0.04979253 0.08169014
## 5 1 2 2 662407 122 225 0.07377049 0.11111111
## 6 1 2 3 658465 118 180 0.05932203 0.16666667n_idf = number of individuals of a particular group in
2012 in the Paris Region (Ile de France). N_2002 = number
of individuals of a particular group in the Health and Nutrition
Barometer survey in 2002. N_2008 = number of individuals of
a particular group in the Health and Nutrition Barometer in 2008.
conso_5_2002 = proportion of individuals eating 5+ portions
of fruit and vegetables a day in the Health and Nutrition Barometer in
2002. conso_5_2008 = proportion of individuals eating 5+
portions of fruit and vegetables a day in the Health and Nutrition
Barometer in 2008.
Computation of the total proportion of healthy individuals
i.e. eating 5+ portions of fruit and vegetables a day)
data$conso_5_2002 <- as.numeric(data$conso_5_2002)
data$conso_5_2008 <- as.numeric(data$conso_5_2008)
n_tot_idf <- sum(data$n_idf)
data$share_idf <- data$n_idf / n_tot_idf
data$h_idf_2002 <- round(data$n_idf * data$conso_5_2002,0)
data$h_idf_2008 <- round(data$n_idf * data$conso_5_2008,0)
data$h_baro_2002 <- round(data$N_2002 * data$conso_5_2002,0)
data$h_baro_2008 <- round(data$N_2008 * data$conso_5_2008,0)
total_prop <- function(data, year, sample){
if(sample == "idf"){
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data$n_idf)
} else {
tp <-
sum(data[,paste("h", sample, year, sep="_")]) /
sum(data[,paste("N", year, sep="_")])
}
return(tp)
}Results with the Paris region as reference population
total_prop(data=data, year = 2002, sample = "idf")## [1] 0.09568032total_prop(data=data, year = 2008, sample = "idf")## [1] 0.1205141Results with the diet survey as reference population
total_prop(data=data, year = 2002, sample = "baro")## [1] 0.1063931total_prop(data=data, year = 2008, sample = "baro")## [1] 0.1201574Computation of EducIneqIndex
n_tot_baro_2002 <- sum(data$N_2002)
n_tot_baro_2008 <- sum(data$N_2008)
EII <- function(data, year, sample){
if(sample == "idf"){
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, "n_idf"] / n_tot_idf)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, "n_idf"] / n_tot_idf)
)
} else {
if(year == 2002){ n_tot_baro <- n_tot_baro_2002 }
if(year == 2008){ n_tot_baro <- n_tot_baro_2008 }
eii <- (
data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 1 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 1, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 2, paste0("N_", year)] / n_tot_baro)
) + (
data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
paste0("conso_5_", year)] /
data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
paste0("conso_5_", year)]
) * (
sum(data[data$Sex == 2 & data$Age == 3, paste0("N_", year)] / n_tot_baro)
)
}
return(eii)
}Inequality results with the Paris region as reference population
EII(data=data, year = 2002, sample = "idf")## [1] 1.407512EII(data=data, year = 2008, sample = "idf")## [1] 1.89476Inequality results with the dietary survey as reference population
EII(data=data, year = 2002, sample = "baro")## [1] 1.382566EII(data=data, year = 2008, sample = "baro")## [1] 1.905816Evaluation of 5aday ABM model (section 3.4)
Importing calibration results
Pareto front (Figure 4)
pareto$reliability.sample <- ifelse(pareto$evolution.samples>5, 2,1)
pareto$objective.deltaHealth <- as.numeric(gsub(",","",pareto$objective.deltaHealth))
front <- ggplot(data=pareto, aes(x = objective.deltaHealth, y=objective.deltaSocialInequality))
front +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) +
scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
geom_vline(xintercept = 290000000) +
geom_hline(yintercept = 0.36) +
theme_bw() 
Distribution of parameter values along the Pareto front (Appendix 6)
pareto$id <- rownames(pareto)
colnames(pareto)## [1] "evolution.generation" "evolution.evaluated"
## [3] "evolution.samples" "maxProbaToSwitch"
## [5] "constraintsStrength" "inertiaCoefficient"
## [7] "healthyDietReward" "objective.deltaHealth"
## [9] "objective.deltaSocialInequality" "X.Healthy_vs._2002"
## [11] "X.Healthy_vs._2008" "socialIneq_vs._2002"
## [13] "typo" "typo_name"
## [15] "socialInequality" "deltaSocialInequality"
## [17] "opinionMedian" "numberOfHealthy"
## [19] "deltaHealth" "reliability.sample"
## [21] "id"pareto_long <- melt(pareto[,c("maxProbaToSwitch","constraintsStrength" ,
"inertiaCoefficient","healthyDietReward",
"objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample")],
id.vars=c("objective.deltaSocialInequality", "objective.deltaHealth",
"typo_name", "reliability.sample"))
ggplot(pareto_long, aes(objective.deltaHealth, value)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none") +
facet_wrap(~variable, scales = "free_y", ncol = 1)
ggplot(pareto_long, aes(value, objective.deltaSocialInequality)) +
geom_point(aes(group = typo_name, colour = typo_name,
size=reliability.sample)) +
scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
c("#faa911", #orange
"#07b55e", #green
"#1088b0", #blue
"#aae860" #light green
), .5)) +
theme_bw() + theme(legend.position="none")+
facet_wrap(~variable, scales = "free_x", nrow = 1)
Results (section 4)
Sensitivity to random moves (Figure 6)
random_moves <- read.csv("data/random_moves.csv")
head(random_moves)## erreygersE numberOfHealthy openmole.experiment randomMoveRatio seed
## 1 0.03682734 1150567 706 1 -7.331534e+18
## 2 0.03661155 1158102 700 1 -3.395561e+18
## 3 0.03655777 1151786 681 1 8.946059e+18
## 4 0.03665694 1151391 663 1 3.690077e+17
## 5 0.03652460 1151135 672 1 3.693937e+18
## 6 0.03671739 1153199 675 1 4.750891e+18
## socialInequality
## 1 1.429419
## 2 1.418540
## 3 1.422918
## 4 1.423681
## 5 1.421879
## 6 1.424769annotation <- data.frame(
x = c(0,1),
y = c(1.57,1.46),
label = c("Scenario 2C", "Scenario 2B")
)
ggplot(random_moves, aes(x=randomMoveRatio, y=socialInequality)) +
geom_hline(yintercept = 1.424, color = "darkorange2") +
geom_hline(yintercept = 1.536, color = "springgreen3") +
geom_violin(aes(group=randomMoveRatio, fill=randomMoveRatio)) +
geom_point(size=0.05) +
scale_fill_gradient(
low = "springgreen3",
high = "darkorange2"
) +
geom_text(data=annotation, aes(x=x, y=y, label=label),
colour = c("springgreen3","darkorange2"),
size=4 , angle=0, fontface="bold") +
ylim(1.35,1.6) +
guides(fill="none") +
xlab("Proportion of random moves") +
ylab("Social inequality index (EII)")+
theme_bw() 
Proportions of healthy agents for the 5 scenarios (Figure 7)
mean_1A <- mean(sim[sim$numOfScenario == 1,"numberOfHealthy"]/8778898)
mean_1B <- mean(sim[sim$numOfScenario == 2,"numberOfHealthy"]/8778898)
mean_2A <- mean(sim[sim$numOfScenario == 3,"numberOfHealthy"]/8778898)
mean_2B <- mean(sim[sim$numOfScenario == 4,"numberOfHealthy"]/8778898)
mean_2C <- mean(sim[sim$numOfScenario == 5,"numberOfHealthy"]/8778898)
sim %>%
ggplot( aes(x=numberOfHealthy/8778898, group=Scenario, fill=Scenario)) +
geom_vline(xintercept=0.0957, size=1, color="grey30", linetype = "dashed") +
geom_density(alpha=0.6) +
scale_fill_manual(values=c("firebrick2", "mediumorchid4", "dodgerblue2", "darkorange2", "springgreen3"))+
theme_bw() +
geom_vline(xintercept=mean_1A, size=0.4, color="firebrick2") +
geom_vline(xintercept=mean_1B, size=0.4, color="mediumorchid4") +
geom_vline(xintercept=mean_2A, size=0.4, color="dodgerblue2") +
geom_vline(xintercept=mean_2B, size=0.4, color="darkorange2") +
geom_vline(xintercept=mean_2C, size=0.4, color="springgreen3") +
annotate("text", x=mean_1A + 0.001, y= -50, size=3, label=round(mean_1A,3), color="firebrick2") +
annotate("text", x=mean_1B - 0.001 , y= -50, size=3, label=round(mean_1B,3), color="mediumorchid4") +
annotate("text", x=mean_2A + 0.001 , y= -50, size=3, label=round(mean_2A,3), color="dodgerblue2") +
annotate("text", x=mean_2B + 0.001 , y= -50, size=3, label=round(mean_2B,3), color="darkorange2") +
annotate("text", x=mean_2C + 0.0005, y= -50, size=3, label=round(mean_2C,3), color="springgreen3") +
annotate("text", x=0.0957 + 0.005, y= 115, size=3, label="Value at initialisation (0.0957)", color="grey30") +
xlab("Share of agents with a healthy diet") +
ylab("Number of simulations") +
xlim(0.07,0.15) +
guides(fill=guide_legend(title="Type of scenario",nrow=2,byrow=TRUE)) +
theme(legend.position = "bottom") 
Map of healthy agents in their night cells for the 5 scenarios (Figure 8)
s0r <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s0o <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Initial_0", quiet = TRUE)
s1A <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s1B <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2A <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg",
layer="Final", quiet = TRUE)
s2B <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg",
layer="Final", quiet = TRUE)
s2C <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
layer="Final", quiet = TRUE)
bks <- c(0, 0.1, 0.125, 0.15, 0.175, Inf)Maps
Initial Random
tm_shape(s0r) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Initial Observed
tm_shape(s0o) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1A
tm_shape(s1A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 1B
tm_shape(s1B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2A
tm_shape(s2A) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2B
tm_shape(s2B) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Scenario 2C
tm_shape(s2C) + tm_fill(col="X.propHealthy", style="fixed",
breaks = bks, palette="Greens") +
tm_legend(outside=TRUE)
Duncan Dissimilarity Index - night cells (add in Figure 8)
duncan <- data.frame(matrix(ncol = 2, nrow=7))
colnames(duncan) <- c("type", "Duncan" )
dfs0r<- s0r %>% st_drop_geometry()
duncan[1,1] <- "s0r_night"
duncan[1,2]<- ISDuncan(dfs0r[,3:4]) [1]
dfs0o<- s0o %>% st_drop_geometry()
duncan[2,1] <- "s0o_night"
duncan[2,2]<- ISDuncan(dfs0o[,3:4]) [1]
dfs1A<- s1A %>% st_drop_geometry()
duncan[3,1] <- "s1A_final_night"
duncan[3,2]<- ISDuncan(dfs1A[,3:4]) [1]
dfs1B<- s1B %>% st_drop_geometry()
duncan[4,1] <- "s1B_final_night"
duncan[4,2]<- ISDuncan(dfs1B[,3:4]) [1]
dfs2A<- s2A %>% st_drop_geometry()
duncan[5,1] <- "s2A_final_night"
duncan[5,2]<- ISDuncan(dfs2A[,3:4]) [1]
dfs2B<- s2B %>% st_drop_geometry()
duncan[6,1] <- "s2B_final_night"
duncan[6,2]<- ISDuncan(dfs2B[,3:4]) [1]
dfs2C<- s2C %>% st_drop_geometry()
duncan[7,1] <- "s2C_final_night"
duncan[7,2]<- ISDuncan(dfs2C[,3:4]) [1]
duncan_round <- duncan %>% mutate(across(starts_with("Duncan"), round, 3))
duncan_round## type Duncan
## 1 s0r_night 0.040
## 2 s0o_night 0.038
## 3 s1A_final_night 0.084
## 4 s1B_final_night 0.043
## 5 s2A_final_night 0.098
## 6 s2B_final_night 0.045
## 7 s2C_final_night 0.063Moran’s index of spatial autocorrelation - night cells (add in Figure 8)
moran <- data.frame(matrix(ncol = 3, nrow=7))
colnames(moran) <- c("type", "IMoran" , "pvalueMoran")
nbs0r <- poly2nb(s0r, queen=TRUE)
lws0r <- nb2listw(nbs0r, style="W", zero.policy=TRUE)
moran[1,1] <- "s0r_night"
moran[1,2] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [3]
moran[1,3] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [2]
nbs0o <- poly2nb(s0o, queen=TRUE)
lws0o <- nb2listw(nbs0o, style="W", zero.policy=TRUE)
moran[4,1] <- "s0o_night"
moran[4,2] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [3]
moran[4,3] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [2]
nbs1A <- poly2nb(s1A, queen=TRUE)
lws1A <- nb2listw(nbs1A, style="W", zero.policy=TRUE)
moran[2,1] <- "s1A_final_night"
moran[2,2] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [3]
moran[2,3] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [2]
nbs1B <- poly2nb(s1B, queen=TRUE)
lws1B <- nb2listw(nbs1B, style="W", zero.policy=TRUE)
moran[3,1] <- "s1B_final_night"
moran[3,2] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [3]
moran[3,3] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [2]
nbs2A <- poly2nb(s2A, queen=TRUE)
lws2A <- nb2listw(nbs2A, style="W", zero.policy=TRUE)
moran[5,1] <- "s2A_final_night"
moran[5,2] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [3]
moran[5,3] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [2]
nbs2B <- poly2nb(s2B, queen=TRUE)
lws2B <- nb2listw(nbs2B, style="W", zero.policy=TRUE)
moran[6,1] <- "s2B_final_night"
moran[6,2] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [3]
moran[6,3] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [2]
nbs2C <- poly2nb(s2C, queen=TRUE)
lws2C <- nb2listw(nbs2C, style="W", zero.policy=TRUE)
moran[7,1] <- "s2C_final_night"
moran[7,2] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [3]
moran[7,3] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [2]
moran_round <- moran %>% mutate(across(2:3, round, 3))
moran_round ## type IMoran pvalueMoran
## 1 s0r_night 0.010 0.061
## 2 s1A_final_night 0.010 0.057
## 3 s1B_final_night 0.016 0.008
## 4 s0o_night 0.006 0.175
## 5 s2A_final_night 0.033 0.000
## 6 s2B_final_night 0.024 0.000
## 7 s2C_final_night 0.009 0.091Appendices (Supp. Material)
Alternative measure of social inequality in health behaviour (Appendix 4)
Computation of rank-ordered inequality measure
data_edu <- aggregate(data[,c("n_idf", "N_2002", "N_2008",
"h_idf_2002", "h_idf_2008",
"h_baro_2002", "h_baro_2008")],
by = list(data$Edu), FUN = sum)
colnames(data_edu)[1] <- "edu"
data_edu## edu n_idf N_2002 N_2008 h_idf_2002 h_idf_2008 h_baro_2002 h_baro_2008
## 1 1 3800999 1030 1372 336079 359509 122 155
## 2 2 2946028 595 1135 279323 361858 57 132
## 3 3 2021871 471 797 223609 335409 44 110guido <- function(data_edu, year, sample){
if(sample == "idf"){
n <- sum(data_edu[,"n_idf"])
n1 <- data_edu[data_edu$edu == 1,"n_idf"]
n2 <- data_edu[data_edu$edu == 2,"n_idf"]
n3 <- data_edu[data_edu$edu == 3,"n_idf"]
} else {
n <- sum(data_edu[,paste("N", year, sep="_")])
n1 <- data_edu[data_edu$edu == 1,paste("N", year, sep="_")]
n2 <- data_edu[data_edu$edu == 2,paste("N", year, sep="_")]
n3 <- data_edu[data_edu$edu == 3,paste("N", year, sep="_")]
}
h1 <- data_edu[data_edu$edu == 1,paste("h", sample, year, sep="_")]
h2 <- data_edu[data_edu$edu == 2,paste("h", sample, year, sep="_")]
h3 <- data_edu[data_edu$edu == 3,paste("h", sample, year, sep="_")]
t1 <- (1 - n2 - n3) / 2
t2 <- (n1 - n3 + 1) / 2
t3 <- (n1 + n2 + 1) / 2
guido <- ( 8 / n^2 ) * (h1 * t1 + h2 * t2 + h3 * t3)
return(guido)
}
guido(data_edu = data_edu, year = 2002, sample = "idf")## [1] 0.01748087guido(data_edu = data_edu, year = 2008, sample = "idf")## [1] 0.05830406guido(data_edu = data_edu, year = 2002, sample = "baro")## [1] -0.02409715guido(data_edu = data_edu, year = 2008, sample = "baro")## [1] 0.01927629Density distribution of rank-order inequality index per scenario (Appendix 4)
mean_1A <- mean(sim[sim$numOfScenario == 1,"erreygersE"])
mean_1B <- mean(sim[sim$numOfScenario == 2,"erreygersE"])
mean_2A <- mean(sim[sim$numOfScenario == 3,"erreygersE"])
mean_2B <- mean(sim[sim$numOfScenario == 4,"erreygersE"])
mean_2C <- mean(sim[sim$numOfScenario == 5,"erreygersE"])
sim %>%
ggplot( aes(x=erreygersE, group=Scenario, fill=Scenario)) +
geom_vline(xintercept=0.0175, size=1, color="grey30", linetype = "dashed") +
geom_density(alpha=0.6) +
scale_fill_manual(values=c("firebrick2", "mediumorchid4", "dodgerblue2", "darkorange2", "springgreen3"))+
theme_bw() +
geom_vline(xintercept=mean_1A, size=0.4, color="firebrick2") +
geom_vline(xintercept=mean_1B, size=0.4, color="mediumorchid4") +
geom_vline(xintercept=mean_2A, size=0.4, color="dodgerblue2") +
geom_vline(xintercept=mean_2B, size=0.4, color="darkorange2") +
geom_vline(xintercept=mean_2C, size=0.4, color="springgreen3") +
annotate("text", x=mean_1A - 0.0015, y= -50, size=3, label=round(mean_1A,3), color="firebrick2") +
annotate("text", x=mean_1B - 0.0025 , y= -50, size=3, label=round(mean_1B,3), color="mediumorchid4") +
annotate("text", x=mean_2A + 0.001 , y= -50, size=3, label=round(mean_2A,3), color="dodgerblue2") +
annotate("text", x=mean_2B + 0.001 , y= -50, size=3, label=round(mean_2B,3), color="darkorange2") +
annotate("text", x=mean_2C - 0.001, y= -50, size=3, label=round(mean_2C,3), color="springgreen3") +
annotate("text", x= 0.0175 + 0.003, y= 115, size=3, label="Value at initialisation", color="grey30") +
xlab("Rank-Order Inequality Index") +
ylab("Number of simulations") +
guides(fill=guide_legend(title="Type of scenario", nrow=2, byrow=TRUE)) +
theme(legend.position = "bottom") 
Evolution in simulated values over the six series of three time slices (Appendix 8)
Part A & Part B
phase <- read.csv("data/phase_48.csv")
long_table <- data.frame()
i <- 0
j <- 1
time <- data.frame(
day = rep("day",18),
nd = sort(rep(0:5,3)),
step = rep("step",18),
ns = rep(0:2,6)) |>
mutate(time = paste0(day, nd, "_", step, ns)) |>
select(time)
for(i in 1:nrow(phase)){
nsenario <- phase[i, "numOfScenario"]
EII <- t(
str_split(
str_replace(
str_replace( phase[i, "socialInequality"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
healthy <- t(
str_split(
str_replace(
str_replace( phase[i, "numberOfHealthy"] , c("\\["), ""),
c("\\["), ""
),
",", simplify = T
)
)[1:18,]
long_table[j:(j+17),"time"] <- time
long_table[j:(j+17),"scenario"] <- rep(nsenario, 18)
long_table[j:(j+17),"socialInequality"] <- EII
long_table[j:(j+17),"numberOfHealthy"] <- healthy
j <- j+18
}
plot_table <- long_table |>
mutate(prop_healthy = as.numeric(numberOfHealthy) / 8768898,
EII = as.numeric(socialInequality),
scenario = as.factor(ifelse(scenario == 2, "1B",
ifelse(scenario == 5, "2C", NA))
),
group = paste0(time, scenario)
)
ggplot(plot_table,
aes(x = time, y = EII,
colour = scenario, fill = scenario)) +
geom_violin(aes(group=group), alpha = 0.6) +
# geom_point() +
scale_colour_manual(values = c("mediumorchid4", "darkorange2")) +
scale_fill_manual(values = c("mediumorchid4","darkorange2")) +
xlab("Time") +
ylab("B. Social inequality index (EII)") +
theme(legend.position = "bottom") +
ylim(1.35,1.6) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)))
ggplot(plot_table,
aes(x = time, y = prop_healthy,
colour = scenario, fill = scenario)) +
geom_violin(aes(group=group), alpha = 0.6) +
# geom_point() +
scale_colour_manual(values = c("mediumorchid4", "darkorange2")) +
scale_fill_manual(values = c("mediumorchid4","darkorange2")) +
xlab("Time") +
ylab("A. Proportion of healthy agents") +
theme(legend.position = "bottom") +
ylim(0.07,0.15) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)))
Part C & Part D
Simulated proportion of healthy agents and simulated opinion across the 18 sociodemographic categories
specialRead <- function(x) {
read_csv(file=x, col_names = TRUE, col_types= "nnnnnnnnnn")
}
concatFiles <- function (base, folder){
fullPath = paste(base,folder,sep="")
list.files(path = fullPath, full.names = TRUE, pattern = "(\\d_\\d_\\d)") %>% lapply(specialRead) %>% bind_rows
}
recode <- function(x) {
step_to_day <- c (
'00' = 'day 0 step 0',
'01' = 'day 0 step 1',
'02' = 'day 0 step 2',
'10' = 'day 1 step 0',
'11' = 'day 1 step 1',
'12' = 'day 1 step 2',
'20' = 'day 2 step 0',
'21' = 'day 2 step 1',
'22' = 'day 2 step 2',
'30' = 'day 3 step 0',
'31' = 'day 3 step 1',
'32' = 'day 3 step 2',
'40' = 'day 4 step 0',
'41' = 'day 4 step 1',
'42' = 'day 4 step 2',
'50' = 'day 5 step 0',
'51' = 'day 5 step 1',
'52' = 'day 5 step 2'
)
return(step_to_day[x])
}sortie_sc5_higherprop <- concatFiles("data", "/HigherProp")
sortie_sc5_higherprop$propHealthy <- sortie_sc5_higherprop$healthy / sortie_sc5_higherprop$effective
sortieLong_sc5_hp <- gather(sortie_sc5_higherprop, facteur, valeur, c(propHealthy, avgOpinion) )
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, "ds" = gsub(" ", "", paste(paste("[",paste(sortieLong_sc5_hp$day, sortieLong_sc5_hp$slice)),"]")))
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, ds= recode(paste0(day, slice)))
age.l <- c("15-29 yrs", "30-59 yrs", "60-75 yrs")
names(age.l) <- c("1", "2", "3")
sex.l <- c("Male", "Female")
names(sex.l) <- c("1", "2")
p1 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="propHealthy"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p1 <- p1 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p1 <- p1 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, , size = 7))
p1 <- p1 + labs(title="C. Proportion of healthy by gender, age and education level", y = "Value", x = "Time")
p1 <- p1 + geom_line(size=0.5) + facet_grid( sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p1 <- p1 + coord_cartesian( ylim = c(0, 0.3))
p2 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="avgOpinion"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p2 <- p2 + scale_color_discrete(
name = "Education level",
labels = c("low", "middle", "high")
)
p2 <- p2 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 7))
p2 <- p2 + labs(title="D. Average Opinion by gender, age and education level", color = "Education level", y = "Value", x = "Time")
p2 <- p2 + geom_line(size=0.5) + facet_grid(sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p2 <- p2 + coord_cartesian( ylim = c(0.3, 0.7))
propHealtyAndOpinion_higherprop <- plot_grid(
p1, p2,
ncol = 1)
print(propHealtyAndOpinion_higherprop)
##ggsave("out/propHealtyAndOpinion_higherprop.pdf")———UNUSED———-
Spatial distribution of agents with healthy behaviour - in day and evening cells
In day cells
s1A_day <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s1B_day <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2A_day <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2B_day <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_1")
s2C_day <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_1")Duncan Dissimilarity Index
duncan_day <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_day) <- c("type", "Duncan_dayCell" )
dfs1A_day<- s1A_day %>% st_drop_geometry()
duncan_day[1,1] <- "s1A_final_day"
duncan_day[1,2]<- ISDuncan(dfs1A_day[,3:4]) [1]
dfs1B_day<- s1B_day %>% st_drop_geometry()
duncan_day[2,1] <- "s1B_final_day"
duncan_day[2,2]<- ISDuncan(dfs1B_day[,3:4]) [1]
dfs2A_day<- s2A_day %>% st_drop_geometry()
duncan_day[3,1] <- "s2A_final_day"
duncan_day[3,2]<- ISDuncan(dfs2A_day[,3:4]) [1]
dfs2B_day<- s2B_day %>% st_drop_geometry()
duncan_day[4,1] <- "s2B_final_day"
duncan_day[4,2]<- ISDuncan(dfs2B_day[,3:4]) [1]
dfs2C_day<- s2C_day %>% st_drop_geometry()
duncan_day[5,1] <- "s2C_final_day"
duncan_day[5,2]<- ISDuncan(dfs2C_day[,3:4]) [1]
duncan_day <- duncan_day %>% mutate(across(starts_with("Duncan"), round, 3))Moran’s index of spatial autocorrelation
moran_day <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_day) <- c("type", "IMoran_dayCell" , "pvalueMoran_dayCell")
nbs1A_day <- poly2nb(s1A_day, queen=TRUE)
lws1A_day <- nb2listw(nbs1A_day, style="W", zero.policy=TRUE)
moran_day[1,1] <- "s1A_final_day"
moran_day[1,2] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[1,3] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [2]
nbs1B_day <- poly2nb(s1B_day, queen=TRUE)
lws1B_day <- nb2listw(nbs1B_day, style="W", zero.policy=TRUE)
moran_day[2,1] <- "s1B_final_day"
moran_day[2,2] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[2,3] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2A_day <- poly2nb(s2A_day, queen=TRUE)
lws2A_day <- nb2listw(nbs2A_day, style="W", zero.policy=TRUE)
moran_day[3,1] <- "s2A_final_day"
moran_day[3,2] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[3,3] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [2]
nbs2B_day <- poly2nb(s2B_day, queen=TRUE)
lws2B_day <- nb2listw(nbs2B_day, style="W", zero.policy=TRUE)
moran_day[4,1] <- "s2B_final_day"
moran_day[4,2] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[4,3] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [2]
nbs2C_day <- poly2nb(s2C_day, queen=TRUE)
lws2C_day <- nb2listw(nbs2C_day, style="W", zero.policy=TRUE)
moran_day[5,1] <- "s2C_final_day"
moran_day[5,2] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[5,3] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [2]
moran_round_day <- moran_day %>% mutate(across(2:3, round, 3))Combine Duncan and Moran files
Duncan_moran_dayCell <- left_join(duncan_day, moran_round_day, by = "type")
Duncan_moran_dayCell ## type Duncan_dayCell IMoran_dayCell pvalueMoran_dayCell
## 1 s1A_final_day 0.067 0.010 0.073
## 2 s1B_final_day 0.039 -0.014 0.983
## 3 s2A_final_day 0.078 0.030 0.000
## 4 s2B_final_day 0.038 -0.009 0.900
## 5 s2C_final_day 0.088 0.049 0.000In evening cells
s1A_evening <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s1B_evening <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2A_evening <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2B_evening <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
layer="Final_2")
s2C_evening <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
layer="Final_2")Duncan Dissimilarity Index
duncan_evening <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_evening) <- c("type", "Duncan_eveningCell" )
dfs1A_evening<- s1A_evening %>% st_drop_geometry()
duncan_evening[1,1] <- "s1A_final_evening"
duncan_evening[1,2]<- ISDuncan(dfs1A_evening[,3:4]) [1]
dfs1B_evening<- s1B_evening %>% st_drop_geometry()
duncan_evening[2,1] <- "s1B_final_evening"
duncan_evening[2,2]<- ISDuncan(dfs1B_evening[,3:4]) [1]
dfs2A_evening<- s2A_evening %>% st_drop_geometry()
duncan_evening[3,1] <- "s2A_final_evening"
duncan_evening[3,2]<- ISDuncan(dfs2A_evening[,3:4]) [1]
dfs2B_evening<- s2B_evening %>% st_drop_geometry()
duncan_evening[4,1] <- "s2B_final_evening"
duncan_evening[4,2]<- ISDuncan(dfs2B_evening[,3:4]) [1]
dfs2C_evening<- s2C_evening %>% st_drop_geometry()
duncan_evening[5,1] <- "s2C_final_evening"
duncan_evening[5,2]<- ISDuncan(dfs2C_evening[,3:4]) [1]
duncan_evening <- duncan_evening %>% mutate(across(starts_with("Duncan"), round, 3))Moran’s index of spatial autocorrelation
moran_evening <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_evening) <- c("type", "IMoran_eveningCell" , "pvalueMoran_eveningCell")
nbs1A_evening <- poly2nb(s1A_evening, queen=TRUE)
lws1A_evening <- nb2listw(nbs1A_evening, style="W", zero.policy=TRUE)
moran_evening[1,1] <- "s1A_final_evening"
moran_evening[1,2] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[1,3] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs1B_evening <- poly2nb(s1B_evening, queen=TRUE)
lws1B_evening <- nb2listw(nbs1B_evening, style="W", zero.policy=TRUE)
moran_evening[2,1] <- "s1B_final_evening"
moran_evening[2,2] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[2,3] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2A_evening <- poly2nb(s2A_evening, queen=TRUE)
lws2A_evening <- nb2listw(nbs2A_evening, style="W", zero.policy=TRUE)
moran_evening[3,1] <- "s2A_final_evening"
moran_evening[3,2] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[3,3] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2B_evening <- poly2nb(s2B_evening, queen=TRUE)
lws2B_evening <- nb2listw(nbs2B_evening, style="W", zero.policy=TRUE)
moran_evening[4,1] <- "s2B_final_evening"
moran_evening[4,2] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[4,3] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [2]
nbs2C_evening <- poly2nb(s2C_evening, queen=TRUE)
lws2C_evening <- nb2listw(nbs2C_evening, style="W", zero.policy=TRUE)
moran_evening[5,1] <- "s2C_final_evening"
moran_evening[5,2] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[5,3] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [2]
moran_round_evening <- moran_evening %>% mutate(across(2:3, round, 3))Combine Duncan and Moran files
Duncan_moran_eveningCell <- left_join(duncan_evening, moran_round_evening, by = "type")
Duncan_moran_eveningCell## type Duncan_eveningCell IMoran_eveningCell
## 1 s1A_final_evening 0.084 0.010
## 2 s1B_final_evening 0.039 -0.005
## 3 s2A_final_evening 0.098 0.033
## 4 s2B_final_evening 0.038 0.000
## 5 s2C_final_evening 0.089 0.061
## pvalueMoran_eveningCell
## 1 0.057
## 2 0.769
## 3 0.000
## 4 0.472
## 5 0.000
Social inequalities in healthy behaviour for the 5 scenarios (Figure 5)
Importing simulation results
Distribution of simulated values of EducIneqIndex (EII) for the five space-time scenarios (10 000 replications per scenario) (figure 5)