10 March 2015

Exploring the structure of national consumption

This entry is a direct continuation of my:

  1. entry
  2. entry
  3. entry
  4. entry
  5. entry

In the last entry I did a preliminary, but more formal statistical analysis of trends in the structure of national consumption.

In this entry I will explore the consistency of the trends across types of resources.

Also I will look to associated country area data to the dataset.

library(ggplot2);library(Hmisc)
# READING IN DATA

## SETTING DIRECTORY FOR EORA DATA ON LOCAL HARD DRIVE
wd<-"G:/Documents/PostDocKVA/Data/Eora" ### data directory
setwd(wd)

## READING IN DATA
### MATERIAL USE DATA - ENERGY DATASET
energy.df<-read.csv("TradeBalance_I-ENERGY.csv",header=TRUE)
### Reading in .csv file with annual gdp and population sizes
gdppop.df<-read.csv("gdppop.csv",header=TRUE,skip=1) #skipping the first line which includes a description of the file
### Reading in .csv file with various regional memberships
region.df<-read.csv("regionmembership.csv",header=T,skip=1) #skipping the first line which includes a description of the file


## REMOVING NEGATIVE AND ZERO CONSUMPTION ENTRIES
energy.df<-energy.df[which(energy.df[,"Consumption"]>0),]

## REMOVING NEGATIVE AND ZERO CONSUMPTION ENTRIES
#energy.df<-energy.df[which(as.character(energy.df$Country)!="Former USSR"),]


## merging the gdp and population size data onto the energy consumption data frame
energy.df<-merge(energy.df,gdppop.df,by=c("CountryA3","y","Country"),all.x=TRUE)


## To make consumption more comparable let's calculate per capita consumption by associating population data

### calculate per capita consumption and gdp consumption intensity by associating population data
energy.df[,"Consum.pop.int"]<-energy.df[,"Consumption"]/energy.df[,"val"]
energy.df[,"Consum.gdp.int"]<-energy.df[,"Consumption"]/energy.df[,"GDP"]

### calculating scaled consumption intensitities with respect to population size and gdp
energy.df<-energy.df[order(energy.df[,"Country"],energy.df[,"y"]),]
energy.df[,"Consum.pop.int.scale"]<-unlist(by(energy.df,energy.df[,"Country"], function(x) scale(x[,"Consum.pop.int"],center=TRUE,scale=TRUE)))
energy.df[,"Consum.gdp.int.scale"]<-unlist(by(energy.df,energy.df[,"Country"], function(x) scale(x[,"Consum.gdp.int"],center=TRUE,scale=TRUE)))

energy.df<-energy.df[-which( energy.df[,"y"] %in% c(1991,2000,2011)),]


## Calculating CONSUMPTION BALANCE with respect to EXTRACTION and IMPORTS

### Subtracting exports from imports and domestic extraction in proportion to their size.

### Conceptual code
#### consum.extract <- extraction - (export * extraction/(import+extraction))
#### consum.import <- import - (export * import/(import+extraction))

### two ways of estimating domestic extraction

#energy.df[,"Extraction"]<- (energy.df[,"TerritorialEmissions"] + energy.df[,"DirectEmissions"]) ## These two categories both refer to domestically extracted resources, i thik. Check up with Wiedmann dataset.

energy.df[,"Extraction"]<- energy.df[,"TerritorialEmissions"]

### consumed resources that come from domestic extraction

energy.df[,"Consum.extract"]<- (energy.df[,"Extraction"]) - 
    (energy.df[,"Exports"] * ((energy.df[,"Extraction"]) /
       (energy.df[,"Extraction"] + energy.df[,"Imports"])))


### consumed resources that come from imports

energy.df[,"Consum.import"]<- (energy.df[,"Imports"]) - 
    (energy.df[,"Exports"] * ((energy.df[,"Imports"]) /
       (energy.df[,"Extraction"] + energy.df[,"Imports"])))


### consumption balance index

energy.df[,"Consum.balance"]<-(energy.df[,"Consum.extract"]-energy.df[,"Consum.import"])/energy.df[,"Consumption"]
energy.df<-energy.df[-which(energy.df[,"Consum.balance"]>1),] ## removing entries that are above 1 (i.e. all resources being domestically extracted)



### calculating mean gdp and population size

