Exploratory Data Analysis

Binning Baseball Ages

2009-04-08T04:35:00.000-07:00

To motivate the issues on the next EDA topic, I collected the ages of all 651 pitchers who played major league baseball during the 2008 season.

I have a vector age that contains the ages in years for these pitchers.

A histogram is a standard way of graphing this data. Before I graph, I should select reasonable bins; I could use the default selection of bins chosed by the R hist command, but I typically like more control over my graphical displays.

Here I'm interested in the number of players who are each possible age 20, 21, 22, ... etc. So choose cutpoints 19.5, 20.5, ..., 46.5 that cover

the range of the data and so there will be no confusion about data falling on bin boundaries.

cutpoints=seq(19.5,46.5)

Now I can use the hist function using the optional breaks argument.

hist(age,breaks=cutpoints)

What do I see in this display?

The shape of the data looks a bit right-skewed. I'm a little surprised about the number of pitchers who are 40 or older.
The most popular ages are 25 and 27 among MLB pitchers.
Looking more carefully, it might seem a little odd that we have 79 pitchers of ages 25 and 27, but only 66 pitchers who are age 26.

What is causing this odd behavior in the frequencies for popular ages? We don't see this behavior for the bins with small counts.

Actually, this "odd behavior" is just an implication of the basic EDA idea that

LARGE COUNTS HAVE LARGER VARIABILITY THAN SMALL COUNTS

So we typically will see this type of behavior whenever we construct a histogram.

When we plot a histogram, it would seem desirable to remove this "variability problem" so it is easier to make comparisons. For example, when we compare the counts to expected counts assuming a Gaussian model, it will be harder to look at residuals for bins with large counts and bins with small counts since we will have this unequal variability problem.

This discussion motivates the construction of a rootogram and eventually a suspended rootogram to make comparisons with a symmetric curve.

By the way, can we fit a Gaussian curve to our histogram? On R, we

first plot the histogram using the freq=FALSE option -- the vertical scale will be DENSITY rather than COUNTS
use the curve command to add a normal curve where the mean and standard deviation are found by the mean and sd of the ages

Here are the commands and the resulting graph.

hist(age,breaks=cutpoints,freq=FALSE)

curve(dnorm(x,mean(age),sd(age)),add=TRUE,col="red",lwd=3)

It should be clear that a Gaussian curve is not a good model for baseball ages.

Looking at Hits on my Web Page

2009-03-25T08:11:00.000-07:00

The next topic is two-way tables. To illustrate the application of median polish, I took some convenient data, that is, a data that was readily available to me.

A couple of years ago, I wrote a book on Bayesian computing using R. I have a website that gives resource material for the book and I use Google Analytics to monitor hits on this particular website. Each day I observe the number and location of hits; it is interesting data partly since it seems that statisticians from many countries are interested in my book.

Anyway, here is the data -- a number in the table represents the number of hits for a particular day of the week for a particular week.

Week 1 Week 2 Week 3 Week 4
Sunday 22 12 17 15
Monday 23 15 27 17
Tuesday 17 26 21 14
Wednesday 26 13 18 18
Thursday 24 27 28 13
Friday 28 17 17 19
Saturday 14 11 13 13

I am interested in how the website hits vary across days of the week and also how the website hits vary across weeks. I can explore these patterns by means of an additive fit that I do by the median polish algorithm.

Since the data is stored as a matrix, a median polish is done by the medpolish function:

fit=medpolish(web.hits)

The variable "fit" stores the output of medpolish. Let's look at each component of medpolish.

> fit$overall
[1] 18

This tells that the average number of hits (per day) on my website was 18.

> fit$row
Sunday Monday Tuesday Wednesday Thursday Friday
-2.00 0.50 -0.50 0.00 5.50 1.75
Saturday
-5.25

These are the row effects. For Sunday, the row effect is -2 -- this means that on this day, the number of hits tends to be 2 smaller than average. Comparing Sunday and Monday, there tends to be 0.50 - (-2.00) = 2.5 more hits on Monday.

> fit$col
Week 1 Week 2 Week 3 Week 4
4.50 -2.75 1.00 -1.00

These are the column effects. It looks like my website hits across weeks where HIGH, LOW, high, low. On average, there were 4.50 - (2.75) = 7.25 more hits on Week 1 than Week 2.

The remaining component in the additive fit are the residuals. These tell us how the hit values deviate from the fitted values (from the additive model).

> fit$residuals
Week 1 Week 2 Week 3 Week 4
Sunday 1.50 -1.25 0.00 0.00
Monday 0.00 -0.75 7.50 -0.50
Tuesday -5.00 11.25 2.50 -2.50
Wednesday 3.50 -2.25 -1.00 1.00
Thursday -4.00 6.25 3.50 -9.50
Friday 3.75 0.00 -3.75 0.25
Saturday -3.25 1.00 -0.75 1.25

If the residual values are generally small (small compared to the row and column effects), then the additive model is a good description of the patterns in the data. Actually, the residuals look large to me, so I'm not sure I'd get that excited about this additive fit. Specifically, the residual for Tuesday, Week 2 is 11.25 -- for some reason, this particular day had many hits -- many more than one would expect based on its day of the week and week number.

Smoothing on R

2009-03-17T04:57:00.000-07:00

One of you asked how to produce a "3RSSH, twice" smooth on R. It seems that my R notes on the web could be clarified. Here I illustrate a simple function to do the smooth that we want.

Here is a new function that you can use called smooth.3RSSH.twice.

smooth.3RSSH.twice=function(data)

{

SMOOTH=han(smooth(attend,kind="3RSS")) # 3RSSH smooth

ROUGH=data-SMOOTH # computes the rough

SMOOTH+han(smooth(ROUGH,kind="3RSS")) # twicing operation

}

This program does three things:

1. First, one computes a 3RSSH smooth using the smooth command in R and the han function from the LearnEDA package.

2. Then one computes the rough (the residuals) from this smooth.

3. Then one smooths the rough (using the same 3RSSH smooth) and adds this "smoothed rough" to the first smooth to get a "twiced smooth".

Here is an illustration of how it works for the Braves attendance data from the notes.

(I am assuming the above function has been read into R.)

library(LearnEDA)

data(brave.at)

attach(brave.at)

plot(game,attend)

the.smooth=smooth.3RSSH.twice(attend)

lines(game,the.smooth,lwd=3,col="green")

Smoothing Free Throw Percentages

2009-03-16T07:11:00.000-07:00

I hope you had a nice spring break. My son's tennis team is flying to Florida this week -- I wish I could join him.

Since March Madness is starting, I thought it would be appropriate to talk about basketball data. There was an interesting article about free-throw shooting that recently appeared in the New York Times. See the article here.

The main message was that free-throw shooting in professional basketball has hovered about 75% for many years. Unlike other athletic performances such as running and swimming, basketball players don't seem to be getting better in shooting free throws.

Is that really true? Has free-throw shooting accuracy remained constant for all of the years of professional basketball?

