How do you draw a Pirate Hat?: rnorm()

Since today, 19 Sept is talk like a pirate day, I would like to relate a story about stats and pirates. While driving, my wife was in the back seat with our young child asked me how to draw a pirate hat. Since I was driving I could not show her how to do it, so I had to figure out a way to describe it that both of us could understand. My reply, "Draw a normal distribution curve."

She looked at me kind of weird, then drew it, half laughing to her self, and it worked. A simple pirate hat kind of looks like a normal distribution. My child was happy, my wife and I were kind of giggling our selves, because some how I was able to bring stats into a child's drawing. We both were amused. Below is the R code and very crude picture.

x<-rnorm(n=10000, mean=1, sd=.25)
plot(density(x), type='l', main='', ylab='', xlab='')

Created by Pretty R at inside-R.org

Look Familiar? Mapping in R

For those who have been following the R-Bloggers this picture should be yesterdays news, but there was an article on the BBC on it, full story HERE. I find it interesting how it was on Drudge Report (www.drudgereport.com) which gets over 1 million hits a day. The article gives the perspective of 5 different points of view, the 'data visualization' interpretation is pretty interesting, as are the rest. For those who want to do more than just look at a pretty picture here are some links to help recreate a similar, if not better map.

Links:

Mapping the world Airlines

Worlds Largest Airlines

Data for Airlines

Tutorial on how to create the great circles

Down and Dirty Forecasting: Part 1

I wanted to see what I could do in a hurry using the commands found at Forecasting: Principles and Practice . I chose a simple enough data set of Wisconsin Unemployment from 1976 to the present (April 2013). I kept the last 12 months worth of data to test the accuracy of the models. The next blog post will include a multiple regression analysis. The analysis is lacking many important steps, particularly the ARIMA, but this is a down and dirty exercise.

library(forecast)
library(lmtest)
library(caret)
#State Unemployment seasonally adjusted
#http://www.quandl.com/FRED-Federal-Reserve-Economic-Data/WIUR-Unemployment-Rate-in-Wisconsin
 
#Using Quandl data, great little site
wi<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/WIUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date'))
 
#some minor clean up
colnames(wi)<-c('date', 'rate')
wi$date<-as.Date(wi$date)
summary(wi)
 
#base data, 1-436, test data 437-448
 
wi.b<-wi[1:436,]
wi.p<-wi[437:448,]
wi.ts<-ts(wi.b$rate, start=c(1976, 1), frequency=12)
wi.p.ts<-ts(wi.p$rate, start=c(2012, 5), frequency=12)
plot.ts(wi.ts)
 
#Lets test some models
mean<-meanf(wi.ts, 12)
naive<-rwf(wi.ts, 12)
s.naive<-snaive(wi.ts, 12)
drift<-rwf(wi.ts, 12, drift=T)
 
#linear fit
m1<-tslm(wi.ts~trend)
m2<-tslm(wi.ts~trend+season)
 
#checking for autocorrelation
res1 <- residuals(m1)
par(mfrow=c(1,2))
plot(res1, ylab="Residuals",xlab="Year")
Acf(res1, main="ACF of residuals")
 
res2 <- residuals(m2)
par(mfrow=c(1,2))
plot(res2, ylab="Residuals",xlab="Year")
Acf(res2, main="ACF of residuals")
par(mfrow=c(1,1))
 
#Durbin-Watson Test
dwtest(m1, alt="two.sided")
dwtest(m2, alt="two.sided")
#yep autocorrelation city! No surprize here, due to the nature of unemployment
 
#STL ETS Decomposition
m3<-stl(wi.ts, s.window='periodic')
plot(m3)
m4<-ets(wi.ts, model='ZZZ')
plot(m4)
 
#ARIMA
m5<-auto.arima(wi.ts)
plot(forecast(m5, h=12))

#neural networks
m6<-nnetar(wi.ts)
m6
plot(forecast(m6, h=12))
 
#Testing for accuracy the first 4 models
a1<-accuracy(mean, wi.p.ts)
a2<-accuracy(naive, wi.p.ts)
a3<-accuracy(s.naive, wi.p.ts)
a4<-accuracy(drift, wi.p.ts)
 
