couverture_vaccinale.Rmd

---
title: "Couverture vaccinale"
author: "FD"
output: 
  html_document: 
      code_folding: hide
      toc: TRUE
      toc_float: TRUE
      self_contained: no
editor_options:
  chunk_output_type: console
---

<!--

NOTES: 

- I am using Rmd just to organize my code, not to generate an output. 

- Need to manually check age classes from time to time, as they might change in Ameli data

- IMPORTANT : NEED TO PUT THE DATA IN THE SHINY FOLDER

-->

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls()) # I don't care, I am doing it anyway
```

```{r, include = FALSE}
dlData <- FALSE
# Whether to download the data again
```


# Load data

## Vacination par EPCI et commune

Source EPCI <https://datavaccin-covid.ameli.fr/explore/dataset/donnees-de-vaccination-par-epci/information/>  
Source Communes <https://datavaccin-covid.ameli.fr/explore/dataset/donnees-de-vaccination-par-commune/information/>

```{r loadDataAmeli1}

# Download data
if(dlData){
  source("0_loadAmeliData.R")
}

# Load data
dat.EPCI <- read.csv("data/vaccEPCI.csv", sep = ";")
dat.Com <- read.csv("data/vaccCom.csv", sep = ";")
```

## Vaccination par département de résidence

Source <https://datavaccin-covid.ameli.fr/explore/dataset/donnees-vaccination-par-tranche-dage-type-de-vaccin-et-departement/information/?sort=-date_reference>

```{r loadData}
URL.AmeliDep <- "https://datavaccin-covid.ameli.fr/explore/dataset/donnees-vaccination-par-tranche-dage-type-de-vaccin-et-departement/download/?format=csv&timezone=Europe/Berlin&lang=fr&use_labels_for_header=true&csv_separator=%3B"
## does not work

#dataFile <- paste0("data/donnees-vaccination-par-tranche-dage-type-de-vaccin-et-departement_", Sys.Date(),".csv") # name file with today's date
dataFile <- "data/vaccDep.csv"
if(dlData){
  system(paste0("wget -O ", dataFile, " ", URL.AmeliDep))
}

dat.residence <- read.csv(dataFile, sep = ";", stringsAsFactors = FALSE, dec = ",")
head(dat.residence)

# Remove details about vaccines
dat.residence <- dat.residence[dat.residence$type_vaccin == "Tout vaccin", ]
```

## Vaccination par département d'injection

Source <https://www.data.gouv.fr/en/datasets/donnees-relatives-aux-personnes-vaccinees-contre-la-covid-19-1/>  
`vacsi-dep`

```{r}
URL <- "https://www.data.gouv.fr/en/datasets/r/4f39ec91-80d7-4602-befb-4b522804c0af"
URL.SPFadep <- "https://www.data.gouv.fr/en/datasets/r/83cbbdb9-23cb-455e-8231-69fc25d58111"

dataFile <- paste0("data/vacciDeps_injection.csv") # name file with today's date
dataFile2 <- paste0("data/vacciDeps_injection-a.csv") # name file with today's date
if(dlData){
  download.file(URL, dataFile) # download file from repo
  download.file(URL.SPFadep, dataFile2) # download file from repo
}
dat.injection.noage <- read.csv(dataFile, sep = ";", stringsAsFactors = FALSE, dec = ",")
dat.injection <- read.csv(dataFile2, sep = ";", stringsAsFactors = FALSE, dec = ",")
# head(dat.injection)
head(dat.injection)

```

## Vaccination France totale, SPF

```{r}
# SPF, France

URL <- "https://www.data.gouv.fr/fr/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd"

dataFile <- "data/vacciFrance.csv"
if(dlData){
  download.file(URL, dataFile)
}
datFra <- read.csv(dataFile, sep = ";", stringsAsFactors = FALSE, dec = ",")
head(datFra)

# Population sizes are given in another file
# SPF, final, which contains population values
tmp <- read.csv("data/vaccSPF_tot.csv", sep = ";")

# Add this info to the SPF dataset
datFra <- merge(datFra, tmp[, c("clage_vacsi", "pop")], by = "clage_vacsi", all = TRUE)
```



## Clean data

```{r}
# Change "NS" into "NA"
dat.EPCI[dat.EPCI$population_carto == "NS", "population_carto"] <- NA
dat.Com[dat.Com$population_carto == "NS", "population_carto"] <- NA


# Get final date in the data
lastDate <- max(dat.EPCI$date)
lastDate
```

Add population size in SPF

```{r}
tmp <- dat.injection[dat.injection$jour == lastDate, ]
vpop2020 <- tmp$n_cum_dose1 / (as.numeric(tmp$couv_dose1) / 100)
vpop2020[vpop2020 == Inf] <- NA

tmp2 <- cbind(tmp, pop2020 = vpop2020)
tmp2 <- tmp2[, c("clage_vacsi", "dep", "pop2020")]

# Note: pop sizes are actually given in another SPF dataset; 
# they are added later in this code
```

## Age classes 

```{r}
# Compare age classes
# For the moment, we just look at all ages, but note that we could look at the detail as well one day when we find the time
sort(unique(dat.EPCI$classe_age))
sort(unique(dat.Com$classe_age))

sort(unique(dat.residence$classe_age))

sort(unique(dat.injection$clage_vacsi))
```

```{r}
# Define common age classes
agCom <- list("00-39" = c("00-19", "20-39"), 
              "40-64" = c("40-54", "55-64"),
              "65-74" = "65-74",
              "75+" = "75 et +",
              "Tous âges" = "TOUT_AGE")

agEPCI <- agCom 

agDepRes <- list("00-39" = c("12-17", "18-24", "25-39"), 
                 "40-64" = c("40-54", "55-64"),
                 "65-74" = "65-74",
                 "75+" = "75 et +",
                 "Tous âges" = "TOUT_AGE")

