For this exercise, please try to reproduce the results from Study 1 of the associated paper (Maglio & Polman, 2014). The PDF of the paper is included in the same folder as this Rmd file.

Methods summary:

Researchers recruited 202 volunteers at a subway station in Toronto, Ontario, Canada. Half of the sample was traveling East, while the other half was traveling West. In a 2 (orientation: toward, away from) X 4 (station: Spadina, St. George, Bloor-Yonge, Sherbourne) design, each participant was randomly asked to estimate how far one of the four stations felt to them (1= very close, 7= very far). Authors conducted a 2 X 4 ANOVA on distance estimates, and then tested differences in distance estimates between East and West-bound groups for each individual station.

Target outcomes:

Below is the specific result you will attempt to reproduce (quoted directly from the results section of Study 1):

We carried out a 2 (orientation: toward, away from) × 4 (station: Spadina, St. George, Bloor-Yonge, Sherbourne) analysis of variance (ANOVA) on closeness ratings, which revealed no main effect of orientation, F < 1, and a main effect of station, F(3, 194) = 24.10, p < .001, ηp 2 = .27. This main effect was qualified by the predicted interaction between orientation and station, F(3, 194) = 16.28, p < .001, ηp2 = .20. We decomposed this interaction by the subjective-distance ratings between participants traveling east and west for each of the four subway stations. Westbound participants rated the stations to the west of Bay Street as closer than did eastbound participants; this effect was obtained for both the station one stop to the west (St. George, p < .001, ηp2 = .28) and the station two stops to the west (Spadina, p = .001, ηp2 = .20). The opposite pattern held true for stations to the east of Bay Street. Eastbound participants rated the stations to the east of Bay Street as closer than did westbound participants; this effect was obtained for both the station one stop to the east (Bloor-Yonge, p = .053, ηp2 = .08) and the station two stops to the east (Sherbourne, p < .001, ηp2 = .24). Figure 1 summarizes these results.

Step 1: Load packages

library(tidyverse) # for data munging
library(knitr) # for kable table formating
library(haven) # import and export 'SPSS', 'Stata' and 'SAS' Files
library(readxl) # import excel files
library(lsr)
library(ggthemes)

Step 2: Load data

# Just Study 1
d <- read_excel("data/S1_Subway.xlsx")

Step 3: Tidy data

The data are already tidy as provided by the authors.

head(d)
## # A tibble: 6 × 4
##   DIRECTION DISTANCE STN_NUMBER STN_NAME
##   <chr>        <dbl>      <dbl> <chr>   
## 1 EAST             1          3 B-Y     
## 2 EAST             2          3 B-Y     
## 3 EAST             2          3 B-Y     
## 4 EAST             3          3 B-Y     
## 5 EAST             1          3 B-Y     
## 6 EAST             1          3 B-Y

Step 4: Run analysis

Pre-processing

You will notice that R doesn’t recognize some of these variables as factors. So let’s just factorize our variables before we proceed with the analyses.

d = d %>%
  mutate(
    DIRECTION = as.factor(DIRECTION),
    STN_NUMBER = as.factor(STN_NUMBER),
    STN_NAME = factor(STN_NAME, levels = c("SPAD", "STG", "B-Y", "SHER"))
  )
d
## # A tibble: 202 × 4
##    DIRECTION DISTANCE STN_NUMBER STN_NAME
##    <fct>        <dbl> <fct>      <fct>   
##  1 EAST             1 3          B-Y     
##  2 EAST             2 3          B-Y     
##  3 EAST             2 3          B-Y     
##  4 EAST             3 3          B-Y     
##  5 EAST             1 3          B-Y     
##  6 EAST             1 3          B-Y     
##  7 EAST             2 3          B-Y     
##  8 EAST             1 3          B-Y     
##  9 EAST             3 3          B-Y     
## 10 EAST             2 3          B-Y     
## # ℹ 192 more rows

Inferential statistics