energy.df<-energy.df[order(energy.df[,"Country"],energy.df[,"y"]),]
energy.df[,"mean.gdp"]<-unlist(tapply(energy.df[,"GDP"],energy.df[,"Country"],function(x) rep(mean(x,na.rm=TRUE),length(x))))
energy.df[,"mean.val"]<-unlist(tapply(energy.df[,"val"],energy.df[,"Country"],function(x) rep(mean(x,na.rm=TRUE),length(x))))
energy.df[,"mean.gdp.per.cap"]<-energy.df[,"mean.gdp"]/energy.df[,"mean.val"]
energy.df[,"log10.mean.gdp.per.cap"]<-log10(energy.df[,"mean.gdp.per.cap"])

### country membership of grouped population size

val.membership.df<-data.frame("mean.val"=unique(energy.df[,"mean.val"]),"val.membership"=NA)
val.membership.df[,"val.membership"]<-as.numeric(cut2(val.membership.df$mean.val,m=20))

energy.df<-merge(energy.df,val.membership.df,by="mean.val",all.x=TRUE)


#####

## CLEANING OF DATASET

### exclude the two sudan's - Sudan, South Sudan
# exclude Montenegro, Gaza Strip

### include Andorra, Monaco, San Marino, Liechtenstein despite constant -1 index

## Include Czech Republic 1993 - 2010
## Eritrea 1991-2010, Greenland 1991-2007

# Exclude 2008-2010 for countries where it suddently drops to -1


Consum.balance.BA2006.df<-do.call(rbind,by(energy.df,energy.df[,"Country"],function(x) c(mean(x[which(x[,"y"]%in%c(1970:2007)),"Consum.balance"]),
                                                mean(x[which(x[,"y"]%in%c(2008:2010)),"Consum.balance"])))
        )
Remove.0810.vec<-rownames(Consum.balance.BA2006.df[which(Consum.balance.BA2006.df[,1]>-1 & Consum.balance.BA2006.df[,2]==-1),])
Remove.country.vec<-c("Sudan","South Sudan","Montenegro","Gaza Strip","Former USSR","Statistical Discrepancies")
From.1991.vec<-c("Eritrea","Greenland")
From.1993.vec<-c("Czech Republic")

energy.df<-energy.df[-which( energy.df[,"Country"] %in% Remove.0810.vec & energy.df[,"y"] > 2007 ), ]
energy.df<-energy.df[ which( energy.df[,"Country"] %in% setdiff( energy.df[,"Country"], Remove.country.vec) ), ]
energy.df<-energy.df[ -which( energy.df[,"Country"] %in% From.1991.vec &  energy.df[,"y"] < 1991), ]
energy.df<-energy.df[ -which( energy.df[,"Country"] == From.1993.vec &  energy.df[,"y"] < 1993), ]

### We still see a lot of strange dynamics in 2008 possibly related to the financial crisis - so for this first analysis let's onlylook at years up to that time.
energy.df <- energy.df[which(energy.df[,"y"] < 2008), ]

### mean consumption.balance index
energy.df <- energy.df[order(energy.df[,"Country"],energy.df[,"y"]),]
energy.df[,"Consum.balance.avg"]<-unlist(by(energy.df, energy.df[,"Country"], function(x) rep(mean(x[,"Consum.balance"]),nrow(x))))
#### consumption avg - data frame
energy.ConsumAvg.df<-unique(energy.df[,c("Country","Consum.balance.avg","mean.val","mean.gdp")])

### trend in energy consumption balance
energy.df <- energy.df[order(energy.df[,"Country"],energy.df[,"y"]),]
energy.df[,"ConsumBalance.y.trend"]<-unlist(by(energy.df,energy.df[,"Country"],function(x) rep(coef(lm(x[,"Consum.balance"]~x[,"y"]))[2],nrow(x))))
#### consumption balance trend - data frame
energy.consumtrend.df<-unique(energy.df[,c("Consum.balance.avg","mean.gdp.per.cap","ConsumBalance.y.trend","log10.mean.gdp.per.cap"),])

Picking up where we left off.

### 3d-plane and scatterplot of trend in consumption balance as a function of mean consumption balance and log gdp per capita

attach(energy.consumtrend.df)
## The following objects are masked from i:
## 
##     Consum.balance.avg, ConsumBalance.y.trend,
##     log10.mean.gdp.per.cap, mean.gdp.per.cap
model <- lm(ConsumBalance.y.trend ~ log10.mean.gdp.per.cap + Consum.balance.avg)
x <-range(log10.mean.gdp.per.cap)
x <- seq(x[1], x[2], length.out=50)    
y <- range(Consum.balance.avg)
y <- seq(y[1], y[2], length.out=50)


