Showing posts with label data visualization. Show all posts
Showing posts with label data visualization. Show all posts

Monday, November 21, 2011

Revenue and Outlays 1980-1989

As the visualization indicates, revenues grew during the 1980's after cuts in marginal income tax rates, while spending also increased outpacing revenues and driving deficits.
The default charts were produced from default output (html and flash) generated by the R-code that follows.

I added one tweak, as described in the google visualization documentation here:
http://code.google.com/apis/chart/interactive/docs/gallery/motionchart.html which essentially saves the modifications to the settings in the chart. For more info regarding the data format, see my earlier post here.

Note- there are 3 visualization options - bubbles, dynamic bar chart, and a static line graph . These can be selected using the tabs at the top of the graph.

Created Using R- GoogleVis Package
Flash Enable Browser Required!



Data: budget • Chart ID: MotionChartID32de6867
R version 2.13.1 (2011-07-08) • googleVis-0.2.11
Google Terms of UseData Policy

We also saw patterns of increasing revenues during the Bush years following cuts in marginal tax rates, and a period from 2003-2007 where revenues were increasing and deficits were falling.  After the financial crisis, notice how revenues plunged while the deficit exploded.

Revenues and Outlays 2003-2009

(if preloaded labels and settings are not immediately applied, please use direct link for this post here.)


The data source is the CBO Budget/Historical Tables. I'd provide a link but it moves around constantly. A Google search should direct you to the most recent data.

R-Code:(for the first visualization)

#  ----------------------------------------------------------------------------------
# |PROGRAM NAME: budget_vis_R
# |DATE: 11/21/11
# |CREATED BY: MATT BOGARD 
# |PROJECT FILE:              
# |----------------------------------------------------------------------------------
# | PURPOSE: visualization of revenues and outlays in relation to 
# |          cuts in marginal tax rates 1980-89        
# | 
#  ---------------------------------------------------------------------------------
 
# see http://stackoverflow.com/questions/4646779/embedding-googlevis-charts-into-a-web-site/4649753#4649753
# for original R code reference
 
install.packages('googleVis') # install package if first time
 
library(googleVis) # load package
 
#  set R working directory- this is where your data file will go
#  with the script for creating the visualization
 
setwd("C:\\Users\\Documents\\Briefcase\\R Code and Data")
 
# read pre-formated  data
 
