10.1 Load today’s data

Start by loading today’s data. If you need to download the data, you can do so from these links, Save each to your data directory, and you should be set.

• Pizza_Prices.csv
  • Data on the weekly price and sales of pizza in various cities
• mlbCensus2014.csv
  • Basic data on every Major League Baseball player on a 25-man roster as of 2014/May/15th

# Set working directory
# Note: copy the path from your old script to match your directories
setwd("~/path/to/class/scripts")

# Load data
pizza <- read.csv("../data/Pizza_Prices.csv")
mlb   <- read.csv("../data/mlbCensus2014.csv")

10.2 Populations and samples

We have worked a lot to this point to describe data. The tools we have developed are vital to understanding data, but we have thus far missed an important point: Sometimes we don’t have all the data, but we still want to say something about a particular group. This “group” is called the population, and so far we have assumed that all of the data we have are the entirety of the group we are interested in.

For example, imagine that we are interested in the weekly price of pizza in Denver between 1994 and 1996. Then lucky us: we have a dataset that contains the entire population. Here is the top and bottom of the pizza data set:

head(pizza)
tail(pizza)

Note that the row numbers show us that we have 156 weeks worth of data, (this is 3 years times 52 weeks per year). So, from this data set, we can say that the average price of pizza in our population (weeks from 1994 to 1996) was exactly 2.5660256 with certainty.

However, in most cases, collecting such data is much more complicated. For starters, it is often expensive. To get the average price of pizza for a week for an entire city will involve a lot of work: collecting prices, collecting sales, compiling everything. It would not be unreasonable to only collect data in only a subset of the weeks, then figure out the mean from those weeks.

Such a subset would be called a sample, and it is the kind of data we almost always have. This generally gives a very accurate picture of the population as a whole, but only if the sample is collected carefully.

However, it should be obvious that there will be occasional mistakes and problems. Just by chance, we may have sampled in the cheapest weeks (or the most expensive).

The rest of this chapter and course is largely focused on how (and how accurately) we can use a sample to draw conclusions about a population.

10.3 Sampling example

In class, we represented the population of pizza prices as slips of paper in a cup. From these, each student randomly drew one slip of paper, representing a single week for which we have data. Then, we can calculate a mean and see how it compares to the population mean. The values drawn can be seen in this histogram:

Recall from above that the mean of the entire population was 2.5660256. The mean in this first sample was 2.5186667 which is different from the population by $0.047. Does that mean our sample was “wrong”?

Of course not – but it does mean that it wasn’t “right” either. Every time we draw a sample, the exact mean will flucuate a little bit. However, we expect that, most of the time, it will land fairly close to the population mean.

Let’s look at all 4 samples we took in class. Note that I have added a line for both the population and sample means.

10.4 Sampling in R

Now, instead of manually drawing and calculating, let’s try doing it in the computer. We will start by roughly replicating what we did by hand in class. We will draw a sample, then cacluclate the mean from it.

The function sample() draws a sample from a vector (DEN.Price from the pizza data, below) of the size (number of values) that you give in the size = argument. Here, we will use 20, which is similar to a class size.

# Draw a random sample
tempPrice <- sample(pizza$DEN.Price,
                    size = 20)

We can then look at the sample in several different ways.

# Print the values
tempPrice

##  [1] 2.52 2.56 2.56 2.54 2.56 2.72 2.50 2.68 2.85 2.62 2.68 2.31 2.73 2.46
## [15] 2.64 2.54 2.70 2.36 2.48 2.43

# Calculate the mean
mean(tempPrice)

## [1] 2.572

# Visual display of the distribution
hist(tempPrice)

We can then compare these to the population parameters (which we know because we have access to the whole population).

# Display the true population values
mean(pizza$DEN.Price)

## [1] 2.566026

hist(pizza$DEN.Price)

Now, go back and re-run the sample() line and the summary() of it’s outputs a few times. Notice that the output is different because the random sample selected a differnt set each time.

10.5 Repeated sampling using loops

Now, what if we want to see how “wrong” our mean tends to be? Instead of just relying on manually re-running the sample, let’s try making R do it for us. Here, we introduce the concept of a “loop” for the first time. It is just a formal way of telling the computer: “Do this {stuff} a bunch of times.”

