Load data

Datasets are available on the Dryad repository.

social_s<-read.csv("social_thresholds.csv", header=T)
limits_s<-read.csv("planet_boundaries.csv", header = T)

Social and biophysical scores and SEI in 2011

Raw social indicators

# get raw social data from standardised data from O'neill et al. 2018
# y raw = y std * (y threshold - y min) + y min

social_s$Life_Satisfaction<-(social_s$Life_Satisfaction*(6.5-2.93)+2.93)/10
social_s$Healthy_Life_Expect<-(social_s$Healthy_Life_Expect*(65-40)+40)/75
social_s$Nutrition[social_s$Nutrition>1]<-1
social_s$Sanitation<-(social_s$Sanitation*(95-6.2)+6.2)/100
social_s$Sanitation[social_s$Sanitation>1]<-1
social_s$Income<-(social_s$Income*(95-21.5)+21.5)/100
social_s$Income[social_s$Income>1]<-1
social_s$Access_to_Energy<-(social_s$Access_to_Energy*(95-6.5)+6.5)/100
social_s$Access_to_Energy[social_s$Access_to_Energy>1]<-1
social_s$Education<-(social_s$Education*(95-9.6)+9.6)/100
social_s$Education[social_s$Education>1]<-1
social_s$Social_Support<-(social_s$Social_Support*(90-30)+30)/100
social_s$Social_Support[social_s$Social_Support>1]<-1
social_s$Democratic_Quality<-((social_s$Democratic_Quality*(0.8+2.14)-2.14)+2.5)/5
social_s$Democratic_Quality[social_s$Democratic_Quality>1]<-1
social_s$Equality<-(social_s$Equality*(0.7-0.2)+0.2)/100
social_s$Equality[social_s$Equality>1]<-1
social_s$Employment<-(social_s$Employment*(94-68.6)+68.6)/100 
social_s$Employment[social_s$Employment>1]<-1

Social score and pre-treatment for biophysical score

# Social score (average and remove countries with not enough data)

social_s$score<-apply(social_s[,-1], 1, function(x){
  if(length(which(!is.na(x)))>5){1-mean(x,na.rm=T)}else{NA}})


# Biophysical score (prepare qualitative and quantitative information and remove countries with not enough data)

limits_s$score2<-apply(limits_s[,-1], 1, function(x){
  if(length(which(!is.na(x)))>3){length(which(!is.na(x) & x>=1))/length(which(!is.na(x)))}else{NA}})
limits_s$score3<-apply(limits_s[,-c(1,9)], 1, function(x){
  if(length(which(!is.na(x)))>3){if(length(which(x>=1))>0){mean(x[which(x>=1)], na.rm=T)-1}else{0}}else{NA}})

Biophysical score and socio-environmental index

# combining both datasets

soclim_sa<-merge(social_s[,c("Country","score")], limits_s[,c("Country","score2","score3")], by="Country", all=T)
soclim_sa<-na.omit(soclim_sa)

# calculate biophysical score and socio-environmental index

soclim_sa$score4<-rep(NA,nrow(soclim_sa))
soclim_sa$index_sa<-rep(NA,nrow(soclim_sa))

for(i in 1:nrow(soclim_sa)){
  y0=soclim_sa$score[i] # qualitatif social
  x0=soclim_sa$score2[i] # qualitatif limits
  k0=soclim_sa$score3[i] # quantatitif limits
  k0_b=soclim_sa$score3[i]/max(soclim_sa$score3) # normalize quantitative limits to 1
  
  x1=(x0+k0_b)/2 # x-coordinate of country's projection on z
  k1=(k0_b+x0)/2 # k-coordinate of country's projection on z
  l1=sqrt((x1-0)^2+(k1-0)^2) #  # distance between B(0,0,0) and country's projection on z
  lt=sqrt(2) # distance between B(0,0,0) and W(1,1,1)
  z0=l1/lt # normalize z to 1
  soclim_sa$score4[i]<-z0 # save z as the biophysical score
  
  l2=sqrt((y0-0)^2+(z0-0)^2) # distance between the country and B(0,0)
  soclim_sa$index_sa[i]=1-l2/sqrt(2) #  socio-environmental index (divided by sqrt(2) to normalize SEI to 1)
}

Social and biophysical scores and SEI over time

Raw social indicators

For details about raw data, see SM in O’neill et al. (2018) https://static-content.springer.com/esm/art%3A10.1038%2Fs41893-018-0021-4/MediaObjects/41893_2018_21_MOESM1_ESM.pdf

# life satisfaction : https://data.world/alamshukla/world-happiness-report-2019
life_sat<-read.csv("Chapter2OnlineData.csv")

# healthy life expectancy : https://apps.who.int/gho/data/node.main.HALE?lang=en

life_exp<-read.csv("Life_expectancy2.csv")
life_exp<-life_exp[,c("Country","X2005","X2000")]
life_exp<-melt(life_exp, id.vars="Country")
life_exp$Year<-substr(life_exp$variable,2,5)
life_exp$variable<-NULL
life_exp$Country<-as.character(life_exp$Country) # rename countries for consistency 
life_exp$Country[life_exp$Country=="Bolivia (Plurinational State of)"]<-"Bolivia"
life_exp$Country[life_exp$Country=="Democratic Republic of the Congo"]<-"Congo (Kinshasa)"
life_exp$Country[life_exp$Country=="Congo"]<-"Congo (Brazzaville)"
life_exp$Country[life_exp$Country=="Czechia"]<-"Czech Republic"
life_exp$Country[life_exp$Country=="Iran (Islamic Republic of)"]<-"Iran"
life_exp$Country[life_exp$Country=="Côte d'Ivoire"]<-"Ivory Coast"
life_exp$Country[life_exp$Country=="Lao People's Democratic Republic"]<-"Laos"
life_exp$Country[life_exp$Country=="Republic of North Macedonia"]<-"Macedonia"
life_exp$Country[life_exp$Country=="Republic of Moldova"]<-"Moldova"
life_exp$Country[life_exp$Country=="Russian Federation"]<-"Russia"
life_exp$Country[life_exp$Country=="Republic of Korea"]<-"South Korea"
life_exp$Country[life_exp$Country=="Eswatini"]<-"Swaziland"
life_exp$Country[life_exp$Country=="Syrian Arab Republic"]<-"Syria"
life_exp$Country[life_exp$Country=="United Republic of Tanzania"]<-"Tanzania"
life_exp$Country[life_exp$Country=="United Kingdom of Great Britain and Northern Ireland"]<-"United Kingdom"
life_exp$Country[life_exp$Country=="United States of America"]<-"United States"
life_exp$Country[life_exp$Country=="Venezuela (Bolivarian Republic of)"]<-"Venezuela"
life_exp$Country[life_exp$Country=="Viet Nam"]<-"Vietnam"
life_exp$Country<-as.factor(life_exp$Country)
names(life_exp)[1:2]<-c("Country.name","Healthy.life.expectancy.at.birth")
life_exp<-rbind(life_exp,life_sat[,c("Country.name","Healthy.life.expectancy.at.birth","Year")])

# nutrition : http://www.fao.org/faostat/en/#search/food%20supply
nut_old<-read.csv("FAOSTAT_data_old.csv")
nut_rec<-read.csv("FAOSTAT_data.csv")
nut_old<-nut_old[which(nut_old$Item=="Grand Total"),c("Area","Year","Value")]
nut_rec<-nut_rec[which(nut_rec$Item=="Grand Total"),c("Area","Year","Value")]
nutri<-rbind(nut_old,nut_rec)
nutri$Area<-as.character(nutri$Area) # rename countries for consistency 
nutri$Area[nutri$Area=="Bolivia (Plurinational State of)"]<-"Bolivia"
nutri$Area[nutri$Area=="Congo"]<-"Congo (Kinshasa)"
nutri$Area[nutri$Area=="Czechia"]<-"Czech Republic"
nutri$Area[nutri$Area=="China, Hong Kong SAR"]<-"Hong Kong S.A.R. of China"
nutri$Area[nutri$Area=="Iran (Islamic Republic of)"]<-"Iran"
nutri$Area[nutri$Area=="Côte d'Ivoire"]<-"Ivory Coast"
nutri$Area[nutri$Area=="Lao People's Democratic Republic"]<-"Laos"
nutri$Area[nutri$Area=="North Macedonia"]<-"Macedonia"
nutri$Area[nutri$Area=="Republic of Moldova"]<-"Moldova"
nutri$Area[nutri$Area=="Russian Federation"]<-"Russia"
nutri$Area[nutri$Area=="Republic of Korea"]<-"South Korea"
nutri$Area[nutri$Area=="Eswatini"]<-"Swaziland"
nutri$Area[nutri$Area=="China, Taiwan Province of"]<-"Taiwan Province of China"
nutri$Area[nutri$Area=="United Republic of Tanzania"]<-"Tanzania"
nutri$Area[nutri$Area=="United Kingdom of Great Britain and Northern Ireland"]<-"United Kingdom"
nutri$Area[nutri$Area=="United States of America"]<-"United States"
nutri$Area[nutri$Area=="Venezuela (Bolivarian Republic of)"]<-"Venezuela"
nutri$Area[nutri$Area=="Viet Nam"]<-"Vietnam"
nutri$Area<-as.factor(nutri$Area)


