0 · Libraries

library(tidyverse)     # data wrangling + ggplot
## Warning: package 'tidyverse' was built under R version 4.5.3
## Warning: package 'readr' was built under R version 4.5.3
## Warning: package 'forcats' was built under R version 4.5.3
## Warning: package 'lubridate' was built under R version 4.5.3
library(lmtest)        # bptest() — Breusch–Pagan test for heteroscedasticity
library(effectsize)    # eta_squared() — ANOVA effect size (share of variance)
## Warning: package 'effectsize' was built under R version 4.5.3

1 · Load & explore

menu <- read.csv("restaurant_menu.csv")
str(menu)                       # revenue (num), size (chr), item (chr)
## 'data.frame':    48 obs. of  3 variables:
##  $ revenue: num  3140 4564 4641 4348 3629 ...
##  $ size   : chr  "large" "large" "large" "large" ...
##  $ item   : chr  "A" "A" "A" "A" ...
summary(menu$revenue)           # scale of the outcome
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1851    3021    4038    4843    6272   12449
menu$item <- factor(menu$item)  # make the treatment a factor (levels A,B,C,D)
menu$size <- factor(menu$size)  # make the block a factor (small/mid/large)
# Is the experiment balanced? (equal cell sizes -> a clean randomised block design)
table(menu$item, menu$size)     # 4 in every cell -> BALANCED 4x3, 4 replicates each
##    
##     large mid small
##   A     4   4     4
##   B     4   4     4
##   C     4   4     4
##   D     4   4     4
# Group means: which menu looks best on raw revenue?
aggregate(revenue ~ item, data = menu, FUN = mean)   # B highest (~6817), A lowest (~3190)
item revenue
A 3190.073
B 6816.736
C 3976.671
D 5387.246 
aggregate(revenue ~ size, data = menu, FUN = mean)   # revenue rises with restaurant size
size revenue
large 6223.079
mid 5491.730
small 2813.235 
# Visual first look (course always boxplots before ANOVA).
boxplot(revenue ~ item, data = menu, main = "Revenue by menu item", col = "#e3f0ee")

boxplot(revenue ~ size, data = menu, main = "Revenue by restaurant size", col = "#fbeee0")

2 · Week 1–2 drill · transform decision & OLS diagnostics

# DECISION: should we log revenue? Compare fit on the SAME model.
# CAUTION: AIC is only comparable across models with the SAME outcome scale — you may
# NOT compare AIC(raw) with AIC(log) (different units). Here both share 'revenue', so
# we compare specs; to judge logging, look at residual normality instead.
hist(menu$revenue,      main = "revenue", col = "grey85")

hist(log(menu$revenue), main = "log(revenue)", col = "grey85")   # revenue>0 -> log is safe

# Revenue is only mildly skewed and ANOVA is robust; the course keeps the raw scale here.
# ANOVA and OLS are the same model. lm() shows each menu vs the baseline (item A).
lm_block <- lm(revenue ~ item + size, data = menu)
summary(lm_block)               # itemB, itemC, itemD = each menu vs A, holding size fixed
## 
## Call:
## lm(formula = revenue ~ item + size, data = menu)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2517.8  -961.3  -151.7   625.8  4982.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4570.5      557.6   8.196 2.98e-10 ***
## itemB         3626.7      643.9   5.632 1.34e-06 ***
## itemC          786.6      643.9   1.222  0.22867    
## itemD         2197.2      643.9   3.412  0.00144 ** 
## sizemid       -731.3      557.6  -1.312  0.19681    
## sizesmall    -3409.8      557.6  -6.115 2.71e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1577 on 42 degrees of freedom
## Multiple R-squared:  0.6514, Adjusted R-squared:  0.6099 
## F-statistic:  15.7 on 5 and 42 DF,  p-value: 1.055e-08
bptest(lm_block)                # H0: constant error variance. p<0.05 -> heteroscedastic
## 
##  studentized Breusch-Pagan test
## 
## data:  lm_block
## BP = 10.368, df = 5, p-value = 0.06546
                                # (robust SEs are off-syllabus: NOTE the issue, don't patch)
par(mfrow = c(1, 2)); plot(lm_block, which = c(1, 2)); par(mfrow = c(1, 1))  # residual checks

# Residual normality (course habit: hist of residuals + QQ plot). Inference (F/t/Tukey)
# assumes roughly normal residuals; ANOVA is fairly robust, but check and note it.
hist(lm_block$residuals, main = "Residuals", col = "grey85")

qqnorm(lm_block$residuals); qqline(lm_block$residuals)   # points on the line = ~normal

# Explicit ATE of each menu vs the baseline (item A), with 95% confidence intervals.
# Because assignment is RANDOM, these coefficients ARE the causal ATEs (no confounding).
cbind(ATE = coef(lm_block), confint(lm_block))[c("itemB", "itemC", "itemD"), ]
##             ATE     2.5 %   97.5 %
## itemB 3626.6625 2327.2214 4926.104
## itemC  786.5975 -512.8436 2086.039
## itemD 2197.1725  897.7314 3496.614
# Reading: each row is $/week vs menu A; a CI excluding 0 => that menu differs from A.

3 · Two groups · Welch t-test