This is an example of a for loop it travels through the loop once for each value you assign it. Here, we are using i for “index” and assigning it numeric values.

# Example Loop
for(i in c(1,2,3)){
  print(i)
}

## [1] 1
## [1] 2
## [1] 3

The first argument to for() is what will be set in the loop. The part after in tells it what to set i as: here we are saying to use the values 1, 2, and 3. The loop then contains everything betewen the curly brackets “{ }”. We can put anything we want in there. Here is another example.

for(i in c(1,2,3)){
  # Comments are still ignored
  # but tell us that this will print the value of "i"
  # in each loop
  print(i)
  
  # This will multiply i by 3 and print the result
  print(i * 3)
}

## [1] 1
## [1] 3
## [1] 2
## [1] 6
## [1] 3
## [1] 9

We can also save things in loops. Always remember to “initialize” or create the variable first, so that R knows where to put the result. If you ever get an error that reads something like: Error in myResult[i] <- i * 3 : object 'myResult' not found in a loop, the most likely culprit is that you forgot to initialize the variable (or mistyped it in one place). Here is an example of loop where we save something.

# Create an empty numberic value
myResult <- numeric()

# Run a loop and save a value
for(i in c(1,2,3)){
  # save the result in the "ith" position, hence index
  myResult[i] <- i * 3
}

# show what we saved
myResult

## [1] 3 6 9

Now, we can use that to save the means of many samples. Note that I am replacing the c(1,2,3) with 1:3, which is just a short hand for: “give me every digit from 1 to 3.”

# Create an empty numberic value
myMeans <- numeric()

# Run a loop of the sample and save the mean
for(i in 1:3){
  # Copy the sampling from above
  tempPrice <- sample(pizza$DEN.Price,
                      size = 20)
  
  # save the mean in the "ith" position, hence index
  myMeans[i] <- mean(tempPrice)
}

# diplay results
myMeans

## [1] 2.6045 2.5420 2.6235

Up until now, we could have done all of this by hand. The real advantage (and advantage of R over other tools) is that we can easily do this 1000 times (or any large number we want). Simply copy the above and change 3 to 1000 (or some large arbitrary number).

# Create an empty numberic value
myMeans <- numeric()

# Run a loop of the sample and save the mean
for(i in 1:1000){
  # Copy the sampling from above
  tempPrice <- sample(pizza$DEN.Price,
                      size = 20)
  
  # save the mean in the "ith" position, hence index
  myMeans[i] <- mean(tempPrice)
}

R has just calculated the means of 1000 random samples! Now, let’s take a look at those means. Note that we are no longer plotting the sampled values, but are instead plotting the means of each of the samples

# plot results
hist(myMeans)
abline( v = mean(pizza$DEN.Price),
        col = "darkblue", lwd = 3)

Note that the distribution is centered very close to the true population mean. We can even see how close it is:

# Calculate the mean mean:
mean(myMeans)

## [1] 2.564252

# Compare to actual mean of:
mean(pizza$DEN.Price)

## [1] 2.566026

The other thing to note is that the distribution is much narrower, we don’t get mean values nearly as low as the cheapest week ($2.11) or as high as the most expensive week ($2.99). That is because the mean is going to minimize those extreme swings.

Below, and in more detail next lesson, we will begin to explore the details of these distributions.

10.6 Sample size change

What happens if, instead of taking 20 values, we take a sample of 40 values?

# Create an empty numberic value
myMeansForty <- numeric()

# Run a loop of the sample and save the mean
for(i in 1:1000){
  # Copy the sampling from above
  # Change the "size" to 40, to take 40 samples, instead of 20
  tempPrice <- sample(pizza$DEN.Price,
                      size = 40)
  
  # save the mean in the "ith" position, hence index
  myMeansForty[i] <- mean(tempPrice)
}

# plot results
hist(myMeansForty)

# Numeric summary
mean(myMeansForty)

## [1] 2.565133

How does this distribution compare to when we only took a sample of 20 values? Let’s look. Note that I am adding the population mean to each plot, just so we have a reference.

