Lesson 3: Exploratory Data Analysis & Data Cleaning
EOH 635 — Environmental Epidemiology Methods
In our last session (Lesson 1), we read in and manipulated our data set, and merged in NDVI (greenness) data by census tract. This week’s module is about research questions, conceptual frameworks, and descriptive analysis — you’ve already worked through hypothesis components (exposure, outcome, population, expected relationship) and the mediator-vs-modifier distinction in this week’s “Learn” material. This tutorial is the Practice half: hands-on descriptive statistics, data cleaning, and visualization in R.
That’s not busywork. Every study population, exposure, and outcome needs to be characterized before it can be modeled — and this week’s skills are exactly what you’ll need for Project Milestone 1 (research question, hypothesis, and conceptual framework), due this week. Milestone 1 asks you to start characterizing the exposure and outcome you’re considering for your own project, using this same class dataset. Everything below — checking distributions, spotting errors, visualizing a relationship — is directly reusable for that.
The gray boxes below are live R editors — the code actually runs in your browser (no install needed). In the top-left corner of each gray box, you’ll see a green arrow button labeled “Run Code.” Click that green arrow to run the code in that box. The output (whatever the code produces — a number, a table, a plot, etc.) will appear directly underneath the box a moment later.
You can also edit the code first — click inside the box, change something, then click the green Run Code arrow to run your edited version. Nothing you do here can break anything; if you make a mess, just refresh the page to start that box over.
Sections marked 🧠 Your Turn give you a task and a blank (or partially blank) code editor. Try writing the code yourself before revealing the solution, hidden in a collapsed box underneath. Click it open only after a real attempt (or if you’re truly stuck) — the point is practice, not just reading the answer.
A note on browser vs. RStudio code: one chunk in this lesson runs automatically in the background (labeled #| autorun: true) — it fetches the class data file from the web so it’s available in your browser. When you write your own scripts in RStudio, you will not include this fetching code. Immediately after it, you’ll find a green box labeled “In your RStudio script” showing the simpler code you’ll actually use.
This lesson picks up where Lesson 1 left off: reading in class_data_processed.csv, the merged file you (or this tutorial) saved to data_processed/ at the end of Lesson 1. If you worked through Lesson 1 in RStudio and saved your own version there, use that file — it should match. If you skipped straight to this lesson, that’s fine too; a copy is fetched for you below. See the Data Codebook for variable definitions.
As a reminder from Lesson 1: the web-fetch step below is a browser-only workaround. On your own computer, class_data_processed.csv is just a file sitting in your data_processed/ folder — you read it directly, no downloading involved.
Just as in Lesson 1, we first load any packages we need, then read in our data — using file.path() and the same data_raw / data_processed structure from your EOH635_project, not setwd(). (Recall from Lesson 1: on your own computer you’re responsible for install.packages() and library() for every package you use — this tutorial page loads tidyverse for you automatically, which your own script won’t do.)
Let’s read our data into R
As always, make sure the data you’re reading in are where your code expects to find them.
Note that we’re reading in class_data_processed.csv, which already has the NDVI data linked in from Lesson 1.
In RStudio, that’s:
data <- read.csv(file.path("data_processed", "class_data_processed.csv"))In the browser, the chunk below fetches that file for you and puts it in a data_processed/ folder here in webR’s virtual filesystem, matching the real project structure:
Data check
To check that our data loaded correctly, we’ll use the head() function, which displays the first few rows of the dataset.
What are we looking for here? Mainly, that this loaded the way we expected: we should see the same 38 columns from Lesson 1’s class_data_raw.csv, plus the new X2015_NDVI column merged in at the end. If a column looks missing, duplicated, or full of unexpected values already at this stage, that’s worth chasing down before doing anything else — head() is a fast first sanity check, not a substitute for the more careful checks that follow.
Why Explore Data Before Analyzing It?
It’s tempting to skip straight to the “real” analysis — but exploratory data analysis (EDA) isn’t a formality before the real work starts. It is part of the real work. Before you can trust any regression result later in the course, you need to know your data cold: what each variable looks like, where the errors are, and what a “normal” value even is for each one.
This class dataset is a simulated example of something very real in environmental epidemiology: data assembled from multiple sources, each with its own error patterns. Here, individual-level answers came from a household survey, and neighborhood-level measures (tree canopy, traffic, poverty, crime) were linked in by census tract — the same kind of linkage you did yourself in Lesson 1. Every data-collection method has its own way of introducing error:
- Surveys / interviewer-administered questionnaires — recall bias, misunderstood questions, data entry mistakes.
- Monitors and administrative/electronic data — instrument drift, missing readings, records that were never meant for research use.
- Linked data (like the tract-level variables here) — mismatched geographic identifiers, or a value that’s technically correct but describes the neighborhood, not the person.
Broadly, the errors EDA is designed to catch fall into three buckets:
- Missing data
- Logical inconsistencies — combinations of values that can’t both be true
- Outliers — values so extreme they’re more likely to be errors than reality
We’ll work through each of these in turn.
One more framing note before we start: everything in this lesson is descriptive, not inferential. We’re not running hypothesis tests or looking at p-values here — we’re looking at distributions to understand the data and catch problems before any modeling starts. That distinction matters enough that it’s worth stating explicitly: descriptive analysis characterizes what’s in front of you; inferential analysis (coming later in the course) is what lets you draw conclusions and test hypotheses from it.
Checking for Missing Values
Missing values can affect analysis, so it’s essential to check for them.
Questions to consider:
- What variables have missing values?
- How many missing values do you see for each?
- Do you think this is a problem? Why or why not?
The Data Codebook lists the missingness that was deliberately built into this simulated dataset: mental_health (42), current_smoker (20), high_blood_pressure (12), stress (7), has_park (5), age (3), gender (2). If your colSums(is.na(data)) output matches that, your merge from Lesson 1 didn’t introduce (or lose) any missingness — which is exactly the kind of check worth doing after any merge.
A caveat specific to this class exercise: we can check against the codebook here only because this is a teaching dataset, built with a known amount of missingness baked in on purpose. With real-world data, you typically don’t know in advance how much missingness you have, or where — you find out by running exactly the code above and looking at the result. A codebook for a real dataset will usually tell you what a variable means and how it’s coded, but not necessarily how much of it is missing in the specific file you were handed.
For now we will not recode any of our missing variables. In some cases it makes sense to recode missing variables — particularly if we can infer their values from other variables in the data frame — but for now we’ll proceed without doing that. It’s good to be aware of which variables have missing values, since we’ll need to think about this later when conducting analyses.
Since you’ll be working with these data in your project analysis, be sure to take note of the variables that have missing values.
Checking for Logical Inconsistencies
Missing data is easy to spot — R shows you NA directly. Logical inconsistencies are sneakier. When there’s a logical inconsistency, every individual value can look perfectly plausible on its own — it’s only when you look at two variables together that you find a combination that can’t both be true. A classic example from environmental health survey data: someone recorded as biological sex “male” who also reports an age at first menstruation. Neither value looks wrong in isolation — the problem only shows up when you look at both at once.
The tool for catching this is one you already know: table(), applied to two variables instead of one, to see every combination that actually occurs in the data.
As a quick illustration of the technique (not a real inconsistency — just showing the mechanics), let’s break age into two rough groups and cross it with gender. We’ll use the same approach Lesson 1 used to create r_age and num_education — ifelse(), checking a condition and returning one value if it’s TRUE, another if it’s FALSE:
What’s happening here: ifelse(data$age < 40, "Under 40", "40 or older") checks every value of data$age — wherever it’s less than 40, it returns the text "Under 40"; everywhere else, it returns "40 or older". That gives us a new vector of labels the same length as data$age, which we store as age_group_demo. table() then cross-tabulates that new label against gender, one count per combination.
Reading a two-way table like this: each cell is a count of people with that combination of both variables. If a cell that should be logically impossible has a nonzero count, that’s your signal to investigate.
Let’s build a real check, not just an illustration. Think about education: is it possible for someone to have a Bachelor’s degree at age 20? Sure — early admission, an accelerated program, it happens, though it’s uncommon. What about age 19? Now it’s genuinely rare — not a strict logical impossibility the way “male” + “reports menses” is, but rare enough that finding several such cases in one dataset should make you suspicious of a coding error, not assume you’ve found a bunch of prodigies.
Using the same ifelse() pattern as above, build a two-way table crossing an “under 19 / 19 and older” age split against education.
The “Under 19” row should show 0 people with “Bachelor’s degree or higher.” Unlike the strict logical-impossibility example earlier (menses age recorded for someone coded male), this isn’t a constraint that’s guaranteed to hold in every real dataset — it’s a “should be at or near zero; investigate anything that isn’t” check. That distinction matters: not every consistency check has a clean yes/no answer, and part of the skill is judging how surprising a violation would actually be.
Summary Statistics
The summary() function provides a quick look at key statistics (e.g., min, max, median, mean) for the variables in our data frame.
This is a quick review from Lesson 1 — you ran summary(data) there too, on the raw, pre-merge data. It’s worth doing again here, for two reasons: it’s good practice to re-run basic checks like this after any merge (to make sure nothing shifted unexpectedly), and data now has one additional column, X2015_NDVI, that you haven’t summarized yet.
Question to consider: What descriptive statistics do we see for each of the variables? For what types of variables does the summary() function not provide good information? (If this feels familiar, it should — see if your answer matches what you found in Lesson 1.)
Now let’s redo this, looking only at the numeric variables.
Note: is.double() selects only columns stored as double (continuous numeric variables). This avoids including integer columns (like id and year), which are used for identifiers or years rather than quantities we’d summarize with a mean. Use is.numeric() instead if you want to include both integer and double types — is.double() is appropriate here because we want to calculate means, medians, etc. on genuinely continuous measures. It doesn’t make sense to compute a mean or median for a binary yes/no variable — for those, we need a different tool, which is exactly what the next section covers.
Questions to consider:
- What is the median age in this population? What is the minimum age? What is the maximum age?
- For what variables do we not see much difference in value across the population? Why do you think this is?
- What can you tell about percent living in poverty in participants’ neighborhoods, based on the distribution you see here?
Now let’s create some tables to start looking at distributions for the other variables in the data set.
Note about the code above: this selects all columns that are not double-precision numeric (is.double), and also removes id and year. “Double precision” refers to how continuous numbers are stored in R — they use more memory to allow for decimal values with high precision. The minus sign (-) is used inside select() to exclude columns. We’re storing the non-numeric columns in a data frame called non_numeric_data.
Frequency Tables
Now let’s create frequency tables to look at other variables:
Use the same approach (table(), and if you like, prop.table()) to look at the distribution of race_ethnicity in this newly merged dataset. Do the proportions look the same as what you found for race_ethnicity back in Lesson 1? They should — merging in NDVI shouldn’t have changed any of the original columns. If your numbers don’t match, that’s worth investigating!
We can also use a function to run this for all of the non-numeric variables at once:
You’re not expected to write a function like this on your own yet. Save this code (consider building yourself a small code library of “recipes” like this one) and reuse it whenever you need frequency tables for every categorical variable in a data set at once.
Take a look at asthma in the output above. Its prevalence should look like a plausible proportion of the population — a small minority, not nearly everyone. That’s worth pausing on: a variable whose distribution looks implausible (way too common, way too rare, or suspiciously uniform) is itself a signal worth chasing down, the same way a missing-value count or a failed logical-consistency check is.
Visual Exploration of Data
Tables and summary statistics are useful, but it can also be very helpful to visually inspect the distribution of a variable — some patterns (an outlier, an unexpected cluster, a skew) are far easier to spot in a plot than in a column of numbers. Let’s plot some key variables to understand their distributions and relationships.
Example 1: Histogram
Histograms are useful for visualizing the distribution of a single numeric variable. Let’s look at age.
How to read this: the x-axis is age, split into 5-year bins (that’s what binwidth = 5 does); the y-axis is a count — how many participants fall into each bin. Taller bars mean more people at that age range. From this single plot, you can read off roughly where the minimum and maximum ages are (the leftmost and rightmost bars with any height at all), and see at a glance where the bulk of participants’ ages fall — which, here, looks like it’s concentrated in the younger-to-middle adult range, with progressively fewer people as age increases. You should also be able to spot a bar (or bars) sitting well off to the right, separated from the rest of the distribution — that’s a visual signal of an age outlier, the same implausible values we’re about to deal with directly below.
Before we get to age’s outlier problem below, get a look at our primary exposure variable’s distribution. Build a histogram of hrlctree (tree canopy cover) — start by copying the age histogram code above and just swapping in hrlctree.
One thing you’ll need to change: binwidth. A binwidth of 5 made sense for age, which spans roughly 18 to 116 years — but hrlctree is a proportion, ranging only from 0 to 1. A bin 5 units wide would swallow the entire variable into a single bar, telling you nothing. Try a much smaller value, like 0.05, instead.
More generally, whenever you build a histogram of a new variable: first check its rough range (a quick summary() or range() will tell you), then pick a binwidth that divides that range into something like 10–30 bins. Too few bins hides the shape of the distribution; too many makes it look noisy. There’s no single “correct” binwidth — this is a judgment call you refine by trying a value, looking at the plot, and adjusting.
Most tracts should cluster at fairly low canopy values with a long right tail toward a handful of much greener tracts — worth keeping in mind before you characterize this variable for your own project.
Okay, back to age. Some of the values we see for age look suspicious — the ages > 110 are very suspect. This is the third error type from earlier: an outlier. Let’s recode those with age > 105 to be missing. (Another option would be to recode these to some high but plausible number, e.g. 100.)
There’s a trade-off here: we gain more missing values, or we introduce more uncertainty. This is a balancing act you’ll have to consider as you move forward with your own data analysis. Like most things, there isn’t necessarily one “right” answer — just be thoughtful, and document your choice.
Our new recoded age variable will be called r_age:
We can check the number of missing values in the recoded and original age variables:
That is, we now have four additional missing observations, because there were four individuals with ages > 105 — we have set these to missing. (Original missing age count: 3. After recoding: 7. 3 + 4 = 7 — the recode didn’t touch anyone who wasn’t already missing or above the 105 cutoff.)
Recall from Lesson 1 that dplyr’s if_else() (typically used inside mutate()) is another way to do the same recoding ifelse() just did. Recreate r_age using mutate() and if_else() instead, calling the new column r_age2, and confirm it matches r_age (hint: a table of r_age == r_age2, or identical(), can confirm this — remember that NA == NA is NA, not TRUE, so a simple table() of the comparison is often clearer than identical() here).
(Note: if_else() is stricter than base R’s ifelse() about types — it needs NA_real_ rather than a plain NA here, since age is a decimal/“double” column. This is a common if_else() gotcha worth remembering.)
Let’s recreate our histogram with the recoded age variable:
Our histogram no longer has the age outliers we observed earlier.
Example 2: Boxplot
Boxplots are helpful for identifying outliers and comparing distributions across categories.
Let’s look at the distribution of age by gender.
Question to consider: What is the relationship between gender identity and age?
Let’s look at high blood pressure, our primary outcome/dependent variable, by tree canopy cover (our primary exposure/independent variable):
Note the as.character() wrapped around high_blood_pressure in the code above. Why do we need it? (Hint: try running the code without that wrapper. Hint 2: recall from the Data Codebook how high_blood_pressure is coded/stored.)
high_blood_pressure is stored as an integer 0/1 flag, not a category. Without as.character(), ggplot treats it as a continuous numeric axis and tries to draw a single boxplot (or a strange, mostly-empty range of “groups” between 0 and 1) rather than two separate boxes for “0” and “1”. Wrapping it in as.character() (or factor()) tells ggplot to treat it as a categorical variable, so it draws one box per group as intended. This is the same integer-flags-that-are-really-categories issue you first ran into with str(data) back in Lesson 1.
Question to consider: Based on what you see here, what is the relationship between high blood pressure and % tree canopy cover?
Remember this is a crude estimate of the relationship — and, per Week 1, an association (however real it looks here) is not the same thing as causation. Even if this boxplot shows a clear pattern, several other explanations besides “tree canopy directly affects blood pressure” are still on the table, including confounding by some third variable related to both. We’ll return to formal tools for evaluating that (regression models, and later in the course, confounding assessment with causal diagrams) — for now, just describe what you see.
Let’s also look at tree canopy cover by gender:
Question to consider: Based on what you see here, what is the relationship between gender identity and % tree canopy cover? Why might this matter within the context of our study question?
Let’s also look at tree canopy cover by education:
Question to consider: Based on what you see here, what is the relationship between personal educational attainment and % tree canopy cover? Why might this matter within the context of our study question?
Keep these last two plots in mind. You just checked whether two different variables (gender, education) are related to our exposure of interest. That’s exactly the first step in thinking about confounding — a topic the course returns to formally in Week 5. A variable can only confound a relationship if it’s associated with both the exposure and the outcome; here you’ve just checked the “associated with the exposure” half for two candidates. Not every variable you check will turn out to matter — and that’s useful information too.
Let’s check one more candidate variable, this time against a different outcome. Build a boxplot of hrlctree by mental_health status (remember: you’ll need as.character() or factor() around mental_health, for the same reason as high_blood_pressure above).
Example 3: Scatter Plot
Scatter plots are useful for exploring relationships between two numeric variables. Let’s look at the relationship between percentage of the census tract living in poverty and tree canopy cover:
Question to consider: Based on what you see, what do you conclude about the overall relationship between living in a tract with higher % poverty and % tree canopy? Stated another way, what is the general shape and direction of this relationship?
Scatter plots are also a good sanity check when a dataset has two variables that should be related. This dataset has two different greenness measures: hrlctree (tree canopy cover, linked in during data creation) and X2015_NDVI (the NDVI value you yourself merged in during Lesson 1). Build a scatter plot of hrlctree (x-axis) against X2015_NDVI (y-axis). Do they move together the way you’d expect two measures of the same underlying thing to?
A clear positive relationship here is a reassuring sign — two independently sourced measures of “how green is this place” agreeing with each other is exactly the kind of cross-check worth doing before you build an argument around either one.
Example 4: Summary Table of Key Metrics
We can also create summary tables that group data by a categorical variable. Suppose we want to summarize average age by gender:
Question to consider: Compare the results from this summary table to the boxplot you made earlier (age stratified by gender). Do you find one more useful than the other? Why?
Build a similar summary table, but (a) group by education instead of gender, (b) summarize hrlctree instead of age, and (c) add the median alongside the mean, max, and min.
======================================================
🎉 Next Steps
======================================================
You now have a working toolkit for characterizing a dataset: checking for missing values and logical inconsistencies, summarizing distributions, building frequency tables, and visualizing relationships with histograms, boxplots, scatter plots, and grouped summary tables. That’s exactly what Project Milestone 1 (due this week) asks you to start doing for your own chosen exposure and outcome — pick your pair from this dataset, and run them through the same steps you just practiced here.
Next week, we turn to exposure assessment and misclassification — how exposure gets measured and coded, and what happens when that measurement isn’t perfect. Before then, take a closer look at the environmental exposure variables in the Data Codebook (tree canopy, traffic density, crime, air/greenness measures) as you think about which one you want to build your project around.
======================================================
🚀 Quick-Start Script for Your Own Project
======================================================
Everything above walked through each step one at a time, with the browser-only fetching code mixed in. Below is a single reusable script — this is what you’d paste into a new .R or .qmd file inside EOH635_project/scripts/, running on your own computer with class_data_processed.csv already sitting in data_processed/ (from Lesson 1). Swap in your own variable names wherever you see # <- CHANGE THIS, and use it as a starting checklist for Milestone 1.
# ---- Setup ----
library(tidyverse)
# ---- Read in your data ----
data <- read.csv(file.path("data_processed", "class_data_processed.csv"))
# ---- 1. Missing data: check every column at once ----
colSums(is.na(data))
# ---- 2. Logical inconsistencies: cross-tabulate variables that should agree ----
# Example template -- replace both variables with a pair relevant to your own project
age_group <- ifelse(data$age < 19, "Under 19", "19 or older") # <- CHANGE THIS to your own condition
table(age_group, data$education) # <- CHANGE THIS to your own pair
# ---- 3. Outliers + summary statistics ----
summary(data) # all variables
summary(data %>% select(where(is.double))) # numeric variables only
# ---- 4. Frequency tables for categorical variables ----
table(data$education) # <- CHANGE THIS to your own categorical variables
prop.table(table(data$education)) * 100
# Every non-numeric variable at once:
non_numeric_data <- data %>% select(-where(is.double), -id, -year)
lapply(non_numeric_data, function(x) {
if (is.factor(x) || is.character(x) || is.logical(x) || is.integer(x)) {
counts <- table(x, useNA = "ifany")
percentages <- prop.table(counts) * 100
cbind(Count = counts, Percentage = percentages)
} else {
NULL
}
})
# ---- 5. Recode an implausible value once you've found one ----
# Template, following the age > 105 example above:
data$r_age <- ifelse(data$age > 105, NA, data$age) # <- CHANGE THIS to your own variable/cutoff
# ---- 6. Visualize a single continuous variable ----
ggplot(data, aes(x = age)) + # <- CHANGE THIS to your exposure or outcome
geom_histogram(binwidth = 5, fill = "blue", color = "black") +
labs(title = "Distribution of Age", x = "Age", y = "Frequency")
# ---- 7. Visualize your exposure-outcome relationship ----
# If your outcome is a 0/1 flag, wrap it in as.character() so ggplot treats it as a category:
ggplot(data, aes(x = as.character(high_blood_pressure), y = hrlctree)) + # <- CHANGE THIS
geom_boxplot(fill = "orange") +
labs(title = "Tree canopy cover by high blood pressure status",
x = "High blood pressure", y = "Tree canopy cover")
# If both your exposure and outcome are continuous, use a scatter plot instead:
ggplot(data, aes(x = pov_allage, y = hrlctree)) + # <- CHANGE THIS
geom_point(color = "purple") +
labs(title = "% living in poverty vs. tree canopy cover",
x = "% living in poverty", y = "% tree canopy")
# ---- 8. Check a candidate confounder against your exposure ----
ggplot(data, aes(x = education, y = hrlctree)) + # <- CHANGE THIS
geom_boxplot(fill = "orange") +
labs(title = "Tree canopy cover by education", x = "Education", y = "Tree canopy cover")
# ---- 9. Grouped summary table ----
data %>%
group_by(education) %>% # <- CHANGE THIS
summarize(
avg = mean(hrlctree, na.rm = TRUE), # <- CHANGE THIS
median = median(hrlctree, na.rm = TRUE),
max = max(hrlctree, na.rm = TRUE),
min = min(hrlctree, na.rm = TRUE)
)