How to Bin?

One of you asked me a question about the Fathom assignment on constructing histograms.  Let me describe the issues by considering a histogram for a new dataset.

Today I was thinking about snow and wondering if 6 inches is really a lot of snowfall.  So I thought I would look for some data showing levels of snowfall for different locations in the United States.

At the National Climatic Data Center site http://www.ncdc.noaa.gov/ussc/USSCAppController?action=map, one can download a "1.1 MB ASCII text file containing the maximum 1-day, 2-day, and 3-day snowfall for all available stations in the Lower 48 States and Alaska based on data through December 2006."  This sounded interesting so I downloaded this data, created a text datafile, and then read this into R.  (I had to do some data cleaning -- a little tedious.)

data=read.table("snowclim-fema-sep2006.txt",header=T,sep="\t")

This dataset contains the maximum 1-day snowfall, the maximum 2-day snowfall, the maximum 3-day snowfall for 9087 weather stations in 49 states in the U.S.

Let's focus on the max 1-day snowfall for all of the stations in Ohio -- remember, I wanted to figure out the extreme nature of 6 inches.

Using the subset command, I create a new data frame ohio that contains the data for the Ohio stations.

ohio=subset(data,data$State=="OH")

The variable of interest is called Max.1.Day.

Let's say that I wish to construct a histogram by hand -- that is, choose the bins by hand. 

First, I use the summary command to figure out the extremes.  So the variable falls betweeen 6.5 and 21.

  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 

   6.50   11.25   13.00   13.28   15.00   21.00 

Let's try using 7 bins (why not?).  The range of the data is 21-6.50 = 14.5.  If I divide this range by 7, I get 14,5/7 = 2.07 which is approximately 2 (a prettier number than 2.07).

So let's try using the bins

(6, 8),  (8, 10), (10, 12), (12, 14), (14, 16), (16, 18), (18, 20), (20, 22)