agDepInj <- list("00-39" = c(4, 9, 11, 17, 24, 29, 39), 
                 "40-64" = c(49, 59, 64),
                 "65-74" = c(69, 74),
                 "75+" = c(79, 80),
                 "Tous âges" = 0)

ageClasses <- names(agEPCI)
```


## Locations

```{r}
tmp <- dat.EPCI[dat.EPCI$epci == 999999999 & dat.EPCI$classe_age == "TOUT_AGE" & dat.EPCI$libelle_classe_age != "",]
tail(tmp[order(tmp$semaine_injection),])
```


Sanity Check: Comparisons of population sizes in Ameli data

```{r}
# SIREN codes
epci_grandParis <- 200054781
epci_Lyon <- 200046977
epci_Marseille <- 200054807

cols <- c("population_carto", "effectif_cumu_1_inj", "effectif_cumu_termine")

# Compute coverage
cvr <- function(p){
  print(p)
  as.numeric(p$effectif_cumu_1_inj) / as.numeric(p$population_carto)
}

sort(unique(dat.Com$commune_residence))

## PARIS
ParisEPCI <- dat.EPCI[dat.EPCI$epci == epci_grandParis & dat.EPCI$classe_age == "TOUT_AGE", ]

cvr(ParisEPCI[ParisEPCI$date == lastDate, cols])

ParisCom <- dat.Com[(dat.Com$commune_residence >= "75000" & dat.Com$commune_residence != "83120" & dat.Com$commune_residence != "84089") & dat.Com$classe_age == "TOUT_AGE", ]

cvr(as.data.frame(t(apply(ParisCom[ParisCom$date == lastDate, cols], 2, sum))))

## MARSEILLE
MarseilleEPCI <- dat.EPCI[dat.EPCI$epci == epci_Marseille & dat.EPCI$classe_age == "TOUT_AGE", ]

cvr(MarseilleEPCI[MarseilleEPCI$date == lastDate, cols])

MarseilleCom <- dat.Com[(dat.Com$commune_residence <= "14000" | dat.Com$commune_residence >= "83000" & dat.Com$commune_residence <= "84999") & dat.Com$classe_age == "TOUT_AGE", ]

cvr(as.data.frame(t(apply(MarseilleCom[MarseilleCom$date == lastDate, cols], 2, sum))))

## LYON
LyonEPCI <- dat.EPCI[dat.EPCI$epci == epci_Lyon & dat.EPCI$classe_age == "TOUT_AGE", ]

cvr(LyonEPCI[LyonEPCI$date == lastDate, cols])

LyonCom <- dat.Com[dat.Com$commune_residence <= "70000" & dat.Com$commune_residence >= "69000" & dat.Com$classe_age == "TOUT_AGE", ]

cvr(as.data.frame(t(apply(LyonCom[LyonCom$date == lastDate, cols], 2, sum))))

```

Add departments and regions

Departement 

Source <https://www.insee.fr/fr/information/2510634>
```{r}
composition <- read.csv("data/EPCI_composition-communale.csv", encoding = "UFT-8")
```

```{r}
departements <- read.csv("data/departement2020.csv")

dic.depname <- departements$libelle
names(dic.depname) <- departements$dep
```


```{r}
# Information about departement in which the different EPCI are
# Some are across multiple departements: keep the information by collating them with "_"
agg_nbdep <- aggregate(composition$DEP, by = list(composition$EPCI), FUN = function(i) paste(sort(unique(i)), collapse = "_"))
table(agg_nbdep$x, useNA = "ifany")

# Dictionnary of departement(s) associated to EPCI
dic.dep <- agg_nbdep$x
names(dic.dep) <- agg_nbdep$Group.1

# Add the dep information to our data
dat.EPCI$dep <- dic.dep[as.character(dat.EPCI$epci)]

# Add the libelle information
dat.EPCI$dep_libelle <- dic.depname[as.character(dat.EPCI$dep)]

# Add DROM information
dat.EPCI$DROM <- ifelse(nchar(dat.EPCI$dep) == 3, TRUE, FALSE)

# EPCI for which we do not have a dep
unique(dat.EPCI[is.na(dat.EPCI$dep), "epci"])
# This is the code when home is not known

# La répartition par EPCI n’est pas possible pour certains patients dont le département de résidence ou la commune est inconnu(e). Ces patients sont rattachés à l’EPCI « 999999999 » dans le jeu de données.

```

Add region information 

Source region codes <https://www.data.gouv.fr/en/datasets/regions-de-france/>

```{r, results = 'hide'}
# Load region codes 
reg <- read.csv("data/regions-france.csv")
reg

# As dictionnary
reg.dic <- c(reg$nom_region)
names(reg.dic) <- reg$code_region # !! Name has to be in quotes

sort(unique(dat.EPCI$reg_code))

# Regions in the data
sort(unique(dat.EPCI$reg_code))

# missing region
unique(dat.EPCI[dat.EPCI$reg_code == "", "epci"])
```

Add region information to the other datasets
```{r}
dic.depreg <- departements$reg
names(dic.depreg) <- departements$dep

dat.injection$reg <- dic.depreg[dat.injection$dep]

dat.residence[1,]

unique(dat.residence$departement_residence)
unique(dat.residence$region_residence)
```

## Add INSEE 2021 Demographic data

```{r loadDatasets}
# Load demographic data
# Sources: Population par classe d'age
# https://www.insee.fr/fr/statistiques/1893198
#
# Population par année d'âge (par milliers)
# https://www.insee.fr/fr/statistiques/5347620#tableau-figure6_radio1
#
demoAmeli <- read.csv("data/demographic/insee_2021_Ameli.csv")
demoAmeli
demoSPF <- read.csv("data/demographic/insee_2021_SPF.csv")
demoSPF