# sanitation : https://data.worldbank.org/indicator/SH.STA.SMSS.ZS

sanitation<-read.csv("Sanitation2.csv")
sanitation<-melt(sanitation)
sanitation$Year<-substr(sanitation$variable,2,5)
sanitation$variable<-NULL
sanitation$Country.Name<-as.character(sanitation$Country.Name) # rename countries for consistency 
sanitation$Country.Name[sanitation$Country.Name=="Congo, Dem. Rep."]<-"Congo (Kinshasa)"
sanitation$Country.Name[sanitation$Country.Name=="Congo, Rep."]<-"Congo (Brazzaville)"
sanitation$Country.Name[sanitation$Country.Name=="Egypt, Arab Rep."]<-"Egypt"
sanitation$Country.Name[sanitation$Country.Name=="Iran, Islamic Rep."]<-"Iran"
sanitation$Country.Name[sanitation$Country.Name=="Cote d'Ivoire"]<-"Ivory Coast"
sanitation$Country.Name[sanitation$Country.Name=="Gambia, The"]<-"Gambia"
sanitation$Country.Name[sanitation$Country.Name=="Hong Kong SAR, China"]<-"Hong Kong S.A.R. of China"
sanitation$Country.Name[sanitation$Country.Name=="Kyrgyz Republic"]<-"Kyrgyzstan"
sanitation$Country.Name[sanitation$Country.Name=="Lao PDR"]<-"Laos"
sanitation$Country.Name[sanitation$Country.Name=="North Macedonia"]<-"Macedonia"
sanitation$Country.Name[sanitation$Country.Name=="Slovak Republic"]<-"Slovakia"
sanitation$Country.Name[sanitation$Country.Name=="Russian Federation"]<-"Russia"
sanitation$Country.Name[sanitation$Country.Name=="Korea, Rep"]<-"South Korea"
sanitation$Country.Name[sanitation$Country.Name=="Eswatini"]<-"Swaziland"
sanitation$Country.Name[sanitation$Country.Name=="Syrian Arab Republic"]<-"Syria"
sanitation$Country.Name[sanitation$Country.Name=="Yemen, Rep."]<-"Yemen"
sanitation$Country.Name[sanitation$Country.Name=="Bahamas, The"]<-"Bahamas"
sanitation$Country.Name[sanitation$Country.Name=="Korea, Dem. People’s Rep."]<-"Democratic People's Republic of Korea"
sanitation$Country.Name[sanitation$Country.Name=="Venezuela, RB"]<-"Venezuela"
sanitation$Country.Name[sanitation$Country.Name=="Micronesia, Fed. Sts."]<-"Micronesia (Federated States of)"
sanitation$Country.Name[sanitation$Country.Name=="St. Lucia"]<-"Saint Lucia"
sanitation$Country.Name[sanitation$Country.Name=="St. Vincent and the Grenadines"]<-"Saint Vincent and the Grenadines"
sanitation$Country.Name<-as.factor(sanitation$Country.Name)

# income : https://data.worldbank.org/indicator/SI.POV.DDAY?view=chart

income<-read.csv("Poverty2.csv")
income<-melt(income)
income$Year<-substr(income$variable,2,5)
income$variable<-NULL
income$Country.Name<-as.character(income$Country.Name) # rename countries for consistency 
income$Country.Name[income$Country.Name=="Congo, Dem. Rep."]<-"Congo (Kinshasa)"
income$Country.Name[income$Country.Name=="Congo, Rep."]<-"Congo (Brazzaville)"
income$Country.Name[income$Country.Name=="Egypt, Arab Rep."]<-"Egypt"
income$Country.Name[income$Country.Name=="Iran, Islamic Rep."]<-"Iran"
income$Country.Name[income$Country.Name=="Cote d'Ivoire"]<-"Ivory Coast"
income$Country.Name[income$Country.Name=="Gambia, The"]<-"Gambia"
income$Country.Name[income$Country.Name=="Hong Kong SAR, China"]<-"Hong Kong S.A.R. of China"
income$Country.Name[income$Country.Name=="Kyrgyz Republic"]<-"Kyrgyzstan"
income$Country.Name[income$Country.Name=="Lao PDR"]<-"Laos"
income$Country.Name[income$Country.Name=="North Macedonia"]<-"Macedonia"
income$Country.Name[income$Country.Name=="Slovak Republic"]<-"Slovakia"
income$Country.Name[income$Country.Name=="Russian Federation"]<-"Russia"
income$Country.Name[income$Country.Name=="Korea, Rep"]<-"South Korea"
income$Country.Name[income$Country.Name=="Eswatini"]<-"Swaziland"
income$Country.Name[income$Country.Name=="Syrian Arab Republic"]<-"Syria"
income$Country.Name[income$Country.Name=="Yemen, Rep."]<-"Yemen"
income$Country.Name[income$Country.Name=="Bahamas, The"]<-"Bahamas"
income$Country.Name[income$Country.Name=="Korea, Dem. People’s Rep."]<-"Democratic People's Republic of Korea"
income$Country.Name[income$Country.Name=="Venezuela, RB"]<-"Venezuela"
income$Country.Name[income$Country.Name=="Micronesia, Fed. Sts."]<-"Micronesia (Federated States of)"
income$Country.Name[income$Country.Name=="St. Lucia"]<-"Saint Lucia"
income$Country.Name[income$Country.Name=="St. Vincent and the Grenadines"]<-"Saint Vincent and the Grenadines"
income$Country.Name<-as.factor(income$Country.Name)

# energy access : https://data.worldbank.org/indicator/EG.ELC.ACCS.ZS?view=chart

ener_access<-read.csv("Access_electricity2.csv")
ener_access<-melt(ener_access)
ener_access$Year<-substr(ener_access$variable,2,5)
ener_access$variable<-NULL
ener_access$Country.Name<-as.character(ener_access$Country.Name) # rename countries for consistency 
ener_access$Country.Name[ener_access$Country.Name=="Congo, Dem. Rep."]<-"Congo (Kinshasa)"
ener_access$Country.Name[ener_access$Country.Name=="Congo, Rep."]<-"Congo (Brazzaville)"
ener_access$Country.Name[ener_access$Country.Name=="Egypt, Arab Rep."]<-"Egypt"
ener_access$Country.Name[ener_access$Country.Name=="Iran, Islamic Rep."]<-"Iran"
ener_access$Country.Name[ener_access$Country.Name=="Cote d'Ivoire"]<-"Ivory Coast"
ener_access$Country.Name[ener_access$Country.Name=="Gambia, The"]<-"Gambia"
ener_access$Country.Name[ener_access$Country.Name=="Hong Kong SAR, China"]<-"Hong Kong S.A.R. of China"
ener_access$Country.Name[ener_access$Country.Name=="Kyrgyz Republic"]<-"Kyrgyzstan"
ener_access$Country.Name[ener_access$Country.Name=="Lao PDR"]<-"Laos"
ener_access$Country.Name[ener_access$Country.Name=="North Macedonia"]<-"Macedonia"
ener_access$Country.Name[ener_access$Country.Name=="Slovak Republic"]<-"Slovakia"
ener_access$Country.Name[ener_access$Country.Name=="Russian Federation"]<-"Russia"
ener_access$Country.Name[ener_access$Country.Name=="Korea, Rep"]<-"South Korea"
ener_access$Country.Name[ener_access$Country.Name=="Eswatini"]<-"Swaziland"
ener_access$Country.Name[ener_access$Country.Name=="Syrian Arab Republic"]<-"Syria"
ener_access$Country.Name[ener_access$Country.Name=="Yemen, Rep."]<-"Yemen"
ener_access$Country.Name[ener_access$Country.Name=="Bahamas, The"]<-"Bahamas"
ener_access$Country.Name[ener_access$Country.Name=="Korea, Dem. People’s Rep."]<-"Democratic People's Republic of Korea"
ener_access$Country.Name[ener_access$Country.Name=="Venezuela, RB"]<-"Venezuela"
ener_access$Country.Name[ener_access$Country.Name=="Micronesia, Fed. Sts."]<-"Micronesia (Federated States of)"
ener_access$Country.Name[ener_access$Country.Name=="St. Lucia"]<-"Saint Lucia"
ener_access$Country.Name[ener_access$Country.Name=="St. Vincent and the Grenadines"]<-"Saint Vincent and the Grenadines"
ener_access$Country.Name<-as.factor(ener_access$Country.Name)

