Exploratory Data Analysis: Covariation

PSY 410: Data Science for Psychology

Dr. Sara Weston

2026-05-04

From variation to covariation

Psychology is about relationships

Last time, you explored how individual variables behave — distributions, outliers, missing data.

But psychology is about relationships:

• Does treatment predict depression?
• Does age relate to reaction time?
• Does anxiety co-occur with insomnia?

Today we learn to see those relationships in data — before testing them statistically.

Connect back to last session: “Last time we looked at one variable at a time. But the questions you actually care about in psychology involve relationships.” Spend ~2 minutes. The three example questions (treatment/depression, age/RT, anxiety/insomnia) are deliberately chosen to span different variable type pairings — point out that we’ll learn a different plot type for each.

Categorical + Continuous

Example: Mental health by treatment group

# Simulated therapy outcome data
therapy_data <- tibble(
  condition = rep(c("Control", "CBT", "Mindfulness"), each = 50),
  depression_post = c(
    rnorm(50, mean = 18, sd = 5),  # Control
    rnorm(50, mean = 12, sd = 5),  # CBT
    rnorm(50, mean = 14, sd = 5)   # Mindfulness
  )
)

Example: Mental health by treatment group

glimpse(therapy_data)

Rows: 150
Columns: 2
$ condition       <chr> "Control", "Control", "Control", "Control", "Control",…
$ depression_post <dbl> 10.322840, 19.810435, 19.234806, 19.058651, 23.269446,…

Walk through the simulated data structure. Emphasize that this is simulated — we chose the means and SDs on purpose so the CBT group looks better. In real research, you wouldn’t know the answer in advance. The dataset is small (150 rows) and simple so students can focus on the visualization, not the data wrangling. Students may ask about BDI-II — it’s the Beck Depression Inventory, a widely used clinical measure.

Boxplots: The classic choice

ggplot(therapy_data,
       aes(x = condition, y = depression_post, fill = condition)) +
  geom_boxplot(alpha = 0.7, show.legend = FALSE) +
  scale_fill_manual(values = c(
    "Control" = "#0072B2", "CBT" = "#E69F00", "Mindfulness" = "#009E73")) +
  labs(
    title = "CBT shows the lowest depression scores",
    x = NULL,
    y = "Depression score (BDI-II)"
  ) +
  theme_minimal(base_size = 14) +
  theme(panel.grid.major.x = element_blank())

Boxplots: The classic choice

What boxplots show

• Line in middle: median
• Box: 25th–75th percentile (IQR)
• Whiskers: 1.5 × IQR
• Dots: outliers

. . .

Great for comparison — but they hide the actual data.

Spend ~2 minutes on the anatomy of boxplots. Students may have seen boxplots in stats but not know the mechanics. Draw attention to the specific components: “The line in the middle is NOT the mean — it’s the median.” This is a common misconception. Ask students: “Based on this boxplot, which treatment looks most effective?” (CBT has the lowest median.)

The problem with boxplots

Boxplots summarize, but they hide important information:

• The actual distribution shape (is it bimodal? skewed?)
• Individual data points (how many observations are there?)
• The raw data (where do specific values fall?)

We can do better.

This slide motivates the next section on raincloud plots. The key point: boxplots hide bimodality, sample size, and the actual distribution shape. You can show the classic “same boxplot, different distributions” example verbally — two very different distributions can produce identical boxplots. If students push back (“but boxplots are fine for homework”), acknowledge that boxplots are quick and useful, but for publications and presentations, we want more information.

Raincloud plots

The modern psych visualization

Raincloud plots combine three elements:

1. Density curve (distribution shape)
2. Boxplot (summary stats)
3. Jittered points (individual data)

They’re increasingly popular in psychology publications because they show everything.

Use the ggrain package. (More details here.)

Spend ~3 minutes on this conceptual intro. Raincloud plots are increasingly expected in psychology journals — if students go on to do research, they’ll encounter these. The name comes from the shape: the half-density on top is the “cloud” and the jittered points below are the “rain.” (Note: a violin is a mirrored density; rainclouds use just one side, so “density” is the more accurate label.) The ggrain package makes them easy; without it, you’d need to layer density + boxplot + jitter manually. Students do NOT need to memorize the ggrain syntax for the assignment — boxplots are fine — but exposure here is valuable.

Building a raincloud