demoCommon <- read.csv("data/demographic/insee_2021_commonAgeClass.csv")

# Add new population data to SPF
dat.injection <- merge(dat.injection, demoSPF, "clage_vacsi")

# Compute new rates
dat.injection$couv_tot_dose1_2021 <- 100 * dat.injection$n_cum_dose1 / dat.injection$pop2021
dat.injection$couv_tot_complet_2021 <- 100 * dat.injection$n_cum_complet / dat.injection$pop2021

# Add population data to the other Ameli dataset
names(demoAmeli) <- c("classe_age", "pop2021")
#dat.residence2 <- merge(dat.residence, demoAmeli, by = "classe_age")

#unique(dat.residence$classe_age)
#unique(dat.residence2$classe_age)
```


## Population sizes SPF

```{r}
# SPF, final, which contains population values
#tmp <- read.csv("data/vaccSPF_tot.csv", sep = ";")
# 

tmp <- read.csv("data/vaccSPF-dep-tot.csv", sep = ";")
tmp <- tmp[, c("dep", "clage_vacsi", "pop")]
names(tmp) <- c("dep", "clage_vacsi", "pop2020")

# Add this info to the SPF dataset
dat.injection <- merge(dat.injection, tmp, by = c("clage_vacsi", "dep"), all = TRUE)

head(dat.injection)
```

# Useful stuff 

```{r}
is.empty <- function(x){
  is.na(x) | (x == "") | (x == "NS")
}
```

```{r}
trp <- 0.8
colEPCI <- adjustcolor(col = "#1f78b4", alpha.f = trp)
colDep <- adjustcolor(col = "#F1B200", alpha.f = trp)
colSPF <- adjustcolor(col = "#FF581F", alpha.f = trp)
```



# Comparisons, France entiere

```{r}
#thedate <- lastDate
#agcl <- "00-39"
#inj <- "1"
#agcl <- "Tous âges"

getCouv <- function(thedate, agcl, inj){
  
  if(agcl == "Tous âges"){agg <- "TOUT_AGE"}else{agg <- agcl}
  # Population sizes, INSEE 2021
  insee2021 <- demoCommon[demoCommon$classe_age == agg, "population_carto2021"]

  # Define column names depending on the chose type of injection
  # (Ameli and SPF names differ)
  
  if(inj == "1"){
    thecol <- "effectif_cumu_1_inj"
    thecolw <- "effectif_1_inj"
    thecolSPF <- "n_cum_dose1"
    thecolSPFw <- "n_dose1"
  }
  if(inj == "complet"){
    thecol <- "effectif_cumu_termine"
    thecolw <- "effectif_termine"
    thecolSPF <- "n_cum_complet"
    thecolSPFw <- "n_complet"
  }
  
  ## EPCI data
  tmpEPCI <- dat.EPCI[is.element(dat.EPCI$classe_age, agEPCI[[agcl]]) & dat.EPCI$date == thedate, ]
  
  linesOK <- which(!is.empty(tmpEPCI[, thecol]) & !is.empty(tmpEPCI[, "population_carto"]))
  txEPCI <- sum(as.numeric(tmpEPCI[linesOK, thecol])) / sum(as.numeric(tmpEPCI[linesOK, "population_carto"]))
  V1EPCI <- sum(as.numeric(tmpEPCI[, thecol]), na.rm = TRUE)
  N1EPCI <- sum(as.numeric(tmpEPCI[, "population_carto"]), na.rm = TRUE)
  V2EPCI <- sum(as.numeric(tmpEPCI[linesOK, thecol]))
  N2EPCI <- sum(as.numeric(tmpEPCI[linesOK, "population_carto"]))
  
  linesOKw <- which(!is.empty(tmpEPCI[, thecolw]) & !is.empty(tmpEPCI[, "population_carto"]))
  txEPCIw <- sum(as.numeric(tmpEPCI[linesOKw, thecolw])) / sum(as.numeric(tmpEPCI[linesOKw, "population_carto"]))
  V1wEPCI <- sum(as.numeric(tmpEPCI[, thecolw]), na.rm = TRUE)
  V2wEPCI <- sum(as.numeric(tmpEPCI[linesOKw, thecolw]))

  txEPCI
  
  ##----------------------------
  
  ## Departement data, Ameli
  
  tmpResidence <- dat.residence[is.element(dat.residence$classe_age, agDepRes[[agcl]]) & dat.residence$date == thedate & dat.residence$libelle_region == "FRANCE", ]
  
  ## 
  linesOK <- which(!is.empty(tmpResidence[, thecol]) & !is.empty(tmpResidence[, "population_insee"]))
  
  txRes <- sum(as.numeric(tmpResidence[linesOK, thecol])) / sum(as.numeric(tmpResidence[linesOK, "population_insee"]))
  
  V1Res <- sum(as.numeric(tmpResidence[, thecol]), na.rm = TRUE)
  N1Res <- sum(as.numeric(tmpResidence[, "population_insee"]), na.rm = TRUE)
  V2Res <- sum(as.numeric(tmpResidence[linesOK, thecol]))
  N2Res <- sum(as.numeric(tmpResidence[linesOK, "population_insee"]))
  N3 <- insee2021
  
  txRes2 <- sum(as.numeric(tmpResidence[linesOK, thecol])) / insee2021
  
  linesOKw <- which(!is.empty(tmpResidence[, thecolw]) & !is.empty(tmpResidence$population_insee))
  txResw <- sum(as.numeric(tmpResidence[linesOKw, thecolw])) / sum(as.numeric(tmpResidence[linesOKw, "population_insee"]))
  txRes2w <- sum(as.numeric(tmpResidence[linesOKw, thecolw])) / insee2021
  
  V1wRes <- sum(as.numeric(tmpResidence[, thecolw]), na.rm = TRUE)
  V2wRes <- sum(as.numeric(tmpResidence[linesOKw, thecolw]))
  
  #*******************************
  # Departement data, SPF
  
  tmpSPF <- datFra[is.element(datFra$clage_vacsi, agDepInj[[agcl]]) & datFra$jour == thedate, ]
  
  tmpSPFw <- datFra[is.element(datFra$clage_vacsi, agDepInj[[agcl]]) & as.Date(datFra$jour) <= as.Date(thedate) & as.Date(datFra$jour) >= as.Date(thedate) - 6, ]

  
  linesOK <- which(!is.empty(tmpSPF[, thecolSPF]) & !is.empty(tmpSPF$pop))
  
  txSPF <- sum(as.numeric(tmpSPF[linesOK, thecolSPF])) / sum(as.numeric(tmpSPF[linesOK, "pop"]))
  
  txSPF2 <- sum(as.numeric(tmpSPF[linesOK, thecolSPF])) / insee2021
  
  V1SPF <- sum(as.numeric(tmpSPF[, thecolSPF]), na.rm = TRUE) 
  V2SPF <- sum(as.numeric(tmpSPF[linesOK, thecolSPF])) 
  N1SPF <- sum(as.numeric(tmpSPF[, "pop"]), na.rm = TRUE)
  N2SPF <- sum(as.numeric(tmpSPF[linesOK, "pop"]))
  
  linesOKw <- which(!is.empty(tmpSPFw[, thecolSPFw]) & !is.empty(tmpSPFw$pop))
  
  txSPFw <- sum(as.numeric(tmpSPFw[linesOKw, thecolSPFw])) / (sum(as.numeric(tmpSPFw[linesOKw, "pop"])) /7)
  txSPF2w <- sum(as.numeric(tmpSPFw[linesOKw, thecolSPFw])) / insee2021
  
  V1wSPF <- sum(as.numeric(tmpSPFw[, thecolSPFw]), na.rm = TRUE) 
  V2wSPF <- sum(as.numeric(tmpSPFw[linesOKw, thecolSPFw]))   

  out <- c(txEPCI, txRes, txRes2, txSPF, txSPF2, 
           txEPCIw, txResw, txRes2w, txSPFw, txSPF2w, 
           V1EPCI, V2EPCI, N1EPCI, N2EPCI, V1wEPCI, V2wEPCI, 
           V1Res, V2Res, N1Res, N2Res, V1wRes, V2wRes, 
           V1SPF, V2SPF, N1SPF, N2SPF, V1wSPF, V2wSPF, 
           N3
           )
  names(out) <- c("txEPCI", "txRes", "txRes2", "txSPF", "txSPF2", 
                 "txEPCIw", "txResw", "txRes2w", "txSPFw", "txSPF2w", 
                 "V1EPCI", "V2EPCI", "N1EPCI", "N2EPCI", "V1wEPCI", "V2wEPCI", 
                 "V1Res", "V2Res", "N1Res", "N2Res", "V1wRes", "V2wRes", 
                 "V1SPF", "V2SPF", "N1SPF", "N2SPF", "V1wSPF", "V2wSPF", 
                 "N3")
  out
}


