Saturday, April 12, 2014

Working with Data in R -- Import, Manipulate, Export -- Part 1

Goal

Goal

The goals of this post are to demonstrate how to (1) import data; (2) manipulate vectors and data.frames to create and store new variables (from existing data); (3) calculate groupwise summary statistics; and (4) export data.

This post will use a non-trivial example for demonstration purposes. Let's say you are interested in determine the extent to which unemployment rates in Califoria counties are related to the national unemployment rate. Historical unemployment data can be downloaded for California and sub-state geographies here. Using this tool, I downloaded all available unemployment data for all California counties. This will be the primary data used for this post – data for the U.S. writ large will be addressed at a later point in time (trust me, it's much easier to obtain!).

Importing data

The first step is to import data – to tell R to read a structured dataset. This can actually be quite tricky because idiosyncracies of operating systems can get in the way. I exclusively use Windows machines to run R so some of the tricks I show in terms of navigating a system may only apply to Windows. But there are plenty of sources online that can help with specific problems – Google search is your friend! That being said, these problems should only apply to locating files in a system – they will not affect how to program with R.

In order to tell R to read in a file, one must tell the machine where to look for said file. Programmers refer to the location of a file on a computer as the “file path”. Folders are called “directories”. While one can read in files by using a specific file path each time, I find it more efficient to set working directories instead. So instead of telling R to “import file 'A' from path 'x/y/z/½/3', I prefer so tell R that "I'm going to be working in location 'x/y/z/½/3', now import file 'A' or 'B' or 'C'”.

How does one locate and refer to a directory? To locate a directory in windows, simply keep clicking until you find the file that you downloaded. Yes, that is a cheeky answer but this post assumes one can perform basic operations on Windows.

This is what you should see once you're done clicking. See the row at the top of the box? Click that row. This is what you should see now. Copy the highlighted text (CTRL + C or right-click) – that is your working directory. However, note that Windows directories use backslashes ("\") to separate folders, whereas R uses forward slashes ("/), so just change all backslashes to forward slashes and you will be all set. Now we can do some programming!

A side note – most of the learning experience is doing so I will not make life easy and provide the file that I downloaded. Follow the link I provided, use the download tool, download the file, create a well-organized location for your project(s), and move the file there.

# Set working directory
setwd("C:/Users/taylor/Dropbox/democratic-Data/Data/Unemployment/California/County")

# print the working directory to confirm the change
getwd()
[1] "C:/Users/taylor/Dropbox/democratic-Data/Data/Unemployment/California/County"

# Examine the contents of the directory
list.files() # this gives the names of files and folders in the given location. The name of the file that you downloaded will show up here
[1] "DA2014150.txt"

## Read in the data
# The website exports comma delimited .txt files, which are values separated by commas. Comma separated values saved as .txt or .csv are the most common type of files used with R and many other applications. Familiarize yourself with how they look.

# One simply has to tell R where the data is and how the data is formatted in order to read it in
args(read.table)
function (file, header = FALSE, sep = "", quote = "\"'", dec = ".", 
    row.names, col.names, as.is = !stringsAsFactors, na.strings = "NA", 
    colClasses = NA, nrows = -1, skip = 0, check.names = TRUE, 
    fill = !blank.lines.skip, strip.white = FALSE, blank.lines.skip = TRUE, 
    comment.char = "#", allowEscapes = FALSE, flush = FALSE, 
    stringsAsFactors = default.stringsAsFactors(), fileEncoding = "", 
    encoding = "unknown", text) 
NULL

CA_unemployment = read.table("DA2014150.txt", # make sure to include the extension (the .txt part)
                             header = T, # this says that there is a series of column names for the data
                             sep = ",", # this says that the separator for the values is a comma
                             stringsAsFactors = F) # just do this, it will make your life easier. The curious can Google "factors in R" to understand why. 

# What type of object is the data stored as?
class(CA_unemployment)
[1] "data.frame"

# Make sure the data seems to have been read in properly -- compare what R "sees" and what you see when you open the file in windows
head(CA_unemployment, 3)
  Year Period           Area Adjusted Preliminary Labor.Force Employment