# education : https://data.worldbank.org/indicator/SE.SEC.ENRR?view=chart

educ<-read.csv("Education2.csv")
educ<-melt(educ)
educ$Year<-substr(educ$variable,2,5)
educ$variable<-NULL
educ$Country.Name<-as.character(educ$Country.Name) # rename countries for consistency 
educ$Country.Name[educ$Country.Name=="Congo, Dem. Rep."]<-"Congo (Kinshasa)"
educ$Country.Name[educ$Country.Name=="Congo, Rep."]<-"Congo (Brazzaville)"
educ$Country.Name[educ$Country.Name=="Egypt, Arab Rep."]<-"Egypt"
educ$Country.Name[educ$Country.Name=="Iran, Islamic Rep."]<-"Iran"
educ$Country.Name[educ$Country.Name=="Cote d'Ivoire"]<-"Ivory Coast"
educ$Country.Name[educ$Country.Name=="Gambia, The"]<-"Gambia"
educ$Country.Name[educ$Country.Name=="Hong Kong SAR, China"]<-"Hong Kong S.A.R. of China"
educ$Country.Name[educ$Country.Name=="Kyrgyz Republic"]<-"Kyrgyzstan"
educ$Country.Name[educ$Country.Name=="Lao PDR"]<-"Laos"
educ$Country.Name[educ$Country.Name=="North Macedonia"]<-"Macedonia"
educ$Country.Name[educ$Country.Name=="Slovak Republic"]<-"Slovakia"
educ$Country.Name[educ$Country.Name=="Russian Federation"]<-"Russia"
educ$Country.Name[educ$Country.Name=="Korea, Rep"]<-"South Korea"
educ$Country.Name[educ$Country.Name=="Eswatini"]<-"Swaziland"
educ$Country.Name[educ$Country.Name=="Syrian Arab Republic"]<-"Syria"
educ$Country.Name[educ$Country.Name=="Yemen, Rep."]<-"Yemen"
educ$Country.Name[educ$Country.Name=="Bahamas, The"]<-"Bahamas"
educ$Country.Name[educ$Country.Name=="Korea, Dem. People’s Rep."]<-"Democratic People's Republic of Korea"
educ$Country.Name[educ$Country.Name=="Venezuela, RB"]<-"Venezuela"
educ$Country.Name[educ$Country.Name=="Micronesia, Fed. Sts."]<-"Micronesia (Federated States of)"
educ$Country.Name[educ$Country.Name=="St. Lucia"]<-"Saint Lucia"
educ$Country.Name[educ$Country.Name=="St. Vincent and the Grenadines"]<-"Saint Vincent and the Grenadines"
educ$Country.Name<-as.factor(educ$Country.Name)

# social support : https://data.world/alamshukla/world-happiness-report-2019
#social_sup<-read.csv("Chapter2OnlineData.csv")

# democratic quality : https://data.world/alamshukla/world-happiness-report-2019
#dem_qual<-read.csv("Chapter2OnlineData.csv")

# equality : https://data.world/alamshukla/world-happiness-report-2019
#equality<-read.csv("Chapter2OnlineData.csv")

# employment : https://data.worldbank.org/indicator/SL.UEM.TOTL.NE.ZS

employ<-read.csv("Unemployment2.csv")
employ<-melt(employ)
employ$Year<-substr(employ$variable,2,5)
employ$variable<-NULL
employ$value<-100-employ$value
employ$Country.Name<-as.character(employ$Country.Name) # rename countries for consistency 
employ$Country.Name[employ$Country.Name=="Congo, Dem. Rep."]<-"Congo (Kinshasa)"
employ$Country.Name[employ$Country.Name=="Congo, Rep."]<-"Congo (Brazzaville)"
employ$Country.Name[employ$Country.Name=="Egypt, Arab Rep."]<-"Egypt"
employ$Country.Name[employ$Country.Name=="Iran, Islamic Rep."]<-"Iran"
employ$Country.Name[employ$Country.Name=="Cote d'Ivoire"]<-"Ivory Coast"
employ$Country.Name[employ$Country.Name=="Gambia, The"]<-"Gambia"
employ$Country.Name[employ$Country.Name=="Hong Kong SAR, China"]<-"Hong Kong S.A.R. of China"
employ$Country.Name[employ$Country.Name=="Kyrgyz Republic"]<-"Kyrgyzstan"
employ$Country.Name[employ$Country.Name=="Lao PDR"]<-"Laos"
employ$Country.Name[employ$Country.Name=="North Macedonia"]<-"Macedonia"
employ$Country.Name[employ$Country.Name=="Slovak Republic"]<-"Slovakia"
employ$Country.Name[employ$Country.Name=="Russian Federation"]<-"Russia"
employ$Country.Name[employ$Country.Name=="Korea, Rep"]<-"South Korea"
employ$Country.Name[employ$Country.Name=="Eswatini"]<-"Swaziland"
employ$Country.Name[employ$Country.Name=="Syrian Arab Republic"]<-"Syria"
employ$Country.Name[employ$Country.Name=="Yemen, Rep."]<-"Yemen"
employ$Country.Name[employ$Country.Name=="Bahamas, The"]<-"Bahamas"
employ$Country.Name[employ$Country.Name=="Korea, Dem. People’s Rep."]<-"Democratic People's Republic of Korea"
employ$Country.Name[employ$Country.Name=="Venezuela, RB"]<-"Venezuela"
employ$Country.Name[employ$Country.Name=="Micronesia, Fed. Sts."]<-"Micronesia (Federated States of)"
employ$Country.Name[employ$Country.Name=="St. Lucia"]<-"Saint Lucia"
employ$Country.Name[employ$Country.Name=="St. Vincent and the Grenadines"]<-"Saint Vincent and the Grenadines"
employ$Country.Name<-as.factor(employ$Country.Name)


social_time<-life_sat[,c("Country.name","Year","Life.Ladder","Social.support","Democratic.Quality","gini.of.household.income.reported.in.Gallup..by.wp5.year")]
social_time<-merge(social_time,life_exp, by=c("Country.name","Year"),all=T)
social_time$Healthy.life.expectancy.at.birth<-as.numeric(social_time$Healthy.life.expectancy.at.birth)
social_time$gini.of.household.income.reported.in.Gallup..by.wp5.year<-1-social_time$gini.of.household.income.reported.in.Gallup..by.wp5.year
social_time<-merge(social_time,nutri, by.x=c("Country.name","Year"),by.y=c("Area","Year"),all.x=T)
names(social_time)[8]<-"nutri"
social_time<-merge(social_time,sanitation, by.x=c("Country.name","Year"),by.y=c("Country.Name","Year"),all.x=T)
names(social_time)[9]<-"sanitation"
social_time<-merge(social_time,income, by.x=c("Country.name","Year"),by.y=c("Country.Name","Year"),all.x=T)
names(social_time)[10]<-"income"
social_time$income<-100-social_time$income
social_time<-merge(social_time,ener_access, by.x=c("Country.name","Year"),by.y=c("Country.Name","Year"),all.x=T)
names(social_time)[11]<-"energy"
social_time<-merge(social_time,educ, by.x=c("Country.name","Year"),by.y=c("Country.Name","Year"),all.x=T)
names(social_time)[12]<-"educ"
social_time<-merge(social_time,employ, by.x=c("Country.name","Year"),by.y=c("Country.Name","Year"),all.x=T)
names(social_time)[13]<-"employ"

# standardize as in O'neill et al. 2018 to check for consistency

soc_threshold<-data.frame("happy"=c(6.5),"life_exp"=c(65),"nutrition"=c(2700),"sanitation"=c(95),"income"=c(95),"energy"=c(95),
                          "educ"=c(95),"social_supp"=c(0.9),"demo_qual"=c(0.8),"equal"=c(0.7),"employ"=c(94))

