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

  1. Clean environment + setup
  2. Import & understand the data (McGrath experiment)
  3. Comparing means: two-sample t-test
  4. Quick practice

1 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.

2 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:

Variable Meaning
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
treatment book coffee food grandtotal
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
summary(data)
#>    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

3 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
mean_control
#> [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

4 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

5 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=