# If a question compares only TWO menus, use a Welch t-test (course default:
# var.equal = FALSE). H0: the two menu means are equal.
ab <- subset(menu, item %in% c("A", "B"))
ab$item <- droplevels(ab$item)                       # drop unused levels C, D
t.test(revenue ~ item, data = ab, var.equal = FALSE) # Welch (unequal variances)
## 
##  Welch Two Sample t-test
## 
## data:  revenue by item
## t = -3.7339, df = 13.212, p-value = 0.002438
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
##  -5721.539 -1531.786
## sample estimates:
## mean in group A mean in group B 
##        3190.073        6816.736
summary(lm(revenue ~ item, data = ab))                        # identical ATE via regression
## 
## Call:
## lm(formula = revenue ~ item, data = ab)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3877.5 -1113.1  -140.3  1386.0  5631.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3190.1      686.8   4.645 0.000125 ***
## itemB         3626.7      971.3   3.734 0.001151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2379 on 22 degrees of freedom
## Multiple R-squared:  0.3879, Adjusted R-squared:  0.3601 
## F-statistic: 13.94 on 1 and 22 DF,  p-value: 0.001151
# Reading: the coefficient / mean difference is the ATE of menu B vs menu A.

4 · One-way ANOVA · does the menu matter at all?

# 3+ groups -> one-way ANOVA (omnibus F). H0: all four menu means are equal.
m1 <- aov(revenue ~ item
          # + size
          , data = menu)
summary(m1)          # F = 6.51, p = 0.001  -> REJECT H0: the menus do differ
##             Df    Sum Sq  Mean Sq F value  Pr(>F)    
## item         3  92094349 30698116   6.506 0.00097 ***
## Residuals   44 207607388  4718350                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# This ignores restaurant size — size variance is still sitting in the error term.

5 · Randomised BLOCK ANOVA · control for restaurant size

# Add the block factor 'size'. Blocking pulls size-variance OUT of the residual, so the
# item effect is measured more precisely (its F rises) even though nothing about the
# menus changed. This is the whole point of the block design.
m2 <- aov(revenue ~ item + size, data = menu)
summary(m2)          # item F = 12.34 (p<0.001) ; size F = 20.73 (p<0.001)
##             Df    Sum Sq  Mean Sq F value   Pr(>F)    
## item         3  92094349 30698116   12.34 6.35e-06 ***
## size         2 103126646 51563323   20.73 5.47e-07 ***
## Residuals   42 104480742  2487637                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ^^ item's F roughly DOUBLED (6.51 -> 12.34) purely from blocking.
# LIKELY EXAM QUESTION: "why did the item effect become more significant when you added
# size?" ANSWER: blocking removes between-size variance from the ERROR term (smaller MSE),
# not because the menu effect itself grew. Both factors are significant.

6 · Two-way ANOVA with interaction · does the menu effect depend on size?

# H0: no interaction (the menu effect is the same across restaurant sizes).
m3 <- aov(revenue ~ item * size, data = menu)
summary(m3)          # item:size F = 1.87, p = 0.114  -> DO NOT reject: no interaction
##             Df    Sum Sq  Mean Sq F value   Pr(>F)    
## item         3  92094349 30698116  13.869 3.62e-06 ***
## size         2 103126646 51563323  23.295 3.23e-07 ***
## item:size    6  24795634  4132606   1.867    0.114    
## Residuals   36  79685108  2213475                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# So the additive blocked model m2 is the right one to report.
interaction.plot(menu$item, menu$size, menu$revenue,
                 xlab = "Menu item", ylab = "Mean revenue", trace.label = "Size")

# Roughly parallel lines confirm the non-significant interaction (p = 0.114).

7 · Post-hoc · which menus differ (Tukey HSD)

# ONLY after a significant omnibus F. Tukey controls the family-wise error rate across
# all pairwise comparisons. A confidence interval that excludes 0 => that pair differs.
TukeyHSD(m2, which = "item")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = revenue ~ item + size, data = menu)
## 
## $item
##           diff        lwr        upr     p adj
## B-A  3626.6625  1904.2598  5349.0652 0.0000078
## C-A   786.5975  -935.8052  2509.0002 0.6168447
## D-A  2197.1725   474.7698  3919.5752 0.0075514
## C-B -2840.0650 -4562.4677 -1117.6623 0.0003962
## D-B -1429.4900 -3151.8927   292.9127 0.1344209
## D-C  1410.5750  -311.8277  3132.9777 0.1424874
# Reading: look for pairs whose 'p adj' < 0.05 (e.g. B–A, D–A) — those menus differ.

8 · Effect size · how big, not just how significant (eta²)

# Significance != magnitude. eta^2 = share of variance explained by each factor.
eta_squared(m2)      # partial eta^2 for item and for size
Parameter Eta2_partial CI CI_low CI_high
item 0.4684945 0.95 0.2621760 1
size 0.4967388 0.95 0.3055049 1 
# Reading: report BOTH — the F/p (is there an effect?) and eta^2 (how large?).

9 · Summary of results & reasoning

Headline. Menus differ (block ANOVA: item F = 12.34, p < 0.001); revenue also rises with size (F = 20.73); there is no item×size interaction (p = 0.114). Menu B is best, A worst. Because assignment is random, these are ATEs. Blocking on size didn’t change the effect — it sharpened the test by removing block variance from the error (item F went 6.51 → 12.34).

Numbers were verified in Python against the CSV; knit once in R to confirm on your machine.