To answer this question, I collected some data. For each of the seasons 1949-50 through the current season 2008-09, I collected the overall free-throw shooting percentage in the NBA.

Here are the shooting percentages graphed as a function of year.

Actually, although the overall free-throw shooting percentage is approximately 75%, there seems to some interesting patterns in the graph.

To better see the patterns, I use the command

han(smooth(FTP,kind="3RSS"))

to superimpose a 3RSSH smooth on the graph and I get the following:

What patterns do we see?

1. In the early days between 1950-1970, the shooting percentages were relatively low with a valley around 72% in the late 1960's.

2. The shooting percentages increased through the 1970's, had a small valley and hit a peak of about 76% in 1990.

3. Then the percentages decreased again and had a local minimum of 74% around 1995.

4. In recent years, the percentages are increasing. It is interesting that the current free-throw shooting percentage 77.2 is the highest in NBA history.

So in reality, the shooting percentage has not stayed flat across years. But it is surprising that NBA players haven't learned to shoot free throws better in the last 60 years.

Fitting a Line by Eye

2009-03-02T08:11:00.000-08:00

You did fine on the latest Fathom "fitting line" homework. But I sensed a little confusion and I should make a few comments about fitting a line by eye.

Let's return to that homework problem where you are plotting the OBP of the baseball players for two consecutive years.

Here is a plot that many of you produced.

Is this a good graph? Actually, NO since there is too much white space around the points.

You can improve this in Fathom by using the hand tool to fill up the space. Here is a better plot.

Second, it is not easy to fit a movable line. To make this process easier, plot the movable line and then add a residual graph.

A good line will remove any tilt pattern in the residual graph. I still have a downward tilt in this graph - this suggests I have to try a little harder.

This looks better -- I don't see much of an increasing or decreasing pattern in the residuals.

Does my best fit correspond to a least-squares or resistant fit? (By the way, the resistant line is called a median-median line in Fathom.) I show all three lines below. Least-squares is blue, median-median is purple, and my line is brown.

It looks like my line is closer to the resistant fit.

Generally, when one has outliers, I would anticipate that my line would be closer to the resistant median-median line.

Graphical User Interface for R

2009-02-27T07:47:00.000-08:00

Since some of you seem to be struggling with R's interface, where you type commands and you have to know the R functions.

There is an attractive package called R Commander which provides a menu interface for R. You work in a special window and most of the useful R functions are available as options in menus.

On the left, you see a snapshot of R Commander (on a Macintosh).

How do you get R Commander? It is easy -- you just install the package Rcmdr from CRAN. This package uses a number of other R packages. If you don't have these other packages installed yet, then R will automatically install these (it takes a little while -- be patient).

One of my former students who teaches at Youngstown State uses R Commander in his statistics classes. He says "Rcmdr is really cool, and the students eat it up."

Anyway, you might want to try this. You can run the functions in LearnEDA by first loading it and typing commands (like lval) in the top window.

Let me know if you find this helpful, since I haven't done much with it.

Plotting and Straightening on R

2009-02-26T06:15:00.001-08:00

To encourage you to use R for next week's assignment, I put together a new movie showing R in action.

As you probably know, one of the fastest growing companies in the U.S. is Google. I found an interesting graph showing Google's growth (as measured by the number of employees):

Obviously there has been a substantial increase in Google employees over this two year period. But we wish to say more. We'd like to describe the size of this growth. Also, we'd like to look at the residuals that will detect possible interesting patterns beyond the obvious growth.

Here's the plan.

1. We start with plotting the data, fitting a resistant line, and looking at the residuals.

2. By looking at the half-slope ratio and the residuals, we decide if the graph is straight.

3. If we see curvature in the graph, then we try to reexpress one or more of the variables (by a power transformation) to straighten the graph. We use half-slope ratios and residual plots to see if we are successful in straightening the graph.

4. Once the graph is straight, then we summarize the fit (interpret the slope of the resistant line) and examine the residuals.

The key function in the LearnEDA package is rline.

I illustrate using rline in my new movie

http://bayes.bgsu.edu/eda/straightening.google.swf

The dataset can be found at http://bayes.bgsu.edu/eda/data/google.txt

and my script of R commands for this analysis can be found at

http://bayes.bgsu.edu/eda/R/Ch5.google.work.R

By the way, I'm using a new version of the LearnEDA package that you can download at

http://bayes.bgsu.edu/eda/LearnEDA_1.01.zip

Hope this is helpful in your work next week.

Sexy Jobs

2009-02-25T08:36:00.000-08:00

Since many of you will be thinking about jobs soon, here is an interesting posting about desirable jobs. Many of you are in the right area!

In a recent talk by Google's chief economist Hal Varian, he says this:

"I keep saying the sexy job in the next ten years will be statisticians.
People think I’m joking, but who would’ve guessed that computer engineers
would’ve been the sexy job of the 1990s? . . .
The ability to take data—to be able to understand it, to process it,
to extract value from it, to visualize it, to communicate it—that’s going to
be a hugely important skill in the next decades, not only at the
professional level but even at the educational level for elementary school
kids, for high school kids, for college kids. Because now we really do have
essentially free and ubiquitous data. So the complimentary scarce factor is
the ability to understand that data and extract value from it.
I think statisticians are part of it, but it’s just a part. You also
want to be able to visualize the data, communicate the data, and utilize it
effectively. But I do think those skills—of being able to access,
understand, and communicate the insights you get from data analysis—are
going to be extremely important. Managers need to be able to access and
understand the data themselves. "

The full article is at

http://www.mckinseyquarterly.com/Hal_Varian_on_how_the_Web_challenges_managers_2286

Reducing Skewness

2009-02-23T07:15:00.000-08:00

Do you want to see me quickly reduce skewness?

skewness

More seriously, here are some comments after looking at your homework (the grades are posted and you did generally well on this homework).

1. Why is symmetry important?

Some of you may be wondering why it is useful to make datasets symmetric. I don't think making a dataset symmetric is as important as making datasets have equal spreads, but here are a couple of reasons why symmetry is helpful.

-- Symmetric datasets are simpler to summarize. There is an obvious "average". Also if the data is bell-shaped, one can use the empirical rule (2/3 of the data fall within one standard deviation of the mean).

-- Many statistical procedures like ANOVA assume normally distributed data. One might wish to reexpress data before applying these procedures.

2. Monotone increasing reexpressions.

-- In our definition of power transformations, recall that when p is negative, we consider

"minus data raised to a power"

We did that so that all of the transformations are monotone increasing. So when you have a right-skewed dataset and you move from p=1 to p=0 to p=-1, you should be moving toward symmetry and left-skewness.

3. Skewness in the middle and skewness in the tails.

Sometimes it is tough to make a data symmetric, since there will be different behavior in the middle half and in the tails. You can detect this by a symmetry plot. The points to the left may be close to the line and the points to the right may fall off the line. This indicates that the middle of the data is symmetric, but there is a long tail to the right (or the left).

4. When does reexpression work?

