library(lfe)
library(reshape2)
library(htmltab)
library(purrr)
library(httr)
library(plyr)
library(readxl)
library(zoo)
library(RColorBrewer)
library("texreg")
darkcols <- c(brewer.pal(8, "Dark2"),brewer.pal(8, "Paired"))
######################################################################################################

######################################################################################################
# Useful functions
######################################################################################################
simplelag <-function(x,by=NULL,mylag=1,outside=NA)
{
  myend<-length(x)- mylag
  if (!is.null(by)) {
    lby<-c(replicate(mylag,""),as.character(by[1:myend]))
    y0<-c(replicate(mylag,outside),x[1:myend])
    y<-ifelse(as.character(by)==lby,y0,outside)
    
  }
  else {
    y<-c(replicate(mylag,outside),x[1:myend])
  }
}
retain<-function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE), length(x)+1)
  values <- c(outside,x[event==TRUE])
  y<- rep(values, diff(indices)) 
}
######################################################################################################


######################################################################################################
# SOURCES
######################################################################################################
setwd("D:\\Olivier\\ownCloud\\Data\\Covid19")

#dc<-read.csv2("French_Dept/DC_Jan-Avr_2018-2020_det.csv")
dc<-read.csv2("France_dc/DC_jan2018-avr2020_det.csv")
dc_age<-read.csv2("France_dc/2020-04-30_deces_parage_regdept.csv")

#Covid deaths per age and region
dccov<-read.csv2("https://www.data.gouv.fr/fr/datasets/r/08c18e08-6780-452d-9b8c-ae244ad529b3")

#Covid deaths per gender and departement
dccov_2<-read.csv2("https://www.data.gouv.fr/fr/datasets/r/63352e38-d353-4b54-bfd1-f1b3ee1cabd7")

#French regions and departments
france<-read.csv("https://raw.githubusercontent.com/opencovid19-fr/data/master/dist/chiffres-cles.csv",encoding="UTF-8")


#DATA ON AGE AND INCOME PER DEPARTMENT

dep_age<-read_excel("French_Dept/TCRD_021.xls",skip=3)
str(dep_age)
colnames(dep_age)<-c("code","dept","number","p_female","p_male","p_age0_24","p_age25_59",
                          "p_age60_pl","p_age75_pl")
dep_inc<-read_excel("French_Dept/TCRD_022.xls",skip=4)
colnames(dep_inc)<-c("code","dept","number","p_imposed","median_inc","D1_inc","D9_inc")
str(dep_inc)

#SOURCES
france<-read.csv("https://raw.githubusercontent.com/opencovid19-fr/data/master/dist/chiffres-cles.csv",encoding="UTF-8")


# Original source for French population. Downloaded on my web site to increase spead.
# url<-"https://www.insee.fr/fr/statistiques/fichier/3545833/ensemble.xls"
# GET(url, write_disk(tf <- tempfile(fileext = ".xls")))
# france_pop <- read_excel(tf, 2L,skip=7)
# write.csv(france_pop,file="france_pop.csv",fileEncoding="UTF-8")
france_pop<-read.csv("http://olivier.godechot.free.fr/hopfichiers/france_pop.csv",encoding="UTF-8")
colnames(france_pop)<-c("X","code_region","nom_region","code_departement","nom_departement","nombre_arrondissements"
                        ,"nombre_cantons","nombre_communes","population_municipale","population_totale")

france_reg_pop<-as.data.frame.table(tapply(france_pop$population_totale,france_pop$code_region,sum))
colnames(france_reg_pop)<-c("code","population")

france_dep_pop<-as.data.frame.table(tapply(france_pop$population_totale,france_pop$code_departement,sum))
colnames(france_dep_pop)<-c("code","population")

france_dep_pop2<-as.data.frame(map2("https://fr.wikipedia.org/wiki/Liste_des_d%C3%A9partements_fran%C3%A7ais_class%C3%A9s_par_population_et_superficie",1, htmltab,rm_nodata_cols = F))
lasctol<-ncol(france_dep_pop2)
colnames(france_dep_pop2)[c(2,lasctol)]<-c("code","area")
str(france_dep_pop2)

#Hospital Beds
france_beds<-read.csv2("http://olivier.godechot.free.fr/hopfichiers/france_hospital_beds.csv",
                       skip=6,header=TRUE,encoding="latin1")
france_beds<-france_beds[,c(1,15,27)]
colnames(france_beds)<-c("unit","beds2000","beds2012")
######################################################################################################



######################################################################################################
# DATA STEPS
######################################################################################################

# dc data
dc$AGE<-dc$ADEC-dc$ANAIS
dc$un<-1
dc$date<-as.Date(paste(dc$ADEC,dc$MDEC,dc$JDEC,sep="-"))
max(dc$date)
dc$day<-as.numeric(dc$date)
dc$day2<-ifelse(format(dc$date,"%Y")=="2018",dc$day+365*2,
                         ifelse(format(dc$date,"%Y")=="2019",dc$day+365,
                                dc$day))
dayconfinement<-as.Date("2020-03-17")
daystart<-as.Date("2020-01-01")
dayend<-max(dc$date)

dc$confinement<-(dc$date>=dayconfinement)*1
dc$confinement_period<-(dc$day2>=as.integer(dayconfinement))*1
dc$day3<-dc$day2-as.numeric(daystart)+1
dc$week<-floor((dc$day2-as.numeric(daystart))/7)+1
lastweeknbdays<-as.numeric(dayend)-as.numeric(daystart)+1-floor((as.numeric(dayend)-as.numeric(daystart))/7)*7


dc$age_cat<-cut(dc$AGE,c(-1,29,49,64,74,84,130))
dc$weekday<-weekdays(dc$date)
dc$un<-1
dc$cl_age90<-10*(floor(dc$AGE/10)+1)-1
dc$cl_age90<-ifelse(dc$cl_age90>89,90,dc$cl_age90)