a.table<-rbind(a1, a2, a3, a4)
 
#Creating the forecast and accuracy for the next 6 models
f1<-forecast(m1, h=12)
f2<-forecast(m2, h=12)
f3<-forecast(m3, h=12)
f4<-forecast(m4, h=12)
f5<-forecast(m5, h=12)
f6<-forecast(m6, h=12)
 
a5<-accuracy(f1, wi.p.ts)
a6<-accuracy(f2, wi.p.ts)
a7<-accuracy(f3, wi.p.ts)
a8<-accuracy(f4, wi.p.ts)
a9<-accuracy(f5, wi.p.ts)
a10<-accuracy(f6, wi.p.ts)
 
#Combining into a table with row names
a.table<-rbind(a.table, a5, a6, a7, a8, a9, a10)
row.names(a.table)<-c('Mean', 'Naive', 'S. Naive', 'Drift', 'Lm~Trend', 
'Lm~Trend+Sea', 'STL', 'ETS', 'ARIMA', 'Neuro')
 
#make into a data frame so the best model is first, according to MAPE
a.table<-as.data.frame(a.table)
a.table<-a.table[order(a.table$MAPE),]
a.table

Created by Pretty R at inside-R.org

Results so far: Looks like the mean like forecasts are doing the best, the fancy models are not doing very well.

Spring Cleaning Data: 6 of 6- Saving the Data

With all the cleaning done, the only thing left to do is save the data to be analyzed, for future use, and I hope by others. The data I thought would be simple, but there were a few interesting twist, like the Primary Credit*, and using ifelse() to edit the districts.

I have included the product as well as the R-code in a single file for people to use and learn from. I would like to thanks all those who made comments, I find all of them helpful. Below are the links to the files generated and used in the series, and the r-code used to exporting and reloading the data.

List of files used and their links

• DiscountWindow.R -the code used to download and clean up discount window data
• Districts.csv -the table used in post 5 of 6 to seperate the districts
• DiscountWindow.csv -a csv of the final product
• DiscountWindow.RData -RData of the final cleaned up product

#Export the data, csv and RData
setwd("C:/Users Defined/")
write.csv(dw, file='DiscountWindow.csv')
save(dw, file='DiscountWindow.RData')
 
#note when loading the data the envir= needs to be defined
#with larger files the RData is definately the way to go
#this file is small enough it does not matter
load('DiscountWindow.RData', envir=.GlobalEnv)
dw<-read.csv(file.choose(), header=T)

Created by Pretty R at inside-R.org

Previous Posts (Part 1, Part 2, Part 3, Part 4, Part 5)

Spring Cleaning Data: 5 of 6- 2 ifelse vs Merge

The blog in the data cleaning series looks at separating out the Federal Reserve Districts. What I wanted was two additional columns, where I had the name of the city and the number for each district. Since I was on a separation kick I thought it would be fun to do this using ifelse() function.

Well, what started out as a fun romp in the fields turned to an exercise in precision and frustration that did end well, but took too much time, and too many lines of code to do what I wanted.

While I was banging my head against the keyboard in frustration, the thought occurred to me. Instead of using the ifelse() function, create a table with the new columns of data then merge the original data with the table just created. Two lines of code for both columns of data, definitely one of those eureka moments.

The lesson in all of this, ifelse() functions are good within a limited use, I would say 5 or less. Unless you really like doing them, then have fun. If there are limited number of occurrences like this example 12 different districts, the table works very well. What took me 2 hours of work using the ifelse() function, took me 15 minutes using the table method. The code is simpler, and easier to understand. Sure, there is the extra table to be imported, but it is small and very manageable. 

I have placed the code below, with the merge code first, followed by the ifelse() code. The table I used can be downloaded from here (District Data). Read the district data in by using the read.csv() then merge the two files using the 'district' as the column they both have in common. The ifelse(logic, true, false), the logic is if the column looks like one of the districts, if true a 1/Boston, at the end there is the 'Error' just in case.

#Merging the data