z <- outer(x,y, 
           function(log10.mean.gdp.per.cap,Consum.balance.avg)
                     predict(model, data.frame(log10.mean.gdp.per.cap,Consum.balance.avg)))

p <- persp(x,y,z, theta=120, phi=10, 
           col="lightblue",expand = 0.5,shade = 0.2,
           xlab="log10 mean gdp per capita", ylab="Mean consumption balance", zlab="Consumption balance trend")
#Project 3d points to 2d, so you can use points and segment:

obs <- trans3d(log10.mean.gdp.per.cap, Consum.balance.avg, ConsumBalance.y.trend, p)
pred <- trans3d(log10.mean.gdp.per.cap, Consum.balance.avg,fitted(model),p)
points(obs, col="green",pch=16,cex=0.5)
segments(obs$x, obs$y, pred$x, pred$y,col="white",lwd=0.25)

detach(energy.consumtrend.df)
rm("obs","pred","model","x","y","z","p")

Now let’s do the same analysis for five related datasets on material footprint covering biomass, fossil fuels, ores, construction materials and their summarized material footprint

## SETTING DIRECTORY FOR Wiedmann - material footprint - EORA DATA ON LOCAL HARD DRIVE
setwd(paste(wd,"/Wiedmann",sep=""))
dir()
##  [1] "TradeBalance_BauxiteaEtc-grossore.csv"                       
##  [2] "TradeBalance_BauxiteaEtc-grossore.xlsx"                      
##  [3] "TradeBalance_BIOMASS.csv"                                    
##  [4] "TradeBalance_BIOMASS.xlsx"                                   
##  [5] "TradeBalance_Browncoal.csv"                                  
##  [6] "TradeBalance_Browncoal.xlsx"                                 
##  [7] "TradeBalance_Cereals.csv"                                    
##  [8] "TradeBalance_Cereals.xlsx"                                   
##  [9] "TradeBalance_Chalk.csv"                                      
## [10] "TradeBalance_Chalk.xlsx"                                     
## [11] "TradeBalance_Chemicalandfertilizerminerals (1).csv"          
## [12] "TradeBalance_Chemicalandfertilizerminerals.csv"              
## [13] "TradeBalance_Chemicalandfertilizerminerals.xlsx"             
## [14] "TradeBalance_CONSTMA.csv"                                    
## [15] "TradeBalance_CONSTMA.xlsx"                                   
## [16] "TradeBalance_CopperOres-grossore.csv"                        
## [17] "TradeBalance_CopperOres-grossore.xlsx"                       
## [18] "TradeBalance_Cropresidues_used.csv"                          
## [19] "TradeBalance_Cropresidues_used.xlsx"                         
## [20] "TradeBalance_Crudeoilandnaturalgasliquids (1).csv"           
## [21] "TradeBalance_Crudeoilandnaturalgasliquids.csv"               
## [22] "TradeBalance_Crudeoilandnaturalgasliquids.xlsx"              
## [23] "TradeBalance_FFUEL.csv"                                      
## [24] "TradeBalance_FFUEL.xlsx"                                     
## [25] "TradeBalance_Fibres.csv"                                     
## [26] "TradeBalance_Fibres.xlsx"                                    
## [27] "TradeBalance_Fruits.csv"                                     
## [28] "TradeBalance_Fruits.xlsx"                                    
## [29] "TradeBalance_GoldSilverEtc-grossore.csv"                     
## [30] "TradeBalance_GoldSilverEtc-grossore.xlsx"                    
## [31] "TradeBalance_Grazedbiomass.csv"                              
## [32] "TradeBalance_Grazedbiomass.xlsx"                             
## [33] "TradeBalance_Hardcoal.csv"                                   
## [34] "TradeBalance_Hardcoal.xlsx"                                  
## [35] "TradeBalance_I-CMFA-TOT.csv"                                 
## [36] "TradeBalance_I-CMFA-TOT.xlsx"                                
## [37] "TradeBalance_IronOres.csv"                                   
## [38] "TradeBalance_IronOres.xlsx"                                  
## [39] "TradeBalance_LeadOres-grossore.csv"                          
## [40] "TradeBalance_LeadOres-grossore.xlsx"                         
## [41] "TradeBalance_Naturalgas.csv"                                 
## [42] "TradeBalance_Naturalgas.xlsx"                                
## [43] "TradeBalance_NickelOres-grossore.csv"                        
## [44] "TradeBalance_NickelOres-grossore.xlsx"                       
## [45] "TradeBalance_Non-Metallicminerals-primarilyconstruction.csv" 
## [46] "TradeBalance_Non-Metallicminerals-primarilyconstruction.xlsx"
## [47] "TradeBalance_Nuts.csv"                                       
## [48] "TradeBalance_Nuts.xlsx"                                      
## [49] "TradeBalance_Oilbearingcrops.csv"                            
## [50] "TradeBalance_Oilbearingcrops.xlsx"                           
## [51] "TradeBalance_ORES.csv"                                       
## [52] "TradeBalance_ORES.xlsx"                                      
## [53] "TradeBalance_OrnamentalStone.csv"                            
## [54] "TradeBalance_OrnamentalStone.xlsx"                           
## [55] "TradeBalance_Othercrops.csv"                                 
## [56] "TradeBalance_Othercrops.xlsx"                                
## [57] "TradeBalance_Otherminingandquarryingproductsnec.csv"         
## [58] "TradeBalance_Otherminingandquarryingproductsnec.xlsx"        
## [59] "TradeBalance_OtherOres-grossore.csv"                         
## [60] "TradeBalance_OtherOres-grossore.xlsx"                        
## [61] "TradeBalance_Peat.csv"                                       
## [62] "TradeBalance_Peat.xlsx"                                      
## [63] "TradeBalance_Pulses.csv"                                     
## [64] "TradeBalance_Pulses.xlsx"                                    
## [65] "TradeBalance_Rootsandtubers.csv"                             
## [66] "TradeBalance_Rootsandtubers.xlsx"                            
## [67] "TradeBalance_Salt.csv"                                       
## [68] "TradeBalance_Salt.xlsx"                                      
## [69] "TradeBalance_Sugarcrops.csv"                                 
## [70] "TradeBalance_Sugarcrops.xlsx"                                
## [71] "TradeBalance_Timber_Industrialroundwood.csv"                 
## [72] "TradeBalance_Timber_Industrialroundwood.xlsx"                
## [73] "TradeBalance_TinOres-grossore.csv"                           
## [74] "TradeBalance_TinOres-grossore.xlsx"                          
## [75] "TradeBalance_UraniumEtc-grossore.csv"                        
## [76] "TradeBalance_UraniumEtc-grossore.xlsx"                       
## [77] "TradeBalance_Vegetables.csv"                                 
## [78] "TradeBalance_Vegetables.xlsx"                                
## [79] "TradeBalance_Woodfuelandotherextraction.csv"                 
## [80] "TradeBalance_Woodfuelandotherextraction.xlsx"                
## [81] "TradeBalance_ZinOrres-grossore.csv"                          
## [82] "TradeBalance_ZinOrres-grossore.xlsx"
#### We are interested in the following files