# dccov data
dccov$date<-as.Date(dccov$jour)
dccov$day<-as.numeric(dccov$date)
dccov$week<-floor((dccov$day-as.numeric(daystart))/7)+1
dayend2<-max(dccov$day)


dccov<-dccov[order(dccov$reg,dccov$cl_age90,dccov$day),]
dccov$new_deaths<-dccov$dc-simplelag(dccov$dc,paste(dccov$reg,dccov$cl_age90,sep="#"),outside=0)

dc<-dc[dc$day2<=as.integer(dayend2),]
lastweeknbdays2<-as.numeric(dayend2)-as.numeric(daystart)+1-floor((as.numeric(dayend2)-as.numeric(daystart))/7)*7



#french_beds
france_beds$code_dep<-ifelse(substr(france_beds$unit,1,1)=="D",substr(france_beds$unit,2,4),NA)
france_beds$code_reg<-ifelse(substr(france_beds$unit,1,1)=="R",substr(france_beds$unit,2,4),NA)


#France
france$location<-france$maille_nom
france$total_deaths<-france$deces
france$date<-as.Date(france$date)
france$day<-as.numeric(france$date)
france$week<-floor((france$day-as.numeric(daystart))/7)+1

france_fr<-france[france$maille_code=="FRA",]
france_fr<-france_fr[order(france_fr$day),]


#French all
france_fr<-france[france$maille_code %in% "FRA",]

france_fr<-france_fr[order(france_fr$location,france_fr$date,-france_fr$total_deaths),]
france_fr$sel<-ave((france_fr$location !=""),paste(france_fr$location,france_fr$date),FUN=function(x) cumsum(x))
france_fr<-france_fr[france_fr$sel %in% 1,]
row.names(france_fr)<-NULL