dist<-read.csv(file.choose(), header=T)
dw<-merge(dw, dist, by='district')
 
 
#re-coding the district data to numerical
tmp1<-ifelse(dw$district=='Boston (1)', 1,
   ifelse(dw$district=='New York (2)', 2,
   ifelse(dw$district=='Philadelphia (3)', 3,
   ifelse(dw$district=='Cleveland (4)', 4,
   ifelse(dw$district=='Richmond (5)', 5,
   ifelse(dw$district=='Atlanta (6)', 6,
   ifelse(dw$district=='Chicago (7)', 7,
   ifelse(dw$district=='St. Louis (8)', 8,
   ifelse(dw$district=='Minneapolis (9)', 9,
   ifelse(dw$district=='Kansas City (10)', 10,
   ifelse(dw$district=='Dallas (11)', 11,
   ifelse(dw$district=='San Francisco (12)', 12,
   'Error'))))))))))))
 
dw$dist.no<-as.numeric(tmp1)
 
 
#Isolating the names, making to factor
tmp2<-ifelse(dw$district=='Boston (1)', 'Boston',
   ifelse(dw$district=='New York (2)', 'New York',
   ifelse(dw$district=='Philadelphia (3)', 'Philadelphia',
   ifelse(dw$district=='Cleveland (4)', 'Cleveland',
   ifelse(dw$district=='Richmond (5)', 'Richmond',
   ifelse(dw$district=='Atlanta (6)', 'Atlanta',
   ifelse(dw$district=='Chicago (7)', 'Chicago',
   ifelse(dw$district=='St. Louis (8)', 'St. Louis',
   ifelse(dw$district=='Minneapolis (9)', 'Minneapolis',
   ifelse(dw$district=='Kansas City (10)', 'Kansas City',
   ifelse(dw$district=='Dallas (11)', 'Dallas',
   ifelse(dw$district=='San Francisco (12)', 'San Francisco',
   'Error'))))))))))))
 
dw$dist.city<-as.factor(tmp2)

Created by Pretty R at inside-R.org

Previous Posts (Part 1, Part 2, Part 3, Part 4)

Spring Cleaning Data: 4 of 6- Combining the files & Changing the Dates/Credit Type

So far the individual files have been left on their own, it is now time to combine using the rbind function, simple enough after all we have done so far, then the quick check with summary.

Now that we have one data frame, time to make larger changes to the data. The first is to get the dates into a format that R can understand. The as.Date() function does this by defining the variable, then the pattern for the date. At this point, I had a hard time figuring out what each one meant; basically you are defining what the date looks like now in the data frame, not in the future.

For this data set the '%b %d %Y' or in other words Feb 01 2011, if the date looked like Feb-01-2011, then the code would be '%b-%d-%Y', or if the date was 02-02-2011, then '%m-%d-%Y'. For a more comprehensive tutorial, see the post on Quick-R.

#Changing the date variables, then 
#isolating the year variable for alter use
library(stringr)
dw$loan.date<-as.Date(dw$loan.date, '%b %d %Y')
dw$mat.date<-as.Date(dw$mat.date, '%b %d %Y')
dw$repay.date<-as.Date(dw$repay.date, '%b %d %Y')

At this point, I like to have two extra variables so I can aggregate the data later for some nice results, in particular the year and the month. The reason is I want to know if there is a difference in the years.  I know there are only 2 years so far, but every quarter new data will be released so I am setting up the code for it now. The month I want to know if there is any seasonality to it. If I choose to I can isolate the day, but this gets messy because February has 28/29 days, then the rest of the months fluctuate between 30 and 31. The data is scattered and blotchy as is, making the day too small of a unit to be useful.

The code assumes the date has been changed to the R default of YYYY-MM-DD, for the year I selected the first 4 numbers using the str_sub() function, while making it a numerical value- as.numeric(). The year and date variable I made it a factor for easier sorting and categorizing, with a similar process as above except I want both.

#Create a year variable
dw$year<-as.numeric(str_sub(dw$loan.date, start=1, 
   end=4))
 
#Create a year and month variable
dw$year.month<-as.factor(str_sub(dw$loan.date, 
   start=1, end=7))

The next step is to change the credit type to something simpler for tables and graphs. I used the gsub, one of the most interesting and fun functions I never knew existed until I did this. Basically it will take a string then replace it with another. For this data I wanted to replace the "Primary Credit" with "primary" because it make things so much easier for graphs and tables. Then I changed it to a factor.

 
#Changing the type of credit to one word
dw$type.credit<-with(dw, 
   gsub("Primary Credit", 'primary', type.credit))