(Actually this is 8 bins, but that's ok -- it is close to 7 bins.)

In R, I create a vector that gives these bin endpoints -- we call these my.breaks.

my.breaks=seq(6,22,2)

I use the hist command with the breaks option.  Also I want to show "6 inches" on the graph.  I do this by adding a vertical line to the graph (abline command) and label it by the text command. 

hist(ohio$Max.1.Day,breaks=my.breaks)

abline(v=6,lwd=3,col="red")

text(8,45,"6 inches",col="red")

Is this a good looking histogram?  In other words, did I make a reasonable choice of bins?  

This is a pretty smooth and I have a clear idea of the shape of the data (a little right-skewed).  But maybe I need more bins to see finer structure in the histogram.

Let's try using bins of length 1.

my.breaks=seq(6,22,1)

hist(ohio$Max.1.Day,breaks=my.breaks)

abline(v=6,lwd=3,col="red")

text(8,30,"6 inches",col="red")

Well, I think we went too far.  The histogram is less smooth -- I see some bumps in the histogram that may be a byproduct of random noise.  I definitely prefer the first histogram to this one.

I am guessing that an "optimal" bin width is between 1 and 2.

If you use the hist command without a breaks option, it will use the "Sturges" method to determine the optimal bin width for this particular sample size.  This method is similar to the one described in the Fathom lab.

I tried it.  Suprisingly (really), it turns out that R used the same breaks as I did in the first histogram.  

Wrap-up:  How does this discussion relate to the Fathom homework?  In Fathom, it is very easy to play with the bin width and see the effect of the choice of bins on the histograms.  In the famous kid's tale "Goldilocks and the Three Bears", remember that Goldilocks tried several bowls of porridge -- one bowl was "too hot", a second bowl was "too cold", and a third bowl was "just right".  Similarly, you'll get bad histograms when you choose too few bins or too many bins, and there will be a middle number of bins that will be "just right".  By playing with different histograms, you will decide (by eye) what choice of bin width is "just right", and then you will see if your choice of best bin width corresponds to the bin choice using the optimal histogram rule.

By the way, how does the "optimal rule" 

  bin width = 2 (fourth spread) / n^(1/3)

perform here?  Here the quartiles are 11.25 and 15, the fourth spread is 15-11.25 = 3.75, so the optimal bin width is 2 (3.75)/179^(1/3) = 1,33 which is a bit smaller than the bin width of 2 used in the hist command.

Final thought -- "6 inches" is really not a lot of snow, although it seemed to cripple the city where I live.

President Ages

As you know, yesterday was a big day in the history of the US as Barack Obama became our 44th president.  In honor of that event, I collected the ages at inauguration and lifespans of all 44 presidents and read this dataset into R.

data=read.table("president.ages.txt",header=T,sep="\t")

I attach the dataframe to make the variables visible.

attach(data)

Then using the stem.leaf function in the aplstats package, I construct a stemplot of the president ages:

> stem.leaf(Age)

1 | 2: represents 12

 leaf unit: 1

            n: 44

   2     t | 23

         f | 

   6     s | 6677

   9    4. | 899

  15    5* | 011111

  17     t | 22

  (9)    f | 444445555

  18     s | 6667777

  11    5. | 8

  10    6* | 0111

   6     t | 2

   5     f | 445

         s | 

   2    6. | 89

What do I see?  Here are some things I notice:

• The ages at inauguration seem pretty symmetric shaped about the values of 54 or 55.
• President ages range from 42 to 69.  Actually, one has to be a particular age to be president, so I believe that 42 is close to the rule.
• Barrack Obama is one of the youngest presidents in history at 48

I didn't play with the options for stemplot -- I just used the default settings in the stem.leaf function.  Could I produce a better stemplot?

Let's break between the tens and ones with five leaves per stem.

> stem.leaf(Age,unit=1,m=2)

1 | 2: represents 12

 leaf unit: 1

            n: 44

    2    4* | 23

    9    4. | 6677899

   22    5* | 0111112244444

  (12)   5. | 555566677778

   10    6* | 0111244

    3    6. | 589

> 

Or I could go the other way and split between the ones and tenths with 10 leaves per stem.

> stem.leaf(Age,unit=.1,m=1)

1 | 2: represents 1.2

 leaf unit: 0.1

            n: 44

   1    42 | 0

   2    43 | 0

        44 | 

        45 | 

   4    46 | 00

   6    47 | 00

   7    48 | 0

   9    49 | 00

  10    50 | 0

  15    51 | 00000

  17    52 | 00

        53 | 

  22    54 | 00000

  (4)   55 | 0000

  18    56 | 000

  15    57 | 0000

  11    58 | 0

        59 | 

  10    60 | 0

   9    61 | 000

   6    62 | 0

        63 | 

   5    64 | 00

   3    65 | 0

HI: 68 69

> 

I think it is pretty obvious that the first stemplot is the best.   If I have too few lines, then I lose some of the structure of the distribution; with too many lines, I don't see any structure at all.

Learning R by Viewing Movies

As you are learning R, you should try out the movies that I have made that are posted at

http://bayes.bgsu.edu/eda

(I call them R Encounters.) You'll see movies that show

-- how to create a dataset and read it into R
-- how to install the LearnEDA package in R
-- illlustrating a simple data analysis on R

I'm assuming that most of you are working in Windows. Let me know if you work on a Macintosh -- things work on R a little differently on mac.

Students have struggled in the past on creating datafiles and reading them in R -- let me know if you are experiencing problems with this.

Also, let me know if you like the R movies. They are easy to make and I can make more of them if you find them useful.

Welcome to Exploratory Data Analysis

Hi EDA students.

I will be using this blog to provide advice on using R, give you new examples and help on the homework.

This week, you'll be learning R and Fathom and creating a couple of datasets that you'll be reading into R.

I have provided help on R on my website http://bayes.bgsu.edu/eda

You'll find guidance on installing R and the LearnEDA package, and several R lessons that you can work through.  Also I made some movies that illustrate some basic R stuff.

I know that learning R can be difficult, so let me know if you have specific concerns.

Why do we flog?

I just put an example of using flogs to compare scores on two years of placement scores.  But some of you may be confused why we take flogs in the first place.  Here are the main points.

1.  We want to compare two proportions, say p1 and p2.  There are two problems with a direct comparision like p1/p2.

a)  This measure depends on whether we consider the proportions p1, p2 or the proportions 1-p1, 1-p2, and that's a problem.  The choice of p1, p2 or 1-p1, 1-p2 shouldn't change our comparision.

b)  Small proportions tend to have smaller variation than large proportions (close to .5).  The ratio p1/p2 = 1.5 is more meaningful (significant) if the p's are close to zero, than if the p's are close to 0.5.

So we want to transform a proportion p so that 

-- it doesn't matter if we consider p or 1-p

-- the variance of the transformed p will be roughly the same for p near 0 and p near 0.5.

By using the flog reexpression log(p/(1-p)) we achieve these two goals.

By the way, the flog reexpression is the basis for logistic regression models.  You will likely be using logistic regression in your statistical life, but we are trying here to motivate why we use the flog reexpression.

Flogging placement data

Here's an illustration on transforming data for proportion data.

Every summer, a math placement test is given to over 3000 entering freshmen.  The score on this test is used to determine which math course they are allowed to take in the fall.

The score on the placement test is a course number that indicates that the student can take that course in the fall.  Here are the counts of the placement scores for the freshmen in the last six years:

     2003 2004 2005 2006 2007 2008

131H   33   54   51   57   44   41

131   196  364  342  361  320  251

130   192  245  236  248  211  208

128   557  707  647  700  603  618

126   428  518  489  480  442  428

122   501  612  580  565  464  498

215   661  912  888  838  792  747

112   207  216  230  208  212  180

095   418  524  545  419  407  335

090    46   78   89   58   62   46

How can we compare the scores for different years?  A first step to compute percentages of each column.

     2003 2004 2005 2006 2007 2008

131H  1.0  1.3  1.2  1.4  1.2  1.2

131   6.1  8.6  8.3  9.2  9.0  7.5

