#########################################################################################################

#                                         European DataLab
#                         La confiance dans la police en progression depuis 10 ans

#########################################################################################################
#########################################################################################################
####David Marguerit
###Mars 2021



#=========================================
# Set direction & manage packages
#=========================================

  #---------------
  # Set direction
  #---------------

setwd("...")


  #---------------
  # Manage packages
  #---------------

    #......................
    # Updates library
    #......................

update.packages(ask = FALSE)

    #......................
    # Load library
    #......................

library (plm)
library (openxlsx)
library (eurostat)
library (survey)
library (haven)
library (lubridate)
library (tidyverse)


#=========================================
# Arrange contextual variables
#=========================================

  #--------------
  # Eurostat data
  #--------------

    #...............................
    # Police officer par inhabitant
    #...............................

# Police officers: 2008-2018
polof <- get_eurostat ("crim_just_job")
polof <- filter (polof, sex == "T", isco08 =="OC5412", unit == "NR")
polof$year <- year (polof$time)
polof <- select (polof, iso2 = geo, police = values, year)


# Police officers: 1993-2007
polof2 <- get_eurostat ("crim_plce")
polof2$year <- year (polof2$time)
polof2 <- select (polof2, iso2 = geo, police = values, year)
context <- bind_rows (polof, polof2)
rm (polof2, polof)

# Population
pop <- get_eurostat ("demo_pjan")
pop <- filter (pop, sex == "T", age =="TOTAL", unit == "NR")
pop$year <- year (pop$time)
pop <- select (pop, iso2 = geo, pop = values, year)

# Create indicators
context <- full_join (context, pop, by = c('iso2', 'year'))
rm (pop)

context$police.inhab <- context$police / context$pop * 1000


    #...................................
    # Public spending by police officer
    #...................................

# Public spending
cofog <- get_eurostat ("gov_10a_exp")
cofog <- filter (cofog, sector == "S13", unit =="MIO_NAC" | unit =="PC_GDP", na_item == "TE", cofog99 == "GF0301")
cofog$year <- year (cofog$time)
cofog <- select (cofog, iso2 = geo, values, year, unit)
cofog <- pivot_wider (cofog, names_from = "unit", values_from = "values")
cofog <- rename (cofog, spend = "MIO_NAC", spend.gdp = "PC_GDP")

# HICP
hicp <- get_eurostat ("prc_hicp_aind")
hicp <- filter (hicp, coicop == "CP00", unit =="INX_A_AVG")
hicp$year <- year (hicp$time)
hicp <- select (hicp, iso2 = geo, hicp = values, year)

# PPP
ppp <- get_eurostat ("prc_ppp_ind")
ppp$year <- year (ppp$time)
ppp <- filter (ppp, ppp_cat == "GDP", na_item =="PPP_EU27_2020", year == 2015)
ppp <- select (ppp, iso2 = geo, ppp = values)

# Create indicators
context <- full_join (context, cofog, by = c('iso2', 'year'))
context$spend.police <- (context$spend * 1000000) / context$police
rm (cofog)

context <- full_join (context, ppp, by = c('iso2'))
context$spend.police.ppp <- context$spend.police / context$ppp
rm (ppp)

context <- full_join (context, hicp, by = c('iso2', 'year'))
context$spend.police.ppp <- (context$spend.police.ppp / context$hicp) * 100
rm (hicp)


    #...................................
    # Number of crime
    #...................................

# Crime 1993-2007
crime <- get_eurostat ("crim_gen")
crime <- filter (crime, iccs == "ICCS0101" | iccs == "ICCS0401" | iccs == "ICCS05012")
crime$year <- year (crime$time)
crime <- select (crime, geo, values, iccs,  year)

# Police officers: 2008-2018
crime2 <- get_eurostat ("crim_off_cat")
crime2 <- filter (crime2, iccs == "ICCS0101" | iccs == "ICCS0401" | iccs == "ICCS05012", unit == "NR")
crime2$year <- year (crime2$time)
crime2 <- select (crime2, geo, values, iccs,  year)