social_time2<-social_time[,c("Country.name","Year")]
social_time2$Life_Satisfaction<-(social_time$Life.Ladder-min(social_time$Life.Ladder, na.rm=T))/(soc_threshold[[1]]-min(social_time$Life.Ladder, na.rm=T))
social_time2$Healthy_Life_Expect<-(social_time$Healthy.life.expectancy.at.birth-min(social_time$Healthy.life.expectancy.at.birth , na.rm=T))/(soc_threshold[[2]]-min(social_time$Healthy.life.expectancy.at.birth , na.rm=T))
social_time2$Nutrition<-(social_time$nutri-min(social_time$nutri, na.rm=T))/(soc_threshold[[3]]-min(social_time$nutri, na.rm=T))
social_time2$Sanitation<-(social_time$sanitation-min(social_time$sanitation, na.rm=T))/(soc_threshold[[4]]-min(social_time$sanitation, na.rm=T))
social_time2$Income<-(social_time$income-min(social_time$income, na.rm=T))/(soc_threshold[[5]]-min(social_time$income, na.rm=T))
social_time2$Access_to_Energy<-(social_time$energy-min(social_time$energy, na.rm=T))/(soc_threshold[[6]]-min(social_time$energy, na.rm=T))
social_time2$Education<-(social_time$educ-min(social_time$educ, na.rm=T))/(soc_threshold[[7]]-min(social_time$educ, na.rm=T))
social_time2$Social_Support<-(social_time$Social.support-min(social_time$Social.support, na.rm=T))/(soc_threshold[[8]]-min(social_time$Social.support, na.rm=T))
social_time2$Democratic_Quality<-(social_time$Democratic.Quality-min(social_time$Democratic.Quality, na.rm=T))/(soc_threshold[[9]]-min(social_time$Democratic.Quality, na.rm=T))
social_time2$Equality<-(social_time$gini.of.household.income.reported.in.Gallup..by.wp5.year-min(social_time$gini.of.household.income.reported.in.Gallup..by.wp5.year, na.rm=T))/(soc_threshold[[10]]-min(social_time$gini.of.household.income.reported.in.Gallup..by.wp5.year, na.rm=T))
social_time2$Employment<-(social_time$employ-min(social_time$employ, na.rm=T))/(soc_threshold[[11]]-min(social_time$employ, na.rm=T))

social_time3<-droplevels(subset(social_time2, Year!=2018))
social_time3$group<-rep(NA,nrow(social_time3))
social_time3$Year<-as.numeric(social_time3$Year)
social_time3$group[which(social_time3$Year %in% c(2003:2007))]<-2005
social_time3$group[which(social_time3$Year %in% c(2008:2012))]<-2010
social_time3$group[which(social_time3$Year %in% c(2013:2017))]<-2015
social_time3<-as.data.frame(social_time3 %>% group_by(Country.name,group) %>% summarize(Life_Satisfaction=mean(Life_Satisfaction,na.rm=T), Healthy_Life_Expect=mean(Healthy_Life_Expect,na.rm=T), Nutrition=mean(Nutrition,na.rm=T), Sanitation=mean(Sanitation,na.rm=T), Income=mean(Income,na.rm=T), Access_to_Energy=mean(Access_to_Energy,na.rm=T), Education=mean(Education,na.rm=T), Social_Support=mean(Social_Support,na.rm=T), Democratic_Quality=mean(Democratic_Quality,na.rm=T), Equality=mean(Equality,na.rm=T), Employment=mean(Employment,na.rm=T)))

# use raw data and maxima to compute the social score

social_time_init<-social_time3
social_s2<-social_time_init
social_s2$Life_Satisfaction<-(social_s2$Life_Satisfaction*(6.5-2.93)+2.93)/10
social_s2$Healthy_Life_Expect<-(social_s2$Healthy_Life_Expect*(65-40)+40)/75
social_s2$Nutrition[social_s2$Nutrition>1]<-1
social_s2$Sanitation<-(social_s2$Sanitation*(95-6.2)+6.2)/100
social_s2$Sanitation[social_s2$Sanitation>1]<-1
social_s2$Income<-(social_s2$Income*(95-21.5)+21.5)/100
social_s2$Income[social_s2$Income>1]<-1
social_s2$Access_to_Energy<-(social_s2$Access_to_Energy*(95-6.5)+6.5)/100
social_s2$Access_to_Energy[social_s2$Access_to_Energy>1]<-1
social_s2$Education<-(social_s2$Education*(95-9.6)+9.6)/100
social_s2$Education[social_s2$Education>1]<-1
social_s2$Social_Support<-(social_s2$Social_Support*(90-30)+30)/100
social_s2$Social_Support[social_s2$Social_Support>1]<-1
social_s2$Democratic_Quality<-((social_s2$Democratic_Quality*(0.8+2.14)-2.14)+2.5)/5
social_s2$Democratic_Quality[social_s2$Democratic_Quality>1]<-1
social_s2$Equality<-(social_s2$Equality*(0.7-0.2)+0.2)
social_s2$Equality[social_s2$Equality>1]<-1
social_s2$Employment<-(social_s2$Employment*(94-68.6)+68.6)/100 
social_s2$Employment[social_s2$Employment>1]<-1

Raw biophysical indicators

# climate change : https://worldmrio.com/footprints/carbon/

clim_chg<-read.delim("cba.txt")
clim_chg[,c("X")]<-NULL
clim_chg_long<-melt(clim_chg)
clim_chg_long$year<-substr(clim_chg_long$variable,2,5)
clim_chg_long$variable<-NULL

# land-system change : https://www.pnas.org/content/110/25/10324 and https://onlinelibrary.wiley.com/doi/full/10.1111/jiec.12238

