Multiple comparisons of group-level means is a tricky problem in statistical inference. A standard practice is to adjust the threshold for statistical significance according to the number of pairwise tests performed. For example, according to the widely-known Bonferonni method, if we have 3 different groups for which we want to compare the means of a given variable, we would divide the standard significance level (.05 by convention) by the number of tests performed (3 in this case), and only conclude that a given comparison is statistically significant if its p-value falls below .017 (e.g. .05/3).In this post, we'll employ an approach often recommended by Andrew Gelman, and use a multi-level model to examine pairwise comparisons of group means. In this method, we use the results of a model built with all the available data to draw simulations for each group from the posterior distribution estimated from our model. We use these simulations to compare group means, without resorting to the expensive corrections typical of standard multiple comparison analyses (e.g. the Bonferroni method referenced above).
The Data
We will return to a data source that we have examined previously on this blog: walking data from the Accupedo app on my phone. Accupedo is a free step-counting app that I've been using for more than 3 years to track how many steps I walk each day. It's relatively straightforward to import these data into R and to extract the cumulative number of steps taken per hour.In the data used for this post, I exclude all measurements taken before 6 AM (because the step counts are cumulative we don't lose any steps by excluding these observations). This results in 18,117 observations of cumulative hourly step counts for 1113 days. The first observation is on 3 March 2015 and the last on 22 March 2018, resulting in just over 3 years of data!A sample of the dataset, called aggregate_day_hour, is shown below: oneday
dow
hour
steps
week_weekend
1
2015-03-07
Sat
16
14942
Weekend
2
2015-03-07
Sat
17
14942
Weekend
3
2015-03-07
Sat
18
14991
Weekend
4
2015-03-07
Sat
19
15011
Weekend
5
2015-03-07
Sat
20
15011
Weekend
6
2015-03-07
Sat
21
15011
Weekend
7
2015-03-07
Sat
22
15011
Weekend
8
2015-03-07
Sat
23
15011
Weekend
9
2015-03-08
Sun
11
2181
Weekend
10
2015-03-08
Sun
12
2428
Weekend
Oneday represents the calendar year/month/day the observation was recorded, dow contains the name of the day of the week, hour is the hour of the day (in 24-hour or &quot;military&quot; style), steps gives the cumulative step count and week_weekend indicates whether the day in question is a weekday or a weekend.Analytical Goal Our ultimate purpose here is to compare step counts across a large number of days, without having to adjust our threshold for statistical significance. In this example, we will compare the step counts for 50 different days. The standard Bonferonni correction would require us to divide the significance level by the number of possible pairwise comparisons, 1225 in this case. The resulting significance level is .05/1225, or .00004, which imposes a very strict threshold, making it more difficult to conclude that a given comparison is statistically significant.As mentioned above, the method adopted in this post is one that Andrew Gelman often advocates on his blog. This advice often comes in the form of: &quot;fit a multilevel model!,&quot; with a link to this paper by Gelman, Hill and Yajima (2012). I have to admit that I never fully understood the details of the method from the blog posts. However, by looking at the paper and the explanations and code presented in Chapter 12 of Gelman and Hill's (2007) book, I think I've figured out what he means and how to do it. This post represents my best understanding of this method.*In brief (see the paper for more details), we calculate a multilevel model using flat priors (e.g. the standard out-of-the-box model in the lme4 package). We extract the coefficients from the model and the standard deviation of the residual error of estimated step counts per day (also known as sigma hat), and use these values to simulate step count totals for each day (based on the estimated total step count and sigma hat). We then compare the simulations for each day against each other in order to make pairwise comparisons across days.Sound confusing? You're not alone! I'm still grappling with how this works, but I'll gently walk through the steps below, explaining the basic ideas and the code along the way.Step 1: Calculate the Model
In order to perform our pairwise comparisons of 50 days of step counts, we will construct a multi-level model using all of the available data. We'll employ a model similar to one we used in a previous post on this blog. The model predicts step count from A) the hour of the day and B) whether the day is a weekday or a weekend. We specify fixed effects (meaning regression coefficients that are the same across all days) for hour of the day and time period (weekday vs. weekend). We specify a random effect for the variable that indicates the day each measurement was taken (e.g. a random intercept per day, allowing the average number of estimated steps to be different on each day). We center the hour variable at 6 PM (more on this below).We can calculate this model using the lme4 package. We will also load the arm package (written to accompany Gelman and Hill's (2007) book) to extract model parameters we'll need to compute our simulations below.# load the lme4 packagelibrary(lme4)# load the arm package# to extract the estimate of the # standard deviation of our residuals # (sigma hat)library(arm)# center the hour variable for 6 pmaggregate_day_hour$hour_centered Column')# sort the matrix by the dates# first turn the date columns# into date types in R using# the lubridate packagelibrary(lubridate)melted_cormat$x melted_cormat Our final formatted matrix for plotting, called melted_cormat, looks like this (first 10 rows shown):
x
y
count
meaningful_diff
1
2015-03-18
2015-03-18
0
No Difference
2
2015-03-18
2015-03-29
222
No Difference
3
2015-03-18
2015-05-12
737
No Difference
4
2015-03-18
2015-06-29
222
No Difference
5
2015-03-18
2015-07-19
515
No Difference
6
2015-03-18
2015-09-15
768
No Difference
7
2015-03-18
2015-09-22
884
No Difference
8
2015-03-18
2015-10-14
683
No Difference
9
2015-03-18
2015-10-18
837
No Difference
10
2015-03-18
2015-10-22
728
No Difference
The count variable here gives the number of times that the y date simulations were greater than the x date simulations. So the second row above indicates that the simulations on '2015-03-29' were greater than those for '2015-03-18' a total of 222 times. We can make the heatmap with the following code:# make the heatmapggplot(melted_cormat, aes(as.factor(x), as.factor(y), fill = meaningful_diff)) + # tile plot geom_tile() +# remove x and y axis titles# rotate x axis dates 90 degrees theme(axis.title.x=element_blank(),axis.title.y=element_blank(), axis.text.x = element_text(angle = 90, vjust = .5)) + # choose the colors for the plot # and specify the legend title scale_fill_manual(name = &quot;Comparison&quot;, values=c(&quot;#999999&quot;, &quot;#800000&quot;)) + # add a title to the plot labs(title = &quot;Pairwise Comparisons of Model-Estimated Step Counts at 6 PM&quot;)
Which gives us the following plot:
How can we read this plot? All of our 50 dates are displayed on the x and y axes. When the row step count is meaningfully larger than the column step count (according to our method, e.g. if more than 950 of the simulations for the row date are greater than those for the column date), the cell is colored in maroon. Otherwise, cells are colored in grey.Rows that have a lot of maroon values are days that have higher step counts. To take a concrete example, the row representing '2016-08-08' has many red values. If we look at the estimated step count for that day on the first graph above, we can see that it has the highest predicted step count in the sampled dates - over 25,000! It therefore makes sense that the estimated step count on this day is &quot;significantly&quot; larger than those from most of the other dates.Columns that have many red values are days that have especially low step counts. For example, the column representing '2015-12-25' is mostly red. As we saw above, this date has the lowest predicted step count in our sample - under 5,000! It is no wonder then that this date has a &quot;significantly&quot; lower step count than most of the other days. How Baysian Am I?
I'm a newbie to Baysian thinking, but I get the sense that Baysian statistics comes in many different flavors. The approach above strikes me as being a little bit, but not completely, Baysian.That's So BaysianThis approach is Baysian in that it uses posterior distributions of daily step counts to make the comparisons between days. The approach recognizes that observed daily step counts are just observations from a larger posterior distribution of possible step counts, and we make explicit use of this posterior distribution in the data analysis. In contrast, frequentist methods basically only make use of point estimates of coefficients and the standard errors of those estimates to draw conclusions from data. Not So Fast!There are some things about this perspective that aren't so Baysian. First, we use flat priors in the model; a &quot;more Baysian&quot; approach would assign prior distributions to the intercept and the coefficients, which would then be updated during the model computation. Second, we use the point estimates of our coefficients to compute the estimated daily step counts. A &quot;more Baysian&quot; approach would recognize that the coefficient estimates also have posterior distributions, and would incorporate this uncertainty in the simulations of daily step counts.What Do I Know?This is my current best understanding of the differences in Baysian and frequentist perspectives as they apply to what we've done here. I'm reading Statistical Rethinking (and watching the accompanying course videos) by Richard McElreath (love it), and these differences are what I've understood from the book and videos thus far. The method used in this blog post is a nice compromise for someone comfortable with frequentist use of multi-level models (like myself). It requires a little bit more work than just interpreting standard multilevel model output, but it's not a tremendous stretch. Indeed, the more &quot;exotic&quot; procedures (for a Baysian newbie) like assigning priors for the model parameters are not necessary here.Summary and Conclusion
In this post, we used multilevel modeling to construct pairwise comparisons of estimated step counts for 50 different days. We computed the multilevel model, extracted the coefficients and sigma hat value, and computed 1000 simulations from each day's posterior distribution. We then conducted pairwise comparisons for each day's simulations, concluding that a day's step count was &quot;significantly&quot; larger if 950 of the 1000 simulations were greater than the comparison day. We used a heat map to visualize the pairwise comparisons; this heat map highlighted days with particularly high and low step counts as being &quot;significantly&quot; different from most of the other days. This allowed us to conduct these comparisons without the expensive adjustments that come with lowering our significance threshold for multiple comparisons!  Coming Up NextIn the next post, we will move away from classical statistics and talk about machine learning techniques. Specifically, we will use a type of deep learning model to automatically generate band names.Stay tuned!-------* Please correct me if I'm wrong!