france_fr$cas_confirmes_approx<-ave(france_fr$cas_confirmes,france_fr$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_fr$death_approx<-ave(france_fr$total_deaths,france_fr$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_fr$hospitalises_approx<-ave(france_fr$hospitalises,france_fr$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_fr$gueris_approx<-ave(france_fr$gueris,france_fr$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_fr$reanimation_approx<-ave(france_fr$reanimation,france_fr$location,FUN=function(x) na.approx(x,na.rm=FALSE))

france_fr$mild_cases_approx<-france_fr$cas_confirmes_approx-france_fr$death_approx-france_fr$hospitalises_approx-france_fr$gueris_approx


france_fr$hospitalises_evol<-france_fr$hospitalises_approx-simplelag(france_fr$hospitalises_approx,france_fr$location,outside=0)
france_fr$reanimation_evol<-france_fr$reanimation_approx-simplelag(france_fr$reanimation_approx,france_fr$location,outside=0)

france_fr$lcas_confirmes_approx<-simplelag(france_fr$cas_confirmes_approx,outside=0)
france_fr$lmild_cases_approx<-simplelag(france_fr$mild_cases_approx,outside=0)
france_fr$llocation<-simplelag(as.character(france_fr$location),outside = "")

france_fr$continuity_approx0<-ifelse(as.character(france_fr$location) != france_fr$llocation,
                                      0,
                                      ifelse(is.na(france_fr$lmild_cases_approx)==TRUE,
                                             france_fr$lcas_confirmes_approx-france_fr$hospitalises_approx+france_fr$gueris_approx+france_fr$death_approx,
                                             france_fr$lmild_cases_approx))

france_fr$continuity_approx<-retain(france_fr$continuity_approx0,
                                     (is.na(france_fr$cas_confirmes_approx)==TRUE & is.na(france_fr$lcas_confirmes_approx)==FALSE) 
                                     | as.character(france_fr$location) != france_fr$llocation
                                     ,0)

france_fr$total_cases<-ifelse(is.na(france_fr$cas_confirmes_approx) %in% FALSE,france_fr$cas_confirmes_approx,
                               france_fr$hospitalises_approx+france_fr$gueris_approx+france_fr$death_approx+france_fr$continuity_approx)

france_fr$new_deaths<-france_fr$total_deaths-simplelag(france_fr$total_deaths,france_fr$location,outside=0)
france_fr$new_cases<-france_fr$total_cases-simplelag(france_fr$total_cases,france_fr$location,outside=0)


france_fr$code<-substr(france_fr$maille_code,5,6)
#france_fr<-merge(france_fr,france_fr_pop,by="code",all.x=TRUE)

france_fr$cases_2_pop<-10000*france_fr$total_cases/sum(france_pop$population_totale)
france_fr$deaths_2_pop<-100000*france_fr$death_approx/sum(france_pop$population_totale)

france_fr<-france_fr[order(france_fr$location,france_fr$date,-france_fr$total_deaths),]


#French departments
france_dep<-france[france$granularite %in% "departement",]
france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]
france_dep$sel<-ave((france_dep$location !=""),paste(france_dep$location,france_dep$date),FUN=function(x) cumsum(x))
france_dep<-france_dep[france_dep$sel %in% 1,]

france_dep$cas_confirmes_approx<-ave(france_dep$cas_confirmes,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$death_approx<-ave(france_dep$total_deaths,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$hospitalises_approx<-ave(france_dep$hospitalises,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$gueris_approx<-ave(france_dep$gueris,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$reanimation_approx<-ave(france_dep$reanimation,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))

france_dep$hospitalises_evol<-france_dep$hospitalises_approx-simplelag(france_dep$hospitalises_approx,france_dep$location,outside=0)
france_dep$reanimation_evol<-france_dep$reanimation_approx-simplelag(france_dep$reanimation_approx,france_dep$location,outside=0)


france_dep$mild_cases_approx<-france_dep$cas_confirmes_approx-france_dep$death_approx-france_dep$hospitalises_approx-france_dep$gueris_approx

france_dep$lcas_confirmes_approx<-simplelag(france_dep$cas_confirmes_approx,outside=0)
france_dep$lmild_cases_approx<-simplelag(france_dep$mild_cases_approx,outside=0)
france_dep$llocation<-simplelag(as.character(france_dep$location),outside = "")

france_dep$continuity_approx0<-ifelse(as.character(france_dep$location) != france_dep$llocation,
                                      0,
                                      ifelse(is.na(france_dep$lmild_cases_approx)==TRUE,
                                             france_dep$lcas_confirmes_approx-france_dep$hospitalises_approx+france_dep$gueris_approx+france_dep$death_approx,
                                             france_dep$lmild_cases_approx))

france_dep$continuity_approx<-retain(france_dep$continuity_approx0,
                                     (is.na(france_dep$cas_confirmes_approx)==TRUE & is.na(france_dep$lcas_confirmes_approx)==FALSE) 
                                     | as.character(france_dep$location) != france_dep$llocation
                                     ,0)
france_dep$total_cases<-ifelse(is.na(france_dep$cas_confirmes_approx) %in% FALSE,france_dep$cas_confirmes_approx,
                               france_dep$hospitalises_approx+france_dep$gueris_approx+france_dep$death_approx+france_dep$continuity_approx)


france_dep$new_deaths<-france_dep$total_deaths-simplelag(france_dep$total_deaths,france_dep$location,outside=0)
france_dep$new_cases<-france_dep$total_cases-simplelag(france_dep$total_cases,france_dep$location,outside=0)


france_dep$code<-substr(france_dep$maille_code,5,7)
france_dep<-merge(france_dep,france_dep_pop,by="code",all.x=TRUE)
france_dep_pop2$code<-ifelse(nchar(france_dep_pop2$code)==1,paste("0",france_dep_pop2$code,sep=""),france_dep_pop2$code)
france_dep<-merge(france_dep,france_dep_pop2[,c("code","area")],by="code",all.x=TRUE)
france_dep$area<-iconv(france_dep$area,"UTF-8","latin1")
france_dep$area2<-gsub("\u00A0","",france_dep$area)



france_dep$density<-france_dep$population/as.numeric(france_dep$area2)

france_dep$cases_2_pop<-10000*france_dep$total_cases/france_dep$population
france_dep$deaths_2_pop<-100000*france_dep$death_approx/france_dep$population

france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep<-merge(france_dep,france_beds,by.x="code",by.y="code_dep",all.x=TRUE)
france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep$nb_beds<-france_dep$beds2012*france_dep$population/100000
france_dep$hospitalized_per_beds<-france_dep$hospitalises/france_dep$nb_beds
france_dep$reanimation_per_beds<-france_dep$reanimation/france_dep$nb_beds


france_dep$calendar_day<-as.numeric(france_dep$date)-min(as.numeric(france_dep$date))+1
france_dep$confinement_day<-pmax(as.numeric(france_dep$date)-as.numeric(as.Date("2020-03-17"))+1,0)

france_dep$threshold_day<-ave(ifelse(france_dep$death_approx>10,as.numeric(france_dep$date),999999999),france_dep$location,FUN=function(x) min(x,na.rm=TRUE))
france_dep$day_from_threshold<-ifelse(france_dep$threshold_day<999999999,as.numeric(france_dep$date)-france_dep$threshold_day+1,NA)
france_dep$day_from_threshold<-ifelse(france_dep$day_from_threshold<=0,NA,france_dep$day_from_threshold)
france_dep$time<-factor(france_dep$day_from_threshold)
france_dep$nb_beds_2012<-france_dep$nb_beds
france_dep$beds_2_pop_2012<-france_dep$beds2012
france_dep$beds_decline<-france_dep$beds2000-france_dep$beds2012

france_dep<-merge(france_dep,dep_age,by="code",all.x=TRUE)
france_dep<-merge(france_dep,dep_inc,by="code",all.x=TRUE)

france_dep<-france_dep[order(france_dep$location,france_dep$date),]

france_dep$total_cases_approx<-ave(france_dep$total_cases,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$total_deaths_approx<-ave(france_dep$total_deaths,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$new_cases_approx<-france_dep$total_cases_approx-simplelag(france_dep$total_cases_approx,france_dep$location,outside=0)
france_dep$new_deaths_approx<-france_dep$total_deaths_approx-simplelag(france_dep$total_deaths_approx,france_dep$location,outside=0)

#French regions
france_reg<-france[france$granularite %in% "region",]

france_reg<-france_reg[order(france_reg$location,france_reg$date,-france_reg$total_deaths),]
france_reg$sel<-ave((france_reg$location !=""),paste(france_reg$location,france_reg$date),FUN=function(x) cumsum(x))
france_reg<-france_reg[france_reg$sel %in% 1,]
row.names(france_reg)<-NULL

france_reg$cas_confirmes_approx<-ave(france_reg$cas_confirmes,france_reg$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_reg$death_approx<-ave(france_reg$total_deaths,france_reg$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_reg$hospitalises_approx<-ave(france_reg$hospitalises,france_reg$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_reg$gueris_approx<-ave(france_reg$gueris,france_reg$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_reg$reanimation_approx<-ave(france_reg$reanimation,france_reg$location,FUN=function(x) na.approx(x,na.rm=FALSE))

france_reg$mild_cases_approx<-france_reg$cas_confirmes_approx-france_reg$death_approx-france_reg$hospitalises_approx-france_reg$gueris_approx


france_reg$hospitalises_evol<-france_reg$hospitalises_approx-simplelag(france_reg$hospitalises_approx,france_reg$location,outside=0)
france_reg$reanimation_evol<-france_reg$reanimation_approx-simplelag(france_reg$reanimation_approx,france_reg$location,outside=0)

france_reg$lcas_confirmes_approx<-simplelag(france_reg$cas_confirmes_approx,outside=0)
france_reg$lmild_cases_approx<-simplelag(france_reg$mild_cases_approx,outside=0)
france_reg$llocation<-simplelag(as.character(france_reg$location),outside = "")

france_reg$continuity_approx0<-ifelse(as.character(france_reg$location) != france_reg$llocation,
                                      0,
                                      ifelse(is.na(france_reg$lmild_cases_approx)==TRUE,
                                             france_reg$lcas_confirmes_approx-france_reg$hospitalises_approx+france_reg$gueris_approx+france_reg$death_approx,
                                             france_reg$lmild_cases_approx))

france_reg$continuity_approx<-retain(france_reg$continuity_approx0,
                                     (is.na(france_reg$cas_confirmes_approx)==TRUE & is.na(france_reg$lcas_confirmes_approx)==FALSE) 
                                     | as.character(france_reg$location) != france_reg$llocation
                                     ,0)

france_reg$total_cases<-ifelse(is.na(france_reg$cas_confirmes_approx) %in% FALSE,france_reg$cas_confirmes_approx,
                               france_reg$hospitalises_approx+france_reg$gueris_approx+france_reg$death_approx+france_reg$continuity_approx)

france_reg$new_deaths<-france_reg$total_deaths-simplelag(france_reg$total_deaths,france_reg$location,outside=0)
france_reg$new_cases<-france_reg$total_cases-simplelag(france_reg$total_cases,france_reg$location,outside=0)


france_reg$code<-substr(france_reg$maille_code,5,6)
france_reg<-merge(france_reg,france_reg_pop,by="code",all.x=TRUE)

france_reg$cases_2_pop<-10000*france_reg$total_cases/france_reg$population
france_reg$deaths_2_pop<-100000*france_reg$death_approx/france_reg$population

france_reg<-france_reg[order(france_reg$location,france_reg$date,-france_reg$total_deaths),]


#French departments
france_dep<-france[france$granularite %in% "departement",]
france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]
france_dep$sel<-ave((france_dep$location !=""),paste(france_dep$location,france_dep$date),FUN=function(x) cumsum(x))
france_dep<-france_dep[france_dep$sel %in% 1,]

france_dep$cas_confirmes_approx<-ave(france_dep$cas_confirmes,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$death_approx<-ave(france_dep$total_deaths,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$hospitalises_approx<-ave(france_dep$hospitalises,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$gueris_approx<-ave(france_dep$gueris,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$reanimation_approx<-ave(france_dep$reanimation,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))

france_dep$hospitalises_evol<-france_dep$hospitalises_approx-simplelag(france_dep$hospitalises_approx,france_dep$location,outside=0)
france_dep$reanimation_evol<-france_dep$reanimation_approx-simplelag(france_dep$reanimation_approx,france_dep$location,outside=0)


france_dep$mild_cases_approx<-france_dep$cas_confirmes_approx-france_dep$death_approx-france_dep$hospitalises_approx-france_dep$gueris_approx

france_dep$lcas_confirmes_approx<-simplelag(france_dep$cas_confirmes_approx,outside=0)
france_dep$lmild_cases_approx<-simplelag(france_dep$mild_cases_approx,outside=0)
france_dep$llocation<-simplelag(as.character(france_dep$location),outside = "")

france_dep$continuity_approx0<-ifelse(as.character(france_dep$location) != france_dep$llocation,
                                      0,
                                      ifelse(is.na(france_dep$lmild_cases_approx)==TRUE,
                                             france_dep$lcas_confirmes_approx-france_dep$hospitalises_approx+france_dep$gueris_approx+france_dep$death_approx,
                                             france_dep$lmild_cases_approx))

france_dep$continuity_approx<-retain(france_dep$continuity_approx0,
                                     (is.na(france_dep$cas_confirmes_approx)==TRUE & is.na(france_dep$lcas_confirmes_approx)==FALSE) 
                                     | as.character(france_dep$location) != france_dep$llocation
                                     ,0)
france_dep$total_cases<-ifelse(is.na(france_dep$cas_confirmes_approx) %in% FALSE,france_dep$cas_confirmes_approx,
                               france_dep$hospitalises_approx+france_dep$gueris_approx+france_dep$death_approx+france_dep$continuity_approx)


france_dep$new_deaths<-france_dep$total_deaths-simplelag(france_dep$total_deaths,france_dep$location,outside=0)
france_dep$new_cases<-france_dep$total_cases-simplelag(france_dep$total_cases,france_dep$location,outside=0)


france_dep$code<-substr(france_dep$maille_code,5,7)
france_dep<-merge(france_dep,france_dep_pop,by="code",all.x=TRUE)
france_dep_pop2$code<-ifelse(nchar(france_dep_pop2$code)==1,paste("0",france_dep_pop2$code,sep=""),france_dep_pop2$code)
france_dep<-merge(france_dep,france_dep_pop2[,c("code","area")],by="code",all.x=TRUE)
france_dep$area<-iconv(france_dep$area,"UTF-8","latin1")
france_dep$area2<-gsub("\u00A0","",france_dep$area)



france_dep$density<-france_dep$population/as.numeric(france_dep$area2)

france_dep$cases_2_pop<-10000*france_dep$total_cases/france_dep$population
france_dep$deaths_2_pop<-100000*france_dep$death_approx/france_dep$population

france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep<-merge(france_dep,france_beds,by.x="code",by.y="code_dep",all.x=TRUE)
france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep<-france_dep[order(france_dep$location,france_dep$date,-france_dep$total_deaths),]

france_dep$nb_beds<-france_dep$beds2012*france_dep$population/100000
france_dep$hospitalized_per_beds<-france_dep$hospitalises/france_dep$nb_beds
france_dep$reanimation_per_beds<-france_dep$reanimation/france_dep$nb_beds


france_dep$calendar_day<-as.numeric(france_dep$date)-min(as.numeric(france_dep$date))+1
france_dep$confinement_day<-pmax(as.numeric(france_dep$date)-as.numeric(as.Date("2020-03-17"))+1,0)

france_dep$threshold_day<-ave(ifelse(france_dep$death_approx>10,as.numeric(france_dep$date),999999999),france_dep$location,FUN=function(x) min(x,na.rm=TRUE))
france_dep$day_from_threshold<-ifelse(france_dep$threshold_day<999999999,as.numeric(france_dep$date)-france_dep$threshold_day+1,NA)
france_dep$day_from_threshold<-ifelse(france_dep$day_from_threshold<=0,NA,france_dep$day_from_threshold)
france_dep$time<-factor(france_dep$day_from_threshold)
france_dep$nb_beds_2012<-france_dep$nb_beds
france_dep$beds_2_pop_2012<-france_dep$beds2012
france_dep$beds_decline<-france_dep$beds2000-france_dep$beds2012

france_dep<-merge(france_dep,dep_age,by="code",all.x=TRUE)
france_dep<-merge(france_dep,dep_inc,by="code",all.x=TRUE)

france_dep<-france_dep[order(france_dep$location,france_dep$date),]

france_dep$total_cases_approx<-ave(france_dep$total_cases,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$total_deaths_approx<-ave(france_dep$total_deaths,france_dep$location,FUN=function(x) na.approx(x,na.rm=FALSE))
france_dep$new_cases_approx<-france_dep$total_cases_approx-simplelag(france_dep$total_cases_approx,france_dep$location,outside=0)
france_dep$new_deaths_approx<-france_dep$total_deaths_approx-simplelag(france_dep$total_deaths_approx,france_dep$location,outside=0)



france_reg$new_deaths<-france_fr$total_deaths-simplelag(france_reg$total_deaths,france_reg$location,outside=0)

tapply(france_fr$new_deaths,france_fr$week,sum,na.rm=TRUE)[1:9]
sum(france_fr$new_deaths,na.rm=TRUE)
str(france_fr)

nursing_home_deaths<-france_fr$deces_ehpad[france_fr$date=="2020-04-14"]
france_fr$deces[france_fr$date=="2020-04-14"]

###############################################
# Graphs 2020 over-mortality per age
###############################################

mycol<-c("gray60","gray70","black","red","darkorange1","darkviolet","deepskyblue")
mypch<-c(2,6,17,16,20,15,26)
mylty<-c(2,2,1,1,1,1,2)

results<-NULL
gdata<-function(agecat,data=dc)
{
  # data<-dc2
  # agecat<-0
   agecat2<-paste("Age: ",floor(agecat/10)*10,"-",agecat,sep="")
  if(agecat==0) {agecat2<-"Age: all"}
  if(agecat==90) {agecat2<-"Age: 90 +"}
  
  
  mydc<-data
  w<-as.data.frame.matrix(table(mydc$week[mydc$cl_age90==agecat],mydc$ADEC[mydc$cl_age90==agecat]))
  colnames(w)<-c("y2018","y2019","y2020")
  w$week<-as.numeric(rownames(w))
  w$y2020[w$week==16]<-NA #w$y2020[w$week==16]*7/lastweeknbdays
  w$y2020[w$week==17]<-NA
  w$y2020[w$week==18]<-NA
  w$y2019[w$week==18]<-w$y2019[w$week==18]*7/lastweeknbdays2
  w$y2018[w$week==18]<-w$y2018[w$week==18]*7/lastweeknbdays2
  w$y2018_19<-0.5*(w$y2018+w$y2019)
  
  
  wcov<-data.frame(tapply(dccov$new_deaths[dccov$cl_age90==agecat],dccov$week[dccov$cl_age90==agecat],sum))
  colnames(wcov)<-"covid_hospital_deaths"
  wcov$week<-as.numeric(rownames(wcov))
  
  w<-merge(w,wcov,by="week",all.x=TRUE)
  w<-w[order(w$week),]
  
  correction_week12<-sum(tapply(france_fr$new_deaths,france_fr$week,sum,na.rm=TRUE)[9])/sum(tapply(france_fr$new_deaths,france_fr$week,sum,na.rm=TRUE)[1:9])
  correction_week11<-sum(tapply(france_fr$new_deaths,france_fr$week,sum,na.rm=TRUE)[8])/sum(tapply(france_fr$new_deaths,france_fr$week,sum,na.rm=TRUE)[1:9])
  correction_week10<-1-correction_week11-correction_week12
  
  w$covid_hospital_deaths[w$week==9]<-0
  w$covid_hospital_deaths[w$week==10]<-correction_week10*w$covid_hospital_deaths[w$week==12]
  w$covid_hospital_deaths[w$week==11]<-correction_week11*w$covid_hospital_deaths[w$week==12]
  w$covid_hospital_deaths[w$week==12]<-correction_week12*w$covid_hospital_deaths[w$week==12]
  w$y2020_without_cov<-w$y2020-w$covid_hospital_deaths
  w$covid_hospital_deaths_scaled<-retain(w$y2020,w$week==9)+w$covid_hospital_deaths
  w
  
  
  rightmax<-max(
    c(max(w$covid_hospital_deaths_scaled[w$week>9],na.rm=TRUE),
    max(w$y2020_without_cov[w$week>9],na.rm=TRUE),
    max(w$y2020[w$week>9],na.rm=TRUE))-w$covid_hospital_deaths_scaled[w$week==9])
  if (rightmax>10)  {
    roundUpNice <- function(x, nice=c(1,2,4,5,6,8,10)) {
      if(length(x) != 1) stop("'x' must be of length 1")
      10^floor(log10(x)) * nice[[which(x <= 10^floor(log10(x)) * nice)[[1]]]]
    }
    rightmax2<-roundUpNice(rightmax)
    ticks<-c(0,rightmax2/4,2*rightmax2/4,3*rightmax2/4,4*rightmax2/4,5*rightmax2/4)
    
  } else {ticks<-c(0,rightmax) }
    
  posticks<-ticks+w$covid_hospital_deaths_scaled[w$week==9]    
  
  g<-w[,c("y2018","y2019","y2018_19","y2020","y2020_without_cov","covid_hospital_deaths_scaled")]
  
  matplot(g,pch=mypch,
          ylab=list("Number of deaths per week",cex=0.8),xlab="",col=mycol,xaxt="n",cex.axis=0.8,las=2)
  
  min_a_2 <- pmin(w$y2018_19[w$week>8 & w$week<16]#,w$y2020_without_cov[w$week>9 & w$week<16]
                  )
  max_a_2 <- pmax(#w$y2018_19[w$week>9 & w$week<16],
                  w$y2020_without_cov[w$week>8 & w$week<16])
  polygon(c(w$week[w$week>8 & w$week<16], rev(w$week[w$week>8 & w$week<16])), c(max_a_2 ,rev(min_a_2)), 
          col = adjustcolor("goldenrod1",alpha.f=0.5),border=NA)
  
  min_a <- pmin(w$y2020[w$week>8 & w$week<16],w$y2020_without_cov[w$week>8 & w$week<16])
  max_a <- pmax(w$y2020[w$week>8 & w$week<16],w$y2020_without_cov[w$week>8 & w$week<16])
  polygon(c(w$week[w$week>8 & w$week<16], rev(w$week[w$week>8 & w$week<16])), c(max_a ,rev(min_a)), 
          col = adjustcolor("plum2",alpha.f=0.5),border=NA)
  
  title(agecat2)
  matlines(g,col=mycol,lty=mylty,pch=mypch)
  matpoints(g,col=mycol,pch=mypch)
  abline(v=9,lty=2,col="darkviolet")
  abline(h=w$y2020[w$week==9],lty=2,col="darkviolet")
  axis(1, at=c(1:18),c("Jan 1-7","Jan 8-14","Jan 15-21","Jan 22-28","Jan 29-Feb 4",
                       "Feb 5-11","Feb 12-18","Feb 19-25","Feb 26-Mar 3","Mar 4-10",
                       "Mar 11-17","Mar 18-24","Mar 25-31","Apr 1-7","Apr 8-14",
                       "Apr 15-21","Apr 22-28","Apr 29-May 3"),cex.axis=0.8,las=2) 
  axis(4, at=posticks, 
       labels=ticks,
       col.ticks = "purple",col="purple",col.axis="purple",las=1,cex.axis=0.8)
  abline(v=11.5,lty=2,col="deepskyblue")
  
  m<-as.data.frame.table(cbind(tapply(w$y2018[w$week<16],w$week[w$week<16]>9,sum),
    tapply(w$y2019[w$week<16],w$week[w$week<16]>9,sum),
    tapply(w$y2020[w$week<16],w$week[w$week<16]>9,sum),
    tapply(w$y2020_without_cov[w$week<16],w$week[w$week<16]>9,sum)))
  colnames(m)<-c("period","year","deaths")
  levels(m$year)<-c("y2018","y2019","y2020","y2020_without_covid")
  m$deaths[7]<-m$deaths[5]
  
  mod1<-glm(deaths~I(period=="TRUE")*I(year=="y2020"), data=m[1:6,],family="poisson")
  summary(mod1)
  
  c1<-mod1$coefficients[4]
  se1<-summary(mod1)$coefficients[, 2][4]
  p1<-summary(mod1)$coefficients[, 4][4]
  
  r1<-exp(mod1$coefficients[4])-1
  r1l<-exp(mod1$coefficients[4]-1.96*se1)-1
  r1u<-exp(mod1$coefficients[4]+1.96*se1)-1
  
  d1<-m$deaths[6]-m$deaths[6]/(1+r1)
  d1l<-m$deaths[6]-m$deaths[6]/(1+r1l)
  d1u<-m$deaths[6]-m$deaths[6]/(1+r1u)
  
  mod2<-glm(deaths~I(period=="TRUE")*I(year=="y2020_without_covid"), data=m[c(1:4,7:8),],family="poisson")
  summary(mod2)
  
  c2<-mod2$coefficients[4]
  se2<-summary(mod2)$coefficients[,2][4]
  p2<-summary(mod2)$coefficients[,4][4]
  
  r2<-exp(mod2$coefficients[4])-1
  r2l<-exp(mod2$coefficients[4]-1.96*se2)-1
  r2u<-exp(mod2$coefficients[4]+1.96*se2)-1
  
  d2<-m$deaths[6]-m$deaths[6]/(1+r2)
  d2l<-m$deaths[6]-m$deaths[6]/(1+r2l)
  d2u<-m$deaths[6]-m$deaths[6]/(1+r2u)

  
  mod3<-glm(deaths~I(year=="y2020"), data=m[c(2,4,6),],family="poisson")
  summary(mod3)
  
  c3<-mod3$coefficients[2]
  se3<-summary(mod3)$coefficients[,2][2]
  p3<-summary(mod3)$coefficients[,4][2]
  
  r3<-exp(mod3$coefficients[2])-1
  r3l<-exp(mod3$coefficients[2]-1.96*se3)-1
  r3u<-exp(mod3$coefficients[2]+1.96*se3)-1
  
  d3<-m$deaths[6]-m$deaths[6]/(1+r3)
  d3l<-m$deaths[6]-m$deaths[6]/(1+r3l)
  d3u<-m$deaths[6]-m$deaths[6]/(1+r3u)
  
  d18_19<-0.5*(m$deaths[2]+m$deaths[4])
  d20<-m$deaths[6]
  dcovhos<-m$deaths[6]-m$deaths[8]
  
  
  mod4<-glm(deaths~I(year=="y2020_without_covid"), data=m[c(2,4,8),],family="poisson")
  summary(mod4)
  
  c4<-mod4$coefficients[2]
  se4<-summary(mod4)$coefficients[,2][2]
  p4<-summary(mod4)$coefficients[,4][2]
  
  r4<-exp(mod4$coefficients[2])-1
  r4l<-exp(mod4$coefficients[2]-1.96*se4)-1
  r4u<-exp(mod4$coefficients[2]+1.96*se4)-1
  
  d4<-m$deaths[6]-m$deaths[6]/(1+r4)
  d4l<-m$deaths[6]-m$deaths[6]/(1+r4l)
  d4u<-m$deaths[6]-m$deaths[6]/(1+r4u)
  
  d18_19<-0.5*(m$deaths[2]+m$deaths[4])
  d20<-m$deaths[6]
  dcovhos<-m$deaths[6]-m$deaths[8]
  
  
  mod5<-glm(deaths~I(year=="y2020"), data=m[c(6,8),],family="poisson")
  summary(mod5)
  
  c5<-mod5$coefficients[2]
  se5<-summary(mod5)$coefficients[,2][2]
  p5<-summary(mod5)$coefficients[,4][2]
  
  r5<-exp(mod5$coefficients[2])-1
  r5l<-exp(mod5$coefficients[2]-1.96*se5)-1
  r5u<-exp(mod5$coefficients[2]+1.96*se5)-1
  
  d5<-m$deaths[6]-m$deaths[6]/(1+r5)
  d5l<-m$deaths[6]-m$deaths[6]/(1+r5l)
  d5u<-m$deaths[6]-m$deaths[6]/(1+r5u)
  
  d18_19<-0.5*(m$deaths[2]+m$deaths[4])
  d20<-m$deaths[6]
  dcovhos<-m$deaths[6]-m$deaths[8]
  
  
  results<-rbind.fill(results,data.frame(agecat,agecat2,d18_19,d20,dcovhos,c1,se1,p1,r1,r1l,r1u,d1,d1l,d1u,
                             c2,se2,p2,r2,r2l,r2u,d2,d2l,d2u,
                             c3,se3,p3,r3,r3l,r3u,d3,d3l,d3u,
                             c4,se4,p4,r4,r4l,r4u,d4,d4l,d4u,
                             c5,se5,p5,r5,r5l,r5u,d5,d5l,d5u))
  row.names(results)<-NULL
  
  return(results)
  
}



dev.off()

  
  pdf(file="France_Mortality.pdf",paper="a4",width=8,height=11,family="Bookman")
  results<-NULL
  #par(oma=c(0,0,0,1))
  par(mar=c(5.1, 4.1, 2.1, 3.1))
  
      par(mfrow=c(4,2))
      
      
      #Legend
      
      plot(0:30, type = "n", xaxt="n", yaxt="n", bty="n", xlab = "", ylab = "")
      title("Covid-19 epidemic and mortality\n in France by age category")
      legend("bottomleft",title="Legend",c("Total deaths 2018","Total deaths 2019","Average deaths 2018 and 2019",
                         "Total deaths 2020 (provisional)",
                         "Total deaths 2020 without Covid-19 hospital deaths",
                         "Covid-19 hospital deaths (right scale)",
                         "Lockdown"),
      pch=mypch,lty=mylty,
      col=mycol,bty = "n",
      text.col=mycol,cex=0.9,title.col = "black")
      
      
      dc2<-dc
      dc2$cl_age90<-0
      results<-gdata(0,data=dc2)
      
      results<-gdata(90)
      
      results<-gdata(89)
      
      results<-gdata(79)
      
      results<-gdata(69)
      
      results<-gdata(59)
      
      results<-gdata(49)
      
      
      #################################################
      # Page 2
      par(mfrow=c(4,2))
      
      results<-gdata(39)
      
      results<-gdata(29)
      
      results<-gdata(19)
      
      results<-gdata(9)
      
      results$agecat3<-gsub("Age: ","",results$agecat2)
      
      #dev.off()
      
      
      bar1<-barplot(t(cbind(c(results$d3[1],results$d1[1]),
                            c(results$d5[1],results$d5[1]),
                            c(results$d4[1],results$d2[1]),
                            c(nursing_home_deaths,nursing_home_deaths))),
                    beside=TRUE,
                    names.arg=c("Simple difference","Difference-in-difference"),
                    col=c(adjustcolor("red",alpha.f=0.5),
                          adjustcolor("darkviolet",alpha.f=0.5),
                          adjustcolor("darkorange",alpha.f=0.5),
                          adjustcolor("darkgreen",alpha.f=0.5)),
                    ylim=c(0,max(1.2*results$d1u[1])),
                    las=1,xlab=list("Age: all",cex=0.8),
                    ylab=list("Number of deaths in excess",cex=0.8),cex.axis =0.8,cex.names = 0.8)
      abline(h=0)
      title("Excess deaths between March 4th and April 14th",cex.main=1)
      segments(bar1, t(cbind(c(results$d3l[1],results$d1l[1]),c(NA,NA),c(results$d4l[1],results$d2l[1]),c(NA,NA))),bar1,
               t(cbind(c(results$d3u[1],results$d1u[1]),c(NA,NA),c(results$d4u[1],results$d2u[1]),c(NA,NA))), lwd = 0.8,col="gray30")
      arrows(bar1, t(cbind(c(results$d3l[1],results$d1l[1]),c(NA,NA),c(results$d4l[1],results$d2l[1]),c(NA,NA))),bar1,
             t(cbind(c(results$d3u[1],results$d1u[1]),c(NA,NA),c(results$d4u[1],results$d2u[1]),c(NA,NA))), lwd = 0.8, angle = 90,
             code = 3, length = 0.05,col="gray30")
      text(c(1.7,3.7,6.7,8.7),c(results$d3[1],results$d4[1],results$d1[1],results$d2[1]),
                            paste("+",100*round(c(results$r3[1],results$r4[1],results$r1[1],results$r2[1]),2),"%")
           ,pos=1,cex=0.8)
      legend("topleft",  
             c("Excess deaths",
               "Covid-19 hospital deaths",
               "Excess deaths without\n Covid-19 hospital deaths",
               "Covid-19 nursing home deaths"), 
             fill=c(adjustcolor("red",alpha.f=0.5),
                    c(adjustcolor("darkviolet",alpha.f=0.5)),
                    adjustcolor("darkorange",alpha.f=0.5),
                    adjustcolor("darkgreen",alpha.f=0.5)),
             ncol=1, cex=0.7,bty="n",
             border=1,bg=NA)
      
      
      bar1<-barplot(t(cbind(results$d1[-1],results$d5[-1],results$d2[-1])),beside=TRUE,names.arg=results$agecat3[-1],
                    col=c(adjustcolor("red",alpha.f=0.5),
                          adjustcolor("darkviolet",alpha.f=0.5),
                          adjustcolor("darkorange",alpha.f=0.5)),
                    ylim=c(1.01*min(results$d2l[-1]),max(1.01*results$d1u[-1])),las=2,xlab=list("Age",cex=0.8),
                    ylab=list("Number of deaths in excess (diff-in-diff)",cex=0.8),cex.axis =0.8,cex.names = 0.8)
      abline(h=0)
      title("Excess deaths between March 4th and April 14th",cex=0.7)
      segments(bar1, t(cbind(results$d1l[-1],rep(NA,nrow(results)-1),results$d2l[-1])),bar1,
               t(cbind(results$d1u[-1],rep(NA,nrow(results)-1),results$d2u[-1])), lwd = 0.8,col="gray50")
      arrows(bar1, t(cbind(results$d1l[-1],rep(NA,nrow(results)-1),results$d2l[-1])),bar1,
             t(cbind(results$d1u[-1],rep(NA,nrow(results)-1),results$d2u[-1])), lwd = 0.8,col="gray50", angle = 90,
             code = 3, length = 0.01)
      legend("topright", 
             c("Excess deaths","Covid-19 hospital deaths","Excess deaths without\n Covid-19 hospital deaths"), 
             fill=c(adjustcolor("red",alpha.f=0.5),
                    adjustcolor("darkviolet",alpha.f=0.5),
                    adjustcolor("darkorange",alpha.f=0.5)), ncol=1, cex=0.8,bty="n",
             border=1,bg=NA)
      
      bar1<-barplot(t(cbind(results$r1[-1],results$r1[-1]-results$r2[-1],results$r2[-1])),beside=TRUE,names.arg=results$agecat3[-1],
                    col=c(adjustcolor("red",alpha.f=0.5),
                          adjustcolor("darkviolet",alpha.f=0.5),
                          adjustcolor("darkorange",alpha.f=0.5)),
                    ylim=c(1.3*min(results$r2l[-1]),max(1.3*results$r1u[-1])),las=2,xlab=list("Age",cex=0.8),
                    ylab=list("Increase in deaths (diff-in-diff)",cex=0.8),cex.axis =0.8,cex.names = 0.8,yaxt="n")
      abline(h=0)
      axis(2,las=2,at=pretty(t(cbind(results$r1[-1],results$r2[-1]))), 
                             lab=paste(pretty(t(cbind(results$r1[-1],results$r2[-1])) * 100),"%"),cex.axis=0.8)
      title("Increase in deaths between March 4th and April 14th",cex.main=0.9)
      segments(bar1, t(cbind(results$r1l[-1],rep(NA,nrow(results)-1),results$r2l[-1])),bar1,
               t(cbind(results$r1u[-1],rep(NA,nrow(results)-1),results$r2u[-1])), lwd = 0.8,col="gray50")
      arrows(bar1, t(cbind(results$r1l[-1],rep(NA,nrow(results)-1),results$r2l[-1])),bar1,
             t(cbind(results$r1u[-1],rep(NA,nrow(results)-1),results$r2u[-1])), lwd = 0.8,col="gray50", angle = 90,
             code = 3, length = 0.01)
      legend("topright", 
             c("Increase in deaths","Increase due to Covid-19 hospital deaths",
               "Increase in deaths without\n Covid-19 hospital deaths"), 
             fill=c(adjustcolor("red",alpha.f=0.5),
                    adjustcolor("darkviolet",alpha.f=0.5),
                    adjustcolor("darkorange",alpha.f=0.5)), ncol=1, cex=0.8,bty="n",
             border=1,bg=NA)
      
          plot(0:30, type = "n", xaxt="n", yaxt="n", bty="n", xlab = "", ylab = "")
      title("Sources")
      text(1, 23, "Deaths in France",pos=4,cex=1)
      text(1, 21, "Insee",pos=4,cex=0.8)
      text(1, 19, "https://www.insee.fr/fr/information/4470857",pos=4,cex=0.6)
      text(1, 17, "File: Fichier individuel comportant des informations sur chaque décès",pos=4,cex=0.6)
      
      text(1, 13, "Covid deaths per age and region",pos=4,cex=1)
      text(1, 11, "Santé Public France - data.gouv.fr",pos=4,cex=0.8)
      text(1, 9, "https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/",pos=4,cex=0.6)
      text(1, 7, "File: donnees-hospitalieres-etablissements-covid19-2020-##-##.csv",pos=4,cex=0.6)
  
      text(1, 3, "R Script",pos=4,cex=1)
      text(1, 1, "http://olivier.godechot.free.fr/hopfichiers/French_Overmortality.R",pos=4,cex=0.6)
  
  dev.off()