getCouv(lastDate, "75+", "1")

getCouv(lastDate, "75+", "complet")

vv <- getCouv(lastDate, "Tous âges", "1")
vv

vv['N1Res'] - vv['N1SPF']

vv['V1EPCI'] / vv['N1EPCI']
vv['V2EPCI'] / vv['N2EPCI']

vv['V2SPF'] / vv['N2SPF']
vv['V2Res'] / vv['N2Res']

#

names(vv) <- c("txEPCI", "txRes", "txRes2", "txSPF", "txSPF2", 
                 "txEPCIw", "txResw", "txRes2w", "txSPFw", "txSPF2w", 
                 "V1EPCI", "V2EPCI", "N1EPCI", "N2EPCI", "V1wEPCI", "V2wEPCI", 
                 "V1Res", "V2Res", "N1Res", "N2Res", "V1wRes", "V2wRes", 
                 "V1SPF", "V2SPF", "N1SPF", "N2SPF", "V1wSPF", "V2wSPF", 
                 "N3")
```


```{r}
parms <- expand.grid(thedate = sort(unique(dat.EPCI$date)), agcl = names(agCom), inj = c("1", "complet"))

nrow(parms)

res <- lapply(seq_len(nrow(parms)), function(i) do.call(getCouv, parms[i, ]))

mres <- as.data.frame(matrix(unlist(res), byrow = TRUE, nrow = nrow(parms)))
names(mres) <- c("txEPCI", "txRes", "txRes2", "txSPF", "txSPF2", 
                 "txEPCIw", "txResw", "txRes2w", "txSPFw", "txSPF2w", 
                 "V1EPCI", "V2EPCI", "N1EPCI", "N2EPCI", "V1wEPCI", "V2wEPCI", 
                 "V1Res", "V2Res", "N1Res", "N2Res", "V1wRes", "V2wRes", 
                 "V1SPF", "V2SPF", "N1SPF", "N2SPF", "V1wSPF", "V2wSPF", 
                 "N3")

mres <- cbind(parms, mres)

mres$date <- as.Date(mres$thedate)
```

### Plot tests

```{r, eval = FALSE}
#thedate <- lastDate
#agcl <- "00-39"
inj <- "1"
agcl <- "Tous âges"

mres1 <- mres[mres$inj == inj, ]
mm <- mres1[mres1$agcl == agcl, ]
names(mm)

nend <- nrow(mm)

plot(mm$date, mm$N1EPCI, ylim = c(60*10^6, 70*10^6))
points(mm$date, mm$V1EPCI)

