Defaults, Lists and Classes: A Functional Post

In this post, I demonstrate a couple of useful tricks to writing functions in R. The context is a function I wrote called called samp() that allows for an easy demonstration of sampling distribution properties.



Defaults

By default, this function draws K = 500 samples of size N from a normal distribution with a mean of mu = 0 and a standard deviation of std = 1. As you can see, the syntax for specifying a default for an argument to the function is to set it equal to its default value inside the function statement. Arguments can be given no default by naming them as an argument to the function without specifying their default value.

Sensibly chosen defaults can save a lot of time for people who use your function. For example, although this samp() function has four arguments, I don't need to type samp(100, 0, 1, 500) every time I wish to take 500 samples of size N=100 from a Normal(0,1) random variable. If sampling from a Normal(0,1) is especially common, setting the defaults this way can save a lot of keystrokes.

On the other hand, using defaults instead of hard coding the Normal(0,1) choice has the distinct advantage that users get flexibility from your function giving them the option to specify arguments other than the default. The right use of defaults strikes the perfect balance between ease of use and flexibility.

Returning a List... with Names

If your function warrants the special treatment, you have probably performed several useful calculations that you would like to return. Maybe you have conducted several related hypothesis tests that you would like to reference later (report and/or use in another calculation). A great way to do this is to store the results of your function into a list.

For example, in the samp() function, I wanted to store vectors of the sample means, sample standard deviations and standard errors (500, one for each replication). For good measure, I wrote the function to store the means and standard deviations of these vectors as well as the true parameter values that went into the function.

Although it is often a good idea to store the results in a list, it is a better idea to name the elements of the list for easy extraction. If I ran the command

mysamp = samp(100)

The object mysamp would contain a list of the results from the function samp(). Suppose I want to extract the vector of means (for example, to compute a histogram). Without running the names command, I would have to remember that the vector of means was in the first position in the list and use the command

hist(mysamp[1])

to plot a histogram of the means vector from my sampling object. As it is easy to forget the precise order in which the elements of a function are stored, this syntax can lead to too much thinking. It is better to name the first list element something like means. Naming the list elements allows the user to use the $ extractor on objects created with the function. That is, the syntax becomes

hist(mysamp$means)

which is much easier to remember and easier to read. This latter fact is often unappreciated, but if you get in the habit of naming elements in a list, you will be a better collaborative coder. Even if you're not into collaboration, when you return to that project after two weeks of doing something else, it is much easier to remember where you left off.

A Benefit of Keeping It Classy

Another coding practice that often goes unappreciated is the use of classes to make your life easier (and your console cleaner). At this point, if you type mysamp into your console window, R will bombard you with more than a screen's worth of output. Most of this is output that I don't want to see printed out. I might want to save it for later, but I don't want to save it.

This is where the class(result) = "samp" command comes in handy. By "classing" my object this way and writing a short a print method for this class (named print.samp()), I can cut down on the amount of output I have to see when I inspect my object. Here's code to define a print method.



Now, try typing mysamp into R. You will see much less output (only the output from print.samp), but the object mysamp still has all of the information you want to keep around. If you don't believe me, just type mysamp$means to be sure.

Just as we could "class and print," we can also "class and plot." This can be handy if you want there to be a standard meaning for the command plot(mysamp). In our working example, let's make this command mean plot three histograms, one for means, one for standard deviations and one for standard errors. Here's some code that imposes this meaning on objects with class "samp"



In case you were wondering what this plotting out looks like, here is the picture that is produced from the command plot(mysamp).

Describing Data: Frequently Used Commands

Obtaining a coherent numerical summary of data is a common task, and it is common to want to port these summary statistics into a table of results. When I am in interactive mode with my data, I use the summary() command applied to my data frame. For example, the following code loads and summarizes a data frame on Yogurt advertising and prices:

library(Ecdat) ## Econometrics Data (useful!)
data(Yogurt) ## Loads Yogurt from Ecdat
summary(Yogurt) ## Summarizes Yogurt

For each quantitative variable, the summary() command provides a five-number summary (min, max, Q1, Q3, median) plus the mean. For categorical variables, the counts of each level are provided. This provides an excellent summary measure of each variable, but you may prefer a richer set of information (especially when it comes to typing up tables).

I recently discovered a great way to obtain a richer set of information on a data frame. This method involves using the psych library, which contains functions describe() and describe.by(). Continuing with the code from above, here is the basic syntax:

library(psych)
describe(Yogurt) ## Describes in more detail the Yogurt data frame

Suppose you also want to break your summary statistics into two (or four) tables for comparison sake (perhaps to illustrate stark differences across select subsets of your data). The describe.by() command is a convenient technique to break the data down by the levels of a factor. Here's an example with on the Yogurt data.

describe.by(Yogurt, Yogurt$choice)