We carried out a 2 (orientation: toward, away from) × 4 (station: Spadina, St. George, Bloor-Yonge, Sherbourne) analysis of variance (ANOVA) on closeness ratings, which revealed no main effect of orientation, F < 1, and a main effect of station, F(3, 194) = 24.10, p < .001, ηp 2 = .27. This main effect was qualified by the predicted interaction between orientation and station, F(3, 194) = 16.28, p < .001, ηp2 = .20.

Note from teaching team: ηp2 = “partial eta squared”

First, check the number of observations for each condition

d %>%
  group_by(DIRECTION, STN_NAME) %>%
  summarize(
    num_obs = n(),
    .groups = "keep"
  ) %>%
  ungroup()
## # A tibble: 8 × 3
##   DIRECTION STN_NAME num_obs
##   <fct>     <fct>      <int>
## 1 EAST      SPAD          26
## 2 EAST      STG           26
## 3 EAST      B-Y           23
## 4 EAST      SHER          26
## 5 WEST      SPAD          25
## 6 WEST      STG           25
## 7 WEST      B-Y           26
## 8 WEST      SHER          25

Looks pretty even to me!

Now run the model

When I see “analysis of variance” in the paper, my first thought is to use the aov function in R. This is a wrapper function that fits a linear model to the data.

Note on statistics: the goal of this assignment is for you to get used to reproducing analyses. It is not intended to teach you about the subtleties of statistics; that will be covered in PSYC201B. We are happy to chat about these details in office hours, but do not worry if you do not understand them right now!