# Change the par(mfrow) to make one plot above another
# Note also that I am setting the x limit to the same
#   range for each plot to capture the top and bottom
par(mfrow=c(2,1))

# plot the set from sample size == 10
hist(myMeans,      xlim = c(min(myMeans,myMeansForty),
                            max(myMeans,myMeansForty)))
# Add a vertical line for the mean
abline(v=mean(pizza$DEN.Price), col = "darkblue",lwd = 3)

# plot the set from sample size == 40
hist(myMeansForty, xlim = c(min(myMeans,myMeansForty),
                            max(myMeans,myMeansForty)))
# Add a vertical line for the mean
abline(v=mean(pizza$DEN.Price), col = "darkblue",lwd = 3)

# Don't forget to set your par back
par( mfrow = c(1,1))

From the plot, which of these seems to give more accurate estimates of the population mean? Why do you think that is?

As sample size goes up, the extreme values affect the means less often. This suggests that any single mean we take from a larger sample size is likely to be closer to the true value than a mean from a smaller sample size.

10.7 Examine accuracy

We can do even better than approximating our accuracy by adding a numeric summary of how consistent we are. A measure of accuracy for samples is called the Standard Error (SE) and is just a special case of the Standard Deviation.

Because of the 95% rule, (in roughly unimodal, symmetric data, 95% of values fall within 2 sd of the mean), we can estimate that roughly 95% of our sample means will fall within 2 SE of the true population mean. This is an important point that will be discussed in more detail in the next chapter. To calculate, we use the same function we use to calculate sd’s of population and samples.

# for our samples of size 20
sd(myMeans)

## [1] 0.03730978

# For our samples of size 40
sd(myMeansForty)

## [1] 0.02485234

Note how much smaller the SE is from the larger sample sizes. This suggests that we should have more confidence in larger samples, and supports our intuition about the differnces.

10.8 The importance of a random sample

Always be careful that you are drawing a representative sample. If your sample is not representative, your estimate of the population is not likely to be accurate. For instance, compare these approaches to sampling Mlb Weights.

# See where we should be
mean(mlb$Weight)

## [1] 209.72

# Generate several samples
mlbSamples <- numeric()
for(i in 1:1000){
  # Copy the sampling from above
  tempWeight <- sample(mlb$Weight,size = 40)
  
  # save the mean in the "ith" position, hence index
  mlbSamples[i] <- mean(tempWeight)
}

# plot results
hist(mlbSamples,
     xlab = "Mean of sample")

# Add lines
abline(v = mean(mlb$Weight), col = "darkblue", lwd = 3)
abline(v = mean(mean(mlbSamples)),
       col = "red3", lwd = 3, lty = 3)

# Add legend
legend("topright",
       legend = c("Mean of sample means","Population Mean"),
       lty = c(3,1), lwd = 3,
       col = c("red3","darkblue")
       )

However, what if, instead, we accidentally sample only Pitchers? Note here that, within sample, I am limiting to only those players where mlb$POS.Broad == "Pitcher".

# Generate several samples
mlbBiasedSamples <- numeric()
for(i in 1:1000){
  # Copy the sampling from above
  tempWeight <- sample(mlb$Weight[mlb$POS.Broad == "Pitcher"],
                       size = 40)
  # save the mean in the "ith" position, hence index
  mlbBiasedSamples[i] <- mean(tempWeight)
}

# plot results
hist(mlbBiasedSamples,
     xlab = "Mean of sample")

# Add lines to show means 
abline(v = mean(mlb$Weight), col = "darkblue", lwd = 3)
abline(v = mean(mean(mlbBiasedSamples)),
       col = "red3", lwd = 3, lty = 3)

# Add a legend
legend("topright",
       legend = c("Mean of sample means","Population Mean"),
       lty = c(3,1), lwd = 3,
       col = c("red3","darkblue")
       )

While, in R, this seems like it would be difficult to “accidentally” do, in the real world, such errors are rather common. Imagine, if you will, that we wanted to sample player weights as early as we could, so we went to Spring Training the day that it opened for each team. Do you know who the first players are to report?

It turns out that pitchers and catchers are the first ones to camp. This would make it relatively easy to mistakenly get the wrong picture of MLB player weights if we weren’t careful about our samples.