budget <- read.csv("budget.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(budget)
 
M <- gvisMotionChart(budget, "ITEM","YEAR")
 
#  look at data object- this includes the script that
#  will be used if you want to publish on your web page/blog
 
print(M) 
 
# open your browser and display the visualization
 
plot(M)
Created by Pretty R at inside-R.org

Tuesday, March 15, 2011

Applied Anaytics with R and Venn Diagrams


For a particular client I developed a predictive model that scored a set of patrons or donors at different points of time, providing the predicted probability that they would stop making contributions. At each point in time, they were more and more experienced with the service and more data about the patron was collected. As a result the model’s predictive accuracy improved with time. The client wanted to know, looking at the same cohort of customers over time, how often were the same customers predicted to stop making donations. In other words, at t=1, when the model is weakest, how many customers predicted to stop contributions were also on the ‘list’ at say t=3 when the model is much more accurate?

To do this I used the 'limma' package from the 'bioconductor' R mirror. (see reference below and R code that follows)

Before attempting to construct the Venn diagram, I had to take the scored donor data set and subset it based on all those patrons ever indicated to be 'high risk.' Then I created a data set with one row per patron and a binary indicator tracking their movement from 'novice' to 'intermediate' to 'experienced.' (t=1,2,3 respectively)



The format for the data set is similar to the layout below:

ID NOVICE INTERMEDIATE EXPERIENCED
1 1 1 1
2 1 1 0
3 1 0 0
.  .  .  .
etc.

The resulting Venn Diagram is below:



Note, there were 50 patrons in this example data set, and only 12 of those were predicted to be 'high risk' every time as they moved across each experience category or time period. 15 were high risk 'novice' patrons but never became part of the 'intermediate' or 'experienced' segments.

References:

How can I generate a Venn diagram in R?
http://www.ats.ucla.edu/stat/r/faq/venn.htm

R code:

# *------------------------------------------------------------------
# | PROGRAM NAME: R_Venn
# | DATE: 3/15/11
# | CREATED BY: MATT BOGARD  
# | PROJECT FILE: stats blog        
# *----------------------------------------------------------------
# | PURPOSE: CREATE VENN DIAGRAMS FOR MEMBERSHIP IN MULTIPLE GROUPS              
# |
# *------------------------------------------------------------------
# | COMMENTS:               
# |
# |  1: REFERENCES: How can I generate a Venn diagram in R? 
# |     http://www.ats.ucla.edu/stat/r/faq/venn.htm
# | 
# |  2: 
# |  3: 
# |*------------------------------------------------------------------
# | DATA USED: data scored by predictive model  
# |
# |*------------------------------------------------------------------
# | CONTENTS:               
# |
# |  PART 1: Run UCLA example code for practice 
# |  PART 2: My data
# |  PART 3: 
# *-----------------------------------------------------------------
# | UPDATES:               
# |
# |
# *------------------------------------------------------------------
 
 
 
 
 rm(list=ls()) # get rid of any existing data 
 ls() # view open data sets
 
 
 
# for 1st time use- get source code for bioconductor limma library
 
 
source("http://www.bioconductor.org/biocLite.R")
 
 
biocLite("limma")
 
ls() # see what data is there
 
library(limma) # load package
 
# *------------------------------------------------------------------
# | Part 1: Run UCLA example code            
# *-----------------------------------------------------------------
 
 
# read data
 
hsb2<-read.table("http://www.ats.ucla.edu/stat/R/notes/hsb2.csv", sep=',', header=T)
 
fix(hsb2) # view data set 
 
# create column vectors to represent the data sets
 
hw<-(hsb2$write>=60)
hm<-(hsb2$math >=60)
hr<-(hsb2$read >=60)
c3<-cbind(hw, hm, hr)
 
# create the matrix that will be used to plot the venn diagram
a <- vennCounts(c3)
a
 
vennDiagram(a) # plot venn diagram
 
# *------------------------------------------------------------------
# | Part 2: My Data           
# *-----------------------------------------------------------------
 
 
setwd('/Users/wkuuser/Desktop/R Data Sets') # set working directory
 
 
list<- read.csv("CUSTOMER_LOYALTY.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8") # read data
 
fix(list) # view data set
 
names(list) # get variable names (for cutting and pasting below)
 
# look at summary statistics for each data group
 
library(Hmisc) # for describe function
 
novice <-list[(list$NOVICE==1),] #subset novice segment
describe(novice) # n =37
 
intermediate <- list[(list$INTERMEDIATE==1),] # subset intermediate segment
describe(intermediate) # n=25
 
experienced <- list[(list$EXPERIENCED==1),] # subset experienced segment
describe(experienced)  # n= 25
 
# format data for use in venn diagram function below
 
l <- list[c("NOVICE","INTERMEDIATE","EXPERIENCED"  )] # keep only indicator variables
 
l3 <- as.matrix(l) # convert to a matrix # save as matrix
 
a <- vennCounts(l3) # create counts for venn digram
a
 
# plot venn digram 
vennDiagram(a, include = "both", names = c("Novice (n =37)", "Intermediate (n=25)", "Experienced (n=25)"), cex = 1, counts.col = "blue")
title("Donors Likely to Stop Contributions by Experience")
Created by Pretty R at inside-R.org

Friday, January 21, 2011

Flexibility of R Graphics

(note scroll all the way down to see 'old code' and 'new more flexible code'

Recall and older post that presented overlapping density plots using R (Visualizing Agricultural Subsidies by KY County) see image below.


The code I used to produce this plot makes use of the rbind and data.frame functions (see below)

library(colorspace) # package for rainbow_hcl function
 
 
ds <- rbind(data.frame(dat=KyCropsAndSubsidies[,][,"LogAcres"], grp="All"),
            data.frame(dat=KyCropsAndSubsidies[,][KyCropsAndSubsidies$subsidy_in_millions > 2.76,"LogAcres"], grp=">median"),
            data.frame(dat=KyCropsAndSubsidies[,][KyCropsAndSubsidies$subsidy_in_millions <= 2.76,"LogAcres"], grp="<=median"))
 
 
# histogram and density for all ears
hs <- hist(ds[ds$grp=="All",1], main="", xlab="LogAcres", col="grey90", ylim=c(0, 25), breaks="fd", border=TRUE)
 
dens <- density(ds[ds$grp=="All",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[1])
 
# density for above median subsidies
dens <- density(ds[ds$grp==">median",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[2])
 
# density for below median subsidies
dens <- density(ds[ds$grp=="<=median",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[3])
 
# Add a rug to illustrate density.
 
rug(ds[ds$grp==">median", 1], col=rainbow_hcl(3)[2])
rug(ds[ds$grp=="<=median", 1], col=rainbow_hcl(3)[3])
 
# Add a legend to the plot.
 
legend("topright", c("All", ">median", "<=media"), bty="n", fill=rainbow_hcl(3))
 
# Add a title to the plot. 
 
title(main="Distribution of Acres Planted by Subsidies Recieved Above or Below Median", sub=paste("Created Using R Statistical Package"))
Created by Pretty R at inside-R.org

I really don't understand the ins and outs of the rbind or data.frame functions, and in another project, when I tried to repeat a similar analysis, it wouldn't work. I could not figure out what my error was, but I new enough about R to create the plots with an alternative implementation. It is not as compact, but more general, and it worked. (see code below, although it references a new data set with new vars and produces 4 density curves vs. 3)

# histogram and density estimates for all data 
hs <- hist(trade_by_yr$logTrade,main="", xlab="trade", col="grey90", ylim=c(0, 95), breaks="fd", border=TRUE)  # histogram 
 
dens <- density(trade_by_yr$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density 
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[1]) # plot densiy 
 
# density estimates for year 2000 trade data
 
y2000 <- trade_by_yr[trade_by_yr$year==2000,] # subset data for year
dens <- density(y2000$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density  
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[2]) # plot densiy 
 
# density estimates for year 2004 trade data
 
y2004 <- trade_by_yr[trade_by_yr$year==2004,]  # subset data for year 
dens <- density(y2004$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density 
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[3]) # plot densiy 
 
# densty estimates for year 2008 trade data
 
y2008 <- trade_by_yr[trade_by_yr$year==2008,] # subset data for year
dens <- density(y2008$logTrade) # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density  
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[4]) # plot densiy 
 
# Add a legend to the plot.
 
legend("topright", c("All", "2000", "2004", "2008"), bty="n", fill=rainbow_hcl(4))
 
# Add a title to the plot. 
 
title(main="Distribution of Total World Trade Volume by Country by Year", sub=paste("Created Using R Statistical Package"))
Created by Pretty R at inside-R.org

See graph below:

Thursday, January 13, 2011

R GoogleVis Visualizing Taxes and Deficits

(flash enabled browser required)



For the simplest visualization, deselect 'trails' and select (checkbox) the variables DEFICIT, INCOME_TAX, TOTAL_REVENUE.

For the best visualization, deselect 'trails' under color select 'unique colors' for Size select 'IN_BILLIONS' and change the X-axis to PCT_GDP and the y-axis to IN_BILIONS. Select variables DEFICIT, INCOME_TAX, TOTAL_REVENUE

In any case, notice how early on, in the years following cuts in marginal income taxes, total revenues are increasing, revenues from income taxes are increasing, and the DEFICIT IS PLUNGING. All along spending is steadily incresing. Then about 2008 the deficit literally takes off running, both in billions of dollars and as a percentage of GDP as tax revenues start to plunge. Spending finds itself trapped as far as it can go in the northeast corner of the graph.

This is consistent with my previous graphical visualization here.

For more info on the code and references for this visualization see here.

The data source is the CBO Budget/Historical Tables. I'd provide a link but it moves around constantly. Just Google it and did for it and you can find the data. (or see below)

This is the format required for the R googleVis package. (I saved it as a csv file)

BUDGET_ITEM YEAR IN_BILLIONS PCT_GDP
CORP_TAX 2003 131.8 34.90604764
CORP_TAX 2004 189.4 45.88989817
CORP_TAX 2005 278.3 87.42060525
CORP_TAX 2006 353.9 142.5975397
CORP_TAX 2007 370.2 230.3657102
CORP_TAX 2008 304.3 667.9837559
CORP_TAX 2009 138.2 9.782782586
DEFICIT 2003 377.585 3.381560093
DEFICIT 2004 412.727 3.507495538
DEFICIT 2005 318.346 2.556790619
DEFICIT 2006 248.181 1.881156674
DEFICIT 2007 160.701 1.15429536
DEFICIT 2008 45.555 0.311849671
DEFICIT 2009 1412.686 9.313594409
INCOME_TAX 2003 793.7 7.108185563
INCOME_TAX 2004 809 6.875159344
INCOME_TAX 2005 927.2 7.446791422
INCOME_TAX 2006 1043.9 7.912529372
INCOME_TAX 2007 1163.5 8.357276253
INCOME_TAX 2008 1145.7 7.84296276
INCOME_TAX 2009 915.3 6.034414557
SPENDING 2003 2159.906 19.34359663
SPENDING 2004 2292.853 19.48545084
SPENDING 2005 2471.971 19.85359409
SPENDING 2006 2655.057 20.12474039
SPENDING 2007 2728.702 19.59992817
SPENDING 2008 2982.554 20.41726451
SPENDING 2009 3517.681 23.19146229
TOTAL_REVENUE 2003 131.8 1.180368977
TOTAL_REVENUE 2004 189.4 1.609586131
TOTAL_REVENUE 2005 278.3 2.235161834
TOTAL_REVENUE 2006 353.9 2.682483135
TOTAL_REVENUE 2007 370.2 2.659100704
TOTAL_REVENUE 2008 304.3 2.083105148
TOTAL_REVENUE 2009 138.2 0.911128692

Wednesday, January 12, 2011

Using R to Generate Motion Charts via Google's Visualization API

Just press the play button to get the full effect.(requires a flash enabled browser) Notice the ability to select variables and follow particular variables, the label options etc. are pretty robust. This utilizes Google's visualization tools but my local data.

The R code I used can be found here.






I figured out how to do this as a result of a recent the Revolution Analytics blog post demonstrating how to create motion charts with the new GoogleVis package. They link to an easy to follow tutorial on stack-overflow that walks you through the code.

Apparently the GoogleVis package generates script that accesses google’s Visualization API and allows you to use their data and visualizations. Or you can visualize your own data using their motion charts as I have done using the well known iris data set (adding a time variable as did the blogger kohske).

I found the Hans Rosling data visualization video inspiring and back in December I tried to create a similar visualization of agricultural subsidies, but it was a one time period static view. See below

or see here for the original post)

My next project will be to revisit the agricultural subsidy data, which is available across time, and convert my original bubble chart to a motion chart.

Sunday, December 26, 2010

Visualizing Taxes and Deficits

There has been a lot of debate about the impact of the early decade tax cuts on economic activity and deficits.

As the chart below depicts, from 2000-2009, we saw drastic increases in revenues (nearly 30% from 2001-2007) in the face of marginal tax cuts. Any deficit that resulted would have to be attributed to expenditures or outlays, and could not be attributed to cuts in marginal tax rates. As the graph shows, outlays also increased during this period, but even more drastically by 46%!


 As the next graphic shows, early on we saw a fairly rapid increase in the budget deficit from 2002-2003, a tapering off from 2003-2004 and  a rapid decline from 2004-2007, by as much as 61%! This is very impressive given the large amounts of spending increases depicted above. If it were not for the large influx of tax revenues during this period (in the face of marginal tax cuts) the deficit likely would have been on the increase vs. the precipitous fall depicted below.



However, on the heals of the financial crisis, going into 2008 & 2009, we start to see declining revenues, and unprecedented increases in spending and the deficit. From 2007 - 2009 we saw an increase in spending by about 28%, and an 88% increase over 2001 levels.  (indicated by the drastic upturn in outlays in the first graph)

But the impacts on the deficit were even more dramatic. From 2007-2008 we saw a 185% increase in the deficit, from 2008-2009 the deficit increased by 208%! Overall, compared to the 2002 levels that was an increase in the deficit of almost 800% over 7 years. If you compare to the 2007 low, considering the drastic reductions in the deficit after the tax cuts,  that is nearly an 800% increase in the deficit in just 3 years!

From 2004-2007 there was a rapid decline followed by a spike in the deficit during 2008 & 2009

Looking at the data, it appears that the reduction in marginal tax rates in the 2000's did not coincide with the rapid increase in the budget deficit that occurred at the end of the decade, but in fact were in step with the very rapid reduction in the budget deficit through 2007.

Most likely the deficit was the result from decreased revenues and increased expenditures associated with the financial crisis, not cuts in marginal tax rates. The real question becomes what was the cause of the financial crisis? There is no macroeconomic theory that I am aware of that links tax rates to business cycles, but many competing theories on business cycles as they relate to monetary policy or shocks to the production function. Pinning the blame for the deficit on one administration or the other during this transition period seems difficult, and a political subject that I am not interested in pursuing on this blog.

This adhoc analysis however does not prove that the effect of marginal tax cuts on the economy as a whole were positive or nrgative. Looking at one or two variables at a time (I haven't even included GDP or unemployment data) leaves one subject to mistakes. Only by building and testing models that specify multiple relationships among variables can you truly gauge the impact of the tax cuts on the deficit and economic output.  Lawrence Lindsey did this in 1987, looking specifically at revenue from income taxes paid by those earning over $200,000. Others have looked at the impact of tax cuts on economic activity, in terms of multipliers, and other research has been done relating taxes, spending, and unemployment (see references below). That is the proper context to view the impact of tax cuts, and as of yet, I am not aware of any empirical work that has been done to formally evaluate the true impact of the latest cuts in marginal income taxes.

References:

Lindsey, Lawrence B. 1987. “Individual Taxpayer Response to Taxcuts, 1982-1984.” J. of Public Economics 33 (July) 173-206

WHY DO EUROPEANS WORK (MUCH) LESS? IT IS TAXES AND GOVERNMENT SPENDING
Economic Inquiry, 2008, vol. 46, issue 2, pages 197-207

Economist Greg Mankiw gives a great review of the empirical work related to tax cuts and spending multipliers here on his blog:  http://gregmankiw.blogspot.com/2008/12/spending-and-tax-multipliers.html

Data Used: U.S. Budget Historical Tables http://www.whitehouse.gov/omb/budget/fy2009/hist.html (accessed Feb 2, 2009)


RECEIPTS  OUTLAYS      DEFICIT
2000 ............................................................................... 2,025,198  1,788,957 236,241
2001 ............................................................................... 1,991,142 1,862,906 128,236
2002 ............................................................................... 1,853,149 2,010,907 –157,758
2003 ............................................................................... 1,782,321 2,159,906 –377,585
2004 ............................................................................... 1,880,126 2,292,853 –412,727
2005 ............................................................................... 2,153,625 2,471,971 –318,346
2006 ............................................................................... 2,406,876 2,655,057 –248,181
2007 ............................................................................... 2,568,001 2,728,702 –160,701
2008 ............................................................................... 2,523,999 2,982,554 –458,555
2009 ............................................................................... 2,104,995 3,517,681 –1,412,686






































































Sunday, December 12, 2010

Visualizing Agricultural Subsidies by Kentucky County


In this post,  I provide results from my first full blown application of R to read, merge, clean, subset, manipulate, analyze and visualize data related to agricultural subsidies by Kentucky counties. This is very similar to the work I do on a daily basis, and was a great test of the capabilities of doing these tasks open source with R.  This is not breakthrough research, and doesn’t really even provide any new insights, but does demonstrate that farm subsidies accrue to those that grow the most crops.  It does demonstrate R’s usefulness in manipulating data sets and visualization tools.

Data

I extracted CSV files from the USDA Agricultural Statistics service related to the production (2008 data) of three major staple commodities grown in Kentucky and number of farm operations by county in Kentucky. I merged these data sets (via a series of left joins) with data from the Environmental Working Group related to subsidy accounts by Kentucky county.  No apparent way was available to extract this data to csv so it was copied, pasted, and cleaned up in excel to produce a csv file.  Additional transformations (conversion from character to numeric etc.) were conducted in R.

The data sets were combined through a series of left joins made possible from the ‘merge’ function in R, creating a data set called KyCropsAndSubsidies.

Overlapping Density Charts

The following is a histogram plotting the log acres planted by county with kernal density estimates for three populations: the entire KY county sample, those  counties receiving above median subsidies, and those receiving below median subsidies. It can be seen that the density curve for those receiving above median subsidies is to the right of the distributions of the population as a whole, and those receiving below median subsidy amounts (indicating more acres planted by this group of producers)

 

Bubble Charts

The following is essentially a scatter plot, plotting  acres planted in each county by the number of farms in each county. The size of each data point reflects the relative amount of subsidies received by all producers in each county. (note some obvious counties like Warren, Barren, Logan are missing due to missing data). It can be seen that the large circles (those counties receiving the most subsidies) float to the top indicating more acres planted. The number of farms per county seems to be unrelated to acres planed and subsidy amounts. 

 
Maps

The data set KyCropsAndSubsidies was merged with a spatial data frame (by county) which contains data related to county boundaries necessary for creating maps and plotting data by county.  The relationship between acres planted and subsidies by county can be seen below. Those counties planting the most acres tend to receive the most in subsidies. It is well known in Kentucky, that most of the row crops are grown in the western part of the state.



 

Code and References:

#   -------------------------------------------------------------
#  | PROGRAM NAME: Bubble_Crops
#  | DATE: 12-4-2010
#  | CREATED BY: Matt Bogard
#  | PROJECT FILE: /Users/user/Desktop/R Programs            
#  |-------------------------------------------------------------
#  | PURPOSE: Intitially to create bubble charts to demonstrate  
#  |  allocation             
#  | of farm subsides to producers by KY county, but expanded to 
#  | include
#  | numerous other visualization tools available using R
#  |
#  |-------------------------------------------------------------
#  | COMMENTS:              
#  |
#  |     1: Final data set for analysis is read in in the section
#  |     labeled 'basic statistical analysis'
#  | 2: In analysis section data for acres planted and subsidies
#  | are rescaled for
#  |     maps generated later in the program
#  |     3: Code for bubble charts adapted from:
#  | http://flowingdata.com/2010/11/23/how-to-make-bubble-charts/
#  |     4: Code for spatial analysis and mapping adapted from 
#  |  Harvard Applied Spatial Statistics workshop at:
#  |        
#  |http://www.people.fas.harvard.edu/~zhukov/spatial.html
#  |     5:
#  |
#  |
#  |-------------------------------------------------------------
#  | DATA USED:               
#  |     
#  |      1) various data sets downloaded from USDA NASS
#  |http://www.nass.usda.gov/#  |
#  |          2) data copied and pasted and formated in excel from EWG 
#  |subsidy data base: http://www.nass.usda.gov/
#  |      3) spatial data frame data used in the Harvard Applied 
#  |Spatial Statistics workshop at:
#  |       
#  | http://www.people.fas.harvard.edu/~zhukov/spatial.html
#  |      4) final data set located in project file:
#  |crops_and_subsidies.csv
#  |
#  |-------------------------------------------------------------
#  | CONTENTS:              
#  |
#  |          PART 1: get data    
#  |          PART 2: merge data sets
#  |          PART 3: recode and aggregate data sets
#  |          PART 4: basic statisticl analysis
#  |          PART 5: densisty plots
#  |          PART 6: bubble charts
#  |    PART 7: spatial analysis- Maps
#  |    
#  |
#  |-------------------------------------------------------------
#  | UPDATES:              
#  |
#  |
#  |
#   -------------------------------------------------------------


setwd('/Users/wkuuser/Desktop/R Data Sets')
rm(list=ls()) # get rid of any existing data
     ls() # view open data sets
    
#  *-------------------------------------------------------------
#  |
#  |
#  |
#  | get data
#  |
#  |
#  |
#  *-------------------------------------------------------------    
#  *-------------------------------------------------------------
#  | get farm # data
#  *-------------------------------------------------------------

farms <- read.csv("KyFarms.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(farms)

# get only the most recentl available year (2007)
farms07 <- farms[farms$Year =="2007",]
dim(farms07) # n=120
print(farms07)
names(farms07)

# keep only variables you want, convert 'Value' to numeric 'Operations' and rename county code for common merge key
farms07 <-farms07[c("Value","County.Code")]
farms07<-transform(farms07, Operations=as.numeric(paste(farms07$Value)))
farms07 <- rename(farms07, c(County.Code="CoFips"))
dim(farms07)
names(farms07)
farms07



#  *-------------------------------------------------------------
#  | get corn data
#  *-------------------------------------------------------------

corn <- read.csv("CR08.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(corn)

# subset- get only KY counties data

kycorn <- corn[corn$State =="Kentucky",]
dim(kycorn) # n=90
print(kycorn$County)
names(kycorn)

# keep and rename only relevant varaibles
# library(reshape)

kycorn <-kycorn[c("County","Yield","Harvested","CoFips","Planted.All.Purposes")]
names(kycorn)
kycorn <- rename(kycorn, c(Yield="cornYield",Harvested="cornHarvested",Planted.All.Purposes="cornPlanted"))
names(kycorn)
kycorn

#  *-------------------------------------------------------------
#  | get soybean data
#  *-------------------------------------------------------------

soybeans <- read.csv("SB08.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(soybeans)


# subset- get only KY counties data

kysoybeans <- soybeans[soybeans$State =="Kentucky",]
dim(kysoybeans) # n=79
print(kysoybeans$County)
names(kysoybeans)

# keep and rename only relevant varaibles

library(reshape)

kysoybeans <-kysoybeans[c("County","Yield","Harvested","Planted.All.Purposes")]
names(kysoybeans)
kysoybeans <- rename(kysoybeans, c(Yield="soybeanYield",Harvested="soybeanHarvested",Planted.All.Purposes="soybeanPlanted"))
names(kysoybeans)
kysoybeans

#  *-------------------------------------------------------------
#  | get wheat data
#  *-------------------------------------------------------------

wheat <- read.csv("AW08.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(wheat)

# subset- get only KY counties data

kywheat <- wheat[wheat$State =="Kentucky",]
dim(kywheat) # n= 70
print(kywheat$County)
names(kywheat)

# keep and rename only relevant varaibles
# library(reshape)

kywheat <-kywheat[c("County","Yield","Harvested","Planted.All.Purposes")]
names(kywheat)
kywheat <- rename(kywheat, c(Yield="wheatYield",Harvested="wheatHarvested",Planted.All.Purposes="wheatPlanted"))
names(kywheat)

#  *-------------------------------------------------------------
#  | get EWG subsidy data
#  *-------------------------------------------------------------

kysubsidies <- read.csv("EWGKYSubsidies.csv", na.strings=c(".", "NA", "", "?"), encoding="UTF-8")
names(kysubsidies)
kysubsidies


#  *-------------------------------------------------------------
#  |
#  |
#  |
#  | merge data sets
#  |
#  |
#  |
#  *-------------------------------------------------------------

# corn and soybeans left join
# note this will keep all corn growing counties (n=90) and add info about soybeans
# this will be the driver file for all subsequent merges, as a result, the staring base
# data is corn growing counties, i.e. the end data set will be based only on counties that
# grow corn in addition to other crops, but as a result counties that do not grow corn will be excluded
# from the analysis

corn_and_soybeans <- merge(kycorn,kysoybeans, by.kycorn=County,by.kysoybeans=County, all=FALSE, all.x=TRUE, all.y=FALSE)
dim(corn_and_soybeans)
names(corn_and_soybeans)
corn_and_soybeans

# all 3 crops- left join with wheat data
allKyCrops <- merge(corn_and_soybeans,kywheat, by.corn_and_soybeans=County,by.kywheat=County, all=FALSE, all.x=TRUE, all.y=FALSE)
dim(allKyCrops)
names(allKyCrops)

# left join with farm operations data on CoFips as key
# since County won't match due to case differences
farmAndCrops <- merge(allKyCrops,farms07, by.allKyCrops=CoFips,by.farms07=CoFips, all=FALSE, all.x=TRUE, all.y=FALSE)
dim(farmAndCrops)
names(farmAndCrops)
farmAndCrops

# crops and subsidies- left join with subsidy data
KyCropsAndSubsidies <- merge(farmAndCrops,kysubsidies, by.farmAndCrops=County,by.kysubsidies=County, all=FALSE, all.x=TRUE, all.y=FALSE)
dim(KyCropsAndSubsidies)
names(KyCropsAndSubsidies)

# quick report
KyCropsAndSubsidies[ c("County","cornPlanted", "soybeanPlanted","wheatPlanted","cornHarvested","soybeanHarvested","wheatHarvested","Operations","Subsidy")]



#  *-------------------------------------------------------------
#  |
#  |
#  |
#  | recode and aggregate variables
#  |
#  |
#  |
#  *-------------------------------------------------------------

# recode missing values-acres planted

KyCropsAndSubsidies$SoybeanAcresPlanted <- ifelse (is.na(KyCropsAndSubsidies$soybeanPlanted) == 'TRUE', (KyCropsAndSubsidies$SoybeanAcresPlanted <- 0),(KyCropsAndSubsidies$SoybeanAcresPlanted <- KyCropsAndSubsidies$soybeanPlanted))

KyCropsAndSubsidies[ c("County","SoybeanAcresPlanted")]

KyCropsAndSubsidies$WheatAcresPlanted <- ifelse (is.na(KyCropsAndSubsidies$wheatPlanted) == 'TRUE', (KyCropsAndSubsidies$WheatAcresPlanted <- 0),(KyCropsAndSubsidies$WheatAcresPlanted <- KyCropsAndSubsidies$wheatPlanted))

KyCropsAndSubsidies[ c("County","WheatAcresPlanted")]

# recode missing values - acres harvested

KyCropsAndSubsidies$SoybeanAcresHarvested <- ifelse (is.na(KyCropsAndSubsidies$soybeanHarvested) == 'TRUE', (KyCropsAndSubsidies$SoybeanAcresHarvested <- 0),(KyCropsAndSubsidies$SoybeanAcresHarvested <- KyCropsAndSubsidies$soybeanHarvested))

KyCropsAndSubsidies[ c("County","SoybeanAcresHarvested")]

KyCropsAndSubsidies$WheatAcresHarvested <- ifelse (is.na(KyCropsAndSubsidies$wheatHarvested) == 'TRUE', (KyCropsAndSubsidies$WheatAcresHarvested <- 0),(KyCropsAndSubsidies$WheatAcresHarvested <- KyCropsAndSubsidies$wheatHarvested))

KyCropsAndSubsidies[ c("County","WheatAcresHarvested")]

# aggregations

KyCropsAndSubsidies <- transform(KyCropsAndSubsidies, AcresPlanted = cornPlanted + SoybeanAcresPlanted + WheatAcresPlanted,
                                             AcresHarvested = cornHarvested + SoybeanAcresHarvested + WheatAcresHarvested)

KyCropsAndSubsidies <- transform(KyCropsAndSubsidies, LogAcres =log(AcresPlanted))