30.3.1 Entering data

Count data like these are usually entered in the form of a matrix, though we could do similar things with data.frames. Here, we use rbind() which binds its arguments into rows (row bind). By naming the two values, we automatically create the row names. In addition, we probably want to set column names so that we can keep track of what colors we are talking about. For that, we use the function colnames(), which (as its name suggests) gives or sets the column names.

# enter data as matrix
mmCounts <- rbind(
  Mini    = c(34, 31, 36, 41, 30, 31),
  Regular = c(35, 19, 39, 37, 20, 26)
  )
# Add descriptive column names
colnames(mmCounts) <-
  c("Blue", "Brown", "Green",
    "Orange", "Red", "Yellow")

# Look at the data
mmCounts

If we wanted to enter the data as column names instead, we could use cbind() (for column bind), and set the row names separately with row.names(). You will notice that this looks a lot like the count data we have gotten from the function table() in the past. Because, like from table() all of the values are numbers (with the descriptors in the column and row names), we can do lot’s of things working with the data directly.

30.4.0.1 Play with colors

Those colors aren’t bad, but are a little garrish. I prefer to pick better colors, and this is a nice place to remind you how to do that. For more information, refer to the chapter on functions, grep, and colors. This is just to play with some colors and would not be necesary. First, we need to add the function we used for exploring colors.

# To look at colors, from previous
showpanel <- function(Colors){
  image(matrix(1:length(Colors), ncol=1),
        col=Colors, xaxt="n", yaxt="n" )
}

I like the blue and orange, and know that I like “green3”" and “red3”“, however, feel free to look for other colors to use there. I will stick to the R builtin colors for now, but will also show you how to use hex colors you find elsewhere below. I am going to focus on picking a better brown and yellow. Staring with brown, I want to look at all of the options.

# Plot all of the browns
showpanel(grep("brown",colors(), value = TRUE))

# Show the names
grep("brown",colors(), value = TRUE)

##  [1] "brown"       "brown1"      "brown2"      "brown3"      "brown4"     
##  [6] "rosybrown"   "rosybrown1"  "rosybrown2"  "rosybrown3"  "rosybrown4" 
## [11] "saddlebrown" "sandybrown"

Of those options, I like “saddlebrown” though you may prefer to dig some more. For yellow, we then do the same thing.

# Look for better yellow
showpanel(grep("yellow",colors(), value = TRUE))

grep("yellow",colors(), value = TRUE)

##  [1] "greenyellow"          "lightgoldenrodyellow" "lightyellow"         
##  [4] "lightyellow1"         "lightyellow2"         "lightyellow3"        
##  [7] "lightyellow4"         "yellow"               "yellow1"             
## [10] "yellow2"              "yellow3"              "yellow4"             
## [13] "yellowgreen"

I like “yellow2” best on my screen, though all are a bit ugly. Finally, we set the values of colors to use. Make sure that they match the order of the color names in mmCounts.

# Set some new colors
# based on ones I like and what I found above
myColors <- c("blue", "saddlebrown","green3",
              "orange","red3","yellow2")

Then, we can plot with these colors:

# Plot the data
barplot(mmCounts["Mini",],
        col = myColors,
        main = "Distribution of Mini M&M colors")

Alternatively, if you found colors that you liked online or elsewhere, you can set them directly from the hex values (or from rgb values). For example, this link pulled up some suitable hex colors. I have them named below for clarity, though what matters is that they are in the same order as our columns.

mmColors <- c(
  Blue   = "#0168A3",
  Brown  = "#492004",
  Green  = "#2FB44A",
  Orange = "#FF6E17",
  Red    = "#C0012E",
  Yellow = "#FFF200"
  )

And they (or colors of your choosing) can then be used in the plot instead:

# Plot the data
barplot(mmCounts["Mini",],
        col = mmColors,
        main = "Distribution of Mini M&M colors")

and for the Regular M&M’s:

# Plot the data
barplot(mmCounts["Regular",],
        col = mmColors,
        main = "Distribution of Mini M&M colors")