You need sufficient spread in the data as measured by HI/LO. If this ratio is not much different from 1, reexpressions won't help much.

5. Are some of you "R resistant"?

Many of you seem to prefer using Fathom or Minitab. That's okay, but the best way to get comfortable using R is to practice using it. I made the package LearnEDA to make R easier to use, but some of you aren't taking advantage of the special functions.

Last Homework and Dreamtowns

2009-02-16T04:34:00.000-08:00

We've been talking about taking power reexpressions of data to achieve particular objectives, such as equalizing spreads between batches or making a batch symmetric.

When do these reexpressions work?

1. First you need to have data with sufficient spread, that we measure by the ratio HI/LO for the power transformations to have much effect.

2. You have to be taking a "significant" power. In the last Fathom homework, most of you tried p = .82 (from a starting value of p = 1) which would typically have little effect on the data. Typically, we try reexpressions in steps of 0.5, so we move from raw to roots to logs, and so on.

This week, we are reexpressing data to achieve approximate symmetry. Here's an example.

Last summer, there was an interesting article posted on bizjournals.com that ranked 140 "dreamtowns" -- these towns offer refuge from big cities and conjested traffice. I got interested in the article since my town, Findlay, made the list.

Anyway, they collected a number of variables from each city, including the percentage of adults (25 and over) who hold college degrees.

Here's a histogram of the these percentages from the 140 towns.

This looks right-skewed with one outlier (I never knew that Bozeman, Montana had a lot of highly educated people) and this is a good candidate for reexpression.