# Create indicators
crime <- bind_rows (crime, crime2)
rm (crime2)
crime <- pivot_wider (crime, names_from = "iccs", values_from = "values")
crime <- rename (crime, iso2 = 'geo', homicide = 'ICCS0101', robbery = 'ICCS0401', burglary = 'ICCS05012')
context <- full_join (context, crime, by = c('iso2', 'year'))

context$homicide.inhab <- context$homicide / context$pop * 1000
context$robbery.inhab <- context$robbery / context$pop * 1000
context$burglary.inhab <- context$burglary / context$pop * 1000

rm (crime)

  #----------------------------------------------
  # Corruption (from Transparency International)
  # Dcwnload data from: https://www.transparency.org/en/cpi/2020/index/nzl
  #----------------------------------------------

#2012-2020
cpi <- read.xlsx ("CPI2020_GlobalTablesTS_210125.xlsx", startRow = 3, sheet = 'CPI Timeseries 2012 - 2020')
cpi <- select (cpi, iso3 = ISO3, cpi2020 = CPI.score.2020, cpi2019 = CPI.score.2019, cpi2018 = CPI.score.2018, cpi2017 = CPI.score.2017,
               cpi2016 = CPI.score.2016, cpi2015 = CPI.score.2015, cpi2014 = CPI.score.2014, cpi2013 = CPI.Score.2013,
               cpi2012 = CPI.Score.2012)

# 2010
cpi10 <- read_delim ("CPI-2010-new_200601_105629.csv", delim = ",")
cpi10 <- select (cpi10, iso3 = iso, cpi2010 = score)

# 2008
cpi08 <- read_delim ("CPI-Archive-2008-2.csv", delim = ",")
cpi08 <- select (cpi08, iso3 = iso, cpi2008 = score)

# 2006
cpi06 <- read_delim ("CPI-2006-new_200602_095933.csv", delim = ",")
cpi06 <- select (cpi06, iso3 = iso, cpi2006 = score)

# 2004
cpi04 <- read_delim ("CPI-2004_200602_110140.csv", delim = ",")
cpi04 <- select (cpi04, iso3 = iso, cpi2004 = score)

# 2002
cpi02 <- read_delim ("CPI-2002_200602_115328.csv", delim = ",")
cpi02 <- select (cpi02, iso3 = iso, cpi2002 = score)

# Create indicators
cpi <- full_join (cpi, cpi10, by = "iso3")
cpi <- full_join (cpi, cpi08, by = "iso3")
cpi <- full_join (cpi, cpi06, by = "iso3")
cpi <- full_join (cpi, cpi04, by = "iso3")
cpi <- full_join (cpi, cpi02, by = "iso3")
rm (cpi10, cpi08, cpi06, cpi04, cpi02)

cpi <- pivot_longer (cpi, !iso3, names_to = "year", values_to = "cpi")
cpi$year <- substr (cpi$year, 4, 7)
cpi$year <- as.numeric (cpi$year)

cpi <- cpi %>%
  group_by (year) %>%
  mutate (cpi.rank = min_rank (1/cpi)) %>%
  ungroup ()


  #--------------
  # merge data
  #--------------

country <- read_delim("https://wp.me/aa3WWD-dD", ";", escape_double = FALSE, 
                            col_types = cols(edl = col_number(), euro = col_number(),
                                             numeric = col_number(), oecd = col_number(), 
                                             ue = col_number()), locale = locale(encoding = "WINDOWS-1252"), 
                            trim_ws = TRUE)
country <- country %>%
  filter (edl == 1) %>%
  select (name_fr, iso2, iso3)


context$iso2[context$iso2 == "EL"] <- "GR"
context <- left_join (country, context, by = "iso2")

context <- left_join (context, cpi, by = c('year', 'iso3'))

context <- select (context, -iso3)
rm (cpi, country)


#=========================================
# Arrange European Social survey
# Download data from: https://www.europeansocialsurvey.org/
#=========================================

  #---------------
  # Load data
  #---------------

ess <- read_stata ("ESS1-9e01_1.dta")

  #-------------------
  # Arrange variables
  #-------------------

    #........
    # Survey
    #........