Those two distributions don’t look quite the same. In the next section, we will look at how to compare them directly.

30.6.1 Expected counts

That plot is helpful, but what if we want to know if the two samples have the same distribution? In statistics, we say that we are looking for an association between two variables. Here, we ask if the distribution of colors is associated with the type of M&M’s. As always, we need to start with our hypotheses, and we want our null to be the “status quo” and as boring as possible.

H₀: Color is not associated with M&M type

H_a: Color is associated with M&M type

But, what the heck do we mean by “associated”? Essentially, we want to see if the distribution is different in each category, and we call this an association because it means that what color M&M’s you get is influenced by (associated with) the type of M&M’s you get.

For this, we will again calculate a Chi-square statistic. However, it is not immediately clear what our “Expected” counts should be. Let’s play with this a bit. First, because sample size matters, we need to work with our table with the counts instead of the proportions (the mmCounts matrix).

Next, we want to see what the distribution would look like if we thought that everything came from the same pool. This could mean either the distribution of colors if all M&M’s are the same or the distribution of M&M types if we ignore colors. This means that we need to add the marginal totals with addmargins().

# Display with marginal totals
addmargins(mmCounts)

This tells us how many we would get if pooled, but we can go a step further. Let’s see what proportions (of the total, not the rows) are in each “cell” (position) of the table. We will again use addmargins() because one of the things we are interested in is how what proportion of the total are each color and each type. (Round is used to make the printing a bit prettier.)

# Display proportions
round(addmargins( mmCounts/ sum(mmCounts)),3)

So, this tells us that 53.6% of the M&M’s are Mini M&M’s and that 18.2% are Blue. Now, if there is no relationship between color and type, that means we should expect that 18.2% of the Mini M&M’s are Blue (and that 18.2% percent of Regular M&M’s are blue too). For this, our expected counts could be expressed as the proportion of all M&M’s that are blue times the number of Mini M&M’s (and times the number of Regular M&M’s for that expected count). In mathematical notation, that is:

\[ \text{Expected} = \frac{\text{Column Total}}{\text{Overall Total}}*\text{Row Total} \]

We could also say that we expect 53.6% of the Blue M&M’s to be Mini M&M’s (and the same for each color). In mathematical notation, that is:

\[ \text{Expected} = \frac{\text{Row Total}}{\text{Overall Total}}*\text{Column Total} \]

These are equivalent (check the algebra), and can be written as:

\[ \text{Expected} = \frac{\text{Row Total}*\text{Column Total}}{\text{Overall Total}} \]

By making these calculations for each cell, we can figure out our Expected counts for each Color-Type cell. So, we could go through, do each of these calculations by hand[*] To do this, we would have lots of options in R, including loops, apply() and (most likely) crossprod(). , then calclate the chi-sqaure statistic from these. However, as usual we have a shortcut, and in this case, it is worth more to us to jump straight to it.

30.6.2 Running chi-square test for association

In R, that shortcut is the same as the shortcut for the chi-square test for goodness of fit: the chisq.test() function. Here, instead of a vector and (optional) expected proportion, we pass in our table and let R calculate all of those expected values.

# Run chisq.test for association
chisq.test(mmCounts)

## 
##  Pearson's Chi-squared test
## 
## data:  mmCounts
## X-squared = 3.7538, df = 5, p-value = 0.5854

Here, the p-value tells us how likely we are to get data that are this different from our expected values.

Incidentally, assuming that the p-value doesn’t allow us to reject, now is a good point to remind you again that not rejecting the null does not mean that the null is true, only that we don’t have enough evidence to reject it. In this case, the Mars Company has said that thier target distributions are different for the two types of M&M’s. The differences are modest, but real. They are small enough that we are unlikely to detect them in a sample this small, but that doesn’t mean they don’t exist (though it may mean that they don’t matter).

30.6.3 What does the difference mean?

This tests tells us if the variables are associated, but it doesn’t tell us where those differences are. For that, we need to compare the observed and expected counts. Does that mean we need to calculate the expected after all? Didn’t R have to do that? Can’t we use that somehow?