library(ggrain)
therapy_data |>
  ggplot(aes(
    x = condition, 
    y = depression_post, 
    fill = condition, 
    color = condition)) +
  # geom_rain creates all parts of your raincloud
  geom_rain(
    alpha = .6,
    # change just the boxplot part
    boxplot.args = list(color = "black")) + 
  # The mean
  stat_summary(fun = mean, geom = "point", shape = 18, size = 5, color = "black") +
  scale_fill_manual(
    values = c("Control" = "#0072B2", "CBT" = "#E69F00", "Mindfulness" = "#009E73")) +
  scale_color_manual(
    values = c("Control" = "#0072B2", "CBT" = "#E69F00", "Mindfulness" = "#009E73")) +
  labs(
    title = "Post-Treatment Depression by Condition",
    subtitle = "Diamond = mean. Density = distribution. Box = IQR. Dots = individual scores.",
    x = "Treatment condition",
    y = "Depression score (BDI-II)"
  ) +
  coord_flip() +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "none",
    panel.grid = element_blank())

Walk through the code on this slide since the previous slide just showed the output. Point out the key pieces: geom_rain() does the heavy lifting, stat_summary adds the mean diamond, and coord_flip() makes it horizontal. Don’t get bogged down in the color scale code. The important takeaway is that geom_rain() gives you violin + box + points in one function call.

When to use which

Boxplot

• Median, IQR, outliers
• Fast to read
• Best for quick EDA and large datasets where individual points would clutter

Raincloud

• Density + summary + raw data
• Shows shape, sample size, and outliers all at once
• Best for publications and presentations

Quick reference slide — ~30 seconds. Emphasize that boxplots are still perfectly fine for EDA and assignments. Rainclouds are the ideal for final presentations and publications.

Pair coding break

Your turn: Build a raincloud

Using therapy_data, build a raincloud plot of depression scores by condition.

1. Load ggrain
2. Map x = condition, y = depression_post, and fill = condition
3. Add geom_rain()
4. Give it an informative title and labeled axes

Time: 10 minutes

Circulate during this exercise. The main thing students will trip on is forgetting to load ggrain or forgetting the fill aesthetic. If pairs finish early, ask them to add coord_flip() for a horizontal layout, or layer in stat_summary() to add the mean. Common question: “do I need a separate scale_fill_manual?” — no, they can let ggplot pick defaults; we’re not optimizing colors here.

Before we move on

📤 Submit your code on Canvas for participation credit. Paste what you have — it doesn’t need to work perfectly.

Categorical + Categorical

Example: Diagnosis by referral source

# Simulated diagnostic data
diagnosis_data <- tibble(
  referral = sample(c("Self", "Physician", "School"), 200, replace = TRUE),
  diagnosis = sample(c("Depression", "Anxiety", "Both", "Other"),
                     200, replace = TRUE)
)

head(diagnosis_data)

# A tibble: 6 × 2
  referral  diagnosis 
  <chr>     <chr>     
1 School    Both      
2 Physician Depression
3 Physician Depression
4 Physician Anxiety   
5 Physician Other     
6 Self      Both      

Option 1: Stacked bar chart

The next three slides show three options for the same data. Move through them at a moderate pace (~1 minute each). For stacked bars, point out the problem: it’s hard to compare the sizes of non-bottom segments across groups because they don’t share a common baseline. This is a perceptual limitation, not a coding one.

ggplot(diagnosis_data, aes(x = referral, fill = diagnosis)) +
  geom_bar() +
  labs(
    title = "Diagnosis by referral source",
    x = "Referral source",
    y = "Count",
    fill = "Diagnosis"
  ) +
  theme_minimal()

Option 1: Stacked bar chart

Option 2: Side-by-side bars

ggplot(diagnosis_data, aes(x = referral, fill = diagnosis)) +
  geom_bar(position = "dodge") +
  labs(
    title = "Diagnosis by referral source",
    x = "Referral source",
    y = "Count",
    fill = "Diagnosis"
  ) +
  theme_minimal()

Option 2: Side-by-side bars

Option 3: Proportions

Highlight the difference between position = “dodge” (previous slide) and position = “fill” (this slide). “fill” normalizes each bar to 100%, which makes it easy to compare proportions across groups but hides differences in group sizes. Ask: “If one referral source has 10 people and another has 100, this chart makes them look equally important. Is that a problem?” This gets students thinking about when proportions vs. counts are appropriate.

ggplot(diagnosis_data, aes(x = referral, fill = diagnosis)) +
  geom_bar(position = "fill") +
  labs(
    title = "Diagnosis distribution by referral source",
    x = "Referral source",
    y = "Proportion",
    fill = "Diagnosis"
  ) +
  theme_minimal()

Option 3: Proportions

Continuous + Continuous

The classic: Scatterplots