# Year
ess$year <- NA
ess$year [ess$essround == 1] <- 2002
ess$year [ess$essround == 2] <- 2004
ess$year [ess$essround == 3] <- 2006
ess$year [ess$essround == 4] <- 2008
ess$year [ess$essround == 5] <- 2010
ess$year [ess$essround == 6] <- 2012
ess$year [ess$essround == 7] <- 2014
ess$year [ess$essround == 8] <- 2016
ess$year [ess$essround == 9] <- 2018

# Weights
ess$wgt <- ess$anweight



    #................................................
    # confiance police, aggression & discrimination
    #................................................

# Confiance police
ess$police <- ess$trstplc

# Respondent or household member victim of burglary/assault last 5 years
ess$agression <- NA
ess$agression [ess$crmvct == 1] <- 1
ess$agression [ess$crmvct == 2] <- 2
ess$agression <- factor(ess$agression,
                     levels = c(1, 2),
                     labels = c("Agressé", "Pas agressé"))

# Feeling of safety of walking alone in local area after dark
ess$securite <- NA
ess$securite [ess$aesfdrk == 1] <- 1
ess$securite [ess$aesfdrk == 2] <- 2
ess$securite [ess$aesfdrk == 3] <- 3
ess$securite [ess$aesfdrk == 4] <- 4
ess$securite <- factor(ess$securite,
                        levels = c(1, 2, 3, 4),
                        labels = c("Très sur", "Sur", "peu sur", "pas sur"))
 
# Discrimination of respondent's group: colour or race
ess$disc.race <- NA
ess$disc.race [ess$dscrrce == 0] <- 2
ess$disc.race [ess$dscrrce == 1] <- 1
ess$disc.race <- factor(ess$disc.race,
                       levels = c(1, 2),
                       labels = c("Oui", "Non"))

# Discrimination of respondent's group: nationality
ess$disc.nat <- NA
ess$disc.nat [ess$dscrntn == 0] <- 2
ess$disc.nat [ess$dscrntn == 1] <- 1
ess$disc.nat <- factor(ess$disc.nat,
                        levels = c(1, 2),
                        labels = c("Oui", "Non"))

# Discrimination of respondent's group: religion
ess$disc.rel <- NA
ess$disc.rel [ess$dscrrlg == 0] <- 2
ess$disc.rel [ess$dscrrlg == 1] <- 1
ess$disc.rel <- factor(ess$disc.rel,
                       levels = c(1, 2),
                       labels = c("Oui", "Non"))



    #...........................
    # characteristics sociodemo
    #...........................

# Age
ess <- ess %>%
  mutate (age = case_when(
    agea < 25 ~ 1,
    agea >= 25 & agea < 65 ~ 2,
    agea >= 65 ~ 3))
ess$age <- factor(ess$age,
                     levels = c(1, 2, 3),
                     labels = c("Moins de 25 ans", "25-64", "65 ans et plus"))


# Sex
ess$sex <- NA
ess$sex[ess$gndr == 1 ] <- 1
ess$sex[ess$gndr == 2] <- 2
ess$sex <- factor(ess$sex,
                  levels = c(1, 2),
                  labels = c("Homme", "Femme"))


# Education
ess$edu <- NA
ess$edu[ess$edulvla == 1 | ess$edulvla == 2] <- 1
ess$edu[ess$edulvla == 3 | ess$edulvla == 4] <- 2
ess$edu[ess$edulvla == 5] <- 3
ess$edu <- factor(ess$edu,
                    levels = c(1, 2, 3),
                    labels = c("Faible", "Moyen", "Eleve"))


# Place of birth
ess$birth <- NA
ess$birth[ess$brncntr == 1] <- 1
ess$birth[ess$brncntr == 2] <- 2
ess$birth <- factor(ess$birth,
                  levels = c(1, 2),
                  labels = c("Autochtone", "Allochtone"))


# Belong to minority ethnic group in country
ess$minorite <- NA
ess$minorite[ess$blgetmg == 2] <- 1
ess$minorite[ess$blgetmg == 1] <- 2
ess$minorite <- factor(ess$minorite,
                    levels = c(1, 2),
                    labels = c("Non", "Oui"))


# couple (missing value for 2018)
ess$couple <- NA
ess$couple [ess$partner == 1] <- 1
ess$couple [ess$partner == 2] <- 2
ess$couple <- factor(ess$couple,
                         levels = c(1, 2),
                         labels = c("En couple", "Célibataire"))