Finally, you may want to port your data into LaTeX format and/or select particular summary statistics from the list. I wrote a function that serves as a convenience interface to describe.by() and toLatex(). As toLatex() does not work directly on objects created using describe.by(), you might find this helpful.



If you do not like knowing about the kurtosis of your data, you could read up on the options of describe.by() to learn about how to shut it down. If you're going to port it into a LaTeX table anyway, you could also just modify the code I wrote here to eliminate the summary statistics you don't want and produce LaTeX output.

FYI: Quick R has a nice summary of some other methods for summarizing data. Of the methods at Quick R that I didn't describe, pastecs looks most like a method I would use.

Merging Data Video Tutorial

Here's a video tutorial where I walk through some code that does what the previous post describes.



The FRED data is used extensively for macroeconomics. You might these data useful for joining in graph fights in the blogosphere.

Common Data Creation Commands

Here is a video tutorial where I go through some of the most commonly used commands in creating and manipulating data. As soon as I want to do more than just running a single regression, I use these commands more than any other set of commands (in some of the other videos, you may have seen these).



Here is the code I use in the video if you would like to try out these commands for yourself.

A Population Regression

Here's a video on some of the theory behind simple linear regression.



There's no R involved with this video, but the video provides some theory behind what it is that R's lm() command estimates.

A Creative Use of R

Update (5/18/2011): Looks like Freakonomics approves as well. Let the record show that I approved first :)

I approve: "I use the open-source program R to create the patterns."



But, I'm not sure I approve of calling these distributions "evil."

In case you were wondering, here are the commands to create the density plots for these "evil" distributions.

Poisson: dpois
Weibull: dweibull
Gumbel: dgumbel
Cauchy: dcauchy
Erlang: dgamma. It's a special case of the Gamma distribution. A special case of Erlang is the Chi-squared distribution.

How did I make this plot?

To make this plot, I used R's plot(), points() and lines() commands. If you have been wanting to learn how to plot in R, watch it unfold in this video tutorial:



Also, here's the code I used:

How to Bootstrap in R: The Case of R-squared

I recorded a new video tutorial whose original intent was to demonstrate how to write a for loop. As I wanted to make the for loop count for something, I decided that my application would be to write some code that computes the bootstrap approximation to the sampling distribution. This is a common econometric application.

Here is the video where I describe the code and the method.



As with most of the videos here, I am merely describing the application. You should understand why you want to bootstrap before you consider running this code (unless you just want to see a for loop in action).

Here is the code I used:



I am not going to post the data on this one, but the code should be easy to modify to work with any data set you use.

Basic Plots in R

Here's a tutorial I recorded on producing basic plots in R.



I lost the script file I used to create the video to a horrifying black screen of death, but I used the data from the previous post (available here). Hopefully, the video is clear enough that you can follow along with your own window of R.

Testing Multiple Linear Restrictions in R

This video shows a couple of things:

1. How to test (multiple) linear restrictions in R.
2. How to install a package in R.



Here is the code I use in the video (with some extra comments).



Here is the data set that I use in the video.

How Does R Deal with Categorical Explanatory Variables?

Here's a video tutorial where I demonstrate an answer to this question.



This video is designed to instill a baseline level of practical knowledge. There is more to how R treats factors in regression models. An interested reader should Google "R contrasts" for more information.

Here is the code I used in the video:



Here is a link to the data set.

Running a Regression in R

I created another video tutorial on R. This time, I discuss R's lm() command and how to use it for a variety of standard applications.



Here is the code that goes with the video:

Enjoy!

Tutorial on Distributions in R

Here's a video tutorial I put together to go over how to generate a random sample from one of the commonly known parametric distributions in R.



Along the way, I also discuss how some of the properties of estimators are reflected in the computations I perform.

How to read and write Stata data (.dta) files into R

Here's an R tutorial where I explain how to read Stata data files into R (even if you don't own the program Stata). I also offer some other basic tips.



Of note, you can also write Stata .dta files from R (if your coauthors or journals insist on having ".dta" data).

Suppose your R data frame has the name fuzzybunny and you want to save the file to the C: drive as myfuzzydata.dta. After loading the foreign library,

library(foreign)

Here's the syntax:


write.dta(fuzzybunny, "C://myfuzzydata.dta")

Try it out and impress your friends (you know, the friends who only know how to use Stata).

Summary Statistics in R

There are plenty of useful techniques for manipulating data frame objects in R. This post summarizes some of these useful techniques, and gives code to implement them. Before running the code in this post, you should read in the data set (which is available for download here). For tips on how to read Stata ".dta" files into an R data frame, see my previous post.

## Print the entire data frame to the screen by typing its name ##

caschool.df

## Look just as the variable names ##

names(caschool.df)

## Examine a particular variable using the $ extractor ##

caschool.df$testscr

## You can also think of the data frame as a matrix ##

caschool.df[,4]       ## Extracts the 4th column
caschool.df[4,]       ## Extracts the 4th observation
caschool.df[4,4:18]