dw$type.credit<-with(dw, 
   gsub("Seasonal Credit", 'seasonal', type.credit))
dw$type.credit<-with(dw, 
   gsub("Secondary Credit", 'secondary', type.credit))
 
#change to factor
dw$type.credit<-as.factor(dw$type.credit)
summary(dw)

Created by Pretty R at inside-R.org

Links to the previous posts (post 1, post 2, post 3)

Spring Cleaning Data: 3 of 6- The Little but Big Correction

Building on the previous posts (post 1 & post 2) I found there were 12 instances with the type of credit where there was a "Primary*" which means the lender borrowed twice in the same day, in the 2010 q4 data. It would seem simple enough in Excel, use the filter, find the 12 instances, delete the "*" and be done. For R this turned out to be much more difficult to do. Mainly because the "*" is an operator used in multiplication, so it is hard to get rid of. What seemed to be a minor correction turned into something much bigger. Why would I go through all the trouble?

Mainly because I know that there will be cases where I will be dealing with data much larger than this one, that will not fit into Excel. Where I will have to be able to make corrections in R because I will not be able to do it otherwise. Plus, it is a good way to hone my data manipulation skills.

I used the following, but failed doing it (suggestions for improvement are welcomed):

• Gsub()- did not work because of the '*'
• Variable[condition]<-replacement- looking for a numerical value, which there is none
• Ifelse()- did not find the difference between Primary and Primary*
• Used the stringer package to try and delete the last character, but did not work

So what I did was use the subset() function to flush the problem variables out. Knowing there was 12 I replaced the wrong variables with the correction, using the rep() function. Then I made sure it was a factor using the as.factor().

#Correcting the Primary Credit* 
#Step 1 isolate the rows of data
tmp<-subset(x=dw.2010.q4, 
   subset=type.credit=='Primary Credit*')
 
#Step 2 change the column of data to the correct label
tmp$type.credit<-rep('Primary Credit', 12)
 
#Step 3 make sure the data is a factor
tmp$type.credit<-as.factor(tmp$type.credit)
summary(tmp)

I then erased the problem variables from the original data using the file[-which(file$variable=='Primary Credit*'),]. Note the use of the '-' before the which to remove those rows. Then I added the cleaned up rows to the data. I am not particularly interested in keeping the data in any particular order as the date will keep things in chronological order.

#Step 4 remove the Primary Credit* 

#rows of data from original dataframe
dw.2010.q4<-dw.2010.q4[-which(
   dw.2010.q4$type.credit=='Primary Credit*'),]
 
#Step 5 add the corrected data back into dataframe
dw.2010.q4<-rbind(dw.2010.q4, tmp)
 
#Check to see if correct
summary(dw.2010.q4)

Created by Pretty R at inside-R.org

Spring Cleaning Data: 2 of 6- Changing Column Names and Adding a Column

The first post (found here) we downloaded the data and imported it to R using the gdata package. This post we will be changing the column names to make them more reasonable, and adding a quarter variable. The reason for changing the column names is because the dw.2010.q1 file column names are messed up due to the formatting done in Excel. So if I was going to have to change one, just as well change them all, so i did.

The first chunk of code defines the labels I am going to use as c.label. Then I used the colnames() function to rename each file.

#Defining the new labels
c.label<-c('loan.date', 'mat.date', 'term',
   'repay.date', 'district', 'borrower', 'city',
   'state', 'ABA', 'type.credit', 'i.rate',
   'amount', 'outstanding.credit',
   'total.outstanding', 'collateral',
   'commercial', 'residential.morg',
   'comm.real', 'consumer', 'treasury',
   'municipal', 'corp', 'mbs.cmo',
   'mbs.cmo.other', 'asset.backed',
   'internat', 'tdfd')
 
#Changing the column names
colnames(dw.2010.q3)<-c.label
colnames(dw.2010.q4)<-c.label
colnames(dw.2011.q1)<-c.label