country_land_use<-read.csv("country_hanpp.csv") # reconstituted from Krausmann suppl
country_land_use_2005<-read.csv("planet_boundaries.csv", header = T)
land_use<-merge(country_land_use, country_land_use_2005[,c("Country","eHANPP")], by="Country", all.y=T)
year_land_use<-read.csv("year_hanpp.csv")
land_use$y2000<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y2000[i]<-land_use[i,3]*year_land_use$y2000[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1990<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1990[i]<-land_use[i,3]*year_land_use$y1990[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1980<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1980[i]<-land_use[i,3]*year_land_use$y1980[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1970<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1970[i]<-land_use[i,3]*year_land_use$y1970[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1960<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1960[i]<-land_use[i,3]*year_land_use$y1960[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1950<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1950[i]<-land_use[i,3]*year_land_use$y1950[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1930<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1930[i]<-land_use[i,3]*year_land_use$y1930[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
land_use$y1910<-rep(NA, nrow(land_use))
for(i in 1:nrow(land_use)){
  land_use$y1910[i]<-land_use[i,3]*year_land_use$y1910[which(year_land_use$Region %in% land_use[i,2])]/year_land_use$y2005[which(year_land_use$Region %in% land_use[i,2])]}
names(land_use)[3]<-"y2007"
land_use$y1991<-1*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1992<-2*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1993<-3*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1994<-4*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1995<-5*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1996<-6*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1997<-7*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1998<-8*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y1999<-9*(land_use$y2000-land_use$y1990)/10+land_use$y1990
land_use$y2001<-1*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2002<-2*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2003<-3*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2004<-4*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2005<-5*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2006<-6*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2008<-8*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2009<-9*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2010<-10*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2011<-11*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2012<-12*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2013<-13*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2014<-14*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2015<-15*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2016<-16*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2017<-17*(land_use$y2007-land_use$y2000)/7+land_use$y2000
land_use$y2018<-18*(land_use$y2007-land_use$y2000)/7+land_use$y2000

land_use_long<-melt(land_use)
land_use_long$year<-substr(land_use_long$variable,2,5)
land_use_long$variable<-NULL

land_use_long$Country<-as.character(land_use_long$Country) # rename countries for consistency 
land_use_long$Country<-gsub("_", " ", land_use_long$Country)
land_use_long$Country[land_use_long$Country=="Côte d'Ivoire"]<-"Cote dIvoire"
land_use_long$Country[land_use_long$Country=="West Bank and Gaza"]<-"Gaza Strip"
land_use_long$Country[land_use_long$Country=="Kyrgyz Republic"]<-"Kyrgyzstan"
land_use_long$Country[land_use_long$Country=="Lao PDR"]<-"Laos"
land_use_long$Country[land_use_long$Country=="Slovak Republic"]<-"Slovakia"
land_use_long$Country[land_use_long$Country=="Korea"]<-"South Korea"
land_use_long$Country[land_use_long$Country=="Syrian Arab Republic"]<-"Syria"
land_use_long$Country[land_use_long$Country=="Macedonia"]<-"TFYR Macedonia"
land_use_long$Country[land_use_long$Country=="United Kingdom"]<-"UK"
land_use_long$Country[land_use_long$Country=="United States"]<-"USA"
land_use_long$Country[land_use_long$Country=="United Arab Emirates"]<-"UAE"
land_use_long$Country[land_use_long$Country=="Vietnam"]<-"Viet Nam"
land_use_long$Country<-as.factor(land_use_long$Country)

# ecological footprint : https://data.world/footprint/nfa-2019-edition

ecol_foot<-read.csv("NFA_2019_public_data.csv")
ecol_foot<-droplevels(subset(ecol_foot[which(ecol_foot$record=="EFConsPerCap"),]))
ecol_foot$country<-as.character(ecol_foot$country) # rename countries for consistency 
ecol_foot$country[ecol_foot$country=="Antigua and Barbuda"]<-"Antigua"
ecol_foot$country[ecol_foot$country=="Brunei Darussalam"]<-"Brunei"
ecol_foot$country[ecol_foot$country=="Congo, Democratic Republic of"]<-"DR Congo"
ecol_foot$country[ecol_foot$country=="Iran, Islamic Republic of"]<-"Iran"
ecol_foot$country[ecol_foot$country=="Côte d'Ivoire"]<-"Cote dIvoire"
ecol_foot$country[ecol_foot$country=="Lao People's Democratic Republic"]<-"Laos"
ecol_foot$country[ecol_foot$country=="Libyan Arab Jamahiriya"]<-"Libya"
ecol_foot$country[ecol_foot$country=="Slovak Republic"]<-"Slovakia"
ecol_foot$country[ecol_foot$country=="Russian Federation"]<-"Russia"
ecol_foot$country[ecol_foot$country=="Korea, Republic of"]<-"South Korea"
ecol_foot$country[ecol_foot$country=="Eswatini"]<-"Swaziland"
ecol_foot$country[ecol_foot$country=="Syrian Arab Republic"]<-"Syria"
ecol_foot$country[ecol_foot$country=="Korea, Democratic People's Republic of"]<-"North Korea"
ecol_foot$country[ecol_foot$country=="Venezuela, Bolivarian Republic of"]<-"Venezuela"
ecol_foot$country[ecol_foot$country=="Tanzania, United Republic of"]<-"Tanzania"
ecol_foot$country[ecol_foot$country=="Macedonia TFYR"]<-"TFYR Macedonia"
ecol_foot$country[ecol_foot$country=="United Kingdom"]<-"UK"
ecol_foot$country[ecol_foot$country=="United States of America"]<-"USA"
ecol_foot$country[ecol_foot$country=="United Arab Emirates"]<-"UAE"
ecol_foot$country<-as.factor(ecol_foot$country)

# material footprint : https://worldmrio.com/footprints/

bioch_wat_mat<-rbind(read.delim("tradereport_1970.txt"),read.delim("tradereport_1971.txt"),
                     read.delim("tradereport_1972.txt"),read.delim("tradereport_1973.txt"),
                     read.delim("tradereport_1974.txt"),read.delim("tradereport_1975.txt"),
                     read.delim("tradereport_1976.txt"),read.delim("tradereport_1977.txt"),
                     read.delim("tradereport_1978.txt"),read.delim("tradereport_1979.txt"),
                     read.delim("tradereport_1980.txt"),read.delim("tradereport_1981.txt"),
                     read.delim("tradereport_1982.txt"),read.delim("tradereport_1983.txt"),
                     read.delim("tradereport_1984.txt"),read.delim("tradereport_1985.txt"),
                     read.delim("tradereport_1986.txt"),read.delim("tradereport_1987.txt"),
                     read.delim("tradereport_1988.txt"),read.delim("tradereport_1989.txt"),
                     read.delim("tradereport_1990.txt"),read.delim("tradereport_1991.txt"),
                     read.delim("tradereport_1992.txt"),read.delim("tradereport_1993.txt"),
                     read.delim("tradereport_1994.txt"),read.delim("tradereport_1995.txt"),
                     read.delim("tradereport_1996.txt"),read.delim("tradereport_1997.txt"),
                     read.delim("tradereport_1998.txt"),read.delim("tradereport_1999.txt"),
                     read.delim("tradereport_2000.txt"),read.delim("tradereport_2001.txt"),
                     read.delim("tradereport_2002.txt"),read.delim("tradereport_2003.txt"),
                     read.delim("tradereport_2004.txt"),read.delim("tradereport_2005.txt"),
                     read.delim("tradereport_2006.txt"),read.delim("tradereport_2007.txt"),
                     read.delim("tradereport_2008.txt"),read.delim("tradereport_2009.txt"),
                     read.delim("tradereport_2010.txt"),read.delim("tradereport_2011.txt"),
                     read.delim("tradereport_2012.txt"),read.delim("tradereport_2013.txt"),
                     read.delim("tradereport_2014.txt"),read.delim("tradereport_2015.txt"),
                     read.delim("tradereport_2016.txt"),read.delim("tradereport_2017.txt"))
# e.g. https://worldmrio.com/ComputationsM/Phase199/Loop082/leontief/tradereport_2015.txt for 2015
bioch_mat<-droplevels(subset(bioch_wat_mat[which(bioch_wat_mat$Record=="Footprint"),], 
                                 Indicator.Description %in% c("N emissions from fertilizer and manure (kg)",
                                                              "P emissions from fertilizer and manure (kg)",
                                                              "Raw material inputs, itemized (t)")))

wat_foot<-droplevels(subset(bioch_wat_mat[which(bioch_wat_mat$Record=="Footprint"),], 
                                 Indicator.Description %in% c("WFN: Water footprint of domestic water supply (Mm3/yr) - Blue",
                                                              "WFN: Water footprint of domestic water supply (Mm3/yr) - Grey")))

wat_foot<-as.data.frame(wat_foot %>% group_by(Country.Name,Year) %>% summarize(value=sum(Value)))

biophy_time<-droplevels(subset(clim_chg_long[which(clim_chg_long$Record=="Population"),]))
biophy_time$Record<-NULL
names(biophy_time)[2]<-"population"
biophy_time<-merge(biophy_time,droplevels(subset(clim_chg_long[which(clim_chg_long$Record=="CBA_MtCO2perCap"),]))[,c(1,3,4)], by=c("Country","year"))
names(biophy_time)[4]<-"CBA_MtCO2perCap"
biophy_time<-merge(biophy_time,land_use_long[,c(1,3,4)], by=c("Country","year"),all.x=T)
names(biophy_time)[5]<-"eHANPP"
biophy_time<-merge(biophy_time,ecol_foot[,c("country","year","total")], by.x=c("Country","year"), by.y=c("country","year"),all.x=T)
names(biophy_time)[6]<-"EFConsPerCap"
biophy_time<-merge(biophy_time,droplevels(subset(bioch_mat, Indicator.Description=="N emissions from fertilizer and manure (kg)"))[,c("Country.Name","Year","Value")],
                   by.x=c("Country","year"), by.y=c("Country.Name","Year"),all.x=T)
names(biophy_time)[7]<-"N_emissions"
biophy_time<-merge(biophy_time,droplevels(subset(bioch_mat, Indicator.Description=="P emissions from fertilizer and manure (kg)"))[,c("Country.Name","Year","Value")],
                   by.x=c("Country","year"), by.y=c("Country.Name","Year"),all.x=T)
names(biophy_time)[8]<-"P_emissions"
biophy_time<-merge(biophy_time,wat_foot,by.x=c("Country","year"), by.y=c("Country.Name","Year"),all.x=T)
names(biophy_time)[9]<-"Water_use"
biophy_time<-merge(biophy_time,droplevels(subset(bioch_mat, Indicator.Description=="Raw material inputs, itemized (t)"))[,c("Country.Name","Year","Value")],
                   by.x=c("Country","year"), by.y=c("Country.Name","Year"),all.x=T)
names(biophy_time)[10]<-"Material_footprint"
biophy_time$N_emissions_perCap<-biophy_time$N_emissions/biophy_time$population
biophy_time$P_emissions_perCap<-biophy_time$P_emissions/biophy_time$population
biophy_time$Water_use_perCap<-biophy_time$Water_use/biophy_time$population*1000000
biophy_time$Material_footprint_perCap<-biophy_time$Material_footprint/biophy_time$population

bio_boundaries<-data.frame("co2 (t)"=c(1.61),"P (kg)"=c(0.89),"N (kg)"=c(8.9),
                           "water (m3)"=c(574),"ecol_footprint (ha)"=c(1.72),"material_footprint (t)"=c(7.2))#eHANPP already as a limit
biophy_time2<-biophy_time[,c("Country","year")]
biophy_time2$CO2_Emissions<-biophy_time$CBA_MtCO2perCap/bio_boundaries[[1]]
biophy_time2$Phosphorus<-biophy_time$P_emissions_perCap/bio_boundaries[[2]]
biophy_time2$Nitrogen<-biophy_time$N_emissions_perCap/bio_boundaries[[3]]
biophy_time2$Blue_Water<-biophy_time$Water_use_perCap/bio_boundaries[[4]]
biophy_time2$eHANPP<-biophy_time$eHANPP
biophy_time2$Ecological_Footprint<-biophy_time$EFConsPerCap/bio_boundaries[[5]]
biophy_time2$Material_Footprint<-biophy_time$Material_footprint_perCap/bio_boundaries[[6]]

biophy_time3<-droplevels(subset(biophy_time2, year %in% c(1998:2018)))
biophy_time3$group<-rep(NA,nrow(biophy_time3))
biophy_time3$year<-as.numeric(biophy_time3$year)
biophy_time3$group[which(biophy_time3$year %in% c(1998:2002))]<-2000
biophy_time3$group[which(biophy_time3$year %in% c(2003:2007))]<-2005
biophy_time3$group[which(biophy_time3$year %in% c(2008:2012))]<-2010
biophy_time3$group[which(biophy_time3$year %in% c(2013:2017))]<-2015
biophy_time3<-as.data.frame(biophy_time3 %>% group_by(Country,group) %>% summarize(CO2_Emissions=mean(CO2_Emissions,na.rm=T), Phosphorus=mean(Phosphorus,na.rm=T), Nitrogen=mean(Nitrogen,na.rm=T), Blue_Water=mean(Blue_Water,na.rm=T), eHANPP=mean(eHANPP,na.rm=T), Ecological_Footprint=mean(Ecological_Footprint,na.rm=T), Material_Footprint=mean(Material_Footprint,na.rm=T)))

biophy_time_init<-biophy_time3
biophy_time_init[biophy_time_init==Inf]<-NA
biophy_time_init$Country<-as.character(biophy_time_init$Country) # rename countries for consistency
biophy_time_init$Country[which(biophy_time_init$Country=="Congo")]<-"Congo (Brazzaville)"
biophy_time_init$Country[which(biophy_time_init$Country=="DR Congo")]<-"Congo (Kinshasa)"
biophy_time_init$Country[which(biophy_time_init$Country=="Hong Kong")]<-"Hong Kong S.A.R. of China"
biophy_time_init$Country[which(biophy_time_init$Country=="Cote dIvoire")]<-"Ivory Coast"
biophy_time_init$Country[which(biophy_time_init$Country=="TFYR Macedonia")]<-"Macedonia"
biophy_time_init$Country[which(biophy_time_init$Country=="Taiwan")]<-"Taiwan Province of China"
biophy_time_init$Country[which(biophy_time_init$Country=="UAE")]<-"United Arab Emirates"
biophy_time_init$Country[which(biophy_time_init$Country=="UK")]<-"United Kingdom"
biophy_time_init$Country[which(biophy_time_init$Country=="USA")]<-"United States"
biophy_time_init$Country[which(biophy_time_init$Country=="Viet Nam")]<-"Vietnam"
biophy_time_init$Country[which(biophy_time_init$Country=="Antigua")]<-"Antigua and Barbuda"
biophy_time_init$Country[which(biophy_time_init$Country=="Brunei")]<-"Brunei Darussalam"
biophy_time_init$Country[which(biophy_time_init$Country=="Cape Verde")]<-"Cabo Verde"
biophy_time_init$Country[which(biophy_time_init$Country=="North Korea")]<-"Democratic People's Republic of Korea"
biophy_time_init<-droplevels(subset(biophy_time_init, Country!="Sudan"))
biophy_time_init$Country<-as.factor(biophy_time_init$Country)

Social and biophysical scores over time

social_s2$score<-apply(social_s2[,-c(1:2)], 1, function(x){
  if(length(which(!is.na(x)))>5){1-mean(x,na.rm=T)}else{NA}})

limits_s2<-biophy_time_init
limits_s2$score2<-apply(limits_s2[,-c(1:2)], 1, function(x){
  if(length(which(!is.na(x)))>3){length(which(!is.na(x) & x>=1))/length(which(!is.na(x)))}else{NA}})
limits_s2$score3<-apply(limits_s2[,-c(1,2,10)], 1, function(x){
  if(length(which(!is.na(x)))>3){if(length(which(x>=1))>0){mean(x[which(x>=1)], na.rm=T)-1}else{0}}else{NA}})

SEI over time

soclim_sa2<-merge(social_s2[,c("Country.name","group","score")],
                 limits_s2[,c("Country","group","score2","score3")], by.x=c("Country.name","group"),
                 by.y=c("Country","group"), all=T)
soclim_sa2<-na.omit(soclim_sa2)

soclim_sa2$score4<-rep(NA,nrow(soclim_sa2))
soclim_sa2$index_sa<-rep(NA,nrow(soclim_sa2))
for(i in 1:nrow(soclim_sa2)){
  y0=soclim_sa2$score[i] #qualitatif social
  x0=soclim_sa2$score2[i] #qualitatif limits
  k0=soclim_sa2$score3[i] #quantatitif limits
  k0_b=soclim_sa2$score3[i]/max(soclim_sa2$score3)
  
  x1=(x0+k0_b)/2
  k1=(k0_b+x0)/2
  l1=sqrt((x1-0)^2+(k1-0)^2)
  lt=sqrt(2)
  z0=l1/lt
  soclim_sa2$score4[i]<-z0
  
  l2=sqrt((y0-0)^2+(z0-0)^2) 
  soclim_sa2$index_sa[i]=1-l2/sqrt(2)
}

External socio-economic data

GDP

# from http://data.un.org/Data.aspx?d=WDI&f=Indicator_Code%3aNY.GDP.PCAP.PP.KD
gdp_time<-read.csv("UNdata_Export_20210427_172302587.csv") # GDP per capita, constant 2017 international $
gdp_time<-droplevels(subset(gdp_time, Year %in% c(2003:2017)))
gdp_time$group<-rep(NA, nrow(gdp_time))
gdp_time$group[which(gdp_time$Year %in% c(2003:2007))]<-2005
gdp_time$group[which(gdp_time$Year %in% c(2008:2012))]<-2010
gdp_time$group[which(gdp_time$Year %in% c(2013:2017))]<-2015
gdp_time<-as.data.frame(gdp_time %>% group_by(Country.or.Area,group) %>% summarize(gdp=mean(Value,na.rm=T)))

Population density

# population https://data.un.org/Data.aspx?q=population&d=PopDiv&f=variableID%3a12

pop_time<-read.csv("UNdata_Export_20200716_162002498.csv")

# country sizes https://worldpopulationreview.com/country-rankings/countries-by-density

country_size<-read.csv("country_size.csv")

pop_time<-droplevels(subset(pop_time, Variant=="Medium"))
names(pop_time)[2]<-"Year"
pop_time<-droplevels(subset(pop_time, Year %in% c(2003:2017)))
pop_time$group<-rep(NA, nrow(pop_time))
pop_time$group[which(pop_time$Year %in% c(2003:2007))]<-2005
pop_time$group[which(pop_time$Year %in% c(2008:2012))]<-2010
pop_time$group[which(pop_time$Year %in% c(2013:2017))]<-2015
pop_time<-as.data.frame(pop_time %>% group_by(Country.or.Area,group) %>% summarize(pop=mean(Value,na.rm=T)))
pop_time$pop<-1000*pop_time$pop

# rename countries for consistency 

country_size$country<-as.character(country_size$country)
country_size$country[which(country_size$country=="Bolivia")]<-"Bolivia (Plurinational State of)"
country_size$country[which(country_size$country=="Brunei")]<-"Brunei Darussalam"
country_size$country[which(country_size$country=="Cape Verde")]<-"Cabo Verde"
country_size$country[which(country_size$country=="Hong Kong")]<-"China, Hong Kong SAR"
country_size$country[which(country_size$country=="Macau")]<-"China, Macao SAR"
country_size$country[which(country_size$country=="Republic of the Congo")]<-"Congo"
country_size$country[which(country_size$country=="Ivory Coast")]<-"Côte d'Ivoire"
country_size$country[which(country_size$country=="Czech Republic")]<-"Czechia"
country_size$country[which(country_size$country=="DR Congo")]<-"Democratic Republic of the Congo"
country_size$country[which(country_size$country=="North Korea")]<-"Dem. People's Republic of Korea"
country_size$country[which(country_size$country=="Swaziland")]<-"Eswatini"
country_size$country[which(country_size$country=="Iran")]<-"Iran (Islamic Republic of)"
country_size$country[which(country_size$country=="Laos")]<-"Lao People's Democratic Republic"
country_size$country[which(country_size$country=="Micronesia")]<-"Micronesia (Fed. States of)"
country_size$country[which(country_size$country=="Macedonia")]<-"North Macedonia"
country_size$country[which(country_size$country=="South Korea")]<-"Republic of Korea"
country_size$country[which(country_size$country=="Moldova")]<-"Republic of Moldova"
country_size$country[which(country_size$country=="Russia")]<-"Russian Federation"
country_size$country[which(country_size$country=="Palestine")]<-"State of Palestine"
country_size$country[which(country_size$country=="Syria")]<-"Syrian Arab Republic"
country_size$country[which(country_size$country=="Tanzania")]<-"United Republic of Tanzania"
country_size$country[which(country_size$country=="United States")]<-"United States of America"
country_size$country[which(country_size$country=="Venezuela")]<-"Venezuela (Bolivarian Republic of)"
country_size$country[which(country_size$country=="Vietnam")]<-"Viet Nam"
country_size$country<-as.factor(country_size$country)


pop_time<-merge(pop_time, country_size[, c("country","area")],by.x="Country.or.Area", by.y="country", all.x=T)

# calculate population density

pop_time$density<-pop_time$pop/pop_time$area

Sustainable development index

# SDI from Hickel 2020 https://www.sciencedirect.com/science/article/abs/pii/S0921800919303386 available at https://www.sustainabledevelopmentindex.org/time-series

sdi<-read.csv("SDI.csv")

sdi_time<-melt(sdi[-1])
sdi_time$Year<-substr(sdi_time$variable,2,5)
sdi_time$variable<-NULL

sdi_time$Year<-as.numeric(sdi_time$Year)
sdi_time<-droplevels(subset(sdi_time, Year %in% c(2003:2017)))
sdi_time$group<-rep(NA, nrow(sdi_time))
sdi_time$group[which(sdi_time$Year %in% c(2003:2007))]<-2005
sdi_time$group[which(sdi_time$Year %in% c(2008:2012))]<-2010
sdi_time$group[which(sdi_time$Year %in% c(2013:2017))]<-2015
sdi_time<-as.data.frame(sdi_time %>% group_by(country,group) %>% summarize(sdi=mean(value,na.rm=T)))

Analyses in space and time

Space

Combine all data in space

# combined dataset
comb_soclim_2011b<-soclim_sa
comb_soclim_2011b$Country<-gsub("_"," ",comb_soclim_2011b$Country)
comb_soclim_2011b<-merge(comb_soclim_2011b,gdp_time[gdp_time$group=="2010",-2], by.x=c("Country"), by.y=c("Country.or.Area"),all.x=T)
comb_soclim_2011b<-merge(comb_soclim_2011b,pop_time[pop_time$group=="2010",-2], by.x=c("Country"), by.y=c("Country.or.Area"),all.x=T)
comb_soclim_2011b<-merge(comb_soclim_2011b,sdi_time[pop_time$group=="2010",-2], by.x=c("Country"), by.y=c("country"),all.x=T)
comb_soclim_2011b$gdp_l<-log(comb_soclim_2011b$gdp)
comb_soclim_2011b$density_l<-log(comb_soclim_2011b$density)

Relationship in space

# relationship with GDP

W_GDP_c<-14275.462 # in GDP per capita based on purchasing power parity (PPP) https://data.worldbank.org/indicator/NY.GDP.PCAP.PP.KD


gdp_lm_sa1<-lm(index_sa~gdp_l,comb_soclim_2011b)
gdp_sa1<-segmented(gdp_lm_sa1)
gdp_p_sa1<-summary(gdp_sa1) 

ggplot(comb_soclim_2011b, aes(x = gdp, y = index_sa))+
  geom_point() + theme_modern()+
  labs(x ="GDP (log scale)", y = "SEI")+geom_vline(xintercept=exp(gdp_sa1$psi[2]), linetype="dashed", size=1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]+1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]-1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=W_GDP_c, linetype="dashed", size=1,col="red") +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cs", k=4))+
  scale_x_continuous(trans="log10")

ggplot(comb_soclim_2011b, aes(x = gdp, y = 1-score))+
  geom_point() + theme_modern()+
  labs(x ="GDP (log scale)", y = "Social score")+geom_vline(xintercept=exp(gdp_sa1$psi[2]), linetype="dashed", size=1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]+1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]-1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=W_GDP_c, linetype="dashed", size=1,col="red") +
  scale_x_continuous(trans="log10")

ggplot(comb_soclim_2011b, aes(x = gdp, y = score4))+
  geom_point() + theme_modern()+
  labs(x ="GDP (log scale)", y = "Biophysical score")+geom_vline(xintercept=exp(gdp_sa1$psi[2]), linetype="dashed", size=1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]+1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=exp(gdp_sa1$psi[2]-1.96*gdp_sa1$psi[3]), linetype="dashed", size=.1) +
  geom_vline(xintercept=W_GDP_c, linetype="dashed", size=1,col="red") +
  scale_x_continuous(trans="log10")


# relationship with population density

W_POP_d<-58.423 # https://data.worldbank.org/indicator/EN.POP.DNST world population density in 2018 

density_lm_sa3<-lm(index_sa~density_l,comb_soclim_2011b[comb_soclim_2011b$density_l<7,])
density_sa3<-segmented(density_lm_sa3)
density_p_sa3<-summary(density_sa3)
ggplot(comb_soclim_2011b[comb_soclim_2011b$density_l<7,], aes(x = density, y = index_sa))+
  geom_point() + theme_modern()+
  labs(x ="Density (log scale)", y = "SEI")+
  scale_x_continuous(trans="log10")+
  geom_vline(xintercept=W_POP_d, linetype="dashed", size=1,col="red") +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cs", k=4))

ggplot(comb_soclim_2011b[comb_soclim_2011b$density_l<7,], aes(x = density, y = 1-score))+
  geom_point() + theme_modern()+
  labs(x ="Density (log scale)", y = "Social score")+
  scale_x_continuous(trans="log10")+
  geom_vline(xintercept=W_POP_d, linetype="dashed", size=1,col="red")

ggplot(comb_soclim_2011b[comb_soclim_2011b$density_l<7,], aes(x = density, y = score4))+
  geom_point() + theme_modern()+
  labs(x ="Density (log scale)", y = "Biophysical score")+
  scale_x_continuous(trans="log10")+
  geom_vline(xintercept=W_POP_d, linetype="dashed", size=1,col="red")


# relationship with SDI

W_SDI_c <- mean(sdi$X2011, na.rm=T)

sdi_lm_sa2<-lm(index_sa~sdi,comb_soclim_2011b)
sdi_sa2<-segmented(sdi_lm_sa2)
sdi_p_sa2<-summary(sdi_sa2)
ggplot(comb_soclim_2011b, aes(x = sdi, y = index_sa))+
  geom_point() + theme_modern()+
  labs(x ="SDI", y = "SEI")+
  geom_vline(xintercept=W_SDI_c, linetype="dashed", size=1,col="red") +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cs", k=4))

ggplot(comb_soclim_2011b, aes(x = sdi, y = 1-score))+
  geom_point() + theme_modern()+
  labs(x ="SDI", y = "Social score")+
  geom_vline(xintercept=W_SDI_c, linetype="dashed", size=1,col="red") 

ggplot(comb_soclim_2011b, aes(x = sdi, y = score4))+
  geom_point() + theme_modern()+
  labs(x ="SDI", y = "Biophysical score")+
  geom_vline(xintercept=W_SDI_c, linetype="dashed", size=1,col="red")

Time

Combine all data in time

comb_soclim<-soclim_sa2[,c("Country.name", "group","score","score4", "index_sa")]
comb_soclim<-merge(comb_soclim,gdp_time, by.x=c("Country.name","group"), by.y=c("Country.or.Area","group"),all.x=T)
comb_soclim<-merge(comb_soclim,pop_time, by.x=c("Country.name","group"), by.y=c("Country.or.Area","group"),all.x=T)
comb_soclim<-merge(comb_soclim,sdi_time, by.x=c("Country.name","group"), by.y=c("country","group"),all.x=T)
comb_soclim$gdp_l<-log(comb_soclim$gdp)
comb_soclim$density_l<-log(comb_soclim$density)

comb_soclim_delta<- droplevels(comb_soclim[comb_soclim$group %in% c(2005,2015),])
comb_soclim_delta<-as.data.frame(comb_soclim_delta[order(comb_soclim_delta$Country.name,comb_soclim_delta$group),]
                                 %>% group_by(Country.name)  %>% summarize(d_index_sa=index_sa[2]-index_sa[1],
                                                                          d_gdp=gdp[2]-gdp[1],
                                                                          d_density=density[2]-density[1],
                                                                          d_sdi=sdi[2]-sdi[1],
                                                                          d_score=score[2]-score[1],
                                                                          d_score4=score4[2]-score4[1],
                                                                          d_gdp_l=log(gdp[2])-log(gdp[1]),
                                                                          d_density_l=log(density[2])-log(density[1]),
                                                                          d_score_pc=d_score/score[1],
                                                                          d_score4_pc=d_score4/score4[1],
                                                                          d_gdp_pc=d_gdp/gdp[1],
                                                                          d_density_pc=d_density/density[1],
                                                                          d_gdp_pc_l=d_gdp_l/gdp_l[1],
                                                                          d_density_pc_l=d_density_l/density_l[1],
                                                                          d_sdi_pc=d_sdi/sdi[1],
                                                                          d_index_sa_pc=d_index_sa/index_sa[1],
                                                                          d_year=10))


comb_soclim_delta2<-merge(comb_soclim,comb_soclim_delta, by="Country.name")

Relationships in time

# relationship with GDP

d_gdp_sa_lm<-lm(d_index_sa_pc~gdp_l*d_gdp_pc_l,comb_soclim_delta2[comb_soclim_delta2$group==2015,])

# turning point in time
GDP_reverse_value <- summary(d_gdp_sa_lm)$coef[3,1]/(-summary(d_gdp_sa_lm)$coef[4,1])
# because dSEI = summary(d_gdp_sa_lm)$coef[1,1] + summary(d_gdp_sa_lm)$coef[3,1]+summary(d_gdp_sa_lm)$coef[4,1]*GDP)*dGDP
# then GDP = (-summary(d_gdp_sa_lm)$coef[1,1]/dGDP - summary(d_gdp_sa_lm)$coef[3,1])*1/summary(d_gdp_sa_lm)$coef[4,1]
# this becomes zero for all dGDP when GDP = -summary(d_gdp_sa_lm)$coef[3,1]/summary(d_gdp_sa_lm)$coef[4,1]

sub_data_gdp<-droplevels(comb_soclim_delta2[which(comb_soclim_delta2$group %in% c(2005,2015) &
                                                 !is.na(comb_soclim_delta2$d_index_sa) &
                                                 !is.na(comb_soclim_delta2$d_gdp)),])
ggplot(sub_data_gdp)+
  geom_point(aes(x=gdp,y=index_sa, col=as.factor(group)))+xlab("GDP (log scale)")+ylab("SEI")+
  geom_vline(xintercept=exp(GDP_reverse_value), linetype="dashed", size=1) +
  scale_colour_viridis(discrete = T, end=0.8)+theme_modern()+theme(legend.title = element_blank(), legend.position = c(0.8,0.9))+
  geom_segment(aes(y = index_sa, x = gdp, yend = index_sa+d_index_sa, xend = gdp+d_gdp),
               linetype=4,alpha=0.2,data=sub_data_gdp[which(sub_data_gdp$group==2005 & sub_data_gdp$d_index_sa<0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = gdp, yend = index_sa+d_index_sa, xend = gdp+d_gdp),
               linetype=1,alpha=0.2,data=sub_data_gdp[which(sub_data_gdp$group==2005 & sub_data_gdp$d_index_sa>0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = gdp, yend = index_sa+d_index_sa, xend = gdp+d_gdp),
               linetype=9,alpha=0.2,data=sub_data_gdp[which(sub_data_gdp$group==2005 & sub_data_gdp$d_index_sa==0),],arrow = arrow(length = unit(0.03, "npc")))+
  scale_x_continuous(trans="log10")

# relationship with population density

d_density_sa_lm<-lm(d_index_sa_pc~density_l*d_density_pc_l,comb_soclim_delta2[comb_soclim_delta2$group==2015 & comb_soclim_delta2$density_l<7,])

sub_data_pop<-droplevels(comb_soclim_delta2[which(comb_soclim_delta2$group %in% c(2005,2015) &
                                                 !is.na(comb_soclim_delta2$d_index_sa) &
                                                 !is.na(comb_soclim_delta2$d_density) &
                                                 comb_soclim_delta2$density<1100),])
ggplot(sub_data_pop)+
  geom_point(aes(x=density,y=index_sa, col=as.factor(group)))+xlab("Density (log scale)")+ylab("SEI")+
  scale_colour_viridis(discrete = T, end=0.8)+theme_modern()+theme(legend.title = element_blank(), legend.position = c(0.2,0.9))+
  geom_segment(aes(y = index_sa, x = density, yend = index_sa+d_index_sa, xend = density+d_density),
               linetype=4,alpha=0.2,data=sub_data_pop[which(sub_data_pop$group==2005 & sub_data_pop$d_index_sa<0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = density, yend = index_sa+d_index_sa, xend = density+d_density),
               linetype=1,alpha=0.2,data=sub_data_pop[which(sub_data_pop$group==2005 & sub_data_pop$d_index_sa>0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = density, yend = index_sa+d_index_sa, xend = density+d_density),
               linetype=9,alpha=0.2,data=sub_data_pop[which(sub_data_pop$group==2005 & sub_data_pop$d_index_sa==0),],arrow = arrow(length = unit(0.03, "npc")))+
  scale_x_continuous(trans="log10")


# relationship with SDI

d_sdi_sa_lm<-lm(d_index_sa_pc~sdi*d_sdi_pc,comb_soclim_delta2[comb_soclim_delta2$group==2015,])

  
sub_data_sdi<-droplevels(comb_soclim_delta2[which(comb_soclim_delta2$group %in% c(2005,2015) &
                                                  !is.na(comb_soclim_delta2$d_index_sa) &
                                                  !is.na(comb_soclim_delta2$d_sdi) &
                                                  comb_soclim_delta2$Country.name!="Sudan"),])
ggplot(sub_data_sdi)+
  geom_point(aes(x=sdi,y=index_sa, col=as.factor(group)))+xlab("SDI")+ylab("SEI")+
  scale_colour_viridis(discrete = T, end=0.8)+theme_modern()+theme(legend.title = element_blank(), legend.position = c(0.2,0.9))+
  geom_segment(aes(y = index_sa, x = sdi, yend = index_sa+d_index_sa, xend = sdi+d_sdi),
               linetype=4,alpha=0.2,data=sub_data_sdi[which(sub_data_sdi$group==2005 & sub_data_sdi$d_index_sa<0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = sdi, yend = index_sa+d_index_sa, xend = sdi+d_sdi),
               linetype=1,alpha=0.2,data=sub_data_sdi[which(sub_data_sdi$group==2005 & sub_data_sdi$d_index_sa>0),],arrow = arrow(length = unit(0.03, "npc")))+
  geom_segment(aes(y = index_sa, x = sdi, yend = index_sa+d_index_sa, xend = sdi+d_sdi),
               linetype=9,alpha=0.2,data=sub_data_sdi[which(sub_data_sdi$group==2005 & sub_data_sdi$d_index_sa==0),],arrow = arrow(length = unit(0.03, "npc")))

Save SEI in an extern file

to_save <- comb_soclim_delta2[,c("Country.name","group","score","score4","index_sa","d_index_sa","d_score","d_score4","d_score_pc","d_score4_pc","d_index_sa_pc","d_year")]

names(to_save)<-c("Country","Year","Social_score","Biophysical_score","SEI","delta_SEI","delta_Social_score","delta_Biophysical_score","delta_Social_score_percentage","delta_Biophysical_score_percentage","delta_SEI_percentage","d_year")

write.csv(x= to_save, file = "SEI.csv", row.names = F)

References

O’Neill, D. W., Fanning, A. L., Lamb, W. F., & Steinberger, J. K. (2018). A good life for all within planetary boundaries. Nature sustainability, 1(2), 88-95.