Now that we’re done with the first Data Analysis Report, we’re going to be starting on the second half of the class, in which we move from the basics of how to code for research–learning how R works, learning basic concepts of control flow and data wrangling–to the why of our coding work: analyzing data. This will be the focus of the rest of the class!
In the last few weeks you actually started to learn a lot about how to do exploratory and descriptive data analysis: how to make basic, useful exploratory plots, how to generate summaries of data that tell you something about what you’ve got. You also have spent some time putting together your first data-analysis report, which made use of those skills.
So this week we will start with a summary of what inferential data analysis is, and why we might want to employ it in our work. We’re going to start off very simple with some very basic sensory data–right where I feel at home. I’m going to spend a little time giving you a biography of the data so that you will be able to relate it to your own work.
Then, we’ll step back and talk about how small experiments like this one can be used for inferring knowledge about more general situations. We’ll talk about heuristic (traditional) statistical approaches and resampling approaches to inference. Finally, we’ll talk about the counterexample of how we usually learn and think about statistical inference in science, and how it contrasts with what we’re doing here.
Motivation: inference from experiments
Today we’re going to work with a set of data my undergraduate students in FST 3024 (Principles of Sensory Evaluation) collected for me in 2023. Each year, the class does a service-learning project for a small food producer. In 2023, the students worked with Joell Eifert, our undergraduate-research coordinator, to help Homestead Creamery test a new production method for their Hokie Tracks Ice Cream:
Hokie Tracks Ice Cream
According to Homestead
From farm to market, the collaboration between the nationally ranked Virginia Tech College of Agriculture and Life Sciences’ Department of Food Science and Technology and the award-winning Homestead Creamery, brings innovated food science research to the market in the form of delicious Hokie-inspired ice cream. Hokie Tracks is the perfect sweet & savory combination of vanilla ice cream with salted caramel swirls and chocolate covered pretzel pieces!
The basic issue at hand was that Homestead wanted to scale up production of Hokie Tracks, replacing labor-intensive hand-mixing of the caramel swirls and chocolate pieces with a piece of equipment known as a “variegator”, which would allow continuous instead of batch production. Homestead asked us to test whether people could tell the difference between the same recipe made with the two production methods, and whether subjects preferred one or the other.
We conducted what is known as a “preference test”: consumers were given two samples without labeling information, asked to taste both, and forced to indicate whether they preferred one sample or the other. If, as a group, consumers had no true preference, about 50% should select one sample and the other 50% would select the other. We also collected “free comment” data that allowed consumers to explain their choices.
library(readxl)preference_data<-read_excel("data/Week 9/Hokie Tracks Raw Data.xlsx", range ="A1:C80")product_codes<-read_excel("data/Week 9/Hokie Tracks Raw Data.xlsx", range ="D2:F4")glimpse(preference_data)
I’ve imported both the preference data itself and a little reference box in the workbook, entered by the TA at the time, that gives labels for the codes given to Preferred Sample in the preference_data data frame. Otherwise, we’d have to guess what sample 1 and sample 2 were. How can we join these two tables together so we can have informative labels on our data?
preference_data<-preference_data|>left_join(product_codes, by =c("Preferred Sample"="Sample Number"))
We also might want to tell R that our “numeric” data is in fact just categorical, by changing appropriate columns into factor columns.
preference_data<-preference_data|>mutate(preferred_sample =as.factor(`Sample Name`))|># And we can reorganize and drop unnecessary columnsselect(Panelist, preferred_sample, `Open Comments`)
Now we can use skim() to get us a nice summary (although truly this is overkill with such simple data).
skim(preference_data)
Data summary
Name
preference_data
Number of rows
79
Number of columns
3
_______________________
Column type frequency:
character
1
factor
1
numeric
1
________________________
Group variables
None
Variable type: character
skim_variable
n_missing
complete_rate
min
max
empty
n_unique
whitespace
Open Comments
0
1
13
367
0
79
0
Variable type: factor
skim_variable
n_missing
complete_rate
ordered
n_unique
top_counts
preferred_sample
0
1
FALSE
2
Han: 50, Mac: 29
Variable type: numeric
skim_variable
n_missing
complete_rate
mean
sd
p0
p25
p50
p75
p100
hist
Panelist
0
1
40
22.95
1
20.5
40
59.5
79
▇▇▇▇▇
We are essentially interested in one thing in these data: do we think that there are enough people expressing a preference for a particular sample that we believe it’s a “true” preference? Do we think that overall people are not just selecting samples at random?
preference_data|>count(preferred_sample)
# A tibble: 2 × 2
preferred_sample n
<fct> <int>
1 Hand-mixed 50
2 Machine-fed 29
It certainly looks like more people prefer the Hand-mixed sample. But is this just a random sample? Let’s look into how we might answer this question.
We have two questions we might want to answer from this simple example:
Is the difference between numbers of people who prefer Hand-mixed over Machine-fed large? Are we seeing a real effect?
If we repeated this experiment (say, by putting these two products both on the market at the same time), what kind of range of behavior (proportion of people choosing Hand-fed over Machine-fed) would we observe? What kind of range of proportions is compatible with these data?
Let’s use what we’ve learned about simulation and random sampling to answer each of these questions using a coding approach, and then circle back to a statistical heuristic approach.
“A resampling approach” to statistics
Today we are going to learn some approaches to thinking about inference with data that don’t rely on assumptions about formula or rules. These are called “resampling” approaches to statistics, and they instead rely on us being a bit clever about using computational power to do simulations on our datasets. This sounds fancy and intimidating, and I am pleased to say that while it may be a little bit alien, it is not!
As a caveat, however, I will note that we are going to just learn about the basics of these approaches today (and in this class in general). There are whole books on this topic, and more sophisticated simulationist approaches like Bayesian statistics are going to be the most precise and powerful models (for example see McElreath 2020). My hope is that giving you a taste of a simulationist approach shows you how powerful you are once you know the basics of coding for research!
Confidence intervals around statistics (the “bootstrap”)
The first questions we have is about our sample means, which are “point estimates” for the actual, population mean. What am I saying here? Well, let’s take a minute to unpack the idea of a point estimate, and then apply it to our example.
First, let’s imagine we have the normal distribution:
This is our population distribution. We know that its true population mean is 0, and its true population standard deviation is 1, by definition. Now, let’s pull 3 samples of 100 observations from the normal distribution with mean = 0 and sd = 1, and take a look at where they go on the graph:
We can see what is happening here by plotting these onto our normal distribution:
p+geom_vline(xintercept =mean(normal_samples$x1), color ="red", size =0.75)+geom_vline(xintercept =mean(normal_samples$x2), color ="red", size =0.75, linetype ="dashed")+geom_vline(xintercept =mean(normal_samples$x3), color ="red", size =0.75, linetype ="dotted")
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Even though our samples are coming from the exact same population, they have different means (and standard deviations) because of sampling error. Although, as a side note, I will say that “error” makes it sound like something we can fix–in fact, this variation is fundamental and cannot be avoided.
To conclude: even when we are sampling from a perfectly defined and known population, our samples will only approximate that population. But in the long term, if we draw enough samples they will end up in the limit being unbiased approximations. Armed with that knowledge, let’s return to our example.
These rules apply to all kinds of random variables. In our ice-cream preference example, our data show an observed proportion of \(p_{obs}=50/79\). We could assume that, if there were truly no population preference, our data is a sample coming from a binomial distribution with \(n = 79\). Then we could look at a similar story.
We will use a slightly different plotting strategy because the binomial distribution is a discrete distribution: the number of successes is a positive integer (\(0, 1, 2, ...\)), not a range of real numbers.
p<-tibble(x =0:79, y =dbinom(x, size =79, prob =1/2))|>ggplot(aes(x =x, y =y))+geom_col(fill ="grey")+theme_classic()+labs(x =NULL, y ="density")+scale_y_continuous(expand =expansion(mult =c(0, 0.05)))p
We can see that this is centered right around \(79/2 = 39.5\), which gives equal probability to \(39/79\) or \(40/70\) people saying they prefer the hand-mixed ice cream.
If we take 3 random probability draws from the binomial distribution with these properties, we’ll get results that are close to but not identical to that most likely outcome:
set.seed(9)random_preference<-tibble(draw =1:3, proportion =rbinom(n =draw, size =79, prob =1/2), linetype =as.factor(draw+1))p+geom_vline(aes(xintercept =proportion, linetype =linetype), data =random_preference, show.legend =FALSE, color ="red")
So we have the same situation! Even when we have a “known” population, drawing finite, random samples from that distribution will give us results that only approximate it.
For kicks, we might remember that “large enough” binomial distributions are well-approximated by a normal distribution with \(\mu = p\) and \(\sigma = \sqrt{p(1-p)}\). How well? Let’s see:
p+geom_function(fun =dnorm, n =101, args =list(mean =39.5, sd =sqrt(79/4)), color ="red")
Pretty well, I’d say!
A bootstrapping example
In this example, for our sample we had a specific 79 people taste each ice cream and report which one they preferred. Thus, when we calculate our point-estimate summary statistics, we are basing that estimated group mean proportion (or SD or whatever) on those 79 observations. But how would it change if we had 79 other observations? What we know is our sampling-error biased mean, and what we want to get a better picture of is the population these means came from. In our explanation above, what this is equivalent to is trying to infer the shape of the normal distribution with only 1 of the 100-sample draws (that we showed were full of errors). How can we do this?
The answer here is that we can use an approach called “resampling” to build a simulated population from our observed sample, and then use that population to understand how representative our observed point-estimate statistics are. This sounds complicated, but in practice it is very intuitive. The description from Bruce and Bruce (2017) for the steps to do this is quite clear (paraphrased here in almost pseudocode):
From your sample, draw a sample of the same size with replacement.
Record the mean, SD, or other relevant statistic from your “resampled” sample and record the information.
Repeate steps 1-2 a large number of times (a typical number is \(N=10,000\), but this can vary depending on computational intensivity)
Use the new dataset of statistics to calculate a confidence interval or to draw boxplots or histograms to compare to your point estimate.
The reason this approach works well is the “with replacement” bit. This means that, rather than just shuffling our data, we will be sometimes drawing the same sample multiple times. Over the long run, this creates a population that would plausibly generate our sample, as shown in Bruce and Bruce (2017):
Again, this sounds like we’re doing something mystically complicated, but the steps above are each things we already know how to do. Let’s take a look in practice, with our dataset. We have 4 different means we want to resample, one for each of the different treatments. We’ll write a little for loop to do this for us (although you could check out the boot package for more functionality).
# Here we set a number of iterations for the bootstrap. If your computer is# slower, consider setting this to 100 just for the example.boot_reps<-1000# As a bonus, here's how we make a progress bar to track the taskpb<-progress::progress_bar$new(total =boot_reps)# for each step of the loop, we will generate a new sample for each of our# treatments and take its mean, which is the proportion in a set of 0s and 1s.# We will store the results of the proportion of people who preferred hand-mixed# ice creambooted_proportions<-tibble(boot_id =numeric(), boot_proportion =numeric())for(iin1:boot_reps){preferred_hand_mixed<-preference_data|># by setting prop(ortion) = 1 and asking for replacement, we get a new "resample"slice_sample(prop =1, replace =TRUE)|>count(preferred_sample)|># Here we cheat: we know that the sum of n is 79 each time. We could do# this programatically with a few more steps. How?mutate(proportion =n/79)|># Get the actual cell containing the proportion of hand-mixed ice cream# preference in the new sample_[[1, 3]]# Let's not forget to store our results!booted_proportions<-bind_rows(booted_proportions,tibble(boot_id =i, boot_proportion =preferred_hand_mixed))# Update our progress barpb$tick()}
Now we have a big (1,000 row) data frame, in which each row is a plausible simulation of means we could have gotten if our experimental sample had been different, but drawn from the same population! Let’s go ahead and do a little basic data exploration on this dataset.
p_boot<-booted_proportions|>ggplot(aes(x =boot_proportion))+geom_histogram(bins =25, fill ="grey", color ="white")+theme_classic()+scale_y_continuous(expand =expansion(mult =c(0, 0.05)))p_boot
We can do a couple things with this. We can compare our observed proportion, \(p_{obs}\), to the mean of the bootstrapped proportions to see how representative it is. We can also look at the range of these bootstrapped proportions: what, for example, is an interval that contains 95% of the observations?
p_boot+geom_vline(xintercept =50/79, color ="red")+geom_vline(xintercept =mean(booted_proportions$boot_proportion), color ="red", linetype =2)+annotate(geom ="label", x =51/79, y =25, label ="italic(p)[obs]", parse =TRUE)+annotate(geom ="label", x =mean(booted_proportions$boot_proportion+1/79), y =55, label ="italic(bar(p))[boot]", parse =TRUE)
We can see that the lines are practically overlapping, the actual values are:
Note that there’s nothing magical about the 95% interval–it is just “pretty big”. It tells us that 95% of the bootstrapped proportions fall between 53% and 73% preferring hand-mixed ice cream. For our basic question–do more people prefer hand-mixed ice cream–this is evidence for us to conclude that, if our sampling was an effective representation of the population, it is pretty unlikely that people prefer the machine-mixed recipe. Very few of our “resamples” give us a proportion of 50% or under.
We can actually plot these numbers to make presenting them to stakeholders easier:
p_boot+geom_histogram(data =booted_proportions|>mutate(fill =if_else(boot_proportion<=0.5, "a", "b")), mapping =aes(fill =fill, x =boot_proportion), bins =25, color ="white", show.legend =FALSE)+scale_fill_manual(values =c("black", "grey"))+annotate(geom ="rect", xmin =boot_quantiles[1], xmax =boot_quantiles[2], ymin =0, ymax =Inf, fill ="pink", alpha =1/3)+annotate(geom ="text", label ="95% range of\nbootstrap values", x =0.54, y =95, hjust =0)
This tells us a couple things:
Any values of 0.5 (50%) or less are well outside the lower quantile of our 95% range, which means that–if we think this is a “wide enough” range–they are not very likely to be observed in another test.
This data would be evidence for us that there is some kind of preference in the population for the hand-mixed ice cream.
The advantage of the bootstrap is that it is nonparametric: it makes no assumptions about how the data are generated and works equally well for any kind of statistic or any kind of data. It just requires writing some pseudocode to properly generate the loop, and then you have an empirical distribution for your data. Neat!
Is that difference important?
So far we have used resampling to mostly answer questions about our point-estimates: how reliable and stable are our sample statistics? These approaches work well (conceptually) when we are trying to estimate a single parameter, like a proportion. But sometimes we are interested in comparing two parameters: proportions, means, even variances. The key difference (in my mind) is that in these cases–comparing two statistics or parameters–we have some sort of group ID.
In the case of our ice cream data, we would not apply that: each observation we have is a single choice (hand-mixed or machine extruded). So we’ll look at another class project from the same year: a comparison of Black Walnut Syrup to Maple Syrup.
Some of you may know that maple syrup is made from tapping maple trees for sap and then boiling down the sap. It is not made from corn syrup (former Vermont resident here).
From Wikipedia
Black Walnut Syrup is a novel product made in the same way, but from Black Walnut trees, which are abundant in Appalachia.
From Cornell Extension
Students in FST 3024 in 2023 used an affective test to determine whether consumers in Blacksburg preferred one syrup over the other. Unlike the ice cream data, these students asked research subjects to taste and rate the two syrups on a scale of 1-9, where 1 = “Dislike Very Much”. Therefore, we have a bunch of (quasi-)interval data. Let’s take a quick look.
syrup_data<-read_excel("data/Week 9/Walnut Syrup Test Raw Data.xlsx", sheet =2, range ="A1:C131")syrup_labels<-read_excel("data/Week 9/Walnut Syrup Test Raw Data.xlsx", sheet =2, range ="E2:F4")syrup_data<-left_join(syrup_data,syrup_labels, by =c("Blinding Code"="BC"))glimpse(syrup_data)
# A tibble: 2 × 3
`Sample Name` mean_liking sd
<chr> <dbl> <dbl>
1 Black Walnut Syrup 6.18 1.84
2 Maple Syrup 6.43 1.72
And we could even look at the distributions of liking, perhaps as boxplots:
syrup_data|>ggplot(aes(x =`Sample Name`, y =`Liking Rating`))+geom_boxplot()+theme_classic()+labs(x =NULL)
We’d like to know if this small difference in liking is meaningful–we do have 65 subjects in the study, so perhaps a 1/3 point difference is meaningful? We could bootstrap the data and produce a robust, non-parametric estimate for the range of each mean, etc. And this would actually work great, but I am setting us up for another approach!
Specifically, here we have an \(x\) variable: sample type. We want to know if knowing the sample type (Black Walnut or Maple) gets us any more information to predict liking.
These are the kinds of questions that lead to testable hypotheses. In future classes, we are going to touch on models of statistical analysis like the linear model which give us formulas to deal with these questions. In particular, we might use ANOVA (ANalysis Of VAriance) to answer this kind of question (where \(x\) is categorical), or linear regression (they’re really the same thing). But today we’re going to just use some raw computational power to do so.
Does knowing the type of syrup help us predict liking?
We know we have 130 observations in our data set: half (65) are ratings of Maple Syrup, and half are ratings of Black Walnut Syrup. If that doesn’t tell us anything, then switching labels of observations will not change our data. The term for this is permuting our data. We can do this easily using tools like sample().
# A tibble: 2 × 3
`Sample Name` original_mean_liking permuted_mean_liking
<chr> <dbl> <dbl>
1 Black Walnut Syrup 6.18 6.43
2 Maple Syrup 6.43 6.18
In this case, quite by coincidence, we get reversed mean values. But if we run the same code again (without set.seed(9)) we’ll get slightly different numbers. If we do this enough times, we’ll be able to see whether or not labels are meaningful. How? You guessed it: for() loops!
# A tibble: 20,000 × 3
name mean_liking perm_id
<chr> <dbl> <int>
1 Black Walnut Syrup 6.2 10000
2 Maple Syrup 6.42 10000
3 Black Walnut Syrup 6.38 9999
4 Maple Syrup 6.23 9999
5 Black Walnut Syrup 6.25 9998
6 Maple Syrup 6.37 9998
7 Black Walnut Syrup 6.29 9997
8 Maple Syrup 6.32 9997
9 Black Walnut Syrup 6.12 9996
10 Maple Syrup 6.49 9996
# ℹ 19,990 more rows
We can now look at this in a number of ways. First, we’ll look at the distribution of mean liking ratings, and then we’ll look at our distribution of differences in mean liking compared to our actual results.
permutation_results|>ggplot(aes(x =mean_liking, fill =name))+geom_density(alpha =1/2)+scale_fill_manual(values =c("orange", "maroon"))+theme_classic()+labs(fill =NULL, x ="Liking Rating", y =NULL)
These are so close to identical that it is almost impossible to tell them apart. This implies that our observed difference is probably not very meaningful, but let’s look at it to be sure:
observed_diff<-syrup_data|>select(-`Blinding Code`)|>pivot_wider(names_from =`Sample Name`, values_from =`Liking Rating`)|>mutate(diff =`Black Walnut Syrup`-`Maple Syrup`)|>reframe(observed_diff =mean(diff))|>pull(observed_diff)permuted_diffs<-permutation_results|>pivot_wider(names_from =name, values_from =mean_liking)|>mutate(diff =`Black Walnut Syrup`-`Maple Syrup`)p_perm<-permuted_diffs|>ggplot(aes(x =diff))+geom_histogram(fill ="grey", color ="white", bins =24)+geom_vline(xintercept =observed_diff, color ="red")+theme_classic()+labs(x ="Mean Difference", y =NULL)+annotate(geom ="label", x =-0.3, y =1000, label ="Observed Difference", hjust =1)p_perm
We can see that our observed mean difference, as we expected from our density plots, fall squarely into the middle of our distribution of permuted mean differences. If we wanted to be really sure, we could examine the quantiles, just as we did for the bootstrapped results.
# And for a reminder, our observed differenceobserved_diff
[1] -0.2461538
Our observation falls squarely within this range. If the labels were meaningful–Walnut vs Maple–we’d expect that scrambling them would give us less information than we had when the labels were correct. That doesn’t seem to be the case here!
As a final note, our permutation approach here does not require any of the assumptions that we typically need to make, about normality, independence of errors, etc. It does not matter that (as we could show) our data is very non-normal. Instead, we have just reasoned out a way to examine sources of variation and compare them without making any assumptions of how or why our data is shaped how it is. However, please do note that working out these permutations for complex designs is quite difficult, and often computationally expensive. When your designs become more complex, better-designed models for simulation like Bayesian approaches (McElreath 2020) will work much better!
What do we know about statistics?
Finally, as we move into inference, I want to quickly review some ideas about statistical inference for scientific research that we typically take for granted.
We usually learn about statistics in conjunction with the scientific method in a research-methods or stats class:
(OK, I just can’t stop myself here, but this is one of the problems with our framework, as we’ll see below. The “Conclusion” in this cartoon actually is never reached in our idea of the Null Hypothesis Test, because all we can do is reject the null. This is counterintuitive and also not a necessary part of the scientific method.)
Typically, the basic stats classes we all take in our undergrad and grad careers follow a trajectory that is something like this:
Talk about probability, usually using coin flips and dice rolls.
Introduce the idea of a hypothesis as a falsifiable statement. Concentrate on the counterintuitive but key idea of the null hypothesis.
Talk about how to calculate key statistics about samples: central tendencies (means, medians, etc), variations (variance, deviation).
Introduce probability distributions based on situations, typically:
The binomial distribution for coin flips.
The normal distribution for something like heights of students in the class. Don’t forget about the t-distribution, too…
Maybe some more esoteric distributions like the Poisson distribution or the \(\chi^2\)-distribution.
Introduce a bunch of heuristic calculations for how to take a sample you think comes from some distribution to the theoretical null distribution. These are statistical tests.
This is where you usually learn about a menu of tests that you select based on how you believe your sample is distributed, your specific null hypothesis, etc, sch as the following:
Flowchart for selecting a statistical test, taken from the Nishimura Lab.
I think we have all seen a chart like this one, and we’ve all experienced the problem of asking a question like:
OK, I collected the data from my experiment. Now which is the right statistical test to conduct?
…
Uh…
In order to use this kind of chart, we need to make a lot of assumptions. We have to have a well-defined, falsifiable hypothesis (which we should really have before we start collecting data, but…), we need to understand the various tests enough to understand how to use them, and we need to understand the various assumptions underlying each test (this is the point at which we tend to fail the most often).
I don’t like these charts and this melange of statistical methods for all of these reasons! In particular, I dislike them because the assumptions that underlie the tests are important, but because they are rarely achieved in practice our tests are often failing in unpredictable ways. But even more so, I dislike them because, for many simple problems, we simply do not need to be going to all this trouble.
These types of charts, and the statistical tests they lead to, are the result of the need for heuristic approaches to statistics: as we talked about in Week 3, statistics originated in the 17th centuries with studies of probability in the absence of simulation. It was simply not practical to simulate “null” distributions, so early theorists developed very clever, computationally tractable approaches that avoided the need for simulation. But now we have computers. And we have learned to write simple code to make simulations for ourselves, which can often (much more intelligibly) help us explore our data and test our hypotheses.
So, on that note, and with what we’ve learned today, I look forward to trying out applications of this code with you on Thursday!
BONUS: Bayesian preference test
library(tidyverse)library(brms)library(tidybayes)# We know that 29/79 subjects preferred the machine-fed.m1<-brm(data =list(mf =29), family =binomial(link ="identity"), formula =mf|trials(79)~0+Intercept, iter =5000, chains =4, cores =4, warmup =1000, seed =1)as_draws_df(m1)|>ggplot(aes(x =b_Intercept))+stat_halfeye(point_interval =mode_hdi, fill ="darkorange", color ="maroon")+geom_ribbon(data =tibble(x =c(0.475, 0.525), ymin =-Inf, ymax =Inf),aes(x =x, ymin =ymin, ymax =ymax), inherit.aes =FALSE, fill ="pink", alpha =1/2)+annotate(geom ="text", x =0.5, y =0.5, label ="ROPE:\n0.5 ± .025")+theme_bw()