## The last command extracts the 4th through 18th columns
## of the 4th observation.
##
## Note 4:18 is R shorthand for
##  c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)

## Compute summary statistics on the entire data frame ##
## or just one variable                                ##

mean(caschool.df)          ## Returns NA for categorical
round(mean(caschool.df),2) ## Rounding returns easier
                         ## to read formatting
var(caschool.df)           ## Returns var-cov matrix
diag(var(caschool.df))     ## Returns just the variances
mean(caschool.df$testscr)  ## Returns mean of testscr
var(caschool.df$testscr)   ## Returns scalar variance

## Compute the summary statistics on columns 6 through 18 ##

mean(caschool.df[,6:18])
var(caschool.df[,6:18])    ## Returns var-cov matrix

## If you don't like typing the data frame's name every ##
## time you want to explore the variable, you can use   ##
## the attach() command.

attach(caschool.df)

## Now, the variables in the data frame can be accessed ##
## without extracting them with the $

mean(testscr)


##-------------------------------------------------- ##
## Creating/Storing new variables                    ##
## Just using arithmetic definition of new variables ##
##-------------------------------------------------- ##

math_read_avg = (read_scr+math_scr)/2

## -------------------------------------------##
## You can attach it to the data frame        ##
## In fact, for regression, you want to do so ##
## -------------------------------------------##

caschool.df=cbind(caschool.df,math_read_avg)

## -------------------------------------------------##
## ... and its name will be what you called it when ##
## you defined the variable                         ##
## -------------------------------------------------##

names(caschool.df)
caschool.df$math_read_avg

## ------------------------------------------------- ##
## Too many datasets in R's active memory can cause
## problems math_scr might be a variable name in
## multiple data frames ... especially if you work
## with a lot of data sets      
##                                                
## In case of ambiguity, R will write over the old
## variable name with the new definition. 
##                                        
## To avoid this problem, when you are done with an
## attached data set, you should use detach()
## ------------------------------------------------ ##

detach(caschool.df)

## ------------------------------------------------ ##
## An added note: I like to leave my workspace image
## clutter-free.  After working on some code, I save
## my code in a text file, but I do *not* save my
## workspace image.      
##                       
## Following this practice can help you avoid
## referencing a variable name that you created months
## ago for a separate project
## ------------------------------------------------ ##

Reading in data in R

Reading in data can be a chore in R if you haven't done it before. In this post, I describe how to read data into R's most common data object -- the data frame.

## ----------------------------------------------------##
## The most basic structure in R is a data frame. The  ##
## goal here is to read data from some external source ##
## into an R data frame                                ##
## ----------------------------------------------------##

## --------------------------------------------------- ##
## The standard space-delmited function is read.table()##
##                                                     ##
## To read in the data, you will want to save the file ##
## to a known directory. For me, this is "C:/R/"       ##
##                                                     ##
## Time saving note: these commands sometimes require  ##
## double backspaces in the file path. So, even if I   ##
## saved my file at "C:/R/data.txt" I tend to type     ##
## "C://R//data.txt"                                   ##
## --------------------------------------------------- ##

## No header is default ##
traffic.df = read.table("C://R//traffic.txt")
traffic.df   ## Display the data frame

## How to tell R to store variable names ##
traffic.df = read.table("C://R//traffic.txt",header=T)

## --------------------------------------------------- ##
## A common difficulty is that the basic command does  ##
## not read in Excel Workbook files.  For this, there  ##
## are a couple of workarounds.   My favorite is to    ##
## use read.csv().                                     ##
## --------------------------------------------------- ##

## First, save your Excel file as type "csv"
## Second, use a command like the following:

insulin.df = read.csv("C://R//insulin.csv", header = T)

##----------------------------------------------------##
## Suppose you want to read in an Stata .dta file     ##
## The command is read.dta().                         ##
##                                                    ##
## This requires loading a package.                   ##
## In this case, the "foreign" library.               ##
##                                                    ##
## Because external packages are not part of the base ##
## group of packages, you'll have to tell R to use    ##
## them with the command library()                    ##
##----------------------------------------------------##

library(foreign)

caschool.df = read.dta("C://R//caschool.dta")

------------------------------------------------------

If you came across this page looking for advice on how to read data into R, you likely have a data set that you want to use. If you want to use some of the data referenced in this post, the traffic data set is here (just copy and paste into Notepad, and save to your favorite directory):

density speed
20.4  38.8
27.4  31.5
106.2  10.6
80.4  16.1
141.3   7.7
130.9   8.3
121.7   8.5
106.5  11.1
130.5   8.6
101.1  11.1
123.9   9.8
144.2   7.8
29.5  31.8
30.8  31.6
26.5  34.0
35.7  28.9
30.0  28.8
106.2  10.5
97.0  12.3
90.1  13.2
106.7  11.4
99.3  11.2
107.2  10.3
109.1  11.4