points(mm$date, mm$N1Res, col = colDep)
points(mm$date, mm$N1SPF, col = colSPF)
points(mm$date, mm$N1EPCI, col = colEPCI)

points(mm$date, mm$N2Res, col  = colDep, pch = 2)
points(mm$date, mm$N2SPF, col = colSPF, pch = 2)
points(mm$date, mm$N2EPCI, col = colEPCI, pch = 2)


mm[nend, "N2EPCI"]


plot(mm$date, mm$N1EPCI, ylim = c(0*10^6, 70*10^6))
points(mm$date, mm$V1Res, col = colDep)
points(mm$date, mm$V1SPF, col = colSPF)
points(mm$date, mm$V1EPCI, col = colEPCI)

points(mm$date, mm$V2Res, col  = colDep, pch = 2)
points(mm$date, mm$V2SPF, col = colSPF, pch = 2)
points(mm$date, mm$V2EPCI, col = colEPCI, pch = 2)


plot(mm$date, mm$N1EPCI, ylim = c(0*10^6, 7*10^6))
points(mm$date, mm$V1wRes, col = colDep)
points(mm$date, mm$V1wSPF, col = colSPF)
points(mm$date, mm$V1wEPCI, col = colEPCI)

points(mm$date, mm$V2wRes, col  = colDep, pch = 2)
points(mm$date, mm$V2wSPF, col = colSPF, pch = 2)
points(mm$date, mm$V2wEPCI, col = colEPCI, pch = 2)


mm[nend, ]



```


### Plot France entiere

```{r}
for(inj in c("1", "complet")){
  for(agcl in names(agCom)){
    mres1 <- mres[mres$inj == inj, ]
    mm <- mres1[mres1$agcl == agcl, ]
    
    thetype <- "o"
    thecex <- 0.3
    
    tit1 <- ifelse(inj == "1", "Premières injections", "Vaccinations complètes")
    
    agg <- agcl
    if(agg == "Tous âges") agg <- "tousages"
    if(agg == "75+") agg <- "75-etplus"
    
    fname <- paste0("pics/vaccComparaison_cum_", agg, inj, ".pdf")
    pdf(file = fname, width = 5, height = 4.5)
    par(las = 1, mgp = c(2, 0.5, 0), tck = -0.015, xpd = TRUE, 
        mar = c(3.25, 3, 3, 2.5))
    plot(mm$date, mm$txEPCI, col = colEPCI, ylim = c(0, 1), 
          xaxs = "i", yaxs = "i", 
         xlab = "date", ylab = "proportion", 
         main = paste(tit1, "cumulées, ", agcl), 
         type = thetype, cex = thecex)
    if(!is.element(agcl, c("00-39", "Tous âges"))){
      points(mm$date, mm$txRes, col = colDep, 
         type = thetype, cex = thecex)
    }
    points(mm$date, mm$txSPF, col = colSPF, 
         type = thetype, cex = thecex)
    points(mm$date, mm$txRes2, col = colDep, pch = 2, 
         type = thetype, cex = thecex, lty = 2)
    points(mm$date, mm$txSPF2, col = colSPF, pch = 2, 
         type = thetype, cex = thecex, lty = 2)
    legend("topleft", col = c(colEPCI, colDep, colSPF, colDep, colSPF), lty = c(1, 1, 1, 2, 2), legend = c("Ameli, assurés sociaux", "Ameli, INSEE 2020", "SPF, INSEE 2020", "Ameli, INSEE 2021", "SPF, INSEE 2021"), lwd = 2, bty = "n", cex = 0.5)
    axis(4)
    dev.off()
    
    
    fname2 <- paste0("pics/vaccComparaison_sem_", agg, inj, ".pdf")
    pdf(file = fname2, width = 5, height = 4.5)
    par(las = 1, mgp = c(2, 0.5, 0), tck = -0.015, xpd = TRUE, 
        mar = c(3.25, 3, 3, 2.5))
    plot(mm$date, mm$txEPCIw, col = colEPCI, ylim = c(0, 0.1), type = thetype, 
         xaxs = "i", yaxs = "i", 
         xlab = "date", ylab = "proportion", 
         cex = thecex, 
         main = paste(tit1, "de la semaine, ", agcl))
    if(!is.element(agcl, c("00-39", "Tous âges"))){
      points(mm$date, mm$txResw, col = colDep, type = thetype, cex = thecex)
    }
    points(mm$date, mm$txSPFw, col = colSPF, type = thetype, cex = thecex)
    points(mm$date, mm$txRes2w, col = colDep, pch = 2, type = thetype, cex = thecex, lty = 2)
    points(mm$date, mm$txSPF2w, col = colSPF, pch = 2, type = thetype, cex = thecex, lty = 2)
    
    legend("topleft", col = c(colEPCI, colDep, colSPF, colDep, colSPF), lty = c(1, 1, 1, 2, 2), legend = c("Ameli, assurés sociaux", "Ameli, INSEE 2020", "SPF, INSEE 2020", "Ameli, INSEE 2021", "SPF, INSEE 2021"), lwd = 2, bty = "n", cex = 0.5)
    axis(4)
    dev.off()

    
    system(paste("open", fname))
    system(paste("open", fname2))

  }
}




```

# Comparisons, departements

```{r}
#sort(unique(dat.Com$commune_residence))

#dep <- "75"
#agcl <- "Tous âges"
#thedate <- lastDate

#table(dat.EPCI$dep)
# 83, 84 in Marseille