130   5.9  5.8  5.8  6.3  5.9  6.2

128  17.2 16.7 15.8 17.8 17.0 18.4

126  13.2 12.2 11.9 12.2 12.4 12.8

122  15.5 14.5 14.2 14.4 13.0 14.9

215  20.4 21.6 21.7 21.3 22.3 22.3

112   6.4  5.1  5.6  5.3  6.0  5.4

095  12.9 12.4 13.3 10.7 11.4 10.0

090   1.4  1.8  2.2  1.5  1.7  1.4

We see in 2003 that 6.1% of the students placed in MATH 131 and 20.4% placed in MATH 215.

Let's follow Tukey's strategy for comparing percentage vectors.  To make this simple to explain, let's focus on comparing the percentages of 2006 and 2007.

     2006 2007

131H  1.4  1.2

131   9.2  9.0

130   6.3  5.9

128  17.8 17.0

126  12.2 12.4

122  14.4 13.0

215  21.3 22.3

112   5.3  6.0

095  10.7 11.4

090   1.5  1.7

1.  First we cut the data by some row.  Let's try cutting the data after the second row.

     2006 2007

131H  1.4  1.2

131   9.2  9.0

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

130   6.3  5.9

128  17.8 17.0

126  12.2 12.4

122  14.4 13.0

215  21.3 22.3

112   5.3  6.0

095  10.7 11.4

090   1.5  1.7

2.  We compute a folded log for each year.  For 2006, we see that 1.4 + 9.2 = 10.6% are above the line and 100 - 10.6 = 89.4% are below the line, so the flog is

FLOG for 2006 = log(10.6/89.4) = -2.13

Likewise the flog for 2007 is given by

FLOG for 2007 = log(10.2/89.8) = -2.18

3.  To compare the years 2006 and 2007, we look at the difference in flogs:

Change in FLOG from 2006 to 2007 is -2.18 - (-2.13) = -0.05

The interpretation is that students did 0.05 worse in 2007 (on the flog scale).

What if we cut the table by a different row?  We can repeat the procedure using all possible cuts.

Here is the table of flogs:

       2006  2007

 [1,] -4.22 -4.38

 [2,] -2.13 -2.17

 [3,] -1.59 -1.65

 [4,] -0.63 -0.70

 [5,] -0.12 -0.18

 [6,]  0.46  0.35

 [7,]  1.56  1.44

 [8,]  1.98  1.88

 [9,]  4.20  4.03

To understand these values, -4.22 is the flog if we cut after the first row, -2.13 is the flog if we cut after the second row, etc.

To compare the years 2006 and 2007, we look at the difference in flogs:

     2007 FLOG - 2006 FLOG

 [1,] -0.16

 [2,] -0.04

 [3,] -0.06

 [4,] -0.07

 [5,] -0.06

 [6,] -0.11

 [7,] -0.12

 [8,] -0.10

 [9,] -0.17

What have we learned?  Note that all of the flog differences are negative and the median flog difference is -0.10.  So it is clear the 2007 students did a little worse than the 2006 students.

Binning homework

Some of you appear to be confused on this "binning" homework.  Let's outline what you are supposed to do in this assignment.

1.  First you bin your data and construct a histogram.  This is easy (I hope).  The histogram gives you counts of each bin.

By the way, a rootogram is just like a histogram, but you are graphing the root counts against the bins.  (Why do you do this?  Look at the notes.)

2.  You find a Gaussian curve that fits your data -- you have a mean and a standard deviation.

3.  Now you want to find the expected or fitted counts from the Gaussian curve.  You do this by running the function fit.gaussian -- the inputs to this function are your bins, the vector of raw data, and the mean and standard deviation of the Gaussian curve.   Suppose you save the output of this function into the variable s.  Then

s$counts are the observed bin counts

s$expected are the expected counts

4.  The question at this point is -- does the Gaussian curve provide a good fit to the histogram?  Well, maybe yes, and maybe no.

How can you tell?

You look at the residuals which are essentially the deviations of the counts from the expected counts.

I defined several residuals in the notes.

Simple rootogram residuals are 

r.simple = sqrt(d) - sqrt(e)

These are graphed by use of the rootogram function.

The double root residuals are defined by

DRR = sqrt(2+4d) - sqrt(1+4e).

How do I compute these?  Well, you already have computed the d's and the e's -- you can use R to compute the DRR's.

5.  Once I compute the residuals and graphed them, am I done?  Well, the point of graphing the residuals is to see if they show any systematic pattern -- if there is a pattern, then the Gaussian curve is not good.

6.  Well, we are almost done.  Part 1 (d) asks you to fit a G comparison curve to the root data.  (Boy, your instructor seems to like taking roots.)

All I'm asking here is to FIRST transform the data by a root, and then repeat all of the above with the root data.

7.  Lastly, in #2, I'm asking you to fit a G curve to the women heights.  If you have loaded in the studentdata dataset, then to get the female heights, you type

library(LearnEDA)

data(studentdata)

f.heights=studentdata$Height[studentdata$Gender=="female"]

Good luck!