How to use this document. Run the code along with
me. Comments start with #. Each section mirrors the live
lab; code chunks show the command and its output together.
Today’s plan
- Clean environment + setup
- Import & understand the data (McGrath experiment)
- Comparing means: two-sample t-test
- Quick practice
Clean Environment +
Setup
rm(list = ls()) # clean environment
# Load the packages
library(experimentr) # provides the McGrath dataset
library(dplyr)
library(ggplot2)
library(stats)
In the live session we install experimentr with
install.packages("experimentr") and set the working
directory with setwd(...). Here the dataset ships with the
package, so no download or setwd() is needed.
Import & Understand
Data (McGrath Experiment)
# Load the McGrath dataset (comes with experimentr)
data(mcgrath)
data <- mcgrath
Background: the McGrath experiment. In 2016, Mary C.
McGrath, Peter M. Aronow, and Vivien Shotwell published a randomized
controlled trial (RCT) on the effect of a chocolate scent on
bookstore sales. A bookstore was randomly assigned to have either:
treatment = 1 → chocolate scent present, or
treatment = 0 → control (no scent).
The dataset includes five variables:
treatment |
1 = treatment, 0 = control |
book |
Sales of books |
food |
Sales of foods |
coffee |
Sales of bulk coffee, tea, or spices |
grandtotal |
Total of book, coffee, and food sales |
Goal: estimate the causal effect of the chocolate
scent on total sales (grandtotal) and related outcomes.
# In RStudio you would run View(data); here we preview the first rows.
knitr::kable(head(data), caption = "First rows of the McGrath data")
First rows of the McGrath data
| 0 |
181.71 |
393.52 |
84.20 |
659.43 |
| 0 |
134.35 |
367.57 |
81.15 |
583.07 |
| 1 |
163.55 |
293.35 |
100.05 |
556.95 |
| 0 |
85.65 |
201.92 |
17.05 |
304.62 |
| 1 |
141.90 |
293.21 |
102.71 |
537.82 |
| 0 |
171.00 |
305.03 |
69.19 |
545.22 |
#> treatment book coffee food
#> Min. :0.0000 Min. : 83.9 Min. :201.9 Min. : 17.05
#> 1st Qu.:0.0000 1st Qu.:129.7 1st Qu.:297.8 1st Qu.: 50.12
#> Median :1.0000 Median :165.7 Median :329.1 Median : 78.82
#> Mean :0.5333 Mean :179.6 Mean :331.8 Mean : 83.55
#> 3rd Qu.:1.0000 3rd Qu.:199.1 3rd Qu.:366.3 3rd Qu.:112.62
#> Max. :1.0000 Max. :450.9 Max. :428.4 Max. :180.31
#> grandtotal
#> Min. : 304.6
#> 1st Qu.: 524.8
#> Median : 566.4
#> Mean : 594.9
#> 3rd Qu.: 646.1
#> Max. :1003.1
Comparing Means:
Two-Sample t-test
Goal: estimate the average treatment effect (ATE) on
grandtotal.
Step 1 — Sample means for each group.
mean_treatment <- mean(data$grandtotal[data$treatment == 1])
mean_control <- mean(data$grandtotal[data$treatment == 0])
mean_treatment
#> [1] 612.0019
#> [1] 575.3879
Step 2 — ATE (difference in sample means).
ate <- mean_treatment - mean_control
ate
#> [1] 36.61402
Step 3 — Two-sample t-test (two-sided).
t_test_result <- t.test(grandtotal ~ treatment, data = data)
t_test_result
#>
#> Welch Two Sample t-test
#>
#> data: grandtotal by treatment
#> t = -0.71008, df = 27.992, p-value = 0.4835
#> alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
#> 95 percent confidence interval:
#> -142.23716 69.00913
#> sample estimates:
#> mean in group 0 mean in group 1
#> 575.3879 612.0019
Step 4 — One-sided t-test (alternative:
treatment > control).
t_test_one_sided <- t.test(grandtotal ~ treatment, data = data,
alternative = "greater")
t_test_one_sided
#>
#> Welch Two Sample t-test
#>
#> data: grandtotal by treatment
#> t = -0.71008, df = 27.992, p-value = 0.7582
#> alternative hypothesis: true difference in means between group 0 and group 1 is greater than 0
#> 95 percent confidence interval:
#> -124.3301 Inf
#> sample estimates:
#> mean in group 0 mean in group 1
#> 575.3879 612.0019
Step 5 — Linear-regression version (identical
p-value to the two-sided test).
model <- lm(grandtotal ~ treatment, data = data)
summary(model)
#>
#> Call:
#> lm(formula = grandtotal ~ treatment, data = data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -270.77 -83.19 -26.42 61.56 391.14
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 575.39 38.06 15.119 5.36e-15 ***
#> treatment 36.61 52.11 0.703 0.488
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 142.4 on 28 degrees of freedom
#> Multiple R-squared: 0.01732, Adjusted R-squared: -0.01777
#> F-statistic: 0.4936 on 1 and 28 DF, p-value: 0.4881
Quick Practice
Q1) Recreate the two-sample t-test, but only for
book sales. What do you conclude? Is the effect of the
chocolate scent on book sales significant?
t.test(book ~ treatment, data = data)
#>
#> Welch Two Sample t-test
#>
#> data: book by treatment
#> t = -0.63639, df = 27.884, p-value = 0.5297
#> alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
#> 95 percent confidence interval:
#> -85.94311 45.20579
#> sample estimates:
#> mean in group 0 mean in group 1
#> 168.6907 189.0594
Q2) Run the regression version for book sales.
Compare the coefficient and p-value with your t-test
result.
model_book <- lm(book ~ treatment, data = data)
summary(model_book)
#>
#> Call:
#> lm(formula = book ~ treatment, data = data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -105.16 -55.52 -22.98 15.30 261.89
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 168.69 23.49 7.180 8.16e-08 ***
#> treatment 20.37 32.17 0.633 0.532
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 87.91 on 28 degrees of freedom
#> Multiple R-squared: 0.01412, Adjusted R-squared: -0.02109
#> F-statistic: 0.4009 on 1 and 28 DF, p-value: 0.5318
Q3) Run a one-sided t-test for coffee
sales. Hypothesis: the chocolate scent increases
coffee sales.
t.test(coffee ~ treatment, data = data, alternative = "greater")
#>
#> Welch Two Sample t-test
#>
#> data: coffee by treatment
#> t = 0.26074, df = 27.637, p-value = 0.3981
#> alternative hypothesis: true difference in means between group 0 and group 1 is greater than 0
#> 95 percent confidence interval:
#> -30.45632 Inf
#> sample estimates:
#> mean in group 0 mean in group 1
#> 334.7464 329.2362
Wrap-up
Key takeaways
t.test(y ~ group, data = ...) runs a two-sample
t-test; add alternative = "greater" (or
"less") for a one-sided test.
- The ATE is just the difference in group means; the t-test
tells you whether it is distinguishable from 0.
- Regressing the outcome on the treatment indicator gives the
same estimate and p-value as the two-sample
t-test.
- Always check which direction the difference goes before interpreting
a one-sided test.
LS0tCnRpdGxlOiAiTGFiIDkg4oCUIEluZmVyZW5jZSBmb3IgUHJvcG9ydGlvbnMgJiBUd28tU2FtcGxlIHQtdGVzdCIKc3VidGl0bGU6ICJRdWFudGl0YXRpdmUgUmVhc29uaW5nIMK3IExBQiA0MTMiCmF1dGhvcjogIkluc3RydWN0b3I6IFN1YmluIE5hIgpkYXRlOiAiTm92ZW1iZXIgNSwgMjAyNSIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogZmxhdGx5CiAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICB0b2M6IHRydWUKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKICAgIGRmX3ByaW50OiBwYWdlZAogICAgY29kZV9kb3dubG9hZDogdHJ1ZQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoCiAgZWNobyA9IFRSVUUsIG1lc3NhZ2UgPSBGQUxTRSwgd2FybmluZyA9IEZBTFNFLAogIGZpZy5hbGlnbiA9ICJjZW50ZXIiLCBmaWcud2lkdGggPSA3LCBmaWcuaGVpZ2h0ID0gNC4yLAogIGNvbW1lbnQgPSAiIz4iCikKYGBgCgo+ICoqSG93IHRvIHVzZSB0aGlzIGRvY3VtZW50LioqIFJ1biB0aGUgY29kZSBhbG9uZyB3aXRoIG1lLiBDb21tZW50cyBzdGFydCB3aXRoIGAjYC4KPiBFYWNoIHNlY3Rpb24gbWlycm9ycyB0aGUgbGl2ZSBsYWI7IGNvZGUgY2h1bmtzIHNob3cgdGhlIGNvbW1hbmQgYW5kIGl0cyBvdXRwdXQgdG9nZXRoZXIuCgoqKlRvZGF5J3MgcGxhbioqCgoxLiBDbGVhbiBlbnZpcm9ubWVudCArIHNldHVwCjIuIEltcG9ydCAmIHVuZGVyc3RhbmQgdGhlIGRhdGEgKE1jR3JhdGggZXhwZXJpbWVudCkKMy4gQ29tcGFyaW5nIG1lYW5zOiB0d28tc2FtcGxlICp0Ki10ZXN0CjQuIFF1aWNrIHByYWN0aWNlCgotLS0KCiMgQ2xlYW4gRW52aXJvbm1lbnQgKyBTZXR1cAoKYGBge3IgY2xlYW4tc2V0dXB9CnJtKGxpc3QgPSBscygpKSAgICMgY2xlYW4gZW52aXJvbm1lbnQKCiMgTG9hZCB0aGUgcGFja2FnZXMKbGlicmFyeShleHBlcmltZW50cikgICAjIHByb3ZpZGVzIHRoZSBNY0dyYXRoIGRhdGFzZXQKbGlicmFyeShkcGx5cikKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KHN0YXRzKQpgYGAKCj4gSW4gdGhlIGxpdmUgc2Vzc2lvbiB3ZSBpbnN0YWxsIGBleHBlcmltZW50cmAgd2l0aCBgaW5zdGFsbC5wYWNrYWdlcygiZXhwZXJpbWVudHIiKWAKPiBhbmQgc2V0IHRoZSB3b3JraW5nIGRpcmVjdG9yeSB3aXRoIGBzZXR3ZCguLi4pYC4gSGVyZSB0aGUgZGF0YXNldCBzaGlwcyB3aXRoIHRoZQo+IHBhY2thZ2UsIHNvIG5vIGRvd25sb2FkIG9yIGBzZXR3ZCgpYCBpcyBuZWVkZWQuCgojIEltcG9ydCAmIFVuZGVyc3RhbmQgRGF0YSAoTWNHcmF0aCBFeHBlcmltZW50KQoKYGBge3IgaW1wb3J0fQojIExvYWQgdGhlIE1jR3JhdGggZGF0YXNldCAoY29tZXMgd2l0aCBleHBlcmltZW50cikKZGF0YShtY2dyYXRoKQpkYXRhIDwtIG1jZ3JhdGgKYGBgCgoqKkJhY2tncm91bmQ6IHRoZSBNY0dyYXRoIGV4cGVyaW1lbnQuKiogSW4gMjAxNiwgTWFyeSBDLiBNY0dyYXRoLCBQZXRlciBNLiBBcm9ub3csIGFuZApWaXZpZW4gU2hvdHdlbGwgcHVibGlzaGVkIGEgcmFuZG9taXplZCBjb250cm9sbGVkIHRyaWFsIChSQ1QpIG9uIHRoZSBlZmZlY3Qgb2YgYQoqY2hvY29sYXRlIHNjZW50KiBvbiBib29rc3RvcmUgc2FsZXMuIEEgYm9va3N0b3JlIHdhcyByYW5kb21seSBhc3NpZ25lZCB0byBoYXZlIGVpdGhlcjoKCi0gYHRyZWF0bWVudCA9IDFgIOKGkiBjaG9jb2xhdGUgc2NlbnQgcHJlc2VudCwgb3IKLSBgdHJlYXRtZW50ID0gMGAg4oaSIGNvbnRyb2wgKG5vIHNjZW50KS4KClRoZSBkYXRhc2V0IGluY2x1ZGVzIGZpdmUgdmFyaWFibGVzOgoKfCBWYXJpYWJsZSB8IE1lYW5pbmcgfAp8LS0tfC0tLXwKfCBgdHJlYXRtZW50YCB8IDEgPSB0cmVhdG1lbnQsIDAgPSBjb250cm9sIHwKfCBgYm9va2AgfCBTYWxlcyBvZiBib29rcyB8CnwgYGZvb2RgIHwgU2FsZXMgb2YgZm9vZHMgfAp8IGBjb2ZmZWVgIHwgU2FsZXMgb2YgYnVsayBjb2ZmZWUsIHRlYSwgb3Igc3BpY2VzIHwKfCBgZ3JhbmR0b3RhbGAgfCBUb3RhbCBvZiBib29rLCBjb2ZmZWUsIGFuZCBmb29kIHNhbGVzIHwKCioqR29hbDoqKiBlc3RpbWF0ZSB0aGUgY2F1c2FsIGVmZmVjdCBvZiB0aGUgY2hvY29sYXRlIHNjZW50IG9uIHRvdGFsIHNhbGVzCihgZ3JhbmR0b3RhbGApIGFuZCByZWxhdGVkIG91dGNvbWVzLgoKYGBge3IgaW5zcGVjdH0KIyBJbiBSU3R1ZGlvIHlvdSB3b3VsZCBydW4gVmlldyhkYXRhKTsgaGVyZSB3ZSBwcmV2aWV3IHRoZSBmaXJzdCByb3dzLgprbml0cjo6a2FibGUoaGVhZChkYXRhKSwgY2FwdGlvbiA9ICJGaXJzdCByb3dzIG9mIHRoZSBNY0dyYXRoIGRhdGEiKQoKc3VtbWFyeShkYXRhKQpgYGAKCiMgQ29tcGFyaW5nIE1lYW5zOiBUd28tU2FtcGxlICp0Ki10ZXN0CgoqKkdvYWw6KiogZXN0aW1hdGUgdGhlIGF2ZXJhZ2UgdHJlYXRtZW50IGVmZmVjdCAoQVRFKSBvbiBgZ3JhbmR0b3RhbGAuCgoqKlN0ZXAgMSDigJQgU2FtcGxlIG1lYW5zIGZvciBlYWNoIGdyb3VwLioqCgpgYGB7ciBtZWFuc30KbWVhbl90cmVhdG1lbnQgPC0gbWVhbihkYXRhJGdyYW5kdG90YWxbZGF0YSR0cmVhdG1lbnQgPT0gMV0pCm1lYW5fY29udHJvbCAgIDwtIG1lYW4oZGF0YSRncmFuZHRvdGFsW2RhdGEkdHJlYXRtZW50ID09IDBdKQoKbWVhbl90cmVhdG1lbnQKbWVhbl9jb250cm9sCmBgYAoKKipTdGVwIDIg4oCUIEFURSAoZGlmZmVyZW5jZSBpbiBzYW1wbGUgbWVhbnMpLioqCgpgYGB7ciBhdGV9CmF0ZSA8LSBtZWFuX3RyZWF0bWVudCAtIG1lYW5fY29udHJvbAphdGUKYGBgCgoqKlN0ZXAgMyDigJQgVHdvLXNhbXBsZSAqdCotdGVzdCAodHdvLXNpZGVkKS4qKgoKYGBge3IgdHRlc3QtdHdvfQp0X3Rlc3RfcmVzdWx0IDwtIHQudGVzdChncmFuZHRvdGFsIH4gdHJlYXRtZW50LCBkYXRhID0gZGF0YSkKdF90ZXN0X3Jlc3VsdApgYGAKCioqU3RlcCA0IOKAlCBPbmUtc2lkZWQgKnQqLXRlc3QqKiAoYWx0ZXJuYXRpdmU6IHRyZWF0bWVudCA+IGNvbnRyb2wpLgoKYGBge3IgdHRlc3Qtb25lfQp0X3Rlc3Rfb25lX3NpZGVkIDwtIHQudGVzdChncmFuZHRvdGFsIH4gdHJlYXRtZW50LCBkYXRhID0gZGF0YSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgYWx0ZXJuYXRpdmUgPSAiZ3JlYXRlciIpCnRfdGVzdF9vbmVfc2lkZWQKYGBgCgoqKlN0ZXAgNSDigJQgTGluZWFyLXJlZ3Jlc3Npb24gdmVyc2lvbioqIChpZGVudGljYWwgKnAqLXZhbHVlIHRvIHRoZSB0d28tc2lkZWQgdGVzdCkuCgpgYGB7ciBsbX0KbW9kZWwgPC0gbG0oZ3JhbmR0b3RhbCB+IHRyZWF0bWVudCwgZGF0YSA9IGRhdGEpCnN1bW1hcnkobW9kZWwpCmBgYAoKIyBRdWljayBQcmFjdGljZQoKKipRMSkgUmVjcmVhdGUgdGhlIHR3by1zYW1wbGUgKnQqLXRlc3QsIGJ1dCBvbmx5IGZvciBib29rIHNhbGVzLioqCldoYXQgZG8geW91IGNvbmNsdWRlPyBJcyB0aGUgZWZmZWN0IG9mIHRoZSBjaG9jb2xhdGUgc2NlbnQgb24gYm9vayBzYWxlcyBzaWduaWZpY2FudD8KCmBgYHtyIHExfQp0LnRlc3QoYm9vayB+IHRyZWF0bWVudCwgZGF0YSA9IGRhdGEpCmBgYAoKKipRMikgUnVuIHRoZSByZWdyZXNzaW9uIHZlcnNpb24gZm9yIGJvb2sgc2FsZXMuKioKQ29tcGFyZSB0aGUgY29lZmZpY2llbnQgYW5kICpwKi12YWx1ZSB3aXRoIHlvdXIgKnQqLXRlc3QgcmVzdWx0LgoKYGBge3IgcTJ9Cm1vZGVsX2Jvb2sgPC0gbG0oYm9vayB+IHRyZWF0bWVudCwgZGF0YSA9IGRhdGEpCnN1bW1hcnkobW9kZWxfYm9vaykKYGBgCgoqKlEzKSBSdW4gYSBvbmUtc2lkZWQgKnQqLXRlc3QgZm9yIGNvZmZlZSBzYWxlcy4qKgpIeXBvdGhlc2lzOiB0aGUgY2hvY29sYXRlIHNjZW50ICppbmNyZWFzZXMqIGNvZmZlZSBzYWxlcy4KCmBgYHtyIHEzfQp0LnRlc3QoY29mZmVlIH4gdHJlYXRtZW50LCBkYXRhID0gZGF0YSwgYWx0ZXJuYXRpdmUgPSAiZ3JlYXRlciIpCmBgYAoKIyBXcmFwLXVwCgoqKktleSB0YWtlYXdheXMqKgoKLSBgdC50ZXN0KHkgfiBncm91cCwgZGF0YSA9IC4uLilgIHJ1bnMgYSB0d28tc2FtcGxlICp0Ki10ZXN0OyBhZGQKICBgYWx0ZXJuYXRpdmUgPSAiZ3JlYXRlciJgIChvciBgImxlc3MiYCkgZm9yIGEgb25lLXNpZGVkIHRlc3QuCi0gVGhlIEFURSBpcyBqdXN0IHRoZSBkaWZmZXJlbmNlIGluIGdyb3VwIG1lYW5zOyB0aGUgKnQqLXRlc3QgdGVsbHMgeW91IHdoZXRoZXIgaXQKICBpcyBkaXN0aW5ndWlzaGFibGUgZnJvbSAwLgotIFJlZ3Jlc3NpbmcgdGhlIG91dGNvbWUgb24gdGhlIHRyZWF0bWVudCBpbmRpY2F0b3IgZ2l2ZXMgdGhlICoqc2FtZSoqIGVzdGltYXRlIGFuZAogICpwKi12YWx1ZSBhcyB0aGUgdHdvLXNhbXBsZSAqdCotdGVzdC4KLSBBbHdheXMgY2hlY2sgd2hpY2ggZGlyZWN0aW9uIHRoZSBkaWZmZXJlbmNlIGdvZXMgYmVmb3JlIGludGVycHJldGluZyBhIG9uZS1zaWRlZCB0ZXN0Lgo=