In the notes, I give several methods (plotting the mids, using a symmetry plot, using Hinkley's method) for choosing the "right" rexpression.

I'd suggest that you should try at least 2 of these methods in your homework.

Simple Statistical Comparisons

2009-02-09T08:52:00.000-08:00

I just finished grading your "comparing batches" homework. You did fine in the mechanics (computing the summaries, constructing a spread versus level graph, and deciding on an appropriate reexpression), but it wasn't clear that you understand WHY we are doing all of this work.

Let's review the main points of this section.

1. We wish to compare two batches.

2. What does compare mean? Well, it could mean many things. Batch 1 has a larger spread, batch 2 has three more outliers than batch 2, and so on. You illustrated many type of comparisons in your homework writeups.

Here we wish to compare the general or average locations of the two batches. A simple statistical comparison is something like

"batch 2 is 10 units larger than batch 1"

Maybe we should call this a SSC (statisticians love to use acroynms.)

What this could mean is that if I added 10 units to each value in batch 1, then I would get a new dataset that resembles batch 2.

3. Is it always appropriate to make a SSC?

No.

It won't work if the two batches have unequal spreads. If they have unequal spreads, then adding a number to batch 1 will NOT resemble batch 2.

4. So if the batches have different spreads, we give up?

No.

It is possible that can can reexpress the batches to a new scale, so that the new batches have approximate new scales.

5. So the plan is to (1) try to find a suitable reexpression and (2) do a SSC on the reexpressed data.

The snowfall data example using Fathom is one example where our strategy works. But generally, you did a poor job in making a SCC on the reexpressed data.

A Statistical Salute to the Steelers

2009-02-01T20:01:00.001-08:00

I just finished watching the Super Bowl and I'm happy that the Steelers won. The winning quarterback Ben Roethlisberger is from my home town (Findlay) and it was very exciting seeing him direct his team to the winning touchdown at the end of the game.

As I was watching this game, I wondered:

"Which team has the better offense?"

To answer this question, I collected the yards gained for every single play of the two teams, the Steelers and the Cardinals, for this Superbowl. I created a datafile with two variables, Yards and Team, and entered this data into R.

Using the R function

boxplot(Yards~Team,horizontal=T,xlab="Yards Gained",ylab="Team",col="gold")

I created the following boxplot.

What do we see in this boxplot display?

Both batches of yards gained look a bit right-skewed.
There are four outliers for the Steelers and two for the Cardinals. These correspond to big plays for the teams that gained a lot of yards.

Why did I draw this boxplot display? Several comments:

We wish to make a comparison between the two batches.
What is a comparison? Well, we could say that one batch tends to be larger than the second batch. For this example, if it were true, then I would say that the Cardinals gained more yards (per play) than the Steelers.
But that isn't really saying much. When I say "make a comparison", this means that I want to say that one batch is a particular number larger or smaller than the second batch.
When can we say
"batch 1 is 10 larger than batch 2"?
As you'll read in the notes, we can only make this type of comparison "batch 1 is 10 larger than batch 2" when the two batches have similar spreads.

Returning to my football example, it doesn't appear that the two teams were very different with respect to yards gained per play. Both teams had approximately the same median.

Letter values by group

2009-01-31T06:08:00.000-08:00

In the next topic, we are talking about comparing groups. The first step is to construct some graph (such as a stemplot or boxplot or dotplot) across groups and the next step is to compute summaries (such as letter values) for each group. Here are some comments about using R to do this stuff.

Let's return to the boston marathon data where there are two variables, time and age. We are interested in comparing the completion times for different ages.

data(boston.marathon)

attach(boston.marathon)

By the way, this data frame is organized with two variables, response (here time) and

group (here age).

1. GRAPHS

There is no simple way to construct parallel stemplots in R for different groups. Parallel boxplots are easy to construct by typing, say,

boxplot(time~age,ylab="Time",xlab="Group")

Also, it is easy to construct parallel dotplots by use of the stripchart function.

stripchart(time~age,ylab="Group")

Personally, I prefer dots that are solid black and are stacked:

stripchart(time~age,method="stack",pch=19,ylab="Group")

2. SUMMARIES BY GROUP

If you apply the boxplot option with the plot = FALSE option, it will give you five number summaries (almost) for each group. But the output isn't very descriptive, so I wrote a simple "wrapper" function where the output is easier to follow.

lval.by.group=function(response,group)

{

B=boxplot(response~group, plot=FALSE, range=0)

S=as.matrix(B$stats)

dimnames(S)[[1]]=c("LO","QL","M","QH","HI")

dimnames(S)[[2]]=B$names

}

To use this in R, you just type the function into the R Console window. Or if you store this function in a file called lval.by.group.R, then you read the function into R by typing

source("lval.by.group.R")

Anyway, let me apply this function for the marathon data.

lval.by.group(time,age)

20 30 40 50 60

LO 150 194 163 222 219

QL 222 213 224 251 264

M 231 235 239 262 274

QH 240 259 262 281 279

HI 274 330 346 349 338

attr(,"class")

"integer"

I think the output is pretty clear.

3. WHAT IF YOUR DATA IS ORGANIZED AS (GROUP1, GROUP2, ...)?

Sometimes, it is convenient to read data in a matrix format, where the different groups are in different columns. An example of this is the population densities dataset.

data(pop.densities.1920.2000)

pop.densities.1920.2000[1:4,]

STATE y1920 y1960 y1980 y1990 y2000

1 AL 45.8 64.2 76.6 79.6 87.6

2 AK 0.1 0.4 0.7 1.0 1.1

3 AZ 2.9 11.5 23.9 32.3 45.2

4 AR 33.4 34.2 43.9 45.1 51.3

You see that the 1920 densities are in the 2nd column, the 1960 densities in the second column, etc.

Anyway, you want to put this in the

[Response, Group]

format. You can do this by the stack command. The input is the matrix of data (with the first column removed) and the output is what you want.

d=stack(pop.densities.1920.2000[,-1])

d[1:5,]

values ind

1 45.8 y1920

2 0.1 y1920

3 2.9 y1920

4 33.4 y1920

5 22.0 y1920

Now we can compute the letter values by group (year) by typing

lval.by.group(d$values,d$ind)

y1920 y1960 y1980 y1990 y2000

LO 0.10 0.4 0.7 1.00 1.10

QL 17.30 22.5 28.4 32.05 41.40

M 39.90 67.2 80.8 79.60 88.60

QH 62.35 114.6 157.7 181.65 202.85

HI 7292.90 12523.9 10132.3 9882.80 9378.00

Note that since I didn't attach the data frame d, I'm referring to the response variable as d$values and the grouping variable by d$ind.

How to Bin?

2009-01-28T13:51:00.000-08:00

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

2009-01-21T09:09:00.001-08:00

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

2009-01-14T06:56:00.000-08:00

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

2009-01-13T07:31:00.000-08:00

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?

2008-12-01T05:36:00.000-08:00

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

2008-12-01T05:10:00.000-08:00

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

2008-11-21T11:20:00.000-08:00

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!

Learning about Lengths of Baseball Games

2008-11-18T11:50:00.000-08:00

One complaint about baseball is that it is a relatively long game, compared to American football or basketball.

To understand more about the lengths of baseball games, I collected the lengths (in minutes) of all games played during the 2007 season. Here are my questions:

1. Are the times of baseball games normally distributed?

2. If not, how do the times differ from a normal curve?

I'll illustrate the R work. You'll need a couple of packages: LearnEDA and vcd (for the rootogram function).

1. First I'll read in the baseball game times, construct bins and the rootogram (this is essentially a histogram where you plot the roots of the counts instead of the counts).

library(LearnEDA)

data=read.table("game.time.07.txt",header=T)

attach(data)

# set up bins and construct a rootogram

bins=seq(100,300,by=20)

bin.mids=(bins[-1]+bins[-length(bins)])/2

h=hist(time,breaks=bins)

h$counts=sqrt(h$counts)

plot(h,xlab="TIME",ylab="ROOT FREQUENCY",main="")

2. Now I want to fit a normal comparison curve. I figure out the mean and standard deviation of the normal curve by using letter values.

# find mean and standard deviation of matching Gaussian curve

S=lval(time)

f=as.vector(S[2,2:3])

m=as.integer((f[1]+f[2])/2)

sd=as.integer((f[2]-f[1])/1.349)

I get that the times are approximately N(175, 24).

3. Is this normal approximation ok? The function fit.gaussian computes the expected counts for each bins and computes the residuals sqrt(observed) - sqrt(expected).

# function fit.gaussian.R fits Gaussian curve to the counts

s=fit.gaussian(time,bins,m,sd)

# output observed and expected counts

output=round(data.frame(s$counts, sqrt(s$counts), s$probs, s$expected,

sqrt(s$expected), sqrt(s$counts)-sqrt(s$expected)),2)

names(output)=c("count","root","prob","fit","root fit","residual")

output

count root prob fit root fit residual

1 1 1.00 0.01 1.00 1.00 0.00

2 8 2.83 0.06 6.08 2.47 0.36

3 22 4.69 0.19 19.17 4.38 0.31

4 27 5.20 0.32 31.34 5.60 -0.40

5 20 4.47 0.27 26.60 5.16 -0.69

6 12 3.46 0.12 11.72 3.42 0.04

7 6 2.45 0.03 2.67 1.64 0.81

8 1 1.00 0.00 0.32 0.56 0.44

9 1 1.00 0.00 0.02 0.14 0.86

10 1 1.00 0.00 0.00 0.02 0.98

# place root expected counts on top of rootogram

lines(bin.mids,sqrt(s$expected),lwd=4,col="red")

3. The rootogram function plots the residuals. Both plots below show the same info, but in different ways.

# load library vcd that contains rootogram function

library(vcd)

# again save histogram object

h=hist(time,breaks=bins,plot=F)

# illustrate "hanging" style of rootogram

r=rootogram(h$counts,s$expected,type="hanging") # this is the hanging style

# illustrate "deviation" style of rootogram

r=rootogram(h$counts,s$expected,type="deviation")

What have we learned? Although the lengths of baseball games are roughly normal in shape, we see that there are some large negative residuals on the right tail.

What does this mean? Baseball game lengths have a long right-tail which means one tends to see many long games. If you know anything about baseball, you know that baseball can go into extra-innings when there the game is tied after 9 innings, and these extra-inning games cause the long right tail.

Median Polishing Student Grades

2008-11-09T04:23:00.000-08:00

Our department recently was interested in exploring the relationship between student class attendance and performance in a number of 100-level math classes. We all know that missing math classes will have an adverse effect on grades, but we wanted to learn more about this relationship.

For many students in four 100-level math classes, we collected the student's attendance (percentage of classes attended) and their performance (a percentage where a 90% corresponds to A work, 80-90% a B, etc). We decided to categorize attendance using the categories 0-50%, 50-70%, 70-90%, over 90%) and then we found the mean performance for students in each categorization of attendance and course. We got the following two-way table.

Attendance Pct

COURSE 0-50 50-70 70-90 90-100

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

MATH 112 61.9 68.0 68.8 75.5

MATH 115 54.3 71.8 76.6 83.3

MATH 122 56.1 64.8 73.1 78.1

MATH 126 50.6 66.3 73.9 78.1

I'll demonstrate median polish with this data.

1. First, I have to get this data into a matrix form into R. I'll first use the matrix command to form the matrix (by default, one enters data column by column) and then I'll add row and column labels.

grade=matrix(c(61.9,54.3,56.1,50.6,

68.0,71.8,64.8,66.3,

68.8,76.6,73.1,73.9,

75.5,83.3,78.1,78.1),c(4,4))

dimnames(grade)=list(c("MATH 112","MATH 115","MATH 122","MATH 126"),

c("0-50","50-70","70-90","90-100"))

To check that I've read in the data correctly, I'll display "grade":

> grade

0-50 50-70 70-90 90-100

MATH 112 61.9 68.0 68.8 75.5

MATH 115 54.3 71.8 76.6 83.3

MATH 122 56.1 64.8 73.1 78.1

MATH 126 50.6 66.3 73.9 78.1

2. Now I can implement median polish to get an additive fit. I'll store the output in the variable "my.fit" and then I'll display the different components.

> my.fit=medpolish(grade)

Here's the common value.

> my.fit$overall

[1] 69.7

Here are the row effects.

> my.fit$row

MATH 112 MATH 115 MATH 122 MATH 126

-0.7375 4.5000 0.2875 -0.2875

Here are the column effects.

> my.fit$col

0-50 50-70 70-90 90-100

-16.35000 -2.75625 2.75625 8.40000

I'll interpret this additive fit.

The average performance of these students is 69.7%.
Looking at the row effects, we see that MATH 115 students get grades that are 4.5 - (0.2875) approx 4.2 points higher than MATH 122 students. MATH 122 students tend to be 0.2875 - (-0.2875) approx 0.57 points higher than MATH 126 students, and MATH 112 students are about a half percentage point lower than MATH 126 students.
Looking at the column effects, we see a clear relationship between attendance and performance. The best (90-100%) attenders do 8.4 - 2.75 = 5.65 points better (on average) than the 70-90 attendance group, do 8.4 - (-2.75) = 11.15 points better than the 50-70 attendance group, and a whopping 8.4 - (-16.35) = 24.75 points better than the "no shows" (the under 50% attending group).

3. It might be helpful to plot this fit. I wrote a function plot2way in the LearnEDA package:

> library(LearnEDA)

> plot2way(my.fit$overall+my.fit$row,my.fit$col,dimnames(grade)[[1]],

dimnames(grade)[[2]])

Note that the plot2way function has four arguments: the row part, the column part, the vector of names of the rows, and the vector of names of the columns. Here's the figure.

Actually, this plot would look better if I had figured out how to rotate the figure so that the FIT lines are horizontal. (I'll give extra credit points to anyone who can fix my function to do that.)

From this graph we see that the best performances are the MATH 115 students who attend over 90% of the classes; the worst performances are the MATH 112 students who have under 50% attendance.

4. Are we done? Not quite. We have looked at the fit, but have not looked at the residuals -- the differences between the performance and the fit. They are stored in my.fit$residual -- I will round the values so they are easier to view.

> round(my.fit$residual)

0-50 50-70 70-90 90-100

MATH 112 9 2 -3 -2

MATH 115 -4 0 0 1

MATH 122 2 -2 0 0

MATH 126 -2 0 2 0

I see one large residual that I have highlighted in red. MATH 112 students who don't come to class (under 50% attendance) seem to do 9 points better than one would expect based on the additive fit. This might deserve further study.

My Smoothing Tribute to the Phillies

2008-10-30T05:44:00.000-07:00

Last night was an exciting night. My team, the Philadelphia Phillies, won the World Series for only the second time in their history that started in 1883.

So it seems appropriate to give a statistical tribute to the Phillies.

I collected the Phillies winning percentage (percentage of games won) for each season in their history. Below I plot the winning percentage against the season year.

It is hard to see the general pattern in this graph, so it makes sense to smooth. I could use Tukey's resistant smooth described in the EDA notes, but I'll illustrate the use of an alternative smoothing method called lowess. Here is the R code to graph the scatterplot and overlay the lowess smooth.

plot(Year,Win.Pct,pch=19,col="red")

lines(lowess(Year,Win.Pct,f=.2),lwd=3,col="red")

title("PHILLIES WINNING PERCENTAGES",col="red")

abline(h=50,lwd=2)

I added the horizontal line at WIN.PCT = 50 that corresponds to an average season.

Looking at the graph, you should see why the Philadelphia fans are so excited about the Phillies winning this season.

1. The Phillies teams generally have been crummy, especially between 1920 and 1950. The team hasn't experienced much success.

2. But there have been two recent periods where the Phillies have been successful. One period was in the 1970's and the climax of this success was the Phillies first World Series win in 1980. The second period is in recent history and of course the Phillies won their second World Series in 2008.

Comments on Plotting Homework

2008-10-26T17:59:00.000-07:00

Here are some general and specific comments on your efforts on Homework 4 on plotting and straightening.

First, most of you did well on this homework. You were good in describing the fit and the pattern in the residuals. Also you made reasonable choices at reexpressing the x and/or the y variable so that the graph looked pretty straight.

But there were some things that caused you to lose points.

1. What are you looking for in a residual plot? In this assignment, the focus was looking for nonlinear patterns. For example, is it appropriate to fit a line to the (year, population) data for the England and Wales dataset? We can answer this question by looking for a nonlinear pattern in the plot of residuals against year. There is significant curvature (that is, a quadratic pattern) in the residual plot which tells you that the population growth is not linear.

By the way, some of you plotted log population against log year -- why did you take the log of year? It doesn't make any sense to me.

2. Always talk about the fit and the residuals in the context of the data. Someone plotted x against y without telling me the variables. The fun part of statistics is that you can always talk about the application.

3. Remember that funny problem where the scatterplot shows two group of points with a clear separation? This is one type of nonlinear pattern that you won't be able to straighten with a single choice of power transformation. But since the graph clearly divides into two parts, it makes sense to treat this as two independent problems and try to straighten each part.

4. Should one fit a line by least-squares or by a resistant line? In many situations, it won't make a difference -- either fit will work. But least-squares can give you relatively poor fits when there are outliers.

How can you tell if least squares isn't the best fit? Look at the residual plot. If you still see some increasing or descreasing pattern, then this tells you that least-squares hasn't explained all of the "tilt" pattern in the graph.

Straightening

2008-10-20T09:38:00.000-07:00

This week, the main topic is reexpressing either the x or y variable to straighten a nonlinear pattern that we see in a scatterplot. Although the manipulation may be straightforward, it is possible to miss the main message in this material. Here are some questions that may help in your explanations in the homework.

Question 1: What is a simple description of the pattern in a scatterplot?

Question 2: Why do we prefer to fit lines instead of more complicated curves like quadratic or cubic?

Question 3: Is it possible to straighten all nonlinear patterns by power transformations?

Question 4: Can you think of a situation or an example where it is not possible to straighten by power transformations?

Don't forget to look at your data first -- you may see rightaway that it is not possible to straighten the graph.

Example of Resistant Fitting

2008-10-15T11:06:00.000-07:00

In baseball, the objective is to win games and a team wins a game by scoring more runs than its opponent. An interesting question is "how important is a single run" towards the goal of winning a game? Suppose one collects the runs scored, the runs allowed, the number of wins and the number of losses for a group of teams. Bill James (a famous guy who works on baseball data) discovered the empirical relationship

log(wins/losses) = 2 log(runs scored/runs allowed)

He called this the Pythagorean Relationship.

Let's try to demonstrate this relationship by use of a resistant fit.

1. First, I collected data for all baseball teams in the 2008 season. The dataset teams2008. txt contains for each of the 30 teams ...

Team -- the name of the team
Wins -- the number of wins
Losses -- the number of losses
Runs.Scored -- the total number of runs scored
Runs.Allowed -- the total number of runs allowed

2. I read this dataset into R and compute the variables log.RR and log.WL.

data=read.table("http://bayes.bgsu.edu/eda/data/teams2008.txt",header=T)
attach(data)

log.RR=log(Runs.Scored/Runs.Allowed)
log.WL=log(Wins/Losses)

3. I graph log.RR against log.WL and add team labels to the graph. As we hoped, the relationship looks pretty linear.

plot(log.RR,log.WL,pch=19)
text(log.RR,log.WL,Team,pos=2)

4. I next fit a resistant line using the rline function in the LearnEDA package. I add the
fitted line to the graph.

the.fit=rline(log.RR,log.WL,iter=4)
curve(the.fit$a+the.fit$b*(x-the.fit$xC),add=TRUE)

5. If Bill James' relationship holds, the slope of the resistant line should be close to 2.

the.fit
$a
[1] 0.01006079
$b
[1] 1.801718

It approximately holds since the slope of 1.8 is close to 2.

6. To see if this is a reasonable fit, we compute the fit and the residuals.

FIT=the.fit$a+the.fit$b*(log.RR-the.fit$xC)
RESIDUAL=log.WL-FIT

and then plot the residuals, adding the team labels.

plot(log.RR,RESIDUAL,pch=19)
abline(h=0)
text(log.RR,RESIDUAL,Team,pos=2)

7. What are we looking for in the residual plot? First, we look for general patterns that we didn't see earlier in the first plot. I don't see any trend, so it appears that we removed the tilt by fitting a line.

Also we are looking for unusually small or large residual. Here a "lucky team" corresponds to a team who seemed to win more games than one would expect based on their wins and losses.

Which team was unusually lucky in the 2008 season? A hint: they were a "heavenly" team from the west coast.

Reexpressing for Symmetry

2008-10-06T05:31:00.000-07:00

Your instructor is currently in "Phillies Heaven". His baseball team rarely has a chance to win the World Series (they have won one World Series in over 120 seasons of competing) and they currently in the National League Championship against the Dodgers!

Oh right -- I'm supposed to talk about the EDA class.

I finished the grading on your "symmetry" homework and Fathom assignment. You generally did fine, but I'll explain some issues that may have caused you to lose points.

What was I looking for?

In the notes, we talked about several methods for assessing symmetry of a batch, including looking a the sequence of midsummaries, using a symmetry plot, and Hinkley's quick method. For each dataset you consider, here's a outline of what you should do:

1. Demonstrate (using some method) that your raw data looks nonsymmetric.

2. Experiment with power transformations with different choices of p to try to make the reexpressed data approximately symmetric. Use one of our methods to see if the "p-reexpression" works.

3. Convince me that you have found a reasonable transformation by graphing the reexpressed data (say by a histogram or a stemplot).

It is best if you explain (with words) your process of finding the best reexpression. I'm much more interested in your thought process than your computer output.

Here are a few other pitfalls.

1. Some of you were confused when you considered reexpressions with negative values of p.

p = 1 -- data is really right-skewed (big positive value of Hinkley's d)

p = 0 -- data is slightly right-skewed (smaller positive value of d)

p = -1 -- data is right skewed (positive value of d)

What is going on? The problem is that you were defining a reexpression like p = -1/2 as

data^(-1/2)

when you should have used

-data^(-1/2)

By adding the negative sign, all of your transformations are increasing functions, and then

you can make better sense of the reexpressed graphs and the methods.

2. Reexpression only works when there is sufficient spread in your data. Suppose you have data that ranges from 50 to 60 -- here the range is only 60/50 = 1.2. Reexpression using any value of p won't work -- that is, it won't change the shape of the data.

3. One of you used a normal probability plot to determine if the graph was symmetric. What's wrong with this? Well, it is not one of the methods we discussed in the notes. Second, checking for normality is different (but related) than checking for symmetry. Data can be symmetric but not normally distributed.

Did reexpression work?

2008-09-29T12:32:00.000-07:00

I finished grading your Fathom spread vs level plot assignment. I've posted all of the grades on Blackboard.

When I grade your homework, I'm more interested in your explanation and comments rather than the mechanics. This homework is a good illustration of this.

You started by constructing a spread vs level plot of the weights by the supplements. Here's the graph produced by the spread.level.plot function in the LearnEDA package.

Here are the main questions:

1. Is there a dependence between spread and level?

Yes, but the pattern in the graph is a bit confused with the outlying point in the lower-right section of the plot. If we removed this point, there would appear to be a stronger relationship.

2. Can we improve by a suitable reexpression?

There is a line of slope 0.18 drawn on the graph. This would suggest the use of a power transformation with power p = 1 - 0.18 = 0.82.

3. Is this a reasonable strategy?

If we transform by a 0.82 power, this won't really change things. It is almost equivalant to taking a 1 power which is no change.

But looking more carefully at the graph, we would fit a different line if we ignored that one outlier. Then one would get a line with a smaller slope, like 0.50 and this would suggest the use of a root transformation. This would really be nontrivial and would help the general dependence between spread and level.

If you made some comments that were similar in spirt to the ones I've made above, then you got full credit. You could have lost points if you went through the mechanics without commenting on what you actually did.

Missing the Forest for the Trees

2008-09-23T12:02:00.000-07:00

I just completed grading your "comparison" homework and have posted the grades.

There is a figure of speech called "missing the forest for the trees." This means that one can get caught up in the details of a problem without really understanding the real issues that are involved.

An example of "missing the forest for the trees" is your homework.

Why do we reexpress data?

We reexpress to equalize the spreads between batches.

But why do we care about equalizing spreads?

We care about equalizing spreads since we wish to make a simple comparison between batches.

What is a simple comparison?

A simple comparison is saying that one batch is, say, 5 units large than another batch.

If you don't conclude your analysis with a simple comparison, then all of your work (finding 5-number summaries, constructing a spread vs level plot, reexpressing, etc) is for NOTHING.

In other words, MAKING A USEFUL INTERPRETATION (in this case, a useful comparison) is EVERYTHING.

A few of you were successful in making simple comparisons in your homework and I congratulate you if you got a 30/30 on the comparing batches homework.

Remember, don't forget to look for the forest.

A Exploratory Data Analysis Story

2008-09-23T07:08:00.000-07:00

When I was grading your homework this week, I thought of this story.

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

One day, a boy was interested in taking a course in exploratory data analysis, but he didn't have the money to pay for it. He decided on asking his grandmother for the money for the course. She decided to help, but she said "I hope this is a worthwhile class for you."

Anyway, the boy visited his grandmother recently and told her that the course was going well. The grandmother asked what he was learning and the boy responded:

Last week, we learned how to compare groups of data. We compared the yearly snowfall of Buffalo with Cairo. It was hard to compare the two groups since there was a dependence between level and spread that I learned by constructing a spread versus level graph. But this graph suggested the use of a p = 0.5 power transformation and when I did a spread versus level graph of the transformed data, the dependence between spread and level was reduced.

The grandmother, listening intently, responded "So what did you conclude from your data analysis?"

The boy said proudly "There is more snow in Buffalo than Cairo."

The grandmother then with a heavy sigh said "Can we get our money back?"

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

What is the message in this story? (It relates to the work that you did on your homework.)

Boxplots in R

2008-09-17T07:58:00.000-07:00

Here is a simple example of constructing boxplots and summary stats in R.

I'm interested in comparing the team statistics for baseball teams this year -- I've heard
that American League teams score more runs. Is that true?

Using data from baseball-reference.com, I created a dataset 2008teamstats.txt that contains current statistics for all 30 baseball teams.

Here's my R script. I'll paste in a horizontal-style boxplot display at the end.

> b.data=read.table("http://bayes.bgsu.edu/eda/data/2008teamstats.txt",header=T)
> b.data[1:5,1:5]
Tm League R.G R G
1 TEX American 5.49 835 152
2 BOS American 5.29 799 151
3 MIN American 5.16 779 151
4 DET American 5.06 759 150
5 BAL American 5.03 750 149
> # I am interested in comparing the runs scored per game
> # (variable R.G) for the American and National league teams
>
> attach(b.data)
>
> # Here are the boxplots:
>
> boxplot(R.G~League)
>
> # boxplot has many options -- if you prefer horizontal style ...
>
> boxplot(R.G~League, horizontal=TRUE)
>
> # To get summary stats for each group, just assign boxplot()
> # to a variable, and then display the variable.
>
> b=boxplot(R.G~League, horizontal=TRUE)
>
> b
$stats
[,1] [,2]
[1,] 3.99 3.910
[2,] 4.43 4.295
[3,] 4.85 4.575
[4,] 5.06 4.695
[5,] 5.49 4.910

$n
[1] 14 16

$conf
[,1] [,2]
[1,] 4.583968 4.417
[2,] 5.116032 4.733

$out
[1] 5.34

$group
[1] 2

$names
[1] "American" "National"

Working with subgroups in R

2008-09-16T11:05:00.000-07:00

Since we are comparing groups in EDA, I thought I would give some guidance on how to subset data in R.

Suppose we want to construct stemplots for the areas of the islands in each continent in the homework. Here is some R work for constructing a stemplot of the island areas in the Arctic Ocean. The key command is "subset".

By the way, I don't think there is a simple way of constructing parallel stemplots in R.

> data(island.areas)

> names(island.areas)

[1] "Ocean" "Name" "Area"

> attach(island.areas)

> arctic.areas=subset(Area,Ocean=="Arctic")

> arctic.areas

[1] 16671 195928 27038 6194 21331 75767 16274 12872 9570 15913 83896

[12] 8000 35000 2800 23940

> library(aplpack)

> stem.leaf(arctic.areas)

1 | 2: represents 12000

leaf unit: 1000

n: 15

1 0* | 2

4 0. | 689

5 1* | 2

(3) 1. | 566

7 2* | 13

5 2. | 7

3* |

4 3. | 5

HI: 75767 83896 195928

Bins in a histogram and looking ahead

2008-09-14T11:07:00.000-07:00

Hi EDA folks:

I finished grading your Fathom assignment on the number of bins. Generally, you all did well on this, but there are a couple of things I should mention.

1. The moral of this assignment is that as you have more data (bigger n), you should use a small bin width and have more bins. It seemed that your best histograms by eye were similar to the ones chosen by the "optimal rule" formula.

2. I think the rule wasn't that effective for constructing a histogram of the old faithful data. By using a small number of bins, you didn't see any structure in each of the two humps.

3. If you lost points, it probably was due to some confusion on your calculations or maybe not the best answer to a question -- like the one about the histogram of the old faithful data. If you don't know why you lost points, just email me .

Looking ahead, the next assignment is on EFFECTIVE COMPARISON. You'll learn a specific method for equalizing spreads between batches. Although you might understand the method (spread vs. level plot, reexpressing, etc), it is important not to lose sight of what we are trying to accomplish. We want to make a reasonable comparison between groups.

So, when you do your homework this week, don't forget to think about the BIG PICTURE. Conclude your work by making a comparison.

Last, we'll be using some new R commands. Don't forget to look at the "Chapter 3 work" file that illustrates the use of these new commands.

EDA Grading

2008-09-07T08:20:00.000-07:00

Most of you are doing great on the homework so far. But I thought I should explain how I great and why you may be losing points.

Generally I am more interested in your explanations and how you are answering the main questions of interest. For example, in the graphs and summaries homework, I am not interested as much in your R work and your computation. Most of you are doing ok in getting R to produce stemplots and compute letter values. But the BIG questions are ...

-- what is the best choice of stemplot?

-- what have we learned about the data in terms of shape, average, and spread?

-- are there observations that deviate from the rest and why are these observations unusual?

You should be addressing these BIG questions in the first R homework.

In the Fathom activity, we were looking at the number of outliers one would expect for samples from different population distributions.

For normal data, we don't see many outliers. But if the data comes from a flat-tailed distribution (like the t distribution), outliers are more common. If this general conclusion wasn't obvious from your work, then you may have lost points.

Here is a final quibble (small point). Most of the stemplots you showed me were hard to read and certainly you wouldn't want to use them for any presentation.

Which stemplot do you prefer?

Stemplot A:

1 | 2: represents 1.2

leaf unit: 0.1

n: 50

2 -1. | 55

6 -1* | 0233

11 -0. | 67899

22 -0* | 01111113344

(9) 0* | 112223334

19 0. | 789999

13 1* | 0011233

6 1. | 68

4 2* | 023

1 2. | 9

Stemplot B:

1 | 2: represents 1.2

leaf unit: 0.1

n: 50

2 -1. | 55

6 -1* | 0233

11 -0. | 67899

22 -0* | 01111113344

(9) 0* | 112223334

19 0. | 789999

13 1* | 0011233

6 1. | 68

4 2* | 023

1 2. | 9

The message here is that you should use a monoproportional font where each character takes the same space like Courier.

Using the LearnEDA package

2008-09-06T07:53:00.001-07:00

I'm starting to grade your first R homework. I wrote the LearnEDA package to make it easier for you to read in datasets and do some basic calculations.

If you look at the R folder in the Course Documents section of Blackboard, you'll see the appropriate R commands for each topic.

FOR EACH HOMEWORK, MAKE SURE YOU LOOK AT THE R FOLDER SO YOU KNOW

THE COMMANDS YOU NEED TO USE.

In the first homework, you were supposed to read in the baseball attendance data and compute some letter values.

Here's how you do this in R using the LearnEDA package. (I'm assuming you have already installed this package.)

This loads the package.

> library(LearnEDA)

Read in the dataset:

> data(baseball.attendance)

Attach the data to make the variable names available:

> attach(baseball.attendance)

Compute letter values:

> lval(Home.Attendance)

depth lo hi mids spreads

1 15.5 32783.5 32783.5 32783.50 0.0

2 8.0 23704.0 36164.0 29934.00 12460.0

3 4.5 21614.5 40166.0 30890.25 18551.5

4 2.5 16574.0 41010.0 28792.00 24436.0

5 1.0 8651.0 42067.0 25359.00 33416.0

Some common problems in R and Fathom

2008-09-04T12:09:00.000-07:00

Here are some common questions I've heard recently about R and Fathom.

1. Some of you are having problems reading in datafiles which is a big concern. There are two ways you can mess up.

(a) First, it is important that R can find your files. Put all of your R work in a particular folder, say EDA, and then by choosing menu item File -> Change dir ..., you select the file EDA. To check if the working directory really has changed, type

dir()

and you should see your data files.

(b) A general form to read in a text datafile is

data=read.table(file.name, header=T, sep="\t")

where file.name is in double-quotes. The header option says that the first line in the file contains the variable names and the sep option says that columns are separated by the tab character.

2. How do you plot curves on Fathom?

Suppose you have created a scatterplot and wish to add a curve. You select the graph and choose the menu item Graph -> Plot Function. Then you just type the function (using the variable name on the x axis) in the box.

Example of data analysis on R

2008-08-26T17:35:00.000-07:00

Here is a typical data analysis on R. I'm a baseball fan and I'm interested in the wins and losses for the current Major League teams and how these wins and losses are related to the runs scored and runs allowed.

In Excel, I create a dataset that contains the wins, losses, runs scored, and runs allowed for all 30 teams. Here is the first 4 rows of the dataset.

Team	League	W	L	RS	RA
Tampa Bay	American	79	50	597	515
Boston	American	75	55	670	559
NY Yankees	American	70	60	632	585
Toronto	American	67	63	574	510

I save this dataset as "baseball2008.txt" -- text, tab-delimited format. I save this file in a folder called "eda" and I make this the current R working directory so R will find this file.

Here is my analysis:

# read the datafile into R

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

# attach the data to make the variable names available in R

attach(data)

# compute the winning proportion for all teams

win.prop=W/(W+L)

# what was the largest winning proportion?

max(win.prop)

[1] 0.6183206

# which team had the largest winning proportion?

Team[win.prop==max(win.prop)]

[1] Chicago Cubs

# for each team, compute the number of runs scored per game

runs.game=RS/(W+L)

# construct a stemplot of the runs scored per game
# (I'm assuming you have the aplpack package installed)

library(aplpack)
stem.leaf(runs.game)

1 | 2: represents 1.2
leaf unit: 0.1
n: 30
1 s | 7
4 3. | 889
7 4* | 011
8 t | 2
12 f | 4555
(6) s | 666667
12 4. | 88899
7 5* | 00011
t |
2 f | 45

# (sorry the display got messed up when I put it on this blog)
# you get a similar display by using a histogram

hist(runs.game)

# there are two teams who stand out with respect to runs scored

Team[runs.game>5.3]

[1] Texas Chicago Cubs

# what is the average number of runs per game scored by a team?

mean(runs.game)

[1] 4.621464

# do American League teams score more than National League teams?

mean(runs.game[League=="American"])

[1] 4.756021

mean(runs.game[League=="National"])

[1] 4.503728

# there is an interesting relationship between a team's winning
# proportion and the number of runs scored and runs allowed
#
# the relationship is log(W/L) = k log(RS/RA)
# where k is typically near 2

# we'll compute log(W/L), log(RS/RA)
# and fit a least-squares line -- we'll see if the slope is close to 2

log.WL=log(W/L)
log.RR=log(RS/RA)

lm(log.WL ~ log.RR)

Coefficients:
(Intercept) log.RR
-0.0003097 1.8377065

# here we see the slope is 1.84 which relatively close to 2

accessing Fathom

2008-08-26T13:53:00.000-07:00

As a student that is a couple of hours removed from BGSU is there anyway I can install Fathom without coming to the campus?

Learning R

2008-08-25T08:46:00.000-07:00

R is a great program but it has a relatively steep learning curve. To help you get started, I have three help sessions on R planned this week -- you need only attend one of the sessions.

I have four documents R_INTRO_PART_I, R_INTRO_PART_II, R_INTRO_PART_III, and R_INTRO_PART_IV in the Course Documents section that describe different aspects of R.

1. Manipulating vectors. A basic object in R is a vector. R_INTRO_PART_I discusses how to create and work on vectors.

2. Input and output. One typically wants to read datafiles into R -- read.table is useful for doing this. Data is typically stored in a R object called a data frame. Also you'll want to save R output including graphs and paste this material into a Word document.

3. Matrices. You should be comfortable working and manipulating matrices in R.

4. Plotting. You should be familiar with basic plotting commands and understand how one can add things (like labels and titles) to to graphs.

I hope you have R installed on your laptop. You may find it helpful to bring your laptop to the help session.

Installing R Packages

2008-08-22T08:23:00.000-07:00

Now that you have successfully installed R, you have access to many functions in the R "base" package. But we'll want to add packages that will give us additional functions helpful for the class. I'll illustrate adding the package "aplpack" that we'll need to get the "stem.leaf" function. Also, I'll show you how to add the "LearnEDA" package that I wrote for the class.

Installing the package aplpack from CRAN.

1. In R, choose the menu item Packages -> Install Packages.

2. You'll be asked to choose a CRAN mirror site -- choose one in the United States.

3. Then you'll see a list of all available Packages -- choose aplpack.

4. At this point, the package will be downloaded and installed -- if you see the message

package 'aplpack' successfully unpacked and MD5 sums checked

you're in good shape.

5. The package aplpack is installed but not yet loaded into R. To load the package, you type

library(aplpack)

6. To check to see if you have the aplpack commands, try constructing a stemplot on a random sample of values from a standard normal distribution:

stem.leaf(rnorm(50))

If you get some display, you're set.

Installing the class package LearnEDA.

This is a package that I wrote that has some special functions and all the datasets we'll use. This is not available on CRAN yet, since it is still in development. But it is available from my website.

1. Go to the class web folder http://bayes.bgsu.edu/EDA/ and find the zip file

LearnEDA_1.0.zip

Download this file to the Windows desktop (or some other convenient place).

2. To install this package, select Packages -> Install Package(s) from local zip files

3. Select the LearnEDA_1.0.zip file. In a short time, R will install this package.

4. To see if you have successfully loaded this package, load in the package and check to see if you can read in one of the class datasets.

library(LearnEDA)
data(football)
football

If you see a lot of football scores, then you have succeeded in loading this package.

Installing R

2008-08-22T08:09:00.000-07:00

Welcome to the EDA blog. I thought this would be a useful way of communicating and giving you advice that will help you succeed in this class.

My advice will apply to a Windows user which I think will apply to most to you. R also works fine on a Macintosh, but the precise instructions are a little different.

1. Go to The Comprehensive R Achieve Network at http://cran.r-project.org/

2. First choose a Mirror (click on Mirrors on the left) site somewhere in the United States -- this will shorten your download time.

3. Click on Windows in the Download and Install R section. Click on base and then click on the R-2.7.1-win32.exe link.

4. When this is downloaded on your machine, click on the R-2.7.1-win32.exe file. At this point, you'll just follow the instructions and the most recent version of R (2.7.1) will be installed.

5. To launch R, just double-click on the R icon on your desktop and you'll shortly see the following screen.