# Location
ess$location <- NA
ess$location[ess$domicil == 1] <- 1
ess$location[ess$domicil == 2] <- 2
ess$location[ess$domicil == 3] <- 3
ess$location[ess$domicil == 4 | ess$domicil == 5] <- 4
ess$location <- factor(ess$location,
                    levels = c(1, 2, 3, 4),
                    labels = c("Grande ville", "Banlieue", "Ville ou petite ville", "campagne"))

# Working status
ess$wstatus <- NA
ess$wstatus[ess$mnactic == 1] <- 1
ess$wstatus[ess$mnactic == 2] <- 2
ess$wstatus[ess$mnactic == 3 | ess$mnactic == 4] <- 3
ess$wstatus[ess$mnactic == 5 | ess$mnactic == 6 | ess$mnactic == 8] <- 4
ess$wstatus[ess$mnactic == 7 | ess$mnactic == 9] <- 5
ess$wstatus <- factor(ess$wstatus,
                       levels = c(1, 2, 3, 4, 5),
                       labels = c("Travailleur", "Etudiant", "Chômeur", "Inactif", "Autre"))


# Feeling about household's income nowadays
ess$income <- NA
ess$income[ess$hincfel == 1] <- 1
ess$income[ess$hincfel == 2] <- 2
ess$income[ess$hincfel == 3] <- 3
ess$income[ess$hincfel == 4] <- 4
ess$income <- factor(ess$income,
                      levels = c(1, 2, 3, 4),
                      labels = c("Confortable", "Ca va", "Difficile", "Trèse difficile"))



  #-------------------
  # Add design effect
  #-------------------

def <- read_stata ("ESS9e03_1.dta")

def <- select (def, cntry, idno, stratum, psu)

def$year <- 2018

ess <- left_join (ess, def, by = c('year', 'cntry', 'idno'))

rm (def)


#=================================
# Select countries and variables
#=================================


ess <- select (ess, year, wgt, iso2 = cntry, idno, police, agression, securite, disc.race, disc.nat, disc.rel, age, sex, edu, 
               birth, minorite, location, wstatus, income, stratum, psu)


country <- read_delim("https://wp.me/aa3WWD-dD", ";", escape_double = FALSE, 
                      col_types = cols(edl = col_number(), euro = col_number(),
                                       numeric = col_number(), oecd = col_number(), 
                                       ue = col_number()), locale = locale(encoding = "WINDOWS-1252"), 
                      trim_ws = TRUE)
country <- country %>%
  filter (edl == 1) %>%
  select (country = name_fr, iso2)

ess <- inner_join (ess, country, by = ('iso2'))

rm (country)


#=========================================
# Graphique 1 : Evolution confiance
#=========================================

df <- ess %>%
  select (wgt, country, police, year) %>%
  filter (year == 2008 | year == 2018) %>%
  drop_na ()

design <- svydesign (id=~1, strata=~country, weights=~wgt, data=df)

res <- svyby (~police, ~year+country, design, svymean)

res <- res %>%
  select (-se) %>%
  pivot_wider (names_from = year, values_from = police)



#=========================================
# Graphique 2 : characteristics
#=========================================

df <- ess %>%
  filter (year == 2018) %>%
  drop_na ()

design <- svydesign (id=~psu, strata=~stratum, weights=~wgt, data=df)
options(survey.lonely.psu = "certainty")

model <- svyglm (police ~ age + sex + edu + birth + location + wstatus + income +
                  agression + securite + minorite +
                  country, design=design)

summary(model)
rm (df)

#=========================================
# Graphique 3 : contextual variables
#=========================================

df <- ess %>%
  select (wgt, iso2, police, year) %>%
  drop_na ()

design <- svydesign (id=~1, strata=~iso2, weights=~wgt, data=df)

df <- svyby (~police, ~year+iso2, design, svymean)

df <- select (df, iso2, year, trust = police)

df <- inner_join (df, context, by = c("iso2", "year"))

fixed.time <- plm (trust ~ log(police.inhab) + log(spend.gdp) + log(homicide.inhab) + cpi.rank + factor (year), index=c("name_fr", "year"), model="within", weights = pop, data = df)
summary (fixed.time)