{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Assessing sustainability ambition by response completeness"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"library(tidyverse)\nlibrary(data.table)\nlibrary(FactoMineR)\nlibrary(factoextra)\nlibrary(plotly)\nlibrary(leaflet)\nlibrary(sf)\nlibrary(IRdisplay)\nlibrary(htmltools)\nlibrary(gridExtra)\nlibrary(repr)\nlibrary(scales)\nlibrary(spdep)\nlibrary(ineq)\nlibrary(grid)\nlibrary(gtable) \nlibrary(openxlsx)\n\noptions(warn=-1)\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"# Problem statement\nDevelop a methodology for calculating key performance indicators (KPIs) that relate to the environmental and social issues that are discussed in the CDP survey data. Leverage external data sources and thoroughly discuss the intersection between environmental issues and social issues. Mine information to create automated insight generation demonstrating whether city and corporate ambitions take these factors into account.\n\n\n\n# Summary\nOne property of CDP's questionnaires is that they are very detailed and factual about sustainability development. In order to use all the data and quantify the sustainability ambition of cities or corporates, I use a systematic and flexible indicator based on the completeness of the responses provided to CDP. After describing sustainability ambition KPI and its validation, I use it to discuss the intersection between environmental and social issues in the United States. Main finding is that high ambition cities meet more environmental injustice in intensity and diversity. Also, we show evidence that spatial arrangement of pollution sources and poverty in a city may play an important role in environmental ambition. \nIn the last part, I explore corporates and cities KPI around the water security topic. It reveals that ambitious corporates adapt their activities according the vulnerability level of the population at facilities locations. Some synergy between cities and corporates sustainability ambitions is also highlighted.\n\n# My understanding of the problem\nMeasuring greenhouse gas emissions is difficult and tedious. It requires carrying out a complete and precise inventory, taking into account all the parameters of the operation from a city or a business. These measures are often expensive and require a strong motivation from leaders to implement them. On the other hand, CDP’s questionnaire represents a unique and massive amount of information about cities and corporates actions in sustainability development. In order to exploit and value CDP data, I decided to take an analytical point of view by creating an original indicator based on responses completeness : the **sustainability ambition KPI**. Given the complex nature of the environmental challenge, and the exhaustivity of the questionnaire, this methodology perfectly fits the famous formula: **you can't manage what you don't measure**. \nFirst of all, it is very important to distinguish **ambition** from sustainability **compliance** (i.e. a level of sustainability already achieved). There are KPIs to assess the level of sustainability like United Nations's [SDG index](https://www.sdgindex.org/) or readiness to adapt like [Notre Dame Global Adaptation Initiative](https://gain.nd.edu/). Here, my approach consists in evaluating the level of **effort** implemented by the CDP's respondents. This appreciation can make it easier to identify opportunities for collaboration. It makes it possible to identify the institutions needing support towards sustainable development.\nThis systematic approach allow the use of the questionnaire in all its relevance, by making the lack of response a key element of the method.\nThis indicator is then compared with detailed and complete data relating to environmental and social issues, [EJscreen](https://ejscreen.epa.gov/mapper/) database from EPA.\n\n\n# Table of contents\n1. **Key Performance Indicator : Assessing sustainability ambition by response completeness**\n * Missing data \n * Methodology\n * Cities KPI calculation\n * Cities KPI validation \n * *Official CDP score*\n * *Mayors adhesion within \"Global Covenant of Mayors for Climate & Energy\"*\n * *Responses to question 1.0*\n * KPI evolution between 2018 and 2019 \n * Corporates KPI calculation\n * Corporate KPI validation\n2. **Environmental issues and social inequity : an inextricable link**\n * Characterization of vulnerable population \n * Introducing cities KPI into social - environmental interactions\n * *Population*\n * *Vulnerability*\n * *Pollution intensity*\n * *Environmental Justice index*\n * *Cumulated Environmental Justice Indexes*\n * Focus on Los Angeles city \n3. **Water security: Exploration around cities and corporates KPI** \n * Does corporates care about vulnerable population ?\n * Is there a synergy between cities and corporates about sustainability ambition ?\n4. **Possible improvements & discussion**\n\n\n\n# 1. Key Performance Indicator : Assessing sustainability ambition by response completeness\n## a. Missing values\n\nAny statistical analysis or machine learning project is based on the quantity and quality of available data. Thus, the presence of missing values is an inevitable source of bias due to imputation or sample size reduction.\n \nIn the following figure, I sorted the numerical questions by answer completeness in decreasing order. Bars represent the rate of missing values per question. The curve shows the rate of the sample that has at least one missing value among the set of questions."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\n# read city responses\n\n#### 2019\ndis <- fread(\"../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2019_Cities_Disclosing_to_CDP.csv\",encoding = \"UTF-8\")\ndf <- fread(\"../input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2019_Full_Cities_Dataset.csv\", encoding = \"UTF-8\", stringsAsFactors = FALSE)\ncolnames(df) <- gsub(\" \",\"_\", tolower(colnames(df)))\ndf$column_number <- as.character(df$column_number)\ndf$column_number <- ifelse(nchar(df$column_number) == 1, paste0(\"0\",df$column_number), df$column_number)\ndf$column <- paste(\"col\",df$column_number,df$column_name,sep=\"_\")\ndf$row <- paste(\"row\",df$row_number, df$row_name,sep= \"_\")\n\nliste <- split(df,df$question_number)\nliste <- lapply(liste,function(x){\n suppressMessages(dcast(x, account_number +country + organization + question_number + row ~ column,value.var = \"response_answer\",fun.aggregate=function(x) paste(x, collapse = \";\")))\n})\n# Donnees numeriques\nres <- data.frame(id = unique(df$organization))\n\nnums <- lapply(liste, function(x){\n x <- as.data.frame(x)\n \n cols <- c()\n for (i in 6:ncol(x)){\n x[x == \"\"] <- NA\n orig_na <- sum(is.na(x[,i]))\n if (orig_na == nrow(x)) next\n if (sum(is.na(as.numeric(x[,i]))) == orig_na) cols <- c(cols, colnames(x)[i])\n }\n \n return(cols)\n})\n\nliste2 <- list()\nfor (i in 1:length(nums)){\n x <- nums[[i]]\n if (length(x) == 0) next\n \n temp <- liste[[i]][,c(\"organization\",\"row\",x),with=F]\n temp <- melt(temp, id.vars = c(\"organization\",\"row\"))\n temp <- dcast(temp, organization ~ variable + row)\n colnames(temp)[-1] <- paste0(\"Q\",names(liste)[i],\"_\",colnames(temp)[-1])\n liste2[[names(nums)[i]]] <- temp\n \n}\n\nres <- Reduce(function(dtf1, dtf2) merge(dtf1, dtf2, by = \"organization\", all.x = TRUE),\n liste2)\nfor (i in 2:ncol(res)) res[[i]] <- as.numeric(res[[i]])\n\n\n# courbe cumulée donnees manquantes\nnas <- data.frame(missingness=colSums(is.na(res[,-1]))/nrow(res))\nnas$question <- rownames(nas)\nnas <- nas[order(nas$missingness ),]\n\nnas <- nas[1:100,]\n\nnas$cumul <- NA\nfor (i in 1:nrow(nas)){\n cols <- c(nas[1:i,\"question\"])\n nas[i,\"cumul\"] <- 814 - sum(complete.cases(res[,cols,with=FALSE]))\n}\n\nnas$cumul <- nas$cumul/814\nnas <- nas[1:20,]\nnas$`Data point` <- paste(sapply(strsplit(nas$question, \"_\"), function(x) x[1]),\n sapply(strsplit(nas$question, \"_\"), function(x) x[4]), sep =\" - \")\nnas$question <- factor(nas$question, levels = nas$question) \n\n#Plot \np <- ggplot(nas, aes(x=as.numeric(question), label=`Data point`))+\n geom_bar(aes( y=missingness, fill=\"Missing values \\nrate\"),stat=\"identity\", binwidth=.5)+\n geom_line(aes(y=cumul, group=1, linetype=\"Cumulative rate \\nof uncomplete \\nsample\"), colour=\"red\")+\n geom_point(aes(y=cumul, group=1), colour=\"red\")+\n theme_bw()+ \n labs(x=\"Number of questions\",y=\"Sample rate\",fill =\"\", linetype=\"\")\n \npltly <- ggplotly(p, tooltip=c(\"Data point\")) \n\nhtmlwidgets::saveWidget(pltly, \"plot00.html\")\nIRdisplay::display_html('')\n ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If we just try to extract and exploit numerical responses, more than 75% of the sample become incomplete after 10 questions. "},{"metadata":{},"cell_type":"markdown","source":"## b. Methodology\n\nSustainability ambition KPI is based on the following assumption : \n\n**The more an entity (city or company) knows and provides detailed figures on the environmental problem, the more this entity is committed and likely to manage these problems.** \n\nThanks to CDP questionnaires structure, a blank response means that the respondent doesn't respond. We assume that this also means that he does not know the answer. In other words, these missing values are easily detectable and there is no mismatch with not applicable questions. Then, sustainability ambition KPI is based on the combination of : \n* the completeness of the answers relative to **quantities** (population, GHG emissions, fuel mix, etc.)\n* **overall completeness** by section including all kind of questions (choices, free text, numerical) \n* **lexical analysis** of the free text answers, using a technical lexicon translated in 4 languages\n\nWith completeness defined as the rate of responses provided. \n\nAfter this first step, the variables greater than 1 are standardized by the maximum and then a **PCA** is done to reduce dimensionnality and noise. Sustainability ambition KPI is computed by summing all the principal components per individual.\n\nThe advantages of this method are:\n* Speed of calculation\n* Comparability of individuals\n* Flexibility to calculate KPIs oriented towards specific issues\n* Use of the full sample\n\nThus, if one seeks to evaluate the motivation and the commitment of a city/company only on specific subjects like GHG emission, energy mix or transport, the user can compute the KPI based on the exhaustiveness of the correspondent questions. \n\n\n#### Technical lexicon\nEnvironment related words have been scraped on three websites, corresponding to the three main environmental topics of CDP's questionnaires: \n* [Wikipedia's Climate Change Glossary](https://en.wikipedia.org/wiki/Glossary_of_climate_change)\n* [Wikipedia's Environmental science Glossary](https://en.wikipedia.org/wiki/Glossary_of_environmental_science)\n* [Glossary of Water Management Terms from the Canadian goverment](https://www.aadnc-aandc.gc.ca/eng/1100100028092/1100100028094)\n\nThe total lexicon is composed of about 1000 words, translated in French, Spanish and Portuguese. \nFive text-related features are extracted : \n* averaged number of characters per word\n* total number of words \n* rate of technical words used\n* number of technical words used \n* number of unique technical words used \n\nThe 10 first words of the lexicon : "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# List of technical words\nall_words <- read.csv(\"../input/various-datasets-for-cdp-competition/technical_words.csv\", stringsAsFactors = FALSE)\nall_words <- all_words$mot\ncat(all_words[1:10])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## c. Calculation for cities (2019)\n\n#### Features\nSo, we have created 99 raw features composed of : \n* 65 numerical questions completeness\n* 29 sections overall completeness\n* 5 text-related features \n\n\nHead of features"},{"metadata":{"trusted":true,"_kg_hide-output":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"\n# Sections and parent sections\nsections <- df[,c(\"parent_section\",\"section\")]\nsections <- sections[!duplicated(sections),]\n\n\n# Functions to extract features\nget_kpi_features <- function(df, words_by_section=FALSE){\n \n liste <- split(df,df$question_number)\n liste <- lapply(liste,function(x){\n suppressMessages(dcast(x, account_number +country + organization + question_number + row ~ column,value.var = \"response_answer\",fun.aggregate=function(x) paste(x, collapse = \";\")))\n })\n \n nums <- lapply(liste, function(x){\n x <- as.data.frame(x)\n \n cols <- c()\n for (i in 6:ncol(x)){\n x[x == \"\"] <- NA\n orig_na <- sum(is.na(x[,i]))\n if (orig_na == nrow(x)) next\n if (sum(is.na(as.numeric(x[,i]))) == orig_na) cols <- c(cols, colnames(x)[i])\n }\n \n return(cols)\n })\n # quantities precision\n temp <- df\n temp$qu_col <- paste(temp$question_number,temp$column,sep=\"_\")\n cols <- c()\n for (i in 1:length(nums)){\n if (length(nums[[i]]) == 0) next\n cols <- c(cols, paste(names(nums)[i], nums[[i]], sep=\"_\"))\n }\n \n temp <- temp[temp$qu_col %in% cols,]\n temp$response_answer <- as.numeric(temp$response_answer)\n temp <- temp[,.(num_prec = 1 - sum(is.na(response_answer))/.N),.(section,organization,question_number)]\n temp <- suppressMessages(dcast(temp, organization ~ section + question_number, fill=0))\n colnames(temp)[-1] <- paste0(\"numprec_\",colnames(temp)[-1])\n \n # etape 1 : overall response rate by section\n df2 <- df\n df2$response_answer <- ifelse(df2$response_answer == \"\",NA, df2$response_answer)\n df3 <- df2[,.(nas=sum(is.na(response_answer))/.N),.(organization)]\n df4 <- df2[,.(nas=1-sum(is.na(response_answer))/.N),.(organization,section)]\n df4 <- suppressMessages(dcast(df4, organization ~ section, fill=0))\n colnames(df4)[-1] <- paste(\"comp\",colnames(df4)[-1],sep=\"_\")\n \n # technical lexicon\n test <- df[response_answer != \"\",.(count=length(unique(response_answer))/.N),.(question_number,column_name,row_number)]\n open_text <- paste(test$question_number, test$column_name, test$row_number)\n open_text <- open_text[test$count > 0.75]\n liste2 <- split(df, df$organization)\n \n if (words_by_section){\n env <- lapply(liste2,function(x){\n liste3 <- split(x, x$section)\n env_sect <- lapply(liste3, function(x){\n x$id <- paste(x$question_number, x$column_name, x$row_number)\n x <- x[x$id %in% open_text,]\n x$response_answer <- gsub(\"[[:punct:]]+\", \"\", x$response_answer)\n words <- unlist(strsplit(x$response_answer, \" \"))\n nwords <- length(words)\n \n res <- c()\n for (i in c(\"all_words\")){\n n <- sum(words %in% get(i))\n uniq <- sum(get(i) %in% words)\n rate <- n/nwords\n word_len <- mean(nchar(words))\n \n temp <- c(n,rate)\n names(temp) <- c(\"n\",\"rate\")\n res <- c(res,temp)\n }\n \n \n return(c(res,\"nwords\"= nwords))\n })\n for (i in 1:length(env_sect)) names(env_sect[[i]]) <- paste(names(env_sect)[i],names(env_sect[[i]]),sep=\"_\")\n env_sect <- as.data.frame(t(do.call(c,env_sect)))\n \n \n })\n \n \n env <- data.frame(organization=names(liste2),rbindlist(env, use.names = TRUE, fill=TRUE))\n colnames(env)[-1] <- paste(\"wds\",colnames(env)[-1], sep=\"_\")\n for (i in 2:ncol(env)){\n env[[i]][is.na(env[[i]])] <-0\n env[[i]][is.nan(env[[i]])] <-0\n env[[i]][is.infinite(env[[i]])] <-0\n if(max(env[[i]]) > 1) env[[i]] <- env[[i]]/max(env[[i]])\n } \n } else{\n env <- lapply(liste2,function(x){\n x$id <- paste(x$question_number, x$column_name, x$row_number)\n x <- x[x$id %in% open_text,]\n x$response_answer <- gsub(\"[[:punct:]]+\", \"\", x$response_answer)\n words <- unlist(strsplit(x$response_answer, \" \"))\n nwords <- length(words)\n \n res <- c()\n for (i in c(\"all_words\")){\n n <- sum(words %in% get(i))\n uniq <- sum(get(i) %in% words)\n rate <- n/nwords\n word_len <- mean(nchar(words))\n \n temp <- c(n,uniq,rate,word_len)\n names(temp) <- paste(i,c(\"n\",\"uniq\",\"rate\",\"word_len\"))\n res <- c(res,temp)\n }\n \n \n return(c(res,\"nwords\"= nwords))\n })\n \n \n env <- data.frame(organization=names(liste2),do.call(rbind,env))\n colnames(env)[-1] <- paste(\"wds\",colnames(env)[-1], sep=\"_\")\n for (i in 2:ncol(env)){\n env[[i]][is.na(env[[i]])] <-0\n env[[i]][is.nan(env[[i]])] <-0\n env[[i]][is.infinite(env[[i]])] <-0\n if(max(env[[i]]) > 1) env[[i]] <- env[[i]]/max(env[[i]])\n }\n }\n\n \n #final data frame\n \n res <- Reduce(function(df1,df2) merge(df1,df2,by=\"organization\",all=TRUE),list(temp,df4,env))\n return(as.data.frame(res))\n}\n\n# Functions to aggregate into KPIs\nextract_kpi <- function(x, section=\"all\", plot=TRUE){\n if (section == \"all\"){\n pca <- PCA(x[,-1], graph = plot)\n \n kpi <- rowSums(get_pca_ind(pca)$coord, na.rm=TRUE)\n \n } else{\n cols <- colnames(x)[grepl(paste(section,collapse = \"|\"), colnames(x))]\n pca <- PCA(x[,cols], graph = plot)\n kpi <- rowSums(get_pca_ind(pca)$coord, na.rm=TRUE)\n }\n \n return(kpi)\n}\n\n# Extract features for city responses\nfeat <- get_kpi_features(df, words_by_section = FALSE)\n\n# One global KPI\nkpi <- extract_kpi(feat, plot=FALSE)\n\n# One KPI per parent section\nkpis <- lapply(unique(sections$parent_section),function(x){\n extract_kpi(feat,section = sections$section[sections$parent_section == x], plot=FALSE)\n})\n\nkpis <- as.data.frame(do.call(cbind,kpis))\ncolnames(kpis) <- gsub(\" \",\"_\",tolower(unique(sections$parent_section)))\ncolnames(kpis)[1] <- \"none\"\nkpis$org <- feat$organization\nkpis$total <- kpi\n\n\nhead(feat[,c(grepl(\"numprec\",colnames(feat)))], 10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n"},{"metadata":{},"cell_type":"markdown","source":"#### PCA components\nTo understand the KPI calculation process, we can visualize PCA result before summing all components. "},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"pca <- PCA(feat[,-1], graph=FALSE)\npca_ind <- data.frame(org=feat$organization,\n KPI=kpi,\n get_pca_ind(pca)$coord)\n\np <- ggplot(pca_ind, aes(x= Dim.1, y= Dim.2, colour=kpi ,group=org))+ \n geom_hline(yintercept=0)+\n geom_vline(xintercept=0)+\n geom_point()+\n scale_colour_gradient2(low=\"#F8766D\",high=\"#7CAE00\",mid=\"#fffcbf\")+\n theme_bw()+\n labs(x=paste(\"Dimension 1 :\", round(pca$eig[1,\"percentage of variance\"], digits=2),\"% of variance explained\"),\n y=paste(\"Dimension 2 :\", round(pca$eig[2,\"percentage of variance\"], digits=2), \"% of variance explained\"),\n colour= \"Ambition KPI\")\npltly <- ggplotly(p, tooltip=c(\"org\",\"kpi\"))\n\nhtmlwidgets::saveWidget(pltly, \"plot1.html\")\nIRdisplay::display_html('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's obvious that some cities are far ahead from others in term of sustainability ambition. \n\nLet's plot a map displaying sustainability ambition KPI (a few points are misplaced, maybe due to a switch between latitude and longitude coordinates in disclosure data)"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"# Extract cities coordinates\ndat <- kpis\ndat$geo <- dis$`City Location`[match(dat$org, dis$Organization)]\ndat$geo <- gsub(\"POINT \",\"\",dat$geo)\ndat$geo <- gsub(\"\\\\(\",\"\",dat$geo)\ndat$geo <- gsub(\"\\\\)\",\"\",dat$geo)\ndat$long <- as.numeric(sapply(strsplit(dat$geo,\" \"), function(x) x[1]))\ndat$lat <- as.numeric(sapply(strsplit(dat$geo,\" \"), function(x) x[2]))\npts <- dat\npts <- st_as_sf(pts[!is.na(pts$long),],coords = c(\"long\",\"lat\"),crs=4326)\n \n# Plot interactive map \nmypal <- colorNumeric(c(\"#F8766D\",\"#fffcbf\",\"#7CAE00\"), domain = pts$total)\npts$label <- lapply(paste(\"City :\",pts$org,\"
KPI :\",round(pts$total, digits = 2)), HTML)\n\nm <- leaflet(options=list(\n center = c(20, 0),\n zoom = 2,\n worldCopyJump = FALSE,\n maxBounds = list(\n list(-90, -180),\n list(90, 180)))) %>%\n addProviderTiles(providers$Stamen.Toner) %>%\n addCircleMarkers(data=pts, \n fillColor = ~mypal(total),\n radius = 4,\n stroke = TRUE, \n fillOpacity = 1, \n weight = 1,\n color = \"black\",\n label = pts$label) %>%\n addLegend(\"bottomright\",\n pal = mypal,\n values = pts$total) \n\n\n \n\nhtmlwidgets::saveWidget(m, \"m.html\")\ndisplay_html('') ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This methodology allow us to focus on specific sections or questions. \nFor example, we can look City-wide emissions completeness vs Emissions Reductions section completeness."},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"p <- ggplot(kpis, aes(x= city_wide_emissions, y= emissions_reduction, group=org, colour=total))+\n geom_hline(yintercept = 0)+\n geom_vline(xintercept = 0)+\n geom_point(alpha=0.8)+\n scale_colour_gradient2(low=\"#F8766D\",high=\"#7CAE00\",mid=\"#fffcbf\")+\n theme_bw()+\n labs(colour=\"Ambition KPI\", x=\"City-wide emissions KPI\", y=\"Emissions Reductions KPI\")\n\n\n\npltly <- ggplotly(p, tooltip =c( \"org\",\"total\"))\n\nhtmlwidgets::saveWidget(pltly, \"plot2.html\")\nIRdisplay::display_html('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## d. Validation\n\nUnfortunately, there is no easily computable ground truth for this kind of problematic. But we can compare our results to other data sources relating to **sustainability ambition**. \n\n### i. Official CDP score\n\nFirst, we can simply compare our KPI with the list of cities that have obtained an A score by CDP. Dashed lines represents group mean. "},{"metadata":{"_kg_hide-input":true,"trusted":true,"scrolled":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\ncdp <- read.csv(\"../input/cdp-2019-a-list-cities/2019_A_list_Cities_with_Response_Links.csv\",stringsAsFactors = FALSE)\n\ntemp <- data.frame(kpi= c(kpis$total[!kpis$org %in% cdp$Organization], kpis$total[kpis$org %in% cdp$Organization]),\n group = c(rep(\"Others\",nrow(kpis) - nrow(cdp)), rep(\"A-score\", nrow(cdp))))\ntemp_mean <- as.data.table(temp)[,.(kpi = mean(kpi, na.rm=TRUE)),.(group)]\n\n# Plot\np <- ggplot(temp, aes(x=kpi, fill=group)) + \ngeom_density(alpha = 0.4) + \ngeom_vline(data=temp_mean, aes(xintercept=kpi, fill=group),\n linetype=\"dashed\")+\nxlim(c(-30,30)) + \ntheme_bw()+\ntheme(text = element_text(size=15))+\nlabs(fill=\"\", x= \"Ambition KPI\", y= \"Density\")\n\n\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is a clear gap between the two groups, which is pretty reassuring."},{"metadata":{},"cell_type":"markdown","source":"### ii. Mayors adhesion within \"Global Covenant of Mayors for Climate & Energy\"\n\n“Global Covenant of Mayors for Climate & Energy” is an international coalition aiming to gather mayors who want to fight against global warming by limiting greenhouse gas emissions. Mayors who are members of this association are awarded with several badges according to their level of commitment in the fight against global warming : \n* Mitigation inventory\n* Mitigation target\n* Mitigation plan\n* Adaptation assessment\n* Adaptation goal\n* Adaptation plan\n* Access to energy assessment\n* Access to energy plan\n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"mayors <- read.csv(\"../input/various-datasets-for-cdp-competition/cities_2020-11-10.csv\", stringsAsFactors = FALSE)\nmayors$city <- sapply(strsplit(mayors$name,\",\"), function(x) x[1])\nmayors$score <- rowSums(mayors[,7:15], na.rm = TRUE)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can sum up all badges received by each mayor and construct an ambition indicator. Here is the result : "},{"metadata":{"_kg_hide-input":true,"trusted":true,"scrolled":true},"cell_type":"code","source":"dat <- kpis\n# Format city name\ndat$city <- gsub(\"City of \",\"\", dat$org)\ndat$city <- gsub(\"Town of \",\"\", dat$city)\ndat$city <- sapply(strsplit(dat$city, \",\"), function(x) x[1])\n \n# Merge mayor's score to KPI \ndat$score <- -1 \nfor (i in 1:nrow(mayors)){\n index <- which(grepl(mayors[i,\"city\"], dat$city))\n dat$score[index] <- mayors$score[i]\n \n}\n\ncat(paste(sum(dat$score > -1), \"mayors found among\", nrow(dat), \"cities\")) \n \n# Cut KPI into categories then plot boxplot\noptions(repr.plot.width=10,repr.plot.height=6) \ndat$kpi_class <- cut(dat$total, quantile(dat$total,seq(0,1,0.2)), include.lowest=TRUE)\ndat$score_class <- \"No adhesion\"\ndat$score_class <- ifelse(dat$score == 0, \"No badge\", dat$score_class )\ndat$score_class <- ifelse(dat$score %in% 1:3, \"1 to 3 badges\", dat$score_class )\ndat$score_class <- ifelse(dat$score %in% 4:6, \"4 to 6 badges\", dat$score_class )\ndat$score_class <- factor(dat$score_class, levels=c(\"No adhesion\",\"No badge\",\"1 to 3 badges\",\"4 to 6 badges\"))\n\n\nggplot(dat, aes(x=score_class, y=kpi)) + \n geom_boxplot(outlier.shape = NA) +\n theme_bw() +\n theme(text = element_text(size = 15))+\n labs(x=\"Number of badges\", y= \"Ambition KPI\") ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pretty reassuring again, rewarded mayors have better sustainability ambition KPI than others. Interesting fact, cities with no adhesion at all to the Global Covenant of Mayors for Climate & Energy are not worst than member cities with no badge. \n### iii. Responses to question 1.0\nNaturally, we can compare our KPI to the answer given at question 1.0. \n\nQuestion : **Does your city incorporate sustainability goals and targets (e.g. GHG reductions) into the master planning for the city?**\n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true,"scrolled":true},"cell_type":"code","source":"temp <- df[df$question_number == \"1.0\",]\ndat <- kpis\ndat$ambition <- temp$response_answer[match(dat$org, temp$organization)]\noptions(repr.plot.width=10,repr.plot.height=6)\np <- ggplot(dat, aes(x=ambition, y=total)) + \n geom_hline(yintercept = 0, linetype = \"dashed\", size=.3) +\n geom_boxplot(outlier.shape = NA) +\n coord_flip() +\n theme_bw() +\n theme(text = element_text(size = 15)) +\n labs(x=\"\", y= \"KPI distribution\") +\n ggtitle(\"KPI distribution by response (question 1.0)\")\np \n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Well, *Yes* response is ahead but there is a lot of variance. It means that cities responding Yes doesn't necessarly respond with more rigor than others. Due to the construction of the indicator, blank response have the worst KPI. \n## e. Evolution of sustainability KPI 2018-2019\nWe can visualize sustainability ambition KPI evolution from 2018 to 2019. \n* In green are cities with an increasing KPI \n* in red are cities with a decreasing KPI \n\nNote that we don't visualize real KPI values, but its rank in the distribution of each year (standardized from 1 to 100). "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#### 2018 to 2020\n\ndf18_20 <- lapply(c(\"2018\",\"2019\",\"2020\"), function(x){\n path <- paste0(\"../input/cdp-unlocking-climate-solutions/Cities/Cities Responses/\",x,\"_Full_Cities_Dataset.csv\")\n df <- fread(path, encoding = \"UTF-8\", stringsAsFactors = FALSE)\n colnames(df) <- gsub(\" \",\"_\", tolower(colnames(df)))\n df$column_number <- as.character(df$column_number)\n df$column_number <- ifelse(nchar(df$column_number) == 1, paste0(\"0\",df$column_number), df$column_number)\n df$column <- paste(\"col\",df$column_number,df$column_name,sep=\"_\")\n df$row <- paste(\"row\",df$row_number, df$row_name,sep= \"_\")\n # Extract features\n\n feat <- get_kpi_features(df, words_by_section = FALSE)\n\n # One global KPI\n kpi=extract_kpi(feat, plot=FALSE)\n \n kpi <- data.frame(org=feat$organization[order(kpi)],rank=1:length(kpi),year = x)\n kpi$rank <- 100*kpi$rank/nrow(kpi)\n return(kpi)\n})\n\ndf18_20 <- rbindlist(df18_20)\ndf18_20 <- dcast(df18_20, org ~ year, value.var=\"rank\")\ndf18_20 <- df18_20[!is.na(df18_20$`2018`) & !is.na(df18_20$`2019`),]\ndf18_20$dyn <- ifelse(df18_20$`2019` > df18_20$`2018`, df18_20$`2018`/df18_20$`2019`, 2-df18_20$`2019`/df18_20$`2018`)\n\n\np <- ggplot(df18_20, aes( x= `2018`, y= `2019`, colour=dyn ,label=org))+\ngeom_abline(slope=1, intercept=0)+\ngeom_point(size=2.5)+\nscale_colour_gradient2(high=\"#F8766D\",low=\"#7CAE00\",mid=\"#fffcbf\", midpoint=1)+\ntheme_bw()+\ntheme(legend.position = \"none\")+\nlabs(x= \"KPI rank 2018\", y=\"KPI rank 2019\")\n\npltly <- ggplotly(p, tooltip = c(\"org\",\"2018\",\"2019\"))\n\n\nhtmlwidgets::saveWidget(pltly, \"plot0.html\")\nIRdisplay::display_html('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## e. Calculation for corporates\nWe use exactly the same method to calculate corporates sustainability ambition KPI. \nHere the distributions of the two samples KPI.\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Climate change 2019\n\ndis_c <- fread(\"../input/cdp-unlocking-climate-solutions/Corporations/Corporations Disclosing/Climate Change/2019_Corporates_Disclosing_to_CDP_Climate_Change.csv\",encoding = \"UTF-8\")\ndfc <- fread(\"../input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2019_Full_Climate_Change_Dataset.csv\",encoding = \"UTF-8\")\ndfc$country <- dis_c$country[match(dfc$account_number, dis_c$account_number)]\ndfc$industries <- dis_c$industries[match(dfc$account_number, dis_c$account_number)]\ndfc$sector <- dis_c$primary_questionnaire_sector[match(dfc$account_number, dis_c$account_number)]\n\n\ndfc$column_number <- as.character(dfc$column_number)\ndfc$column_number <- ifelse(nchar(dfc$column_number) == 1, paste0(\"0\",dfc$column_number), dfc$column_number)\ndfc$column <- paste(\"col\",dfc$column_number,dfc$column_name,sep=\"_\")\ndfc$row <- paste(\"row\",dfc$row_number, dfc$row_name,sep= \"_\")\n\n# Change column names to match with cities one\ncolnames(dfc)[colnames(dfc) == \"response_value\"] <- \"response_answer\"\ncolnames(dfc)[colnames(dfc) == \"module_name\"] <- \"section\"\n\n# Extract features for city responses\nfeatc <- get_kpi_features(dfc, words_by_section = FALSE)\n\n# One global KPI\nkpic <- extract_kpi(featc, plot=FALSE)\n\n# One KPI per parent section\nkpis_c <- lapply(unique(dfc$section),function(x){\n extract_kpi(featc,section = x, plot=FALSE)\n})\n\nkpis_c <- as.data.frame(do.call(cbind,kpis_c))\ncolnames(kpis_c) <- gsub(\" \",\"_\",tolower(unique(dfc$section)))\ncolnames(kpis_c)[1] <- \"none\"\nkpis_c$org <- featc$organization\nkpis_c$total <- kpic\n\n# Water security 2019\n\ndis_w <- fread(\"../input/cdp-unlocking-climate-solutions/Corporations/Corporations Disclosing/Water Security/2019_Corporates_Disclosing_to_CDP_Water_Security.csv\",encoding = \"UTF-8\")\ndfw <- fread(\"../input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Water Security/2019_Full_Water_Security_Dataset.csv\",encoding = \"UTF-8\")\ndfw$country <- dis_c$country[match(dfw$account_number, dis_c$account_number)]\ndfw$industries <- dis_c$industries[match(dfw$account_number, dis_c$account_number)]\ndfw$sector <- dis_w$primary_questionnaire_sector[match(dfw$account_number, dis_w$account_number)]\n\n\ndfw$column_number <- as.character(dfw$column_number)\ndfw$column_number <- ifelse(nchar(dfw$column_number) == 1, paste0(\"0\",dfw$column_number), dfw$column_number)\ndfw$column <- paste(\"col\",dfw$column_number,dfw$column_name,sep=\"_\")\ndfw$row <- paste(\"row\",dfw$row_number, dfw$row_name,sep= \"_\")\n\n# Change column names to match with cities one\ncolnames(dfw)[colnames(dfw) == \"response_value\"] <- \"response_answer\"\ncolnames(dfw)[colnames(dfw) == \"module_name\"] <- \"section\"\n\n# Extract features for city responses\nfeatw <- get_kpi_features(dfw, words_by_section = FALSE)\n\n# One global KPI\nkpiw <- extract_kpi(featw, plot=FALSE)\n\n# One KPI per parent section\nkpis_w <- lapply(unique(dfw$section),function(x){\n extract_kpi(featw,section = x, plot=FALSE)\n})\n\nkpis_w <- as.data.frame(do.call(cbind,kpis_w))\ncolnames(kpis_w) <- gsub(\" \",\"_\",tolower(unique(dfw$section)))\ncolnames(kpis_w)[1] <- \"none\"\nkpis_w$org <- featw$organization\nkpis_w$total <- kpiw\n\n# Plot the two distributions\n\ntemp <- data.frame(kpi = c(kpic,kpiw),\n group=c(rep(\"Climate Change\", length(kpic)), rep(\"Water Security\", length(kpiw))),\n org=c(featc$org, featw$org))\n\ntemp$authority[temp$org %in% featc$org] <- dis_c$authority_types[match(temp$org[temp$org %in% featc$org], dis_c$organization)]\ntemp$authority[temp$org %in% featw$org] <- dis_w$authority_types[match(temp$org[temp$org %in% featw$org], dis_w$organization)]\n\np <- ggplot(temp, aes(x=group, y=kpi)) +\n geom_boxplot(width = 0.4, fill = \"white\") +\n geom_jitter(aes(color = group, group=org), \n width = 0.1, size = 1) +\n scale_color_manual(values = c(\"#E7B800\", \"#00AFBB\")) + \n labs(x = NULL, y=\"Sustainability ambition KPI\", color=\"Sample\") +\n theme_bw()\n\npltly <- ggplotly(p, tooltip = c(\"org\", \"kpi\"))\n\nhtmlwidgets::saveWidget(pltly, \"plot3.html\")\nIRdisplay::display_html('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's see the KPI distribution by sector. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"options(repr.plot.width=10,repr.plot.height=8)\n# merge disclosing data\ntemp$sector_c <- dis_c$primary_questionnaire_sector[match(temp$org, dis_c$organization)]\ntemp$sector_w <- dis_w$primary_questionnaire_sector[match(temp$org, dis_w$organization)]\ntemp$sector <- ifelse(is.na(temp$sector_c), temp$sector_w, temp$sector_c)\n\n# plot\np <- ggplot(temp, aes(x=sector, y= kpi, fill=group))+\ngeom_hline(yintercept = 0, linetype = \"dashed\", size=.3) +\ngeom_boxplot(outlier.shape=NA)+\nscale_fill_manual(values = c(\"#E7B800\", \"#00AFBB\"))+\ntheme_bw()+\ntheme(text = element_text(size=15))+\nlabs(x = \"\", y= \"Ambition KPI\", fill = \"\")+\ncoord_flip()\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We notice that some sectors are \"good\" or \"bad\" simultaneously at each of the two questionnaires. For example, \"Transport OEMS\" industry seems unlikely to take seriously these problems whereas \"Electric utilities\" are pretty ambitious.\n\nNow, we can look at the influence of the requester type on sustainability ambition. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"p <- ggplot(temp, aes(x=group, y=kpi, colour=authority)) +\n geom_boxplot(width = 0.4, fill = \"white\") +\n labs(x = NULL, y=\"Ambition KPI\", colour=\"Requester\") +\n theme_bw()+\ntheme(text = element_text(size=15))\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There's no doubt that companies requested by investors have a higher sustainability ambition KPI. Intuitively, we can say that these companies have a clear interest in matching the ambitions of their investors. \n## f. Corporate validation\n\nSince external data on corporates are pretty rare, we can validate KPIs by these specific questions relating to ambition : \n* W3.3 : Does your organization undertake a water-related risk assessment?\n* C4.1 : Did you have an emissions target that was active in the reporting year?\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\ndat <- kpis_c\ntemp <- dfc[dfc$question_number == \"C4.1\",]\ndat$ambition <- temp$response_answer[match(dat$org, temp$organization)]\n\np1 <- ggplot(dat, aes(x=ambition, y=total)) + \n geom_hline(yintercept = 0, linetype = \"dashed\", size=.3) +\n geom_boxplot(outlier.shape = NA) +\n coord_flip() +\n theme_bw() +\n labs(x=\"\", y= \"Climate Change KPI\") +\n ggtitle(\"Question C4.1\")+\n theme(text = element_text(size = 20))\n \n\ndat <- kpis_w\ntemp <- dfw[dfw$question_number == \"W3.3\",]\ndat$ambition <- temp$response_answer[match(dat$org, temp$organization)]\n\np2 <- ggplot(dat, aes(x=ambition, y=total)) + \n geom_hline(yintercept = 0, linetype = \"dashed\", size=.3) +\n geom_boxplot(outlier.shape = NA) +\n coord_flip() +\n theme_bw() +\n labs(x=\"\", y= \"Water Security KPI\") +\n ggtitle(\"Question W3.3\")+\n theme(text = element_text(size = 20))\n\ngrid.arrange(p1,p2,ncol=2)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The gap is clear for both questionnaires. We almost see a linear trend for Climate Change questionnaire with equal KPIs for \"Absolute target\" and \"Intensity target\". \n\n# 2. Environmental issues and social inequity : an inextricable link\n\nIn this section, I will explore the interactions between environmental and social issues in the United States. \nI use various data sources : \n* social and pollution data from [EJscreen](https://www.epa.gov/ejscreen) 2019 database. Census block level\n* disease data from CDC. Census tract level\n* Sustainability ambition KPI for 2019 cities calculated in the previous section. City level\n \nEJscreen is an official database from EPA, mapping pollution and demographic data at census block scale. \nPollution data include : \n* Air toxics cancer risk\n* Air toxics respiratory hazard index\n* Air toxics neurological hazard index\n* Traffic proximity and volume\n* Indicator for major direct dischargers to water\n* Proximity to National Priorities List (NPL) sites\n* Proximity to Risk Management Plan (RMP) facilities\n* Proximity to Treatment Storage and Disposal (TSDF) facilities\n* Ozone level in air\n* PM2.5 level in air\n* % pre-1960 housing (lead paint indicator)\n\nDemographic data :\n\n* % under age 5\n* % over age 64\n* % minority\n* % low-income\n* % less than high school\n* % of households (interpreted as individuals) in linguistic isolation\n* Demographic Index (based on 2 factors, % low-income and % minority)\n\nMore importantly, EJscreen is also mapping **EJ Indexes** (EJ standing for \"Environmental Justice\"), which combines pollution and demographic variables. \nHence, an EJ Index is higher in block groups with large numbers of mainly low-income and/or minority residents with a higher environmental indicator value.\n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Import data\ncensus <- fread(\"../input/cdp-unlocking-climate-solutions/Supplementary Data/CDC Social Vulnerability Index 2018/SVI2018_US.csv\",integer64 = \"character\")\nsvi <- fread(\"../input/cdp-unlocking-climate-solutions/Supplementary Data/CDC 500 Cities Census Tract Data/500_Cities__Census_Tract-level_Data__GIS_Friendly_Format___2019_release.csv\", integer64 = \"character\")\nsvi <- svi[,!grepl(\"5CI\",colnames(svi)),with=FALSE]\nej <- read_sf(\"../input/environmental-justice-by-epa/ej.shp\")\nej_labs <- read.csv(\"../input/environmental-justice-by-epa/2019_EJSCREEN_columns_explained.csv\",stringsAsFactors = FALSE)\nzip <- read_sf(\"../input/2017-zip-code-shapefiles/cb_2017_us_zcta510_500k.shp\")\ncities <- read.csv(\"../input/cdp-unlocking-climate-solutions/Supplementary Data/Simple Maps US Cities Data/uscities.csv\", stringsAsFactors = FALSE)\n\n# Use zip list provided in uscities.csv to merge zip polygons into city polygons\ncities <- split(cities$zips, cities$city)\nzip$city <- NA_character_\nfor (i in 1:length(cities)) zip[zip$ZCTA5CE10 %in% unlist(strsplit(cities[[i]],\" \")),\"city\"] <- names(cities)[i]\nzip <- zip[!is.na(zip$city),]\n\nzip <- suppressMessages(summarise(group_by(zip, city)))\n\nzip <- st_transform(zip, crs=4326)\n\n# recode NA\ncensus$EP_UNEMP[census$EP_UNEMP == -999] <- NA\ncensus$EP_POV[census$EP_POV == -999] <- NA\n\n# Variables of interest from EJ data\n\nvars <- c(\"UNDER5PCT\",\"OVER64PCT\",\"MINORPCT\",\"LOWINCPCT\",\"LESSHSPCT\",\"LINGISOPCT\",\"VULEOPCT\",\"PRE1960PCT\",\n \"ACSTOTPOP\",\"CANCER\",\"RESP\",\"DSLPM\",\"PTRAF\",\"PWDIS\",\"PNPL\",\"PRMP\",\"PTSDF\",\"OZONE\",\"PM25\",\n colnames(ej)[substr(colnames(ej),1,2) == \"D_\"])\n\n# Variables of interest from disease data\nsvi_vars <- colnames(svi)[grepl(\"CrudePrev\", colnames(svi))]\n\n# Merge disease data \nej$ID <- as.character(ej$ID)\nej$FIPS <- substr(ej$ID,1,nchar(ej$ID)-1)\n\nshp <- ej\nst_geometry(ej) <- NULL\nej <- setDT(ej)[,c(\"FIPS\",vars),with=F]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## a - Characterization of vulnerable population\n\nVulnerable population can be defined in multiple ways. Low income, linguistic isolation or disability could be taken in account. \nIn this report, I use a synthetic variable from EJscreen, vulnerability index, which combine **low income** and **minority percentage** by census block. \nDefinitions from EPA :\n* **Percent Minority :** The percent of individuals in a block group who list their racial status as a race other than white alone and/or list their ethnicity as Hispanic or Latino. That is, all people other than non-Hispanic white-alone individuals. The word \"alone\" in this case indicates that the person is of a single race, not multiracial.\n* **Percent Low-Income :** The percent of a block group's population in households where the household income is less than or equal to twice the federal \"poverty level.\"\n\nThe vulnerability index is the percentage of vulnerable people, so it range from 0 to 1 in each census block. \n\nWe explore two facets of social inequity : **health** and **pollution exposure**. \nFirst, we look at the pollution exposure by vulnerability index. The number in parenthesis is the number of census block for each category."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\n\n# Prepare data\nej2 <- ej\n#ej2$pov_class <- cut(ej2$VULEOPCT, quantile(ej2$VULEOPCT,seq(0,1,0.25),na.rm=T))\nej2$pov_class <- cut(ej2$VULEOPCT, seq(0,1,0.25), include.lowest = TRUE)\n\n# Edit labels\nej2[,pov_count := .N,.(pov_class)]\nej2[,pov_class_label := paste0(pov_class,\" (\",pov_count,\")\"),]\npov_class_label <- unique(ej2$pov_class_label)\n\n# Mean by category\npollution <- c(\"CANCER\",\"RESP\",\"DSLPM\",\"PTRAF\",\"PWDIS\",\"PNPL\",\"PRMP\",\"PTSDF\",\"OZONE\",\"PM25\",\"PRE1960PCT\")\n\nej2 <- ej2[,lapply(.SD,mean,na.rm=TRUE),.(pov_class),.SDcols=pollution]\nfor (i in 2:ncol(ej2)) ej2[[i]] <- ej2[[i]]/max(ej2[[i]])\ndfmelt <- melt(ej2, id.vars = \"pov_class\")\ndfmelt <- dfmelt[order(dfmelt$variable,dfmelt$pov_class),] \ndfmelt <- dfmelt[!is.na(dfmelt$pov_class),]\ndfmelt[,q4q1 := value[length(value)]/value[1],.(variable)]\ndfmelt$variable <- reorder(dfmelt$variable, dfmelt$q4q1,unique)\nlevels(dfmelt$variable) <- ej_labs$glossaryfieldname[match(levels(dfmelt$variable), ej_labs$gdbfieldname)]\n\n# Plot\noptions(repr.plot.width=20,repr.plot.height=12)\n\np <-ggplot(dfmelt, aes(x=variable, y=value, fill=pov_class))+\n geom_bar(stat = \"identity\", position = position_dodge(), size= .1)+\n scale_fill_brewer(label = pov_class_label, palette=\"Oranges\") +\n coord_flip()+\ntheme_classic()+\ntheme(text=element_text(size=25))+\nlabs(fill = \"Vulnerability index\",\n x=\"\",\n y=\"Relative pollution intensity\")+\nggtitle(\"Pollution intensity by vulnerability index\")\n\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Pollution exposure is clearly more intense in poorest areas. This must be directly linked with the spatial distribution of these households in the urban fabric, which is characterized by a high proximity to factories and a dense traffic.\n\nNow, we look at disease prevalence by vulnerability index. \nData about disease prevalence comes from the CDC dataset at census tract scale, I selected only diseases and omitted the variables related to care like cholesterol test, mammogram ... etc. Pollution data and disease prevalence have been merged using FIPS code. \n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"#maladie <- colnames(svi)[grepl(\"CrudePrev\", colnames(svi))]\nmaladie <- c(\n\"ARTHRITIS_CrudePrev\",\n\"BPHIGH_CrudePrev\",\n\"CANCER_CrudePrev\",\n\"CASTHMA_CrudePrev\",\n\"CHD_CrudePrev\",\n\"COPD_CrudePrev\",\n\"DIABETES_CrudePrev\",\n\"HIGHCHOL_CrudePrev\",\n\"MHLTH_CrudePrev\",\n\"OBESITY_CrudePrev\",\n\"PHLTH_CrudePrev\",\n\"STROKE_CrudePrev\",\n\"TEETHLOST_CrudePrev\"\n)\n\nmaladie_lab <- c(\n\"Arthritis\",\n\"High blood pressure\",\n\"Cancer\",\n\"Asthma\",\n\"Coronary heart disease\",\n\"Chronic obstructive pulmonary disease\",\n\"Diabete\",\n\"High cholesterol\",\n\"Mental health\",\n\"Obesity\",\n\"Physical health\",\n\"Stroke\",\n\"Teethlost\" \n)\n\n\n# Prepare data\nej2 <- ej\nej2 <- merge(ej2, svi[,c(\"TractFIPS\", svi_vars),with=FALSE],\n by.x=\"FIPS\",by.y=\"TractFIPS\", all.x = TRUE)\nej2$pov_class <- cut(ej2$VULEOPCT, seq(0,1,0.25), include.lowest = TRUE)\n\n# Edit labels\nej2[,pov_count := .N,.(pov_class)]\nej2[,pov_class_label := paste0(pov_class,\" (\",pov_count,\")\"),]\npov_class_label <- unique(ej2$pov_class_label)\n\n# Mean by category\nej2 <- ej2[,lapply(.SD,mean,na.rm=TRUE),.(pov_class),.SDcols=maladie]\nfor (i in 2:ncol(ej2)) ej2[[i]] <- ej2[[i]]/max(ej2[[i]])\ndfmelt <- melt(ej2, id.vars = \"pov_class\")\ndfmelt$variable <- maladie_lab[match(dfmelt$variable, maladie)]\nmaladie_desc <- split(dfmelt, dfmelt$variable)\nmaladie_desc <- maladie_desc[names(maladie_desc) %in% maladie_lab]\nmaladie_desc <- sapply(maladie_desc, function(x){\n \n x <- x[!is.na(x$pov_class),]\n x <- x[order(pov_class),]\n sorted <- x$value\n\n if (sorted[length(sorted)] > sorted[1]) return(TRUE)\n else return(FALSE)\n})\n\n\ndfmelt <- dfmelt[order(dfmelt$variable,dfmelt$pov_class),] \ndfmelt <- dfmelt[!is.na(dfmelt$pov_class),]\ndfmelt[,q4q1 := value[length(value)]/value[1],.(variable)]\ndfmelt$variable <- reorder(dfmelt$variable, dfmelt$q4q1,unique)\n\noptions(repr.plot.width=20,repr.plot.height=15)\np <- ggplot(dfmelt[!is.na(pov_class),,], aes(x=variable, y=value, fill=pov_class))+\n geom_bar(stat = \"identity\", position = position_dodge())+\n scale_fill_brewer(label = pov_class_label, palette=\"Oranges\") +\n coord_flip()+\ntheme_classic()+\ntheme(text=element_text(size=25))+\nlabs(fill = \"Vulnerability index\",\n x=\"\",\n y=\"Relative disease prevalence\")+\nggtitle(\"Disease prevalence by vulnerability index\")\n\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"92% (12 among 13) of recensed diseases have higher prevalence for the most vulnerable people. \nWithout talking about causality, we can safely say that low-income/minority people are :\n* more likely to developp diseases\n* more exposed to pollution\n\nSo here is the question that comes to mind : \n### Is there a link between disease prevalence and pollution exposition ?\n\nInstead of doing correlation matrix, it's more appropriate to do regression analysis. \nIn this way, we care about counfounding factors like age or low income and we measure effects *ceteris paribus*. \nSince pollution and social features have very different ranges, all features are scaled. \nA regression analysis is conducted for each disease : "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Merge disease to EJ\nej2 <- ej\nej2 <- merge(ej2, svi[,c(\"TractFIPS\", svi_vars),with=FALSE],\n by.x=\"FIPS\",by.y=\"TractFIPS\", all.x = TRUE)\n\n\n\n# Scale data because of scales inconsistency\nfor (i in c(\"LOWINCPCT\",\"MINORPCT\",\"OVER64PCT\",\"LESSHSPCT\",\"LINGISOPCT\",pollution, maladie)) ej2[[i]] <- scale(ej2[[i]])\n\n# Regressions analysis\nregs <- list()\nfor (i in maladie){\n form <- as.formula(paste(i,\"~\",paste0(c(\"LOWINCPCT\",\"MINORPCT\",\"OVER64PCT\",\"LESSHSPCT\",\"LINGISOPCT\",pollution), collapse = \"+\")))\n regs[[i]] <- coefficients(summary(lm(form, data = ej2)))\n}\n\n# Reshape for plotting\ncoefs <- do.call(cbind,\n lapply(regs,function(x) x[,1]))\ncoefs <- data.frame(covariate = rownames(coefs),coefs) \ncoefs <- melt(coefs, id.vars=\"covariate\") \n \npval <- do.call(cbind,\n lapply(regs,function(x) x[,4])) \npval <- data.frame(covariate = rownames(pval),pval) \npval <- melt(pval, id.vars=\"covariate\") \ncolnames(pval)[3] <- \"pval\" \n \ncoefs <- merge(coefs, pval, by=c(\"covariate\",\"variable\")) \n \ncoefs$label <- paste0(round(coefs$value, digits=2))\ncoefs$label <- ifelse(coefs$pval < 0.1, paste0(round(coefs$value, digits=2), \"*\"), coefs$label) \ncoefs$label <- ifelse(coefs$pval < 0.05, paste0(round(coefs$value, digits=2), \"**\"), coefs$label) \ncoefs$label <- ifelse(coefs$pval < 0.01, paste0(round(coefs$value, digits=2), \"***\"), coefs$label) \ncoefs <- coefs[coefs$covariate != \"(Intercept)\",]\nbornes <- c(min(coefs$value), max(coefs$value) ) \n# Labels\ncoefs$covariate <- ej_labs$glossaryfieldname[match(coefs$covariate, ej_labs$gdbfieldname)]\ncoefs$variable <- maladie_lab[match(coefs$variable, maladie)] \n\n# Heatmap \noptions(repr.plot.width=18,repr.plot.height=17)\nggheatmap <- ggplot(coefs, aes(variable, covariate, fill = value))+\n geom_tile(color = \"white\")+\n scale_fill_gradient2(low = \"blue\", high = \"red\", mid = \"white\", \n midpoint = 0, limit = bornes, space = \"Lab\",\n name=\"Regression\\nCoefficient\") +\n theme_minimal()+ # minimal theme\n theme(axis.text.x = element_text(angle = 45, vjust = 1, \n size = 15, hjust = 1),\n axis.text.y = element_text(size = 15, hjust = 1),\n text = element_text(size=15))+\n coord_fixed()+\n geom_text(aes(variable, covariate, label = label), color = \"black\", size = 5)+\n labs(x=\"Outcome\", y=\"Covariates\") \nggheatmap ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Other things being equal**, we observe that the largest effects belong to social characteristics like old age or low income. \n\nDiseases like **arthritis, cholesterol, heart disease and cancer** are naturally more influenced by age. \nOn the other hand, **obesity, mental & physical diseases and teethloss** seems to be more influenced by being poor. \n**Diabete and stroke** seems to be equally influenced by age and having a low income. \n\nIn fact, it seems logical since : \n* low income people simply can't easily access to care\n* low income people can't afford quality food which is generally more expensive \n\nHowever, many pollution-related coefficients are highly significant which means that pollution exposure still plays a **significant role in disease prevalence**. For example, high level of ozone and lead paint indicator explains asthma disease pretty well. "},{"metadata":{},"cell_type":"markdown","source":"## b - Introducing sustainability ambition KPI into social - environmental interactions\n\nNow, we can look what is going on when we cross social and environmental data to cities KPI. \nFirst, we have to aggregate social and environmental features by city. For this, I use ZIP lists provided in uscities.csv along with ZIP shapefile to construct cities areas. \nIn order to take advantage of the fine-scale data and capture its distribution in cities area, I compute the following features : \n\n* mean weighted by census block population\n* maximum, to capture extreme values\n* Gini coefficient, to capture the inequity in the distribution :\n\n[Gini coefficient](https://en.wikipedia.org/wiki/Gini_coefficient) is a common indicator for expressing inequality in a distribution. It equals to 0 in the case of perfect equality.\n\n* Moran Index, to capture spatial correlation :\n\n[Moran index](https://en.wikipedia.org/wiki/Moran%27s_I) is a spatial autocorrelation indicator. Spatial autocorrelation measure the spatial dependance between an entity and its neighbors. \nMoran Index range from : \n* -1, perfectly dispersed like (b) figure\n* to 1, perfectly clustered like (a) figure\n* passing by 0, random like (c) figure\n\n\n\n\nFor more readability, I will focus on two pollution features : *Traffic proximity and volume* and *Proximity to Treatment Storage and Disposal (TSDF) facilities*. I choose these ones because they are part of the most widespread pollutions observed in the poorest census blocks, compared to other vulnerability categories. \n\n### i - Population\n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# First, aggregate to census tract\n\nejFIPS <- ej[,lapply(.SD,weighted.mean,w = ACSTOTPOP, na.rm=TRUE),.(FIPS), .SDcols=c(vars)]\npop <- ej[,.(ACSTOTPOP=sum(ACSTOTPOP,na.rm=T)),.(FIPS)]\nejFIPS$ACSTOTPOP <- pop$ACSTOTPOP[match(ejFIPS$FIPS, pop$FIPS)]\n\nejFIPS <- merge(ejFIPS, svi[,c(\"TractFIPS\", svi_vars),with=FALSE],\n by.x=\"FIPS\",by.y=\"TractFIPS\", all.x = TRUE)\n\n\nfor (i in 1:ncol(ejFIPS)) {\n ejFIPS[[i]][is.nan(ejFIPS[[i]])] <- NA \n ejFIPS[[i]][is.infinite(ejFIPS[[i]])] <- NA \n}\n\n\n# Spatial join census block to cities\nshp2 <- suppressMessages(st_join(st_centroid(shp), zip))\n\nshp$city <- shp2$city[match(shp$ID, shp2$ID)]\nshp2 <- shp[!is.na(shp$city),]\n\n# Keep only cities participating in CDP questionnaire to light the whole process\npts_zip <- suppressMessages(st_join(pts,zip))\n\nshp2 <- shp2[shp2$city %in% pts_zip$city[!is.na(pts_zip$city)],]\n\n# Merge polygons to census tract level for L.A map\nshp3 <- suppressMessages(summarise(group_by(shp2, FIPS)))\n\n# Spatial join census tract to cities\nshp4 <- suppressMessages(st_join(st_centroid(shp3), zip))\n\nshp3$city <- shp4$city[match(shp3$FIPS, shp4$FIPS)]\nshp3 <- shp3[!is.na(shp3$city),]\n\n# Moran Index\nej$ID <- shp$ID \nmoran_list <- lapply(unique(shp2$city), function(x){\n\n temp <- shp2[shp2$city == x,]\n \n if (nrow(temp) < 10) return(NULL)\n \n #temp <- merge(temp, ej[,c(\"ID\",vars),with=FALSE], by=\"ID\", all.x=TRUE)\n temp <- as_Spatial(temp)\n tempn <- poly2nb(temp)\n tempw <- nb2listw(tempn, zero.policy = TRUE)\n \n res <- vector(mode = \"numeric\", length=length(c(vars)))\n names(res) <- c(vars)\n for (i in c(vars)){\n mor <- tryCatch(moran.test(temp[[i]], tempw, na.action = na.exclude, zero.policy = TRUE), error=function(e) return(NA))\n res[[i]] <- tryCatch(mor$estimate[1], error=function(e) return(NA))\n }\n return(res)\n})\n\nnames(moran_list) <- unique(shp2$city)\nmoran_list <- moran_list[sapply(moran_list,length) > 0]\nmoran_list <- data.frame(city=names(moran_list), do.call(rbind,moran_list))\n\n \n# Aggregate population, poverty, pollution by city area \nej$city <- shp2$city[match(ej$ID, shp2$ID)]\n#cols <- colnames(ej)[2:(ncol(ejFIPS)-1)] \nej_mean <- ej[,lapply(.SD, weighted.mean, w=ACSTOTPOP,na.rm=TRUE),.(city), .SDcols=vars] \nej_max <- ej[,lapply(.SD, max,na.rm=TRUE),.(city), .SDcols=vars] \nej_gini <- ej[,lapply(.SD, ineq::Gini,na.rm=TRUE),.(city), .SDcols=vars] \n \n# Census block area\nshp_area <- st_area(shp)\nej$area <- as.numeric(shp_area)[match(ej$ID, shp$ID)] \nej_pop <- ej[,.(pop=sum(ACSTOTPOP,na.rm=T),\n density = weighted.mean(ACSTOTPOP/area, ACSTOTPOP, na.rm=TRUE)),.(city)] \n \n# Aggregate population \ntemp <- ej_pop\ntemp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)] \n \noptions(repr.plot.width=14,repr.plot.height=5) \np1 <- ggplot(temp, aes(x= kpi, y=log(pop)))+\n geom_point()+\n geom_smooth(method=\"lm\", formula = y ~ x)+\n theme_bw() +\n theme(text=element_text(size=15))+\n labs(x=\"Ambition KPI\",y=\"log(Population)\")\n \np2 <- ggplot(temp, aes(x= kpi, y=log(density)))+\n geom_point()+\n geom_smooth(method=\"lm\", formula = y ~ x)+\n theme_bw() +\n theme(text=element_text(size=15))+\n labs(x=\"Ambition KPI\",y=\"log(Population density)\")\n \ngrid.arrange(p1,p2, ncol=2) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"High KPI cities tend to be more populated and more densely populated. Hence, these are probably large urbanized cities. \n### ii - Vulnerability"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Plot KPI vs vulnerability \noptions(repr.plot.width=14,repr.plot.height=8)\n\nplots <- list()\nfor (i in c(\"ej_mean\",\"ej_max\",\"ej_gini\",\"moran_list\")){\n if (grepl(\"moran\",i)) title <- \"Moran\"\n else title <- unlist(strsplit(i,\"_\"))[2]\n title <- paste0(toupper(substr(title,1,1)),substr(title,2,nchar(title)))\n \n temp <- get(i)\n temp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n plots[[i]] <- ggplot(temp, aes(x= kpi, y=VULEOPCT))+\n geom_point()+\n geom_smooth(method=\"lm\", formula = y ~ x)+\n theme_bw()+\n labs(x=\"Ambition KPI\",y=\"\")+\n ggtitle(title)\n}\n \ngrid.arrange(grobs=plots,ncol=2) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Well, it's hard to tell anything about these plots. Moran Index may be higher for high KPIS, but it's not really clear. \n\n### ii - Pollution\n\nNow, let's have a look at pollution intensity vs KPI. \nMaximum and mean values have been tranformed to logarithmic scale."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"options(repr.plot.width=18,repr.plot.height=6)\nplots <- list()\n\nfor (i in c(\"PTRAF\",\"PTSDF\")){\n \n for (j in c(\"ej_gini\",\"moran_list\")){\n temp <- as.data.frame(get(j))\n \n temp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n temp <- temp[,c(\"kpi\",i)] \n colnames(temp) <- c(\"kpi\",\"value\")\n temp$indic <- i\n temp$agg <- j \n plots[[length(plots) + 1]] <- temp\n }\n \n for (j in c(\"ej_mean\",\"ej_max\")){\n temp <- as.data.frame(get(j))\n \n temp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n \n \n temp[[i]] <- log(temp[[i]])\n temp <- temp[,c(\"kpi\",i)] \n colnames(temp) <- c(\"kpi\",\"value\")\n temp$indic <- i\n temp$agg <- j \n plots[[length(plots) + 1]] <- temp\n }\n}\n\nplots <- rbindlist(plots)\n\n# Labels\nplots$agg <- ifelse(plots$agg == \"moran_list\", \"Moran Index\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_gini\", \"Gini Index\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_max\", \"Maximum\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_mean\", \"Mean\", plots$agg)\nplots$indic <- ifelse(plots$indic == \"PTRAF\", \"Traffic intensity\", plots$indic)\nplots$indic <- ifelse(plots$indic == \"PTSDF\", \"TSDF proximity\", plots$indic)\n\n# A little hack to correctly display axes\ng1 <- ggplot(plots, aes(x=kpi, y= (value)))+\ngeom_point()+\ngeom_smooth(method=\"lm\", formula = y ~ x)+\ntheme_bw()+\nfacet_wrap(indic ~ agg, scales=\"free\" , ncol=4)+\ntheme(strip.background = element_blank(),strip.text = element_blank())+\nlabs(x= \"City ambition KPI\", y=\"\")\n\ng2 <- ggplot(plots, aes(x=kpi, y= (value)))+\ngeom_point()+\ngeom_smooth(method=\"lm\", formula = y ~ x)+\ntheme_bw()+\ntheme(strip.text.x = element_text(\n size = 16),strip.text.y = element_text(\n size = 16))+\nfacet_grid(indic ~ agg, scales=\"free\" )+\ntheme(strip.background =element_rect(fill=\"#aec9eb\"))\n\ngt1 = ggplot_gtable(ggplot_build(g1))\ngt2 = ggplot_gtable(ggplot_build(g2))\ngt1$grobs[grep('strip-t.+1$', gt1$layout$name)] = gt2$grobs[grep('strip-t', gt2$layout$name)]\n\ngt.side1 = gtable_filter(gt2, 'strip-r-1')\ngt.side2 = gtable_filter(gt2, 'strip-r-2')\ngt.side3 = gtable_filter(gt2, 'strip-r-3')\ngt.side4 = gtable_filter(gt2, 'strip-r-4')\n\ngt1 = gtable_add_cols(gt1, widths=gt.side1$widths[1], pos = -1)\ngt1 = gtable_add_grob(gt1, zeroGrob(), t = 1, l = ncol(gt1), b=nrow(gt1))\n\npanel_id <- gt1$layout[grep('panel-.+1$', gt1$layout$name),]\ngt1 = gtable_add_grob(gt1, gt.side1, t = panel_id$t[1], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side2, t = panel_id$t[2], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side3, t = panel_id$t[3], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side4, t = panel_id$t[4], l = ncol(gt1))\n\ngrid.newpage()\ngrid.draw(gt1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Much more interesting here, Moran Index, mean and maximum values seems to increase along with KPI, especially for traffic proximity. It means that high KPI cities tend to be **more polluted, with more extreme values and these areas are more clustered** in the urban fabric. \n\n### iii - Environmental Justice Indexes\n\nHave a look at the intersection between social and environmental issues. \nMaximum and mean values have been tranformed to logarithmic scale."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"options(repr.plot.width=18,repr.plot.height=6)\nplots <- list()\n\nfor (i in c(\"D_PTRAF_2\",\"D_PTSDF_2\")){\n \n for (j in c(\"ej_gini\",\"moran_list\")){\n temp <- as.data.frame(get(j))\n \n temp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n \n # Remove gini outliers\n if (j == \"ej_gini\"){\n temp <- temp[temp[[i]] > quantile(temp[[i]], 0.1,na.rm=TRUE) & temp[[i]] < quantile(temp[[i]], 0.9,na.rm=TRUE),]\n } \n \n temp <- temp[,c(\"kpi\",i)] \n colnames(temp) <- c(\"kpi\",\"value\")\n temp$indic <- i\n temp$agg <- j \n plots[[length(plots) + 1]] <- temp\n }\n \n for (j in c(\"ej_mean\",\"ej_max\")){\n temp <- as.data.frame(get(j))\n \n temp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n \n temp[[i]] <- log(temp[[i]])\n temp <- temp[,c(\"kpi\",i)] \n colnames(temp) <- c(\"kpi\",\"value\")\n temp$indic <- i\n temp$agg <- j \n plots[[length(plots) + 1]] <- temp\n }\n}\n\nplots <- rbindlist(plots)\n\n# Labels\nplots$agg <- ifelse(plots$agg == \"moran_list\", \"Moran Index\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_gini\", \"Gini Index\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_max\", \"Maximum\", plots$agg)\nplots$agg <- ifelse(plots$agg == \"ej_mean\", \"Mean\", plots$agg)\nplots$indic <- ifelse(plots$indic == \"D_PTRAF_2\", \"Traffic intensity\", plots$indic)\nplots$indic <- ifelse(plots$indic == \"D_PTSDF_2\", \"TSDF proximity\", plots$indic)\n\n# A little hack to correctly display axes\ng1 <- ggplot(plots, aes(x=kpi, y= (value)))+\ngeom_point()+\ngeom_smooth(method=\"lm\", formula = y ~ x)+\ntheme_bw()+\nfacet_wrap(indic ~ agg, scales=\"free\" , ncol=4)+\ntheme(strip.background = element_blank(),strip.text = element_blank())+\nlabs(x= \"City ambition KPI\", y=\"\")\n\ng2 <- ggplot(plots, aes(x=kpi, y= (value)))+\ngeom_point()+\ngeom_smooth(method=\"lm\", formula = y ~ x)+\ntheme_bw()+\ntheme(strip.text.x = element_text(\n size = 16),strip.text.y = element_text(\n size = 16))+\nfacet_grid(indic ~ agg, scales=\"free\" )+\ntheme(strip.background =element_rect(fill=\"#aec9eb\"))\n\ngt1 = ggplot_gtable(ggplot_build(g1))\ngt2 = ggplot_gtable(ggplot_build(g2))\ngt1$grobs[grep('strip-t.+1$', gt1$layout$name)] = gt2$grobs[grep('strip-t', gt2$layout$name)]\n\ngt.side1 = gtable_filter(gt2, 'strip-r-1')\ngt.side2 = gtable_filter(gt2, 'strip-r-2')\ngt.side3 = gtable_filter(gt2, 'strip-r-3')\ngt.side4 = gtable_filter(gt2, 'strip-r-4')\n\ngt1 = gtable_add_cols(gt1, widths=gt.side1$widths[1], pos = -1)\ngt1 = gtable_add_grob(gt1, zeroGrob(), t = 1, l = ncol(gt1), b=nrow(gt1))\n\npanel_id <- gt1$layout[grep('panel-.+1$', gt1$layout$name),]\ngt1 = gtable_add_grob(gt1, gt.side1, t = panel_id$t[1], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side2, t = panel_id$t[2], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side3, t = panel_id$t[3], l = ncol(gt1))\ngt1 = gtable_add_grob(gt1, gt.side4, t = panel_id$t[4], l = ncol(gt1))\n\ngrid.newpage()\ngrid.draw(gt1)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For these two Environmental Justice indexes, there are clearly **more injustices** and these are **more extremes** in high KPI cities. These higher values mean the presence of extremely poor and extremely polluted places simultaneously. \n\n### iv - Cumulated Environmental Justice Indexes\n\nTo not overlook other EJ indexes, we can try to capture the intensity of all EJ indexes (i.e. 11) into one indicator at city scale. I start from maximum EJ Index data at city scale, used below : \n* I cut each EJ Index by decile into 10 classes encoded from 1 to 10\n* I sum all EJ classes by city\n* I cut cities KPI into 5 classes\n\nThus, the higher is the indicator, the more the city meets extreme Environmental Injustice in **multiple** ways. \nHere is the result with one scatterplot and boxplot using KPI classes :"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"EJ <- colnames(ej)[substr(colnames(ej),1,2) == \"D_\"]\n\ntemp <- ej_max\nfor (i in EJ) temp[[paste0(\"class_\",i)]] <- cut(temp[[i]], quantile(temp[[i]], seq(0,1,0.1), na.rm=T), labels=FALSE)\n\ntemp$EJ_class <- rowSums(temp[,paste0(\"class_\",EJ),with=F], na.rm=T)\n#temp$EJ_class <- apply(temp[,paste0(\"class_\",EJ),with=F],1,max)\ntemp$kpi <- pts_zip$total[match(temp$city, pts_zip$city)]\n\n\ntemp$kpi_class <- cut(temp$kpi, 5)\noptions(repr.plot.width=16,repr.plot.height=5) \n\np1=ggplot(temp[!is.na(temp$kpi_class),], aes(x=kpi, y = (EJ_class))) + \ngeom_point() + \ntheme_bw()+ylim(c(-100,200)) + \ntheme(text = element_text(size =15))+\ngeom_smooth(method=\"lm\", formula = y~x) + \nlabs(y=\"Combined EJs\", x=\"Sustainability ambition KPI\")\n\np2=ggplot(temp[!is.na(temp$kpi_class),], aes(x=kpi_class, y = EJ_class)) + \ngeom_boxplot() + \ntheme_bw()+\ntheme(text = element_text(size =15))+\nlabs(y=\"Combined EJs\", x=\"KPI category\")\n\ngrid.arrange(p1,p2, ncol=2)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is pretty clear that sustainability ambition KPI is related to Environmental Injustices intensity, diversity and its spatial concentration. \nHence, **the overlap between environmental issues and social inequity may play a catalytic role in raising awareness.** \n\nThrough the extreme violence and the social pressures this phenomenon generates,this widening of social inequalities can alert and accelerate political decision-making about sustainability developpment."},{"metadata":{},"cell_type":"markdown","source":"## c- Focus on Los Angeles city\n\nWe can display a map just to visualize the overlap between environmental hazards and social inequity in a high KPI city like Los Angeles (KPI = 16.61). \nFor better fluidity, census block polygons have been merged into census tracts. \n\nThe layers shows : \n* Vulnerability Index by EJscreen at decile scale\n* Vulnerability Index LISA \n* Diesel Particule Intensity by EJscreen at decile scale\n* Diesel Particule Intensity LISA\n* Diesel Particule Environmental Justice Index at decile scale\n\nLISA stands for Local Indicator of Spatial Association. Derived from Moran index, it let us localize the significant spatial clusters of a given variable. \nIn red are cluster Low-Low (low value correlated with neighborhood low values). \nIn green are cluster High-High (high value correlated with neighborhood high values). \nCensus tracts are clickable to visualize values. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"shp_la <- shp3[shp3$city == \"Los Angeles\",]\n\nshp_la <- merge(shp_la, ejFIPS, by=\"FIPS\", all.x=TRUE)\n\nshp_la_sp <- as_Spatial(shp_la)\n\n\n# Local moran function\nlisafunction<-function(x){#fonction moran local qui s'arrete pas malgre une erreur\n return(tryCatch(localmoran(x,tempw,zero.policy = T,na.action = na.exclude),error=function(e) NULL))\n} \n\n\nget_lisa_cadran <- function(lisa,x){\n lisa[lisa$Ii>0 & x>0,\"cadran\"]<-1#HH\n lisa[lisa$Ii>0 & x>0,\"color\"]<-\"#2ECC71\"#vert\n lisa[lisa$Ii>0 & x<0,\"cadran\"]<-2#BB\n lisa[lisa$Ii>0 & x<0,\"color\"]<-\"#CD6155\"#rouge\n lisa[lisa$Ii<0 & x>0,\"cadran\"]<-3#HB\n lisa[lisa$Ii<0 & x>0,\"color\"]<-\"blue\"#bleu\n lisa[lisa$Ii<0 & x<0,\"cadran\"]<-4#BH\n lisa[lisa$Ii<0 & x<0,\"color\"]<-\"yellow\"\n lisa[lisa$`Pr(z > 0)`>0.05,\"color\"]<-\"white\"\n return(lisa)\n}\n\n\n\ntempn <- poly2nb(shp_la_sp)\ntempw <- nb2listw(tempn, zero.policy = TRUE) \nfor (j in c(\"VULEOPCT\",\"DSLPM\",\"D_DSLPM_2\")){\n x<-shp_la_sp@data[,j]-mean(shp_la_sp@data[,j],na.rm=T)\n lisa <- lisafunction(shp_la_sp[[j]])\n lisa<-as.data.frame(lisa)\n ind <- complete.cases(lisa)\n x <- x[ind]\n if (length(x) == 0) next\n \n lisa <- lisa[ind,]\n lisa <- get_lisa_cadran(lisa,x)\n \n \n shp_la[[paste0(j,\"_color\")]] <- lisa$color[match(shp_la$FIPS, shp_la_sp@data$FIPS)]\n shp_la[[paste0(j,\"_color\")]][is.na(shp_la[[paste0(j,\"_color\")]])] <- \"white\"\n \n }\n \n \n# Color palettes\npov_pal <- colorQuantile(viridis::viridis(10), n=10, domain = shp_la$VULEOPCT)\n\ndsl_pal <- colorQuantile(viridis::viridis(10), n=10, domain = shp_la$DSLPM)\n\nej_pal <- colorQuantile(viridis::viridis(10), n=10, domain = shp_la$D_DSLPM_2)\n\n# Labels \nshp_la$vul_lisa_lab <- \"N.S\" \nshp_la$vul_lisa_lab <- ifelse(shp_la$VULEOPCT_color == \"#2ECC71\", \"High-High\", shp_la$vul_lisa_lab)\nshp_la$vul_lisa_lab <- ifelse(shp_la$VULEOPCT_color == \"#CD6155\", \"Low-Low\", shp_la$vul_lisa_lab) \nshp_la$dsl_lisa_lab <- \"N.S\" \nshp_la$dsl_lisa_lab <- ifelse(shp_la$DSLPM_color == \"#2ECC71\", \"High-High\", shp_la$dsl_lisa_lab)\nshp_la$dsl_lisa_lab <- ifelse(shp_la$DSLPM_color == \"#CD6155\", \"Low-Low\", shp_la$dsl_lisa_lab) \n \nshp_la$label <- lapply(paste(\"Vulnerability :\",round(shp_la$VULEOPCT, digits=2),\n \"
Vulnerability LISA :\",shp_la$vul_lisa_lab,\n \"
Diesel Particule :\",round(shp_la$DSLPM, digits=2),\n \"
Diesel Particule LISA :\",shp_la$dsl_lisa_lab,\n \"
Diesel Particule EJ :\",round(shp_la$D_DSLPM_2, digits=2)), HTML) \n# Map\nm <- leaflet() %>% addProviderTiles(providers$OpenStreetMap.Mapnik) %>% \naddPolygons(data = shp_la, fillColor= ~pov_pal(VULEOPCT), fillOpacity=1, stroke=TRUE, weight=1.2,color=\"black\", group = \"Vulnerability\",popup=shp_la$label) %>%\naddPolygons(data = shp_la, fillColor= ~VULEOPCT_color, fillOpacity=1,stroke=TRUE, weight=1.2,color=\"black\", group = \"Vulnerability LISA\",popup=shp_la$label) %>% \naddPolygons(data = shp_la, fillColor= ~dsl_pal(DSLPM), fillOpacity=1, stroke=TRUE, weight=1.2,color=\"black\", group = \"Diesel particule\",popup=shp_la$label) %>%\naddPolygons(data = shp_la, fillColor= ~DSLPM_color, fillOpacity=1, stroke=TRUE, weight=1.2,color=\"black\", group = \"Diesel particule LISA\",popup=shp_la$label) %>% \naddPolygons(data = shp_la, fillColor= ~ej_pal(D_DSLPM_2), fillOpacity=1, stroke=TRUE, weight=1.2,color=\"black\", group = \"Diesel particule EJ\",popup=shp_la$label) %>%\n addLayersControl(\n\n baseGroups = c(\"Vulnerability\",\"Vulnerability LISA\", \"Diesel particule\", \"Diesel particule LISA\",\"Diesel particule EJ\"),\n options = layersControlOptions(collapsed = FALSE)\n )\n\nhtmlwidgets::saveWidget(m, \"m2.html\")\ndisplay_html('') ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This map is a nice picture of what we have explored below, high ambition KPI are characterized by : \n* spatial concentration of poverty\n* spatial concentration of pollution\n* overlap of these two clusters\n\n# 3. Water security: Exploration around cities and corporates KPI \nIn this part, I use corporates response at question W5.1 : \n\n\"*For each facility referenced in W4.1c, provide coordinates, total water accounting data and comparisons with the previous reporting year*\" \n\nFirst of all, we can look at the averaged water withdrawal, discharge and consumption versus corporate KPI. \nI use the average instead of the sum to avoid bias due to the way the KPI is constructed. For example, if a company provides information about the consumption of one facility out of 10 and another fills its 10 facilities; then the KPI will potentially be greater for the second. The sum of consumption will also potentially be greater. By using the mean, we get rid of this potential bias. \nNote that water volume has been transformed in logarithmic scale"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Extract facility locations from water questionnaire\nliste <- split(dfw,dfw$question_number)\nliste <- lapply(liste,function(x){\n suppressMessages(dcast(x, account_number +country + sector + organization + question_number + row ~ column,value.var = \"response_answer\",fun.aggregate=function(x) paste(x, collapse = \";\")))\n})\n\ntemp <- liste[[\"W5.1\"]]\ntemp <- temp[col_05_W5.1_C5Latitude != \"\" & `col_03_W5.1_C3Country/Region` == \"United States of America\",,]\ntemp$col_05_W5.1_C5Latitude <- as.numeric(temp$col_05_W5.1_C5Latitude)\ntemp$col_06_W5.1_C6Longitude <- as.numeric(temp$col_06_W5.1_C6Longitude)\ntemp$`col_13_W5.1_C13Total water consumption at this facility (megaliters/year)` <- as.numeric(temp$`col_13_W5.1_C13Total water consumption at this facility (megaliters/year)`)\n\noptions(repr.plot.width=12,repr.plot.height=8)\nfac <- st_as_sf(temp[temp$col_06_W5.1_C6Longitude < -40,],coords = c(\"col_06_W5.1_C6Longitude\",\"col_05_W5.1_C5Latitude\"), crs=4326)\ncolnames(fac)[c(10,13,15,17)] <- c(\"basin\",\"withdrawal\",\"discharge\",\"consumption\")\nfac$withdrawal <- as.numeric(fac$withdrawal)\nfac$discharge <- as.numeric(fac$discharge) \n# Merge KPI\nfac$kpi <- kpis_w$total[match(fac$organization, kpis_w$org)]\n# Aggregate water volumes by organization \nfac2 <- as.data.table(fac)[,.(consumption = mean(consumption, na.rm=TRUE),\n discharge = mean(discharge, na.rm=T),\n withdrawal = mean(withdrawal, na.rm=T),\n kpi=unique(kpi),\n sector=unique(sector)),.(organization)]\nfac2 <- melt(fac2, id.vars=c(\"organization\",\"kpi\",\"sector\"))\nfac2 <- fac2[is.finite(log(fac2$value)),] \n \n#labels \nfac2$variable <- as.character(fac2$variable) \nfac2$variable <- paste0(toupper(substr(fac2$variable,1,1)),\n substr(fac2$variable,2,nchar(fac2$variable))) \n \n \n\np <- ggplot(fac2, aes(y=log(value), x=kpi, org= organization, sector=sector))+\n geom_point(aes(colour=variable))+\n theme_bw()+\n labs(colour=\"\", y = \"Water volume\", x=\"Corporate KPI\") \n\npltly <- ggplotly(p, tooltip = c(\"org\",\"sector\"))\n \nhtmlwidgets::saveWidget(pltly, \"plot4.html\")\nIRdisplay::display_html('') ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see a clear trend between the KPI and the amount of water used. Companies with heavy water use may be more aware of the environmental and social dangers associated with their use intensity. I think it's pretty reassuring and it's going in the right direction. \n\n## a. Does corporates care about vulnerable population ?\n\nWe are going to explore this question by looking at the possible link between discharged waste water quantities and vulnerability index at facility location in the United States. \nFor this, I use EJscreen data for vulnerability index at census block scale and corporates response at question W5.1. \n\n* First, I extract facility coordinates\n* Spatially join census block vulnerability index at facility locations\n\nI will focus on General, Electric utilities and Food, beverage & tobacco sectors which are the most numerous in facility data. \nHere are the plots by sector. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Spatial join of census block\n\nfac_block <- suppressMessages(st_join(fac, shp))\nfac_block$discharge <- as.numeric(fac_block$discharge)\n\nfac_block$discharge <- log(fac_block$discharge)\n\nfac_block$kpi <- kpis_w$total[match(fac_block$organization, kpis_w$org)]\n\noptions(repr.plot.width=10,repr.plot.height=15) \nggplot(fac_block[is.finite(fac_block$discharge) & fac_block$sector %in% c(\"General\",\"Electric utilities\",\"Food, beverage & tobacco\"),], aes(x = discharge, y= VULEOPCT, colour=sector))+\ngeom_point(size = 6)+ \ngeom_smooth(method=\"lm\", formula=y~x)+ \ntheme_bw()+\nylim(c(0,1))+ \ntheme(strip.text.x = element_text(\n size = 16),\ntext = element_text(size=15),\n strip.background =element_rect(fill=\"#aec9eb\"),\n legend.position = \"none\")+\nlabs(y = \"Vunerability Index\", x=\"Corporate discharge\")+\nfacet_wrap(.~sector, scales=\"free\", ncol=1) \n ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There seems to be a decrease for Electric utilities sector, and an increase for General sector but nothing is clear. \nIt would be more interesting to look at company level\n\nI have selected three companies which informs a decent number of facility data: \n* Waste Management, Inc. (General sector)\n* Dominion Energy (Electric utilities)\n* Tyson Foods, Inc. (Food, beverage & tobacco)\n \nCorporates KPI is displayed in each plot title."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"corps <-c(\"Waste Management, Inc.\", # general\n \"Dominion Energy\", # electric\n \"Tyson Foods, Inc.\") # food\n \n \nfac_block <- fac_block[fac_block$organization %in% corps,]\nfac_block$organization <- paste0(fac_block$sector, \" - \", fac_block$organization, paste0(\" (KPI = \", round(fac_block$kpi, digits=2), \")\"))\ncolnames(fac_block)[16] <- \"Discharge_comparison\"\nfac_block$Discharge_comparison <- factor(fac_block$Discharge_comparison, levels = c(\"Much lower\",\"Lower\",\"About the same\",\"Higher\",\"Much higher\"))\n\noptions(repr.plot.width=15,repr.plot.height=15)\nggplot(fac_block[is.finite(fac_block$discharge),], aes(x=discharge, y=VULEOPCT))+\ngeom_point(aes(color=Discharge_comparison), size=6)+\ngeom_smooth(method=\"lm\", formula = y~x)+\nscale_color_brewer(palette = \"Spectral\")+\nylim(c(0,1))+\ntheme_bw()+\ntheme(strip.text.x = element_text(\n size = 16),\ntext = element_text(size=15),\n strip.background =element_rect(fill=\"#aec9eb\"))+\nlabs(y = \"Vunerability Index\", x=\"Corporate discharge\", color= \"Discharge comparison \\nwith previous year\")+\nfacet_wrap(.~organization, scales=\"free\", ncol=1) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This figure confirms the previous one, we see a drop in the volume of water discharged into the most vulnerable census blocks for Electric utilities. However, other sector have not the same trend. \n\nInterestingly, sustainability ambition KPI seems to correctly reflect corporates \"care\" about surrounding population of their facilities, by adapting discharged waste water quantities. A high KPI is related to a discharge decrease, whereas lower KPI corporates are increasing. \n\n## b. Is there a synergy between cities and corporates about sustainability ambition ?\nIn this last part, we're going to look at each facility location again. \nVery few facilities are located in a city which has participated in CDP survey, so we do the following : \n* compute distance matrix from facilities locations to CDP cities\n* assign nearest CDP city KPI to each facility\n* aggregate these KPIs by corporate\n\nAveraged KPI has been computed using this weight for each facility : \n\n$$ \\frac{log(consumption)}{log(distance)} $$ \n\nInspired by the [Gravity model of trade](https://en.wikipedia.org/wiki/Gravity_model_of_trade) ideas, a large water consumption of a company in a city must strengthen the congruence and cooperation between these two entities. On the contrary, the weight naturally decreases with the distance. \n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fac <- st_as_sf(temp[temp$col_06_W5.1_C6Longitude < -40,],coords = c(\"col_06_W5.1_C6Longitude\",\"col_05_W5.1_C5Latitude\"), crs=4326)\ncolnames(fac)[c(10,13,15,17)] <- c(\"basin\",\"withdrawal\",\"discharge\",\"consumption\")\nfac$id <- 1:nrow(fac)\ntest <- suppressMessages(st_join(fac, zip))\nst_geometry(test) <- NULL\nfac_cities <- unique(test$city)\nfac_cities <- fac_cities[!is.na(fac_cities)]\n\n# Distance\ndists <- st_distance(fac, pts_zip)\nfac$city <- test$city[match(fac$id, test$id)]\nfac$city_org <- apply(dists,1, function(x) pts_zip$org[which.min(x)])\nfac$dist <- apply(dists,1, function(x) min(x))\n\n# KPI corporate and cities\nfac$kpi <- kpis_w$total[match(fac$organization, kpis_w$org)]\nfac$kpi_city <- kpis$total[match(fac$city_org, kpis$org)] \n\n# ej\n\nej_city <- shp\nst_geometry(ej_city) <- NULL \nej_city <- setDT(ej_city)[,lapply(.SD, mean, na.rm=TRUE),.(city), .SDcols=vars] \n \n \nfac <- merge(fac, ej_city, by.x=\"city\", by.y=\"city\", all.x=TRUE) \n \noptions(repr.plot.width=12,repr.plot.height=15)\nfac2 <- setDT(fac)[,.(kpi_city = weighted.mean(kpi_city/log(dist) ,log(consumption),na.rm=TRUE),\n kpi = unique(kpi),\n sector=unique(sector)),.(organization)]\n\nsectors <- c(\"Electric utilities\",\"General\", \"Food, beverage & tobacco\")\n\nggplot(fac2[fac2$sector %in% sectors & is.finite(fac2$kpi_city),], aes(x=kpi, y=kpi_city, colour=sector))+\ngeom_point(size = 6)+\ngeom_smooth(method=\"lm\", formula = y~x)+\ntheme_bw()+\ntheme(strip.text.x = element_text(\n size = 16),\n text = element_text(size =15),\n strip.background =element_rect(fill=\"#aec9eb\"),\n legend.position = \"none\")+\nlabs(y = \"Averaged cities KPI\", x=\"Corporate KPI\")+\nfacet_wrap(.~sector, scales=\"free\", ncol=1) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For Electric utilities and General sectors, a synergy in sustainability ambition seems to exist between cities and corporates. \nThere may be sectors, like Food, beverage & tobacco, where public-private cooperation needs to be improved.\n# 4. Possible improvements & discussion\n \n"},{"metadata":{},"cell_type":"markdown","source":"## Possible improvements\n\nIn this report, there are many more features related to quantities completeness than text analysis. Intuitively, and especialy in climate change issues, *a number is worth a thousand words*, so I deliberately did a general lexical analysis, without distinguishing the questions or sections. \n\nHowever, much more features can be extracted, for example : \n \n* counting the number of added rows when the respondent can\n* selecting specific questions of interest\n* defining weights by kind of feature (numerical, text) or by topic\n* clustering after PCA, with K-means algorithm or Hierarchical Ascending Classification\n* doing deeper NLP analysis\n* URL or attached file counting\n \n**Using completeness features in supervised machine learning** \n\nOut of curiosity, we can compute logistic regression using compleness features against A score provided by CDP. \n* \"numprec_\" prefixe means quantity completeness\n* \"comp_\" prefixe means overall section completeness\n* \"wds_\" means lexical analysis \n\nI display only significant features :"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"temp <- data.frame(feat)\ntemp$score <- ifelse(temp$organization %in% cdp$Organization,1,0)\ncovariates <- colnames(temp)[!colnames(temp) %in% c(\"organization\",\"score\")]\n\ncolnames(temp) <- gsub(\"&\",\"\",colnames(temp))\nmodel <- glm(as.formula(paste(\"score ~ \", paste(covariates,collapse=\"+\"))),data= temp)\nresult <- as.data.frame(coefficients(summary(model)))\nresult <- result[result[,4] < 0.05,]\nresult <- result[order(result[,4], decreasing=FALSE),]\nresult","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We find some features related to the [CDP scoring methodology](https://guidance.cdp.net/en/guidance?cid=11&ctype=theme&idtype=ThemeID&incchild=1µsite=0&otype=ScoringModule&page=1) criterions like City-wide emissions attached data or Adaptation planning. \n\n\nNow, we can assess the predicting power of completeness features. I randomly split the data into train and test sets with respective proportions of 2/3 and 1/3, then predict A-score on test set using the model trained on the training set. \nLet's see confusion matrix, with predictions (rows) and labels (columns) : "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"set.seed(123)\nindex <- sample(nrow(temp),round(nrow(temp)/3))\ntrain <- temp[setdiff(1:nrow(temp), index),]\ntest <- temp[index,]\n\n# Logistic model\nmodel <- glm(as.formula(paste(\"score ~ \", paste(covariates,collapse=\"+\"))),data= train)\n\n# Predict \npreds <- predict(model,test)\npreds <- ifelse(preds > 0.5,1,0)\nconf <- as.data.frame.matrix(table(preds, test$score))\nconf\n\nprecision <- conf[2,2]/(conf[2,2] + conf[2,1])\nrecall <- conf[2,2]/(conf[2,2] + conf[1,2])\naccuracy <- (conf[2,2] + conf[1,1])/nrow(test)\n\ncat(paste(\"Precision :\",precision))\ncat(paste(\"\\nRecall :\",round(recall, digits = 2)))\ncat(paste(\"\\nAccuracy :\",round(accuracy, digits = 2)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Well, not so bad ! Of course, this can be improved by using more features or another ML algorithm. \n## Discussion\n\nDespite being a purely statistical indicator, which count missing values instead of using real answers, completeness KPI has proven to be reliable for discussing social and environmental issues for both cities and corporates. To put sustainability ambition KPI in perspective, we can compare **sustainability ambition** to **sustainability compliance** by looking at the global Sustainability Developpment Goal (SDG) Index for cities. SDG Index is a synthetic indicator composed of 17 subindicator, corresponding to the [17 goals](https://en.wikipedia.org/wiki/Sustainable_Development_Goals) set by the United Nations. \nBlack lines are the means on the two variables : \n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nus <- read.xlsx(\"../input/various-datasets-for-cdp-competition/2019USCitiesIndexResults.xlsx\",sheet=2)\neu <- read.xlsx(\"../input/various-datasets-for-cdp-competition/2019EUCitiesIndexResults.xlsx\", sheet = 2)\nsdg <- rbindlist(list(eu[,c(1,2,11:25)], us[,c(3,5,7:21)]),use.names = FALSE)\nkpis$city <- dis$City[match(kpis$org, dis$Organization)]\ntemp <- merge(kpis,sdg, by.x=\"city\", by.y=\"City\")\ntemp$region <- dis$`CDP Region`[match(temp$city, dis$City)]\n\np <-ggplot(temp, aes(x=total, y=(Index.Score), label=city, colour=region))+\n geom_point(size=3)+\n geom_vline(xintercept = mean(temp$total))+\n geom_hline(yintercept = mean(temp$Index.Score))+\n theme_bw()+\n labs(x = \"CDP ambition KPI\", y=\"SDG Index\", colour=\"\")\n\npltly <- ggplotly(p)\n\nhtmlwidgets::saveWidget(pltly, \"plot_sdg.html\")\nIRdisplay::display_html('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As said at the beginning of this report, ambition and compliance are not necessarly correlated. Interestingly, we aknowledge that European cities like Copenhagen, Berlin or Amsterdam are ahead in term of SDG but they are not so ambitious in the CDP questionnaire. One the other hand, US cities are generally below SDG Index mean but they prove to be more motivated in becoming greener. \nWe have to take in account culture difference. For example, USA have grown in a fast way with relative cheap energy in a huge land area. On the other hand, northern European countries like Sweden, has embarked on a sustainable development program since the 1960s, taking full advantage of its natural advantages with for example the exploitation of hydropower. \n\n# Thanks for reading !"}],"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"3.6.3"}},"nbformat":4,"nbformat_minor":4}