In short: yes, we can use what R calculated, we just need to know how to ask for it. Recall that, in R, everything can be saved as a variable, that includes the output of a test. Let’s try it here:

# Save the output of the test
myChiSq <- chisq.test(mmCounts)

# Display it
myChiSq

## 
##  Pearson's Chi-squared test
## 
## data:  mmCounts
## X-squared = 3.7538, df = 5, p-value = 0.5854

Note that nothing (except potential warnings) is displayed when you first run it. Instead, it saves it so you can look at it later. But, it does even more than that. We have used the function names() before to look at the column headers, but the function is even more powerful than that. It will tell us the “names” of any variable that has them. For example, we used it above to access the names of the colors in the count vectors for plotting. Here, we can use it to see what is in the test.

# See what the test has
names(myChiSq)

## [1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
## [7] "expected"  "residuals" "stdres"

This has a lot of useful information that we may want for other things. Each can be accessed using the $ operator, just like for columns of a data.frame. We will be using these values for several other tests in the remainder of this semester. Let’s look at a few. As with other tests (including t.test() and prop.test()), chisq.test() stores the test statistic, in this case the chi-squared vale (and the t-statistic for a t-test).

# Display the test stat
myChiSq$statistic

## X-squared 
##  3.753786

It also makes it easy to access the p-value:

# Display the p-value
myChiSq$p.value

## [1] 0.5853806

Note that this gives us the chi-square stat and p-value from the test (this can be useful, for example if I want to include in this text that the p-value was exactly 0.5853806 even though I wrote this before the class collected the data on Mini M&M’s).

But, the two we are interested in at the moment are the observed and expected values. Let’s look at them.

# Check observed
myChiSq$observed

Note that this is just our input data.

# Check expected
myChiSq$expected

Here, we have the expected counts, which R calculated as we described above. We can see that these match the row/column totals from our inputs as we expected (note, round is used to make the display prettier):

# Add margins to each
addmargins(myChiSq$observed)
round(addmargins(myChiSq$expected), 2)

To see the differences, which will tell us in what way our values are different than the expected values, we can just subtract them

# Display diff from expected
round(myChiSq$observed - myChiSq$expected , 2)

This tells us that there were, for example 2.96 fewer Blue Mini M&M’s than we expected and 3.22 fewer Regular Red M&M’s than we expected[*] I said magically, based on class data I didn’t have, see if you can figure out how, if you are interested .

30.6.4 Try it out

An earlier class collected data on Peanut M&M’s instead of Mini M&M’s. Their data (along with the Regular M&M data) are below. Graph the distributions for both groups in a way that allows direct comparison. Then, run and interpret a chi-square test.

# Enter new data
peanutMMdata <- rbind(
  Peanut  = c(10, 11, 25, 30, 16, 16),
  Regular = c(35, 19, 39, 37, 20, 26)
  )
colnames(peanutMMdata) <-
  c("Blue", "Brown", "Green",
    "Orange", "Red", "Yellow")

You should get a p-value of 0.24 and a plot like this:

30.7.1 Try it out

Make a sampling distribution for the Peanut M&M data and compare it to the results from chisq.test(). They should give very similar values.

30.8.1 Try it out

Is there an association between an individuals perceived change in financial status (column FIN) and zir perception of the direction the state is going in is (column STAT)? Assume that the outlook data are a representative sample of all Montana residents, and run an appropriate test. Don’t forget to make a plot.

You should get a p-value of 0.00975 and a plot that looks something like this (colors, labels etc. are all up to you). No, you do not need to worry abouot all of those tweaks, but you are welcome to try them if you are up for a challenge. I will share my code with anyone that sends me their plot and the code they used to get something similar.

Note that for 2 (and only for 2) groups, I like the use of stacked bars. It visually allows a comparison of either those that said “better” or “worse,” and it wastes less space. However, with 3 (or more) groups, it becomes too hard as the middle groups no longer line up well for easy comparisons (unless the group orders are meaningful, such as a five point ranking of how the state is doing).