getCouvDep <- function(thedate, agcl, inj, dep){
  
  # Define column names depending on the chose type of injection
  # (Ameli and SPF names differ)
  
  if(inj == "1"){
    thecol <- "effectif_cumu_1_inj"
    thecolw <- "effectif_1_inj"
    thecolSPF <- "n_cum_dose1"
    thecolSPFw <- "n_dose1"
  }
  if(inj == "complet"){
    thecol <- "effectif_cumu_termine"
    thecolw <- "effectif_termine"
    thecolSPF <- "n_cum_complet"
    thecolSPFw <- "n_complet"
  }

  #if(is.element(dep, c("13", "75", "91", "92", "93", "94", "95", "69", "83", "84"))){
    # If departement for which we have communes data
    # Actually no need to do it, tmpCom will be empty if not included
  
    ## EPCI data
  # Need to treat separately the case of DROM, 
  # because otherwise "72" can be found in "972"  
    if(nchar(dep) == 2){
      tmpEPCI <- dat.EPCI[which(is.element(dat.EPCI$classe_age, agEPCI[[agcl]]) & dat.EPCI$date == thedate & grepl(dep, dat.EPCI$dep, fixed = TRUE) & nchar(dat.EPCI$dep) != 3), ]
    }else{
      # DROM, dep of lenth 3
      tmpEPCI <- dat.EPCI[which(is.element(dat.EPCI$classe_age, agEPCI[[agcl]]) & dat.EPCI$date == thedate & grepl(dep, dat.EPCI$dep, fixed = TRUE)), ]
    }
    # Remove P, L, M to avoid counting them twice (included in Dep)
    tmpEPCI <- tmpEPCI[!is.element(tmpEPCI$epci, c(epci_grandParis, epci_Lyon, epci_Marseille)), ]
    
      ## Com data for this departement
  tmpCom <- dat.Com[is.element(dat.Com$classe_age, agCom[[agcl]]) & dat.Com$date == thedate & substr(dat.Com$commune_residence, 1, 2) == dep, ]
  
  columns <- c("date", "effectif_1_inj",
               "effectif_termine", "effectif_cumu_1_inj", 
               "effectif_cumu_termine", "population_carto")
  tmp <- rbind(tmpEPCI[, columns], tmpCom[, columns])
  tmp
  nrow(tmp)
  
  linesOK <- which(!is.empty(tmp[, thecol]) & !is.empty(tmp$population_carto))
  
  txAmeli <- sum(as.numeric(tmp[linesOK, thecol])) / sum(as.numeric(tmp[linesOK, "population_carto"]))
  V1EPCI <- sum(as.numeric(tmp[, thecol]), na.rm = TRUE)
  V2EPCI <- sum(as.numeric(tmp[linesOK, thecol]))
  N1EPCI <- sum(as.numeric(tmp[, "population_carto"]), na.rm = TRUE)
  N2EPCI <- sum(as.numeric(tmp[linesOK, "population_carto"]))
  
  linesOKw <- which(!is.empty(tmp[, thecolw]) & !is.empty(tmp$population_carto))
  txAmeliw <- sum(as.numeric(tmp[linesOKw, thecolw])) / sum(as.numeric(tmp[linesOKw, "population_carto"]))
  
  V1wEPCI <- sum(as.numeric(tmp[, thecolw]), na.rm = TRUE)
  V2wEPCI <- sum(as.numeric(tmp[linesOKw, thecolw]))
  
    

  ##----------------------------
  
  ## Departement data, Ameli
  
  tmpResidence <- dat.residence[which(is.element(dat.residence$classe_age, agDepRes[[agcl]]) & dat.residence$date == thedate & dat.residence$departement_residence == dep), ]
  
  linesOK <- which(!is.empty(tmpResidence[, thecol]) & !is.empty(tmpResidence$population_insee))
  
  txRes <- sum(as.numeric(tmpResidence[linesOK, thecol])) / sum(as.numeric(tmpResidence[linesOK, "population_insee"]))
  
  linesOKw <- which(!is.empty(tmpResidence[, thecolw]) & !is.empty(tmpResidence$population_insee))
  txResw <- sum(as.numeric(tmpResidence[linesOKw, thecolw])) / sum(as.numeric(tmpResidence[linesOKw, "population_insee"]))
  
  V1Res <- sum(as.numeric(tmpResidence[, thecol]), na.rm = TRUE)
  V2Res <- sum(as.numeric(tmpResidence[linesOK, thecol]))
  N1Res <- sum(as.numeric(tmpResidence[, "population_insee"]), na.rm = TRUE)
  N2Res <- sum(as.numeric(tmpResidence[linesOK, "population_insee"]))
  
  V1wRes <- sum(as.numeric(tmpResidence[, thecolw]), na.rm = TRUE)
  V2wRes <- sum(as.numeric(tmpResidence[linesOKw, thecolw]))
  
  
  #--------------------------------
  ## Departement data, SPF
  
  tmpSPF <- dat.injection[which(is.element(dat.injection$clage_vacsi, agDepInj[[agcl]]) & dat.injection$jour == thedate & dat.injection$dep == dep), ]
    
  tmpSPFw <- dat.injection[which(dat.injection$dep == dep & is.element(dat.injection$clage_vacsi, agDepInj[[agcl]])), ]
  tmpSPFw <- tmpSPFw[as.Date(tmpSPFw$jour) <= as.Date(thedate) & (as.Date(tmpSPFw$jour) >= as.Date(thedate) - 6), ]
  
  
  linesOK <- which(!is.empty(tmpSPF[, thecolSPF]) & !is.empty(tmpSPF$pop2020))
  
  txSPF <- sum(as.numeric(tmpSPF[linesOK, thecolSPF])) / sum(as.numeric(tmpSPF[linesOK, "pop2020"]))
  
  linesOKw <- which(!is.empty(tmpSPFw[, thecolSPFw]) & !is.empty(tmpSPFw$pop2020))
  txSPFw <- sum(as.numeric(tmpSPFw[linesOKw, thecolSPFw])) / (sum(as.numeric(tmpSPFw[linesOKw, "pop2020"])) /7)
  
  V1SPF <- sum(as.numeric(tmpSPF[, thecolSPF]), na.rm = TRUE)
  V2SPF <- sum(as.numeric(tmpSPF[linesOK, thecolSPF]))
  N1SPF <- sum(as.numeric(tmpSPF[, "pop2020"]), na.rm = TRUE)
  N2SPF <- sum(as.numeric(tmpSPF[linesOK, "pop2020"]))
  
  V1wSPF <- sum(as.numeric(tmpSPFw[, thecolSPFw]), na.rm = TRUE) 
  V2wSPF <- sum(as.numeric(tmpSPFw[linesOKw, thecolSPFw])) 
  

  out <- c(txAmeli, txRes, txSPF, txAmeliw, txResw, txSPFw, 
           V1EPCI, V2EPCI, N1EPCI, N2EPCI, V1wEPCI, V2wEPCI, 
           V1Res, V2Res, N1Res, N2Res, V1wRes, V2wRes, 
           V1SPF, V2SPF, N1SPF, N2SPF, V1wSPF, V2wSPF
  )
  out
  
}