# Simulated reaction time data
rt_data <- tibble(
  age = runif(100, 18, 70),
  reaction_time = 200 + age * 3 + rnorm(100, 0, 40)
)

glimpse(rt_data)

Rows: 100
Columns: 2
$ age           <dbl> 43.53466, 46.08907, 55.75717, 42.43055, 51.37091, 64.599…
$ reaction_time <dbl> 390.1471, 307.7071, 368.6429, 304.2600, 413.9945, 425.64…

Basic scatterplot

Scatterplots should feel familiar from stats classes, but many students haven’t made one in R yet. Walk through the aes() mapping: x goes on the horizontal axis, y goes on the vertical. Ask: “What would you predict about the relationship between age and reaction time?” Most students will guess that older people are slower, which matches the data. Spend ~2 minutes here.

ggplot(rt_data, aes(x = age, y = reaction_time)) +
  geom_point() +
  labs(
    title = "Reaction time by age",
    x = "Age (years)",
    y = "Reaction time (ms)"
  ) +
  theme_minimal()

Basic scatterplot

Add a trend line

geom_smooth(method = “lm”) adds a linear regression line. Students will see lm() again in Session 17 — for now, just explain that it draws the “best-fit” straight line through the data. The gray band is the confidence interval (se = TRUE). Students often ask “what does the gray area mean?” — a simple answer is “where we think the true line probably falls.” Don’t go deep into regression mechanics here.

ggplot(rt_data, aes(x = age, y = reaction_time)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    title = "Reaction time increases with age",
    x = "Age (years)",
    y = "Reaction time (ms)"
  ) +
  theme_minimal()

Add a trend line

The problem: Overplotting

With lots of data, points overlap and hide the true density:

# Lots of data
big_rt_data <- tibble(
  age = runif(5000, 18, 70),
  reaction_time = 200 + age * 3 + rnorm(5000, 0, 40)
)

Overplotting problem demonstrated

ggplot(big_rt_data, aes(x = age, y = reaction_time)) +
  geom_point() +
  labs(title = "Hard to see where the data is dense") +
  theme_minimal()

The plot looks like a solid blob. Ask: “Can you tell where most of the data is?” (No — that’s the problem.) With 5,000 points, geom_point() alone fails. This motivates the next three solutions. In psychology, overplotting is common with large survey datasets or reaction time studies with many trials.

Overplotting problem demonstrated

Solution 1: Transparency (alpha)

ggplot(big_rt_data, aes(x = age, y = reaction_time)) +
  geom_point(alpha = 0.1) +
  labs(title = "Using alpha = 0.1 to show density") +
  theme_minimal()

Alpha is the simplest fix and often the best first attempt. Explain: alpha = 0.1 means each point is 90% transparent, so areas with many overlapping points appear darker. Students should experiment with different alpha values — 0.1 is good for 5000 points, but 0.5 might be better for 200 points. This is the overplotting solution they’ll use most often.

Solution 1: Transparency (alpha)

Solution 2: geom_bin2d()

ggplot(big_rt_data, aes(x = age, y = reaction_time)) +
  geom_bin2d() +
  scale_fill_viridis_c() +
  labs(
    title = "2D bins show density",
    fill = "Count"
  ) +
  theme_minimal()

Solution 2: geom_bin2d()

Correlation coefficient

This is a preview of Session 17 (Correlation and Regression). Keep it brief — ~2 minutes. The goal is just to connect the visual pattern they saw in the scatterplot to a single number. Most students have encountered correlation coefficients in stats, so this should feel familiar. Emphasize the dollar-sign syntax (rt_data$age) as an alternative to piping — cor() expects vectors, not a data frame. Students often ask “is this r significant?” — save that for Session 17.

A single number summary of the linear relationship:

cor(rt_data$age, rt_data$reaction_time)

[1] 0.6894228

Note

• r = 1: perfect positive relationship
• r = 0: no linear relationship
• r = -1: perfect negative relationship

But always plot your data first! (See: Anscombe’s Quartet)

If you have time, briefly describe Anscombe’s Quartet: four datasets with identical means, standard deviations, and correlations that look completely different when plotted. This powerfully reinforces the “always visualize” message. You can also mention the Datasaurus Dozen as a modern, fun version of the same idea.

Patterns and models

What patterns tell us

When you see covariation, ask:

1. Could it be coincidence? (Maybe, especially with small samples)
2. What’s the mechanism? (How are these variables related?)
3. Is there a confound? (Could a third variable explain both?)

Warning

Correlation ≠ Causation

Covariation suggests a relationship, but doesn’t prove one variable causes the other.