Now fill in the aov command with the proper variables (e.g. replace “VAR1” with the first variable of interest. This is a two-way ANOVA that tests for a main effect of direction (east vs west), a main effect of station, and an interaction between direction and station. The interaction term is denoted by the * below. For more information, do a little Googling about two-way ANOVAs.

mod <- aov(data = d,
           formula = DISTANCE ~ DIRECTION + STN_NAME + DIRECTION*STN_NAME)
anova(mod)
## Analysis of Variance Table
## 
## Response: DISTANCE
##                     Df  Sum Sq Mean Sq F value    Pr(>F)    
## DIRECTION            1   0.713  0.7129  0.6644     0.416    
## STN_NAME             3  75.158 25.0525 23.3492 6.011e-13 ***
## DIRECTION:STN_NAME   3  52.413 17.4710 16.2832 1.765e-09 ***
## Residuals          194 208.152  1.0729                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check your results against the instructional team

Note: the values in the parentheses denote the degrees of freedom for the F-test. * Main effect of direction: F(1, 194) = 0.664 * Main effect of station: F(3, 194) = 23.35 * Interaction between direction and station: F(3,194) = 16.28

Note from instructors: We’re unclear as to why the F-values are slightly different from those in the paper. But students in previous versions of this class came to the same answer.

Effect size

Get the partial eta-squared (under the eta.sq.part column). Do your results match the original paper? Hint: they should!

etaSquared(mod)
##                         eta.sq eta.sq.part
## DIRECTION          0.001195964 0.001929304
## STN_NAME           0.223393543 0.265284110
## DIRECTION:STN_NAME 0.155789423 0.201151614

Station-specific analyses

We decomposed this interaction by the subjective-distance ratings between participants traveling east and west for each of the four subway stations. Westbound participants rated the stations to the west of Bay Street as closer than did eastbound participants; this effect was obtained for both the station one stop to the west (St. George, p < .001, ηp2 = .28) and the station two stops to the west (Spadina, p = .001, ηp2 = .20). The opposite pattern held true for stations to the east of Bay Street. Eastbound participants rated the stations to the east of Bay Street as closer than did westbound participants; this effect was obtained for both the station one stop to the east (Bloor-Yonge, p = .053, ηp2 = .08) and the station two stops to the east (Sherbourne, p < .001, ηp2 = .24). Figure 1 summarizes these results.

St. George

# STEP 1: filter the data to get only St. George estimates
st.george = d %>%
  filter(STN_NAME == "STG")

# STEP 2: was the estimated distance to St. George greater for eastbound or westbound travellers?
# no stats needed, just get the mean estimated distance between east and west directions
# Hint: you'll want group_by() and summarize()
st.george %>% 
  group_by(DIRECTION) %>%
  summarize(mean(DISTANCE))
## # A tibble: 2 × 2
##   DIRECTION `mean(DISTANCE)`
##   <fct>                <dbl>
## 1 EAST                  2.77
## 2 WEST                  1.64
# STEP 3: Now run an aov on St. George data with distance as the response variable and direction as the explanatory variable
st.george.aov = aov(data = st.george, formula = DISTANCE ~ DIRECTION)
summary(st.george.aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## DIRECTION    1  16.25  16.252   18.79 7.23e-05 ***
## Residuals   49  42.38   0.865                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# STEP 4: Find the eta squared
etaSquared(st.george.aov)
##              eta.sq eta.sq.part
## DIRECTION 0.2772092   0.2772092

Spadina

## do the same for Spadina
spadina = d %>%
  filter(STN_NAME == "SPAD")

spadina %>% 
  group_by(DIRECTION) %>%
  summarize(mean(DISTANCE))
## # A tibble: 2 × 2
##   DIRECTION `mean(DISTANCE)`
##   <fct>                <dbl>
## 1 EAST                  3.65
## 2 WEST                  2.64
spadina.aov = aov(data = spadina, formula = DISTANCE ~ DIRECTION)
summary(spadina.aov)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## DIRECTION    1  13.10  13.100   11.97 0.00113 **
## Residuals   49  53.64   1.095                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
etaSquared(spadina.aov)
##              eta.sq eta.sq.part
## DIRECTION 0.1962763   0.1962763

Bloor-Yonge

## do the same for Bloor-Yonge
bloor.yonge = d %>%
  filter(STN_NAME == "B-Y")

bloor.yonge %>% 
  group_by(DIRECTION) %>%
  summarize(mean(DISTANCE))
## # A tibble: 2 × 2
##   DIRECTION `mean(DISTANCE)`
##   <fct>                <dbl>
## 1 EAST                  1.61
## 2 WEST                  2.19
bloor.yonge.aov = aov(data = bloor.yonge, formula = DISTANCE ~ DIRECTION)
summary(bloor.yonge.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## DIRECTION    1   4.16   4.157   3.945 0.0528 .
## Residuals   47  49.52   1.054                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
etaSquared(bloor.yonge.aov)
##              eta.sq eta.sq.part
## DIRECTION 0.0774451   0.0774451

Sherbourne

## do the same for Sherbourne
sherbourne = d %>%
  filter(STN_NAME == "SHER")

sherbourne %>% 
  group_by(DIRECTION) %>%
  summarize(mean(DISTANCE))
## # A tibble: 2 × 2
##   DIRECTION `mean(DISTANCE)`
##   <fct>                <dbl>
## 1 EAST                  2.77
## 2 WEST                  4
sherbourne.aov = aov(data = sherbourne, formula = DISTANCE ~ DIRECTION)
summary(sherbourne.aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## DIRECTION    1  19.31  19.306   15.11 0.000305 ***
## Residuals   49  62.62   1.278                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
etaSquared(sherbourne.aov)
##              eta.sq eta.sq.part
## DIRECTION 0.2356667   0.2356667

Step 5: Reflection

Were you able to reproduce the results you attempted to reproduce? If not, what part(s) were you unable to reproduce?

Yes

How difficult was it to reproduce your results?

A little bit

What aspects made it difficult? What aspects made it easy?

Due to being a little rusty…

Bonus: Reproduce Fig 1

summary_data <- d %>%
  group_by(DIRECTION, STN_NAME) %>%
  summarise(
    Mean_Distance = mean(DISTANCE),
    SE = sd(DISTANCE) / sqrt(n())
  )

# Plot the summarized data
ggplot(summary_data, aes(x = STN_NAME, y = Mean_Distance, group = DIRECTION, color = DIRECTION)) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Mean_Distance - SE, ymax = Mean_Distance + SE), width = 0.2) +
  labs(
    x = "Station Name",
    y = "Mean Distance",
    title = "Mean Distance by Station and Direction",
    color = "Direction"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    legend.title = element_blank()
  )