# Check
m <- getCouvDep(lastDate, "Tous âges", "1", "75")
m <- getCouvDep("2021-07-11", "Tous âges", "1", "75")
names(m) <- c("txEPCI", "txRes", "txSPF", 
                 "txEPCIw", "txResw", "txSPFw", 
                 "V1EPCI", "V2EPCI", "N1EPCI", "N2EPCI", "V1wEPCI", "V2wEPCI", 
                 "V1Res", "V2Res", "N1Res", "N2Res", "V1wRes", "V2wRes", 
                 "V1SPF", "V2SPF", "N1SPF", "N2SPF", "V1wSPF", "V2wSPF")
m
getCouvDep(lastDate, "Tous âges", "1", "23")
getCouvDep(lastDate, "Tous âges", "1", "972")
```

```{r}
# All departements
deps <- sort(unique(dat.residence$departement_residence))
deps <- deps[nchar(deps) < 4]
deps

# Initialize output
out <- data.frame(matrix(NA, nrow = 0, ncol = 29))
names(out) <- c("thedate", "agcl", "inj", "dep", 
                "txEPCI", "txRes", "txSPF", 
                "txEPCIw", "txResw", "txSPFw", 
                "V1EPCI", "V2EPCI", "N1EPCI", "N2EPCI", "V1wEPCI", "V2wEPCI", 
                "V1Res", "V2Res", "N1Res", "N2Res", "V1wRes", "V2wRes", 
                "V1SPF", "V2SPF", "N1SPF", "N2SPF", "V1wSPF", "V2wSPF",
                "date")

# Loop on departements
for(dep in deps){

  # Print counter
  print(paste(dep, "en cours"))
  
  # Parameters for this dep
  parms <- expand.grid(thedate = sort(unique(dat.EPCI$date)), agcl = names(agCom), inj = c("1", "complet"), dep = dep, stringsAsFactors = FALSE)
  
  # Compute couvertures
  res <- lapply(seq_len(nrow(parms)), function(i) do.call(getCouvDep, parms[i, ]))
  
  # Format output as data frame
  mres <- as.data.frame(matrix(unlist(res), byrow = TRUE, nrow = nrow(parms)))
  names(mres) <- c("txEPCI", "txRes", "txSPF", 
                   "txEPCIw", "txResw", "txSPFw")
  
  # Add parameters
  mres <- cbind(parms, mres)
  
  # Save date as date
  mres$date <- as.Date(mres$thedate)
  
  # Export
  write.csv(mres, file = paste0("data/calc_vaccDep_", dep, ".csv"), row.names = FALSE)
  
  # Append to the output
  out <- rbind(out, mres)
  
}
# Debut 7h45
# fin vers 9h45

# Debut 15h45
# Fin 17h41
```


```{r}
# Save ouput
write.csv(out, file = "data/calc_vaccDep_alldeps.csv", row.names = FALSE)

save(out, file = "resultsbyDep.RData")
```

## Tests

```{r}
inj <- "1"
agcl <- "Tous âges"
dep <- "75"
mres1 <- out[out$inj == inj & out$dep == dep, ]
mm <- mres1[mres1$agcl == agcl, ]

names(mm)
mm$N1EPCI
mm$V1EPCI
mm$V1wEPCI

plot(mm$V1EPCI, mm$V1Res)
points(mm$V1SPF, mm$V1Res, col = 2)
abline(a = 0, b = 1)

plot(mm$V1wSPF, col = colSPF, type = "o")
points(mm$V1wRes, col = colDep, type = "o")
points(mm$V1wEPCI, col = colEPCI, type = "o")


plot(mm$V2wEPCI / mm$N2EPCI - mm$txEPCIw)
plot(mm$V2EPCI / mm$N2EPCI - mm$txEPCI)
plot(mm$V2SPF / mm$N2SPF - mm$txSPF)


```

## Plot vaccination by dep

```{r}
out <- read.csv("data/calc_vaccDep_alldeps.csv")
head(out)

out$date <- as.Date(out$date)
```


```{r}
addLegendDep <- function(){
    par(family = "mono")
    legend("topleft", 
           col = c(colSPF, colDep, colEPCI), lty = c(1, 1, 1), 
           legend = c(paste0("SPF,   denom. INSEE 2020, couverture finale : ", round(100*mend$txSPF, 1), "%"), 
                      paste0("Ameli, denom. INSEE 2020, couverture finale : ", round(100*mend$txRes, 1), "%"), 
                      paste0("Ameli, assurés sociaux,   couverture finale : ", round(100*mend$txEPCI, 1), "%")), lwd = lwdCurve, bty = "o", cex = 0.8, box.col = "white", box.lwd = -1)
    
    par(family = "sans")
}



inj <- "1"
agcl <- "40-64"#"Tous âges"
dep <- "75"