#"TradeBalance_CONSTMA.csv","TradeBalance_ORES.csv","TradeBalance_FFUEL.csv","TradeBalance_BIOMASS.csv","TradeBalance_I-CMFA-TOT.csv"

### We will read the files we are interested in into a list
dframes.list<-list()

## FOSSIL FUELS
dframes.list[["FFUEL"]]<-read.csv("TradeBalance_FFUEL.csv",header=T)
dframes.list[["FFUEL"]][,"NAME"]<-"FFUEL"
## BIOMASS
dframes.list[["BIOMASS"]]<-read.csv("TradeBalance_BIOMASS.csv",header=T)
dframes.list[["BIOMASS"]][,"NAME"]<-"BIOMASS"

## CONSTRUNCTION MATERIALS
dframes.list[["CONSTMA"]]<-read.csv("TradeBalance_CONSTMA.csv",header=T)
dframes.list[["CONSTMA"]][,"NAME"]<-"CONSTMA"

## ORES
dframes.list[["ORES"]]<-read.csv("TradeBalance_ORES.csv" ,header=T)
dframes.list[["ORES"]][,"NAME"]<-"ORES"

## ALL MATERIALS - FOSSIL FUELS, BIOMASS, CONSTRUCTION MATERIALS, ORES
dframes.list[["TOTAL"]]<-read.csv("TradeBalance_I-CMFA-TOT.csv",header=T)
dframes.list[["TOTAL"]][,"NAME"]<-"TOTAL"

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