1 2014    Jan Alameda County  Not Adj  Not Prelim     780,600    727,600
2 2014    Feb Alameda County  Not Adj      Prelim     782,000    729,900
3 2013 Annual Alameda County  Not Adj  Not Prelim     783,100    725,000
  Unemployment Unemployment.Rate
1       53,000               6.8
2       52,000               6.7
3       58,000               7.4
tail(CA_unemployment, 3)
      Year Period        Area Adjusted Preliminary Labor.Force Employment
16168 1990    Oct Yuba County  Not Adj  Not Prelim      21,700     19,700
16169 1990    Nov Yuba County  Not Adj  Not Prelim      21,600     19,200
16170 1990    Dec Yuba County  Not Adj  Not Prelim      21,500     18,800
      Unemployment Unemployment.Rate
16168        2,000               9.2
16169        2,400              11.3
16170        2,700              12.6

# Everything looks good -- data has been successfully imported!

Summary

That is all for today – part 2 will demonstrate how to manipulate data and how to create/store new variables. As always, please feel free to email me at democraticdata@gmail.com if you have any questions.

Wednesday, March 19, 2014

Education vs. Infant Mortality, Part 1

Goal

Goal

The goal of this post will be to demonstrate how to use R to analyze an open data source and to produce high quality output in a short amount of time. Specifically, this will be a multi-part post that seeks to examine the relationship between education of females and infant mortality on a global scale. I will start with directly with analysis and output and and the end of the series, I will make all of my code available for review.

The Data

Data were obtained from the World Bank, which offers an API or application programmers interface, which greatly reduces the time and effort required to locate and download data series of interest. The WDI package in R, authored by Vincent Arel-Bundockin, R provides a useful interface with World Bank Data. In this case, I simply started with two search terms, “education” and “mortality” and I was able to find the unique identifiers and descriptions for the data I was interested in.

Data were extracted from all available countries from 1980 to 2014. However, it is often the case that global macro data has a significant amount of missingness and this is no exception. Missing data doesn't mean that it has been lost per se – it means that there was no official record of a particular data point for a particular country in a particular year. Missing data can occur for any number of reasons. Some causes of missingness can actually become a problem when trying to estimate relationships between variables if the cause of the missingness is related to the outcome of interest.

In this case, the question to ask is: are countries/years with more missing data more likely to have higher or lower rates of infant mortality? Without having done any analysis to investigate this question, I would assume that the answer is yes. The reasons: (a) data is more likely to be missing in earlier time periods because many countries in the sample were less likely to track and record macro statistics and (b) there is reason to suspect (putting aside my own prior knowledge of the subject) that infant mortality rates have decreased over time. At the very least, this means that estimates of the trend in infant mortality rates over time will be understated.

Comparing Female Eduation to Infant Mortality

Graphical Analysis

A direct comparison between the primary outcome of interest (infant mortality) and the primary predictor (female education) is always a good place to start. A visual inspection of the data allows one to infer a great deal about subsequent analysis. The first things to look for are the type, strength, and direction of a relationship. That being said, I will produce another post on the topic of assessing functional form, strength, and direction of relationships so this post can be more focused on the analysis at hand.

First, note that each plot on the graph represents a country-year, i.e. “the United States in 2000” has one data point on the graph. In addition, the purple line is the “least squares regression line,” which is a common method of estimating and summarizing the relationship between variables. In this case, the relationship between (a lack of) female education and infant mortality is linear, relatively strong, and postive. That is to say, as the proportion of females aged fifteen and over without education increases, so too does the infant mortality rate and the mortality rate of children five and younger. There is some clustering at the bottom left at the plot, which is probably the domain of wealthier, more “developed” countries. Subsequent investigations should determine the extent to which the relationship between female education and mortality differs between developed and developing nations.

plot of chunk main_plot

Estimate the Relationship