colLine <- gray(0.9)
lwdCurve <- 2

mres1 <- out[out$inj == inj & out$dep == dep, ]
mm <- mres1[mres1$agcl == agcl, ]

mend <- mm[mm$date == lastDate, ]

thetype <- "o"
thecex <- 0.3

tit1 <- ifelse(inj == "1", "Premières injections", "Vaccinations complètes")

tit2 <- ifelse(agcl == "Tous âges", "", "ans")

depname <- departements[departements$dep == dep, "libelle"]

agg <- agcl
if(agg == "Tous âges") agg <- "tousages"
if(agg == "75+") agg <- "75-etplus"



#fname <- paste0("pics/vaccComparaison_cum_", agg, inj, ".pdf")
#pdf(file = fname, width = 5, height = 4.5)
par(las = 1, mgp = c(2.25, 0.5, 0), tck = -0.015, xpd = TRUE, 
    mar = c(3.25, 3.25, 3.25, 2.5))

# Initialize plot
plot(mm$date, mm$txEPCI, col = colEPCI, ylim = c(0, 1), 
     xaxs = "i", yaxs = "i", 
     xlab = "Date", ylab = "Proportion estimée de la population", 
     main = paste(tit1, "cumulées, ", agcl, tit2, "\n", dep, "-", depname), 
     type = "n", axes = FALSE)

# Add graduations    
par(xpd = FALSE)
for(i in seq(0, 1, by = 0.1)) abline(h = i, col = gray(0.9), lwd = 1)
for(i in seq(0.05, 0.95, by = 0.1)) abline(h = i, col = colLine, lwd = 0.5)
par(xpd = TRUE)


# Add legend
addLegendDep()

# Add axes
axis(2, at = seq(0, 1, by = 0.1))
axis(4, at = seq(0, 1, by = 0.1))
axis.Date(1, at = seq(as.Date("2021-01-01"), max(mm$date), by="months"))

# Add points
points(mm$date, mm$txEPCI, col = colEPCI, type = thetype, cex = thecex, lwd = lwdCurve)
if(!is.element(agcl, c("00-39"))){
  points(mm$date, mm$txRes, col = colDep, 
         type = thetype, cex = thecex, lwd = lwdCurve)
}
points(mm$date, mm$txSPF, col = colSPF, 
       type = thetype, cex = thecex, lwd = lwdCurve)

#dev.off()


#-----------------------------------------------



#fname2 <- paste0("pics/vaccComparaison_sem_", agg, inj, ".pdf")
#pdf(file = fname2, width = 5, height = 4.5)
par(las = 1, mgp = c(2.25, 0.5, 0), tck = -0.015, xpd = TRUE, 
    mar = c(3.25, 3.25, 3.25, 2.5))
plot(mm$date, mm$txEPCIw, col = colEPCI, ylim = c(0, 0.1), type = "n", 
     xaxs = "i", yaxs = "i", 
     xlab = "Date", ylab = "Proportion estimée de la population", 
     cex = thecex, 
     main = paste(tit1, "de la semaine, ", agcl, "\n Dep", dep, "-", depname), 
     axes = FALSE)

# Add graduations
par(xpd = FALSE)
for(i in seq(0, 0.1, by = 0.01)) abline(h = i, col = colLine)
par(xpd = TRUE)

# Add legend
addLegendDep()

# Add axes
axis.Date(1, at = seq(as.Date("2021-01-01"), max(mm$date), by="months"))
axis(2, at = seq(0, 0.1, by = 0.01))
axis(4, at = seq(0, 0.1, by = 0.01))


# Add points
points(mm$date, mm$txEPCIw, col = colEPCI, type = thetype, cex = thecex, lwd = lwdCurve)
if(!is.element(agcl, c("00-39"))){
  points(mm$date, mm$txResw, col = colDep, type = thetype, cex = thecex, lwd = lwdCurve)
}
points(mm$date, mm$txSPFw, col = colSPF, type = thetype, cex = thecex, lwd = lwdCurve)

#dev.off()
```

# Changing population size?

Check whether population size varies over time for one EPCI x age class (whether e.g. deaths are taken into account)

```{r}
# Aggregate data by EPCI and age class

# Variance in pop sizes
popvar <- aggregate(as.numeric(dat.EPCI$population_carto), by = list(epci = dat.EPCI$epci, clage = dat.EPCI$classe_age), FUN = var)

# Difference between max and min pop size
popvar2 <- aggregate(as.numeric(dat.EPCI$population_carto), by = list(epci = dat.EPCI$epci, libelle_epci = dat.EPCI$libelle_epci, reg = dat.EPCI$reg_code, clage = dat.EPCI$classe_age), FUN = function(x) diff(range(x)))

# Check that there are indeed multiple values when we do this!
unique(aggregate(as.numeric(dat.EPCI$population_carto), by = list(epci = dat.EPCI$epci, clage = dat.EPCI$classe_age), FUN = length)$x)

# EPCI 99999999 is when not residency is found
popvar2[popvar2$epci == 999999999, ]
popvar[popvar$epci == 999999999, ]
# Remove this point
popvar2 <- popvar2[popvar2$epci != 999999999, ]
popvar <- popvar[popvar$epci != 999999999, ]

# Rename last column with more informative name than "x"
names(popvar2)[ncol(popvar2)] <- "diffRange"
names(popvar)[ncol(popvar)] <- "var"

# Plot result 
plot(popvar2$epci, popvar2$diffRange)
plot(popvar$epci, popvar$var)

# Nicer version of the plot
par(las = 1)
plot(popvar2$diffRange, xlab = "epci x age class", ylab = "max difference in pop size")

# Check that it is indeed just 0 like the plots suggest
unique(popvar2$diffRange)
unique(popvar$var)
```