## I will attempt a cleaning function to prepare each data frame for analysi
prepare.dframe.fct<-function(i){
### cleaning data frame
i<-i[which(as.character(i$Country)!="Former USSR"),]#!(df1$id %in% idNums1)
i<-i[-which( i[,"y"] %in% c(1991,2000,2011)),]
i <- i[which(i[,"y"] < 2008), ]
i<-i[ which( i[,"Country"] %in% setdiff( i[,"Country"], Remove.country.vec) ), ]
i<-i[ -which( i[,"Country"] %in% From.1991.vec &  i[,"y"] < 1991), ]
i<-i[ -which( i[,"Country"] == From.1993.vec &  i[,"y"] < 1993), ]

### calculating structure of consumption
#### Territorial emissions as extraction
i[,"Extraction"]<-i[,"TerritorialEmissions"]
#### consumed resources that come from domestic extraction
i[,"Consum.extract"]<- (i[,"Extraction"]) - 
    (i[,"Exports"] * ((i[,"Extraction"]) /
       (i[,"Extraction"] + i[,"Imports"])))
#### consumed resources that come from imports
i[,"Consum.import"]<- (i[,"Imports"]) - 
    (i[,"Exports"] * ((i[,"Imports"]) /
       (i[,"Extraction"] + i[,"Imports"])))
#### consumption balance index
i[,"Consum.balance"]<-(i[,"Consum.extract"]-i[,"Consum.import"])/i[,"Consumption"]
i<-i[which(abs(i[,"Consum.balance"])<1),] ## removing entries that are above 1 (i.e. all resources being domestically extracted or imported)


### associating population size and gdp data with data frame - then calculating basic metrics for later analysis
#### merging the gdp and population size data onto the energy consumption data frame
i<-merge(i,gdppop.df,by=c("CountryA3","y","Country"),all.x=TRUE)

i <- i[order(i[,"Country"],i[,"y"]),]

i[,"mean.gdp"]<-unlist(tapply(i[,"GDP"],i[,"Country"],function(x) rep(mean(x,na.rm=TRUE),length(x))))
i[,"mean.val"]<-unlist(tapply(i[,"val"],i[,"Country"],function(x) rep(mean(x,na.rm=TRUE),length(x))))
i[,"mean.gdp.per.cap"]<-i[,"mean.gdp"]/i[,"mean.val"]
i[,"log10.mean.gdp.per.cap"]<-log10(i[,"mean.gdp.per.cap"])

val.membership.df<-data.frame("mean.val"=unique(i[,"mean.val"]),"val.membership"=NA)
val.membership.df[,"val.membership"]<-as.numeric(cut2(val.membership.df$mean.val,m=20))

i<-merge(i,val.membership.df,by="mean.val",all.x=TRUE)


### Calculating mean consumption balance
### mean consumption.balance index
i <- i[order(i[,"Country"],i[,"y"]),]
i[,"Consum.balance.avg"]<-unlist(by(i, i[,"Country"], function(x) rep(mean(x[,"Consum.balance"]),nrow(x))))

### Calculating trend in consumption balance
i <- i[order(i[,"Country"],i[,"y"]),]
i[,"ConsumBalance.y.trend"]<-unlist(by(i,i[,"Country"],function(x) rep(coef(lm(x[,"Consum.balance"]~x[,"y"]))[2],nrow(x))))

return(i)
}


## apply the data frame preparation function
dframes.list<-lapply(dframes.list,prepare.dframe.fct)

dframes.list[["ENERGY"]]<-energy.df
dframes.list[["ENERGY"]][,"NAME"]<-"ENERGY"


## Get the calculated averages and trends for each data frame
ConsumTrend.df.list<-lapply(dframes.list,function(x) unique(x[,c("Country","Consum.balance.avg","ConsumBalance.y.trend","mean.val","mean.gdp","mean.gdp.per.cap","log10.mean.gdp.per.cap","NAME")]))

                            
## Do the persp plot for all data frames