I also like to add a few additional variables when I see a potential need when I can. At this point the files are individual, and adding the quarter variable might be helpful. Sure I could write a loop to create the new column based on the month of the date, but I like to keep things as simple as possible. Why add complexity when there is no reason. I used the ABA to define the length of the data set because it did not have any missing values, while others did. The new column name is qtr, and the function rep() is used to repeat the quarter number the length of the column ABA.

#defining a quarter variable for future use, so I can 
#isolate quarters to compare and contrast
dw.2010.q3$qtr<-rep(3, length(dw.2010.q3$ABA))
dw.2010.q4$qtr<-rep(4, length(dw.2010.q4$ABA))
dw.2011.q1$qtr<-rep(1, length(dw.2011.q1$ABA))

Created by Pretty R at inside-R.org

Spring Cleaning Data: 1of 6- Downloading the Data & Opening Excel Files

With spring in the air, I thought it would be fun to do a series on (spring) cleaning data. The posts will follow my efforts to to download the data, import into R, cleaned it up, merge the different files, add columns of information created, and then a master file exported. During the process I will be offering at times different ways to do things, this is an attempt to show how there is no one way of doing something, but there are several. When appropriate I will demonstrate as many as I can think of, given the data.

This series of posts will be focusing on the Discount Window of the Federal Reserve. I know I seem to be picking on the Feds, but I am genuinely interested in what they have. The fact that there is data on the discount window is, to be blunt, took legislation from congress to get. The first step in this project was to find the data. The data and additional information can be downloaded here.

The data is in Excel format, with 3 sheets, the first sheet is the information I am interested in, the second is a definition sheet, and then the last is a list of interest rates (which are in the first sheet). At the present date there are 3 excel files, 2010- 3 quarter, 2010- 4 quarter, and 2011- 1 quarter. The data is made available with about a 2 year lag.

There are 2 approaches when it comes to downloading the data.

1. Download the data, clean it up in Excel, save the one sheet as a *.csv and open using read.csv(file, header=T)
2. Download the data using R, import using the gdata package, and clean the data up using R

With a data set where I only have a couple of files to do, and they are relatively small (less than 1,000 rows of data) I will chose option 1. But if there are more than like 10 files, with more than 1,000 rows of data, I will start using R more. With the discount window I would normally use Excel because the data does not come out very often, and the files are not that big. 

But where would the fun be in programming in R?

The first thing to do is download Perl (assuming you don't have it) from this link.

The reason is the gdata package needs to have perl to work. Next install the gdata package, and then find the location for your perl.exe (if windows) and have it handy.

The code below begins with setup the libraries needed, then defining the links to be used. The final section of code is the gdata function to open the files that are in Excel.

# Federal Reserve Data- Discount Window
# http://www.federalreserve.gov/newsevents/reform_discount_window.htm
# Downloading and coding Discount Window data 

# from the Federal Reserve
library(gdata)
 
# Defining the individual links
link1<-"http://www.federalreserve.gov/newsevents/files/dw_data_2010_q3.xls"
link2<-'http://www.federalreserve.gov/newsevents/files/dw_data_2010_q4.xls'
link3<-'http://www.federalreserve.gov/newsevents/files/dw_data_2011_q1.xls'

The dw.2010.q3 is the label for the file, the read.xls is the gdata function, link1 refers to the internet address defined above, sheet=1 means we want the first sheet, then I need to skip the first 3 lines, perl='' is the location of the perl.exe that I mentioned would be needed.

# Downloading the files and opening in R
# Make sure to have Perl installed before using the 

# read.xls() function
dw.2010.q3<-read.xls(link1, sheet=1, skip=3, perl="C:/location/bin/perl.exe")

dw.2010.q4<-read.xls(link2, sheet=1, skip=3, perl="C:/location/bin/perl.exe")
dw.2011.q1<-read.xls(link3, sheet=1, skip=3, perl="C:/location/bin/perl.exe")
 
# checking to see if it worked
summary(dw.2010.q3)
summary(dw.2010.q4)
summary(dw.2011.q1)

Created by Pretty R at inside-R.org

Great Circles, Black Holes, and Community Events Part 3 of 3

The second community event is the Soldier Hollow Junior Olympics (SoHo), again found in the Heber Valley area. Building upon the previous posts (part 1 and part 2) this one will show an event that has more people coming from greater distance. Take the bar charts for the number of participants and the cities they are from. Instead of 2 major cities (Heber Valley Railroad), SoHo has several cities and states with many participants for each.

The histogram for distance shows a similar pattern, where with the railroad it was a nice log looking distribution, this is a little more even. The second histogram is a zoomed in and the bins expanded for greater detail.

While the map is not as great as the railroad, which is why the distance histogram is so important, it does show a good representation of the northwestern states. Unlike the railroad map, each line does represent more participants than 3.

What makes this tool so great is the ability to visually show the interaction between distance, number of participants, and the draw of an event. The numbers are nothing really new, but the charts are what make this analysis shine. When talking to community representatives whose education range from high school graduate to PhD, pictures are critical.

#Soldier Hallow Analysis

soho<-read.csv(file.choose(), header=TRUE)
summary(soho)
table.city<-sort(table(soho$city), decreasing=TRUE)
table.st<-sort(table(soho$state), decreasing=TRUE)

par(mar=c(5, 11, 4, 2), las=2)
barplot(table.city, main=‘SoHo: Cities’, horiz=TRUE, col=‘red’)

par(mar=c(5, 4, 4, 2), las=2)
barplot(table.st, main=‘SoHo: States’, horiz=TRUE, col=‘red’)

heber<-c(-111.33259, 40.511413)
soho.data<-matrix(data=c(soho$long, soho$lat), nrow=373, ncol=2)
soho.ut<-subset(soho, subset=(state==‘UT’))
soho.data.ut<-matrix(data=c(soho.ut$long, soho.ut$lat), nrow=29, ncol=2)
soho.dist<-(distm(heber, soho.data, fun=distVincentyEllipsoid)*0.000621371192)
soho.dist.ut<-(distm(heber, soho.data.ut, fun=distVincentyEllipsoid)*0.000621371192)
dist.soho<-matrix(soho.dist, nrow=373, ncol=1)
dist.soho.ut<-matrix(soho.dist.ut, nrow=29, ncol=1)
summary(dist.soho)
sd(dist.soho)
p.skew.soho<-(3*(mean(dist.soho)-median(dist.soho)))/sd(dist.soho)
hist(dist.soho, main=‘SoHo: Distance Histogram’, col=‘red’)
hist(dist.soho.ut, main=‘SoHo: Distance Histogram Utah’, breaks=20, col=‘red’)

#mapping it out
#US
map("state", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25)
title(main=‘SoHo: US Map’)
for(i in 1:dim(soho.data)[1]){
 inter <- gcIntermediate(heber, soho.data[i, 1:2], n=373, addStartEnd=TRUE)
 lines(inter, col="red")
}
 
#Zoomed into West
par(mfrow=c(1,2), mar=c(5,4,4,2))
map("state", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25, xlim=c(-125, -103), ylim=c(30, 50))
title(main=‘SoHo: Western Region’)
 
for(i in 1:dim(soho.data)[1]){
 inter <- gcIntermediate(heber, soho.data[i, 1:2], n=373, addStartEnd=TRUE)
 lines(inter, col="red")
}

#Utah
map("state", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25, xlim=c(-112.1, -111), ylim=c(40, 42))
title(main=‘SoHo: Utah’)
for(i in 1:dim(soho.data.ut)[1]){
 inter <- gcIntermediate(heber, soho.data.ut[i, 1:2], n=29, addStartEnd=TRUE)
 lines(inter, col="red")
}
par(mfrow=c(1,1))

Created by Pretty R at inside-R.org

Great Circles, Black Holes, and Community Events Part 2 of 3

This post will examine the Heber Valley Railroad, a small town tourist attraction using event gravitational pull. Using the information from part 1 the two factors associated with the events gravity, the number of participants, and the distance they traveled. The number of participants can be shown using bar charts, histograms, and summary tables. The distance traveled can be displayed using bar charts, histograms, summary charts, and most important great circle maps. Below are some of the charts I created when doing the analysis.

The first two are bar charts show where the majority of train riders are coming from. From the graphs Salt Lake City Utah is number one, followed by Ogden, then it goes down hill from there very quickly. The purpose of these first graphs is to show how many people are coming to the event and where.

Where the first set of graphs show the cities, the histograms show where people are coming from in terms of distance. The first histogram shows all the data, the next two are zoomed in to show the majority of people travel less than 100 miles (as the crow flies) to get to the train. While the majority of people who ride the train do so within 100 miles, there are many who travel many miles. But it should be noted, there are only 1-3 customers per line. The map looks really good, but the majority of those traveling long distances, are much fewer.

While the map shows a number of people coming from various locations, there are only 1-3 people per line, nothing to really spend any marketing funds to. The map does show how far people do come.

The next two maps show that within Utah the majority of people riding the train come from Salt Lake, Ogden, and the Wasatch Front area.

require(geosphere)
require(maps) 
 
#HVRR Analysis
#Step 1: basic Stats. Summaries, Histograms, bar charts
#reading the file in
 
hvrr<-read.table(file.choose(), header=TRUE)
 
#summary stats
summary(hvrr)
 
#histograms
par(mfrow=c(1,2))
label.1<-c(‘Utah (424 88%)’, ‘Other(58 12%)’)
state<-c(424, 58)
barplot(state, names.arg=label.1, main=‘HVRR: States’, col=‘blue’)
label.2<-c(‘Salt Lake’, ‘Ogden’, ‘Other’)
cities<-c(173, 102, 149)
barplot(cities, names.arg=label.2, main=‘HVRR: Cities Within Utah’, col=‘blue’)
par(las=2, mar=c(5,12,4,2), mfrow=c(1,1))
city.1<-sort(table(hvrr$city))
city.1<-tail(city.1, n=20)
barplot(city.1, col=‘blue’, hor=TRUE, main=‘HVRR: Utah Cities Top 20’)
par(las=0, mar=c(5,4,4,2))
 
#distance analysis
heber<-c(-111.33259, 40.511413)
data<-matrix(data=c(hvrr$long, hvrr$lat), nrow=482, ncol=2)
ut<-subset(hvrr, subset=(st==‘UT’))
data.ut<-matrix(data=c(ut$long, ut$lat), nrow=424, ncol=2)
dist<-(distm(heber, data, fun=distVincentyEllipsoid)*0.000621371192)
dist.rr<-matrix(dist, nrow=482, ncol=1)
hvrr<-cbind(hvrr, dist.rr)
 
#histograms of various shapes and zooms
summary(dist.rr)
par(mfrow=c(1, 3))
hist(dist.rr, breaks=12, main=‘HVRR Distances: 0-3,000 miles’, xlab=‘Distance in Miles’, col=‘blue’)
hist(dist.rr, breaks=24, main=‘HVRR Distances: 0-500 miles’, xlab=‘Distance in Miles’, xlim=c(0, 500), col=‘blue’)
hist(dist.rr, breaks=50, main=‘HVRR Distances: 0-200 miles’, xlab=‘Distance in Miles’, xlim=c(0, 200), col=‘blue’)
par(mfrow=c(1,1))
 
#mapping it out
 
#US
map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25, xlim=c(-158, -65), ylim=c(15, 50))
title(main=‘HVRR: US Map’)
for(i in 1:dim(data)[1]){
 inter <- gcIntermediate(heber, data[i, 1:2], n=482, addStartEnd=TRUE)
 lines(inter, col="blue")
}
 
#Zoomed into Utah
par(mfrow=c(1,1), mar=c(5,4,4,2))
map("state", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25, xlim=c(-115, -108), ylim=c(37, 42))
title(main=‘HVRR: Utah’)
for(i in 1:dim(data.ut)[1]){
 inter <- gcIntermediate(heber, data.ut[i, 1:2], n=424, addStartEnd=TRUE)
 lines(inter, col="blue")
}
 
#Wasatch Front
map("state", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.25, xlim=c(-112.5, -111), ylim=c(40, 42))
title(main=‘HVRR: Utah- Wasatch Front’)
for(i in 1:dim(data.ut)[1]){
 inter <- gcIntermediate(heber, data.ut[i, 1:2], n=424, addStartEnd=TRUE)
 lines(inter, col="blue")
}
par(mfrow=c(1,1))

Created by Pretty R at inside-R.org