Spend ~3 minutes here — this is conceptually the most important slide in the deck. Walk through each question with the age/reaction-time example: (1) Could it be coincidence? With 100 data points and a clear trend, probably not. (2) What’s the mechanism? Cognitive processing speed declines with age. (3) Is there a confound? Maybe older participants have different levels of education, health, or motivation. The correlation vs. causation callout connects directly to what they learned in stats — reinforce it here.

Your class data: Covariation

Stress by chronotype

ggplot(class_survey, aes(x = chronotype, y = stress, fill = chronotype)) +
  geom_boxplot(alpha = 0.7, show.legend = FALSE) +
  labs(
    title = "Are Night Owls More Stressed?",
    x = "Chronotype",
    y = "Current stress level (1–10)"
  ) +
  theme_minimal(base_size = 14)

This is a categorical + continuous example from your own class. Ask: “Do night owls appear more stressed, or is it about the same?” The answer will vary by class — that’s the fun of using real data. Connect to the earlier lesson: a boxplot shows you medians and spread, but it hides the distribution shape.

Stress by chronotype

Coding excitement vs. coding anxiety

ggplot(class_survey, aes(x = coding_excited, y = coding_anxious)) +
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.6, size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "steelblue") +
  labs(
    title = "Does Excitement Buffer Anxiety About Coding?",
    x = "Excitement about learning to code (1–10)",
    y = "Anxiety about learning to code (1–10)"
  ) +
  theme_minimal(base_size = 14)

This is a continuous + continuous example with a psychologically interesting question. Use geom_jitter because many students gave the same integer responses, creating overplotting. Ask: “What do you notice? Can you be both excited and anxious?” (Yes — ambivalence is real.) The trend line gives a preview of correlation, which they’ll formalize in Session 17.

Coding excitement vs. coding anxiety

Caffeine and sleep

ggplot(class_survey, aes(x = caffeine_per_day, y = sleep_hrs)) +
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.6, size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "steelblue") +
  labs(
    title = "More Caffeine, Less Sleep?",
    x = "Caffeinated drinks per day",
    y = "Average hours of sleep per night"
  ) +
  theme_minimal(base_size = 14)

Another continuous + continuous example that students find relatable. The direction and strength of this relationship will vary by class — that’s fine. Ask: “Does this surprise you? What’s a possible confound?” (Maybe people who sleep less drink more caffeine to compensate, rather than caffeine causing poor sleep.) This is a natural correlation-vs-causation moment.

Caffeine and sleep

Wrapping up

Key takeaways

1. Covariation = relationships between variables
2. Different plot types for different variable combinations:
   • Categorical + continuous: boxplot, raincloud
   • Categorical + categorical: stacked, dodged, or proportional bars
   • Continuous + continuous: scatterplot (with alpha or 2D bins for big data)
3. Rainclouds show distribution, summary stats, and raw data all at once
4. Watch for overplotting — use alpha, jitter, or binning
5. Always visualize first before computing correlations
6. Patterns suggest but don’t prove causation

Quick summary — ~1 minute. If you’re short on time, let students read this slide themselves. The most important takeaway is point 3 (raincloud plots) and point 6 (correlation is not causation). Remind them that Assignment 5 will ask them to use several of these plot types on a real dataset.

Before next class

📖 Read:

• R4DS Ch 12: Logical vectors
• R4DS Ch 13: Numbers

✅ Do:

• Get started on Assignment 5 (due Sun May 10)
• Start thinking about how you’ll compute scale scores — that’s where we’re headed Wednesday

The one thing to remember

Relationships hide in data. Your job is to make them visible — carefully, honestly.

See you Wednesday!

End on this closing line. Remind students: Assignment 5 (EDA on the BFI dataset) was assigned today and is due Sunday May 10. Next session shifts from visualization to data types — logical vectors, numbers, and computing scale scores (like averaging personality items into composite scores). That connects directly to A5, which uses the BFI data.

Get a head start

Assignment 5 preview

Assignment 5 asks you to do EDA on the BFI personality dataset. Two of the tasks line up directly with what we did today:

• Task 2.2: Scatterplot of two Extraversion items with a trend line
• Task 2.3: Boxplots of one Extraversion item by education level

If you have time tonight, take a first pass at either one — the code patterns are the same as today’s class survey examples.

This isn’t an in-class exercise — it’s an “if you have time” head start on A5. Students who want momentum can attempt Task 2.2 (a continuous-by-continuous scatterplot, just like today’s coding-excitement-vs-anxiety plot) or Task 2.3 (a categorical-by-continuous boxplot, just like today’s stress-by-chronotype plot). Remind them they’ll need library(psych) and data(bfi) to load the dataset.