PerspPlot.fct<-function(i){
  
attach(i)
model <- lm(ConsumBalance.y.trend ~ log10.mean.gdp.per.cap + Consum.balance.avg)
x <-range(log10.mean.gdp.per.cap)
x <- seq(x[1], x[2], length.out=50)    
y <- range(Consum.balance.avg)
y <- seq(y[1], y[2], length.out=50)


z <- outer(x,y, 
           function(log10.mean.gdp.per.cap,Consum.balance.avg)
                     predict(model, data.frame(log10.mean.gdp.per.cap,Consum.balance.avg)))

p <- persp(x,y,z, theta=120, phi=10, 
           col="lightblue",expand = 0.5,shade = 0.2,
           xlab="log gdp/cap", ylab="cons avg", zlab="cons trend",main=NAME[1])
#Project 3d points to 2d, so you can use points and segment:

obs <- trans3d(log10.mean.gdp.per.cap, Consum.balance.avg, ConsumBalance.y.trend, p)
pred <- trans3d(log10.mean.gdp.per.cap, Consum.balance.avg,fitted(model),p)
points(obs, col="green",pch=16,cex=0.5)
#segments(obs$x, obs$y, pred$x, pred$y,col="white",lwd=0.25)


detach(i)

rm("obs","pred","model","x","y","z","p")

}
par(mfrow=c(3,2))
plots<-lapply(ConsumTrend.df.list,PerspPlot.fct)

It looks like the negative relationship with log10 gdp per capita gdp is shared among the various resource types. However, that cannot be said for the influence of average consumption, which varies in sign and slope between the various resource types. In particular, ore resources look to have a different dynamic than other types of resources. This is something that requires further exploration and I will return to later.

Let’s just look at log gdp per cap to visualize the shared realtionship in a more graphically simple way.

ConsumTrend.df<-do.call(rbind,ConsumTrend.df.list)


ggplot(ConsumTrend.df,aes(y=ConsumBalance.y.trend,x=log10.mean.gdp.per.cap)) + geom_point(color="grey50",size=1) + ylab("Consumption balance trend") + xlab("log10 gdp/cap") + geom_smooth() + facet_wrap(~NAME, scales="free_y") + xlim(c(-1,2.1))# + ylim(c(-1,1)) 

Ores is still the resource (along with ffuels) that diverges mostly from the hypothesized shared relationship.

How about the average consumption balance?

ConsumTrend.df<-do.call(rbind,ConsumTrend.df.list)


ggplot(ConsumTrend.df,aes(y=Consum.balance.avg,x=log10(mean.val))) + geom_point(color="grey50",size=1) + ylab("Consumption balance") + xlab("log10 Population size") + geom_smooth() + facet_wrap(~NAME) + ylim(c(-1,1))

### Fossil fuels and ores don't follow the hyptohesis

## How about the correlation with gdp per cap?

ggplot(ConsumTrend.df,aes(y=Consum.balance.avg,x=log10.mean.gdp.per.cap)) + geom_point(color="grey50",size=1) + ylab("Consumption balance") + xlab("log10 gdp/cap") + geom_smooth() + facet_wrap(~NAME) + ylim(c(-1,1)) + xlim(c(-1,2.1))

### Only the balance of biomass and total material footprint show negative correlation with log gdp per cap.

Consistency of relationships among continents

Now let’s examine whether the preliminary relationships found are consistent as relationship within each continent separately.

I will read in continent information about countries. I have added this line higher up in this entry. This is what it looks like

setwd(wd)
region.df<-read.csv("regionmembership.csv",header=T,skip=1) #skipping the first line which includes a description of the file
head(region.df)
##   country_id region_id region_type region_type_id region_name  a3
## 1        533      1001 GDP Regions              1 High Income ABW
## 2        660      1001 GDP Regions              1 High Income AIA
## 3        248      1001 GDP Regions              1 High Income ALA
## 4         20      1001 GDP Regions              1 High Income AND
## 5        528      1001 GDP Regions              1 High Income ANT
## 6        804      1001 GDP Regions              1 High Income ARE
##                   name
## 1                Aruba
## 2             Anguilla
## 3        Aland Islands
## 4              Andorra
## 5 Netherlands Antilles
## 6                  UAE
### Get the continent info
continent.df<-region.df[which(region.df[,"region_type"]=="Continents"),c("region_name","name")]

ConsumTrend.df<-merge(ConsumTrend.df,continent.df,by.x="Country",by.y="name",all.z=TRUE)

First, is the trend in resource subsidization with gdp per capita consistent within continents?

ggplot(ConsumTrend.df[which(ConsumTrend.df[,"region_name"] != "Australia"),],aes(y=ConsumBalance.y.trend,x=log10.mean.gdp.per.cap)) + geom_point(color="grey50",size=1) + ylab("Consumption balance trend") + xlab("log10 gdp/cap") + geom_smooth(method="lm") + facet_grid(NAME ~ region_name, scales="free") + xlim(c(-1,2.1)) + ylim(c(-0.02,0.02))