Mmmmmmm

See http://blogs.sas.com/content/iml/2017/02/20/proportion-of-colors-mandms.html for full expository. Any text here is a riff from that. Rick did the hard work.

library(knitr)
library(kableExtra)
library(scimple) # devtools::install_github("hrbrmstr/scimple")
library(hrbrthemes)
library(ggimage)
library(tidyverse)

options(knitr.table.format = "html") 

data_frame(
  color_name = c("Red", "Orange", "Yellow", "Green", "Blue", "Brown"),
  official_color = c("#cb1634", "#eb6624", "#fff10a", "#37b252", "#009edd", "#562f14"), 
  count = c(108, 133, 103, 139, 133, 96),
  prop_2008 = c(0.13, 0.20, 0.14, 0.16, 0.24, 0.13),
  imgs=c("im-red-lentil.png", "im-orange-lentil.png", "im-yellow-lentil.png",
         "im-green-lentil.png", "im-blue-lentil.png", "im-brown-lentil.png")
) %>% 
  mutate(prop = count / sum(count),
         color_name = factor(color_name, levels=color_name)) ->  mms

cap_src <- "Source: Riffed from The DO Loop <http://blogs.sas.com/content/iml/2017/02/20/proportion-of-colors-mandms.html>"

SAS M&M’s Measurements

The breakroom containers at SAS are filled from two-pound bags. So as to not steal all the M&M’s in the breakroom, [Rick] conducted this experiment over many weeks in late 2016 and early 2017, taking one scoop of M&M’s each week. The following data set contains the cumulative counts for each of the six colors in a sample of size N = 712:

ggplot(mms, aes(color_name, prop, fill=official_color)) +
  geom_col(width=0.65) +
  scale_y_percent(limits=c(0, 0.21)) +
  scale_fill_identity(guide = FALSE) +
  labs(x=NULL, y=NULL, title=sprintf("Sample distribution of colors for M&Ms [n=%d]", sum(mms$count)),
       caption=cap_src) +
  theme_ipsum_rc(grid="Y", axis="x")

Data overview

The data set includes the most recent published distribution is from 2008:

mms %>% 
  mutate(prop=scales::percent(prop), 
         prop_2008=scales::percent(prop_2008)) %>% 
  select(Color=color_name, Frequency=count, Percent=prop, `Test Percent`=prop_2008) %>% 
  kable(align="lrrr") %>% 
  kable_styling(bootstrap_options=c("hover", "condensed", "responsive"), full_width = FALSE)

Color	Frequency	Percent	Test Percent
Red	108	15.17%	13%
Orange	133	18.68%	20%
Yellow	103	14.47%	14%
Green	139	19.52%	16%
Blue	133	18.68%	24%
Brown	96	13.48%	13%

Chi-Square Test Results

Let’s test these proportions:

chisq.test(mms$count, p=mms$prop_2008) %>% 
  broom::tidy() %>% 
  kable(align = "rrrl") %>% 
  kable_styling(bootstrap_options=c("hover", "condensed", "responsive"), full_width = FALSE)

statistic	p.value	parameter	method
17.35305	0.0038767	5	Chi-squared test for given probabilities

The chi-square test rejects the test hypothesis at the α = 0.05 significance level (95% confidence). In other words, the distribution of colors for M&M’s in this 2017 sample does NOT appear to be the same as the color distribution from 2008! You can see this visually from the bar chart: the red and green bars are too tall and the blue bar is too short compared with the expected values.

You need a large sample to be confident that this empirical deviation is real. After collecting data for a few weeks, [Rick] did a preliminary analysis that analyzed about 300 candies. With that smaller sample, the difference between the observed and expected proportions could be attributed to sampling variability and so the chi-square test did not reject the null hypothesis. However, while running that test [Rick] noticed that the green and blue colors accounted for the majority of the difference between the observed and theoretical proportions, so [Rick] decided to collect more data.

Simultaneous confidence intervals for the M&M proportions

Rick used some SAS scripts he had written for a previous blog post. We’ve got a few packages for computing simultaneous confidence intervals in R but one of the “better” ones does the computation and prints out the results (literally with print()) vs return data, so I made a new package that does the same computations and returns tidy data frames for the confidence intervals. The original version of this document had the computation in a hidden code block but a package is much better and it includes a function that can compute multiple SCIs and return them in a single data frame, similar to what binom::binom.confint() does.

mms <- bind_cols(mms, scimp_goodman(mms$count, alpha=0.05))

mms %>% 
  select(Color=color_name, Count=count, Proportion=prop, 
                `95% Lower`=lower_limit, `95% Upper`=upper_limit) %>% 
  kable(align=c("lrrrr"), caption="Simultaneous confidence Intervals (Goodman method)") %>% 
  kable_styling(bootstrap_options=c("hover", "condensed", "responsive"), full_width = FALSE)

Simultaneous confidence Intervals (Goodman method)

Color	Count	Proportion	95% Lower	95% Upper
Red	108	0.1516854	0.1196016	0.1905134
Orange	133	0.1867978	0.1513616	0.2282982
Yellow	103	0.1446629	0.1133216	0.1828844
Green	139	0.1952247	0.1590634	0.2372872
Blue	133	0.1867978	0.1513616	0.2282982
Brown	96	0.1348315	0.1045758	0.1721577

The table indicates that the published 2008 proportion for blue (0.24) is far outside the 95% confidence interval, and the proportion for green (0.16) is just barely inside its interval. That by itself does not prove that the 2008 proportion are no longer valid (we might have gotten unlucky during sampling), but combined with the earlier chi-square test, it seems unlikely that the 2008 proportions are applicable to these data.

Cleveland Comparison

We incorporate the Cleveland M&M plant official color distribution values into our data frame and show where their expected distributions fall compared to the sample and where both fall within our computed confidence intervals:

NOTE: Uncomment the commented bits and comment out the geom_image() if you want “normal” points.

url_base <- "http://www.mms.com/Resources/img/"

mms %>% 
  mutate(imgs=sprintf("%s%s", url_base, imgs)) %>% 
  mutate(cleveland_prop=c(0.131, 0.205, 0.135, 0.198, 0.207, 0.124)) %>% 
  ggplot() +
  geom_segment(aes(x=lower_limit, xend=upper_limit, y=color_name, 
                   yend=color_name, color=official_color), size=2) +
  # geom_point(aes(prop, color_name, fill=official_color), 
  #            size=8, shape=21, color="white") +
  geom_image(aes(prop, color_name, image=imgs)) +
  # geom_text(aes(prop, color_name, label="m"), 
  #           color="white", family="Times", fontface="bold") +
  geom_point(aes(cleveland_prop, color_name, color=official_color),
             shape="|", size=6) +
  scale_x_percent(limits=c(0.095, 0.25)) +
  scale_color_identity(guide = FALSE) +
  scale_fill_identity(guide = FALSE) +
  labs(x="Proportion", y=NULL, 
       title="Observed vs 2017 Proportions for M&M Candies",
       subtitle=sprintf("95%% Simultaneous Confidence Intervals, [N=%d]",
                        sum(mms$count)), caption=cap_src) +
  theme_ipsum_rc(grid="X") +
  theme(panel.background=element_rect(color="#5b5b5b", fill="#5b5b5b")) +
  theme(axis.text.x=element_text(hjust=c(0, 0.5, 0.5, 1)))