After taking an initial look at primary variables of interest, it is helpful to estimate the strength of the relationship. In this case, I want to estimate the slope of the purple line that appears in the plot above. The questions is, for every percent increase in females without education, what is the associated increase in the number of infant deaths per 1,000 births? First note that a causal relationship is not assumed here – this pattern can be caused by a number of factors. In order to make causal inferences, one must first consider the plausible causal processes that link x to y and one must do at least a reasonable job of adjusting for other factors that also have an impact on infant mortality. For instance, it may be that countries that are more likely to have a high proportion of uneducated females are also more likely to have poor access to adequate nutrition and health care, which should also have an impact on mortality rates. In that case, it may be that education of males may have a similar relationship to child mortality

That being said, regressions can still be useful investigative tools, even if causal inferences cannot necessarily be made from the results.

Below is a table of regression results:

————————————————–

Outcome Variable Estimate Std. Error t value Pr(>|t|)
Mortality rate, infant (per 1,000 live births) (Intercept) 12.318 1.17286 10.503 1.255e-22
Mortality rate, infant (per 1,000 live births) df$Female_Education 1.064 0.04126 25.792 4.035e-83
Mortality rate, under-5 (per 1,000 live births) (Intercept) 13.469 1.95098 6.904 2.373e-11
Mortality rate, under-5 (per 1,000 live births) df$Female_Education 1.786 0.06863 26.029 4.855e-84

————————————————–

The “Outcome” column indicates which outcome was used in the regression. The “Variable” column indicates which variable was estimated. The “Estimate” column includes the estimates of the slope for a particular variable. Ignore the other columns for now. The “(Intercept)” variables indicate what value of y (when x is zero) best aids in producing the best fit to the data – this doesn't help in interpreting results so it is usually just ignored. The “df$Female_Education” variable is the main predictor, the proportion of females aged fifteen and over with no education. In order to aid in interpretation of estimates, I multiplied this proportion by 100 (to convert the proportion to percent).

The estimate for the first outcome is 1.064. This means, for a 1 unit (percent) increase in x (prorportion of uneducated females), one should expect an increase of 1.064 y (infant mortality rate per 1,000 live births). Therefore, for every 10% increase/decrease in the proportion of uneducated females in a given country-year, one should expect an increase/decrease of 10.64 infant deaths per 1,000 live births. The estimate for the under five mortality rate is 1.786 – this value has the same interpretation.

The “Pr(>|t|)” column indicates the probability that the relationship implied by the estimates reported here are due to chance. In this case, the “P values” are 4 x 10-83 and 4 x 10-84, which are incredibly small numbers. Depending on the subject at hand, the conventional P value for which something is deemed “statistically significant” is 0.05, and these are well below that. This simply means that the association we have observed is unlikely due to chance alone.

Summary

That is all for Part 1. The next step will be to produce summary statistics and associated visualizations to examine underlying characteristics of the data. As always, feel free to email me at democraticdata@gmail.com if you have any questions.

————————————————–

Friday, March 14, 2014

Tutorial on basic R objects and functions

Goal

Goal

The goals of this post are to (1) introduce basic R objects and functions and (2) to show how useful operations can be performed.

vectors, matrices, and data.frames

Below, I will define three types of objects, examine them, and perform operations on them. Note that “#”, without the parentheses, is R's comment character. That means, any text after # on a given line will NOT be interpreted by R. It is important to use comments to explain what the code is doing, both for the programmer's benefit and for anyone that might use the code later.

### Defining a vector
vector_1 = c(1:100)

### How many elements are in the vector?
print(length(vector_1))
[1] 100

### What are the first 5 elements of a vector?  There are multiple ways to
### extract this information, but I will show two methods:

# Method 1
print(vector_1[1:5])
[1] 1 2 3 4 5

# Method 2
print(head(vector_1, 5))
[1] 1 2 3 4 5

### Define objects that equal the mean, median, minimum, maximum, 25th
### percentile, 75th percentile, and sum of the vector
mean_vector_1 = mean(vector_1)
median_vector_1 = median(vector_1)
max_vector_1 = max(vector_1)
min_vector_1 = min(vector_1)
sum_vector_1 = sum(vector_1)
q_25_vector_1 = quantile(vector_1, probs = 0.25)
q_75_vector_1 = quantile(vector_1, probs = 0.75)

### How can this information be extracted more efficiently?  Many R objects
### have summary functions that are quite useful -- it doesn't hurt to try the
### summary() function on an object to see what values are returned.
summary_vector_1 = summary(vector_1)

print(summary_vector_1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0    25.8    50.5    50.5    75.2   100.0 

### What if I wanted to store this information for later use or for export? I
### will define a data.frame for storage
summary_df_vector_1 = data.frame(Mean = mean_vector_1, Median = median_vector_1, 
    Max = max_vector_1, Min = min_vector_1, Sum = sum_vector_1, Q25 = q_25_vector_1, 
    Q75 = q_75_vector_1, stringsAsFactors = F)

row.names(summary_df_vector_1) = NULL

### Let's take a look at the object I just defined (you may have guessed by
### now that the equals sign in 'new thing = some other thing') is used for
### object assignment/definition
print(summary_df_vector_1)
  Mean Median Max Min  Sum   Q25   Q75
1 50.5   50.5 100   1 5050 25.75 75.25

# Object class
print(class(summary_df_vector_1))
[1] "data.frame"

# What are the objects dimensions?
print(dim(summary_df_vector_1))
[1] 1 7

# How many rows?
print(nrow(summary_df_vector_1))
[1] 1

# Columns?
print(ncol(summary_df_vector_1))
[1] 7

### Note that I didn't examine the dimensions of the original vector_1 object
### -- I only examined the length.  This is because, vectors DO NOT HAVE
### dimensions == they only have length. Observe:
print(dim(vector_1))
NULL

# NULL' in R literally means NOTHING or NO VALUE so there are no dimensions

### What values does the Sum variable in the summary_df_vector_1 have? There
### are multiple ways to extract the values of a vector that is held within a
### data.frame

# Method 1
print(summary_df_vector_1$Sum)
[1] 5050

# Method 2
print(summary_df_vector_1[, "Sum"])
[1] 5050

### I can also define a new object using data from summary_df_vector_1
extracted_sum = summary_df_vector_1$Sum

### Is this new vector the same as the original vector?

# One test:
print(identical(extracted_sum, sum_vector_1))
[1] TRUE

# Another:
print(extracted_sum == sum_vector_1)
[1] TRUE

### Note that I used ' == ' and not ' = ' in the instance above.  The ' == '
### is a LOGICAL OPERATOR and tests the TRUTH of some statement, whereas' = '
### assigns a value to some object.

How do functions work?

The first section used several functions but without much explanation. This section demonstrates how one uses R functions. R Functions have names and they have arguments. Some arguments are defined by default and others require additional input from the user. Users can also program their own functions to perform a series of tasks as well as to return output from the work of the function.


### Virtually all R functions have some sort of documentation in help files or
### online that demonstrate how they work, along with examples.  To query the
### help files on one's computer, simply type '?' along with the function name
### into the command line and run the code.

# Example: ?data.frame

### SINCE FUNCTION NAMES AND OBJECT NAMES ARE CASE SENSITIVE IN R, SOMETIMES
### IT IS NECESSARY TO SEARCH THE WEB FOR INFORMATION REGARDING A PARTICULAR
### FUNCTION google: 'data.frame in R' and you'll get some helpful pages as
### the top results.  These help pages show what arguments the function has,
### how they change things, (sometimes) some background on the function, and
### (sometimes) helpful examples the user can run on their own.

### Here is how one can extract the arguments from a function.
print(args(data.frame))
function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, 
    stringsAsFactors = default.stringsAsFactors()) 
NULL

### Arguments can be 'passed' to a particular function 'positionally' or by
### directly assigning values to various arguments.

# Example: rnorm() -- the rnorm() function is one of the most commonly used
# functions in R (and in R help literature!) as it performs random draws
# from a normal distribution with mean and standard deviation parameters
# that are specified by the user.  Let's see how this function works!

print(args(rnorm))
function (n, mean = 0, sd = 1) 
NULL

### The function has arguments of n, mean, and sd. I will define an object
### using different methods of defining object parameters.

# Positional
normal_draw_1 = rnorm(1000, 10, 1)

# Defining each argument
normal_draw_2 = rnorm(n = 1000, mean = 10, sd = 1)

# Both positional and defined
normal_draw_3 = rnorm(1000, mean = 10, sd = 1)

### All three of those objects are the same (bonus points if you test my
### assertion yourself from functions I used in the first code chunk!). Note
### that when the arguments were printed, the mean and sd parameters were
### assigned values of 0 and 1 respectively.  Those are the function defaults,
### which means that, as long as n is assigned, the function can operate
### properly and it will use the default values of mean and sd if none are
### supplied.

Summary

That is all for today. I will end the lesson with a little advanced coding just to give folks something to look forward to (and perhaps cheat ahead a little!). Feel free to email me at democraticdata@gmail.com if you have any questions – I will do my best to respond within a reasonable amount of time. Also, note that there is an excellent listserve that handles everything related to R. Here is the link to sign up to the R help listserve: https://stat.ethz.ch/mailman/listinfo/r-help.


standard_normal_distribution = rnorm(1000)

### Now, I'm going to cut ahead a little and show you some cool things R can
### do -- don't worry about understanding what's going on yet but bonus points
### if you can figure it out!  Loading libraries
require(ggplot2)
## Loading required package: ggplot2
require(reshape2)
## Loading required package: reshape2

# define new data.frame for plotting
combined_normal_draws = data.frame(Index = c(1:length(normal_draw_3)), random_draw = normal_draw_3, 
    standard_normal = standard_normal_distribution)

# take a look
print(head(combined_normal_draws))
##   Index random_draw standard_normal
## 1     1       8.733          0.2096
## 2     2      10.816          0.7517
## 3     3       9.509         -0.7431
## 4     4      10.857         -0.3762
## 5     5      10.560          0.9273
## 6     6       9.657         -0.7291

# Create 'Long' format for ggplot2
long_combined_normal_draws = melt(combined_normal_draws, id = "Index", measure = c("random_draw", 
    "standard_normal"))
# Change variable names for plotting
long_combined_normal_draws$variable = as.character(long_combined_normal_draws$variable)
long_combined_normal_draws$variable[long_combined_normal_draws$variable == "standard_normal"] = "Standard Normal\n(Mean 0, SD 1)"
long_combined_normal_draws$variable[long_combined_normal_draws$variable == "random_draw"] = "Random Draw\n(Mean 10, SD 1)"

# take a look
print(head(long_combined_normal_draws))
##   Index                     variable  value
## 1     1 Random Draw\n(Mean 10, SD 1)  8.733
## 2     2 Random Draw\n(Mean 10, SD 1) 10.816
## 3     3 Random Draw\n(Mean 10, SD 1)  9.509
## 4     4 Random Draw\n(Mean 10, SD 1) 10.857
## 5     5 Random Draw\n(Mean 10, SD 1) 10.560
## 6     6 Random Draw\n(Mean 10, SD 1)  9.657
print(tail(long_combined_normal_draws))
##      Index                        variable   value
## 1995   995 Standard Normal\n(Mean 0, SD 1)  1.1738
## 1996   996 Standard Normal\n(Mean 0, SD 1)  0.4880
## 1997   997 Standard Normal\n(Mean 0, SD 1) -0.8823
## 1998   998 Standard Normal\n(Mean 0, SD 1) -0.6630
## 1999   999 Standard Normal\n(Mean 0, SD 1)  1.1272
## 2000  1000 Standard Normal\n(Mean 0, SD 1)  1.9558

# use a scale parameter
scale = 2
ggplot(long_combined_normal_draws, aes(value)) + facet_wrap(~variable, scales = "free") + 
    theme_bw() + labs(title = "Density plots of random draws from \nthe normal distribution\n", 
    x = "\n", y = "Density\n") + theme(axis.text = element_text(size = 12 * 
    scale, face = "bold"), strip.text = element_text(size = 15, face = "bold"), 
    strip.background = element_rect(fill = "lightblue"), axis.title = element_text(size = 11 * 
        scale, face = "bold"), title = element_text(size = 13 * scale, face = "bold")) + 
    stat_density()

plot of chunk advanced_fun