Learning R for Data Analysis in Journalism

R uses packages to extend its capabilities — think of them like apps you install once and open whenever you need them. The install step only happens once; the library step happens every session.

install.packages("tidyverse")  # data wrangling
install.packages("dplyr")      # R's version of SQL
install.packages("ggplot2")    # graphics
install.packages("psych")      # descriptive stats
install.packages("tidycensus") # census data API
install.packages("janitor")    # data cleaning
install.packages("modelr")     # regression helpers

library(tidyverse)
library(dplyr)
library(ggplot2)
library(psych)
library(tidycensus)
library(janitor)
library(modelr)

Installing a package is like going to the hardware store to buy a tool. Using library() is getting it out of your toolbox. You don't go back to the store every time you need it.

Everything in R lives in a data frame — a table with rows and columns. Use read_csv() to load CSV files. The <- arrow assigns the result to a name you choose.

ind <- read_csv('globalind.csv', col_names = TRUE)
view(ind)   # opens a spreadsheet-like viewer
str(ind)    # shows structure: column types, row count

The global indicator dataset has 222 rows (countries) and 10 columns including GNI, health expenditure, education spending, unemployment rate, and maternal mortality. The first two columns are character strings; the rest are numeric.

The pipe operator %>% chains commands together — read it as "and then." select() picks columns; filter() picks rows; arrange() sorts.

ind %>%
  select(country, code, gni) %>%
  arrange(desc(gni))          # sort by GNI, highest first

filter(ind, country == "United States")   # exact match
filter(ind, gni > 60000)                  # numeric condition

ind %>%
  select(country, code, gni, skulspend) %>%
  filter(gni >= 50000) %>%
  arrange(desc(skulspend))

Shortcut for %>%: Shift + Cmd + M

Before running any statistics, understand your data's shape. Are the mean and median close? Is there missing data? The psych package's describe() gives you everything at once.

mean(ind$gni, na.rm = TRUE)    # na.rm ignores missing values
median(ind$gni, na.rm = TRUE)

describe(ind)                  # mean, sd, min, max, skew...

If you forget na.rm = TRUE and your data has missing values, R will return NA instead of the actual result. Always add it when working with real-world datasets.

A correlation tells you how strongly two variables move together — from -1 (perfect inverse) to +1 (perfect positive). Use pairwise.complete.obs to handle missing data gracefully.

cor(ind$skul, ind$skulspend, use = "pairwise.complete.obs")

round(cor(ind[3:10], use = "pairwise.complete.obs"), 2)

Visualizing before analyzing is essential — it shows you whether data is skewed, whether outliers exist, and whether a relationship actually looks linear.

hist(ind$gni)
hist(ind$healthexp)

# Positively skewed? Try log transformation to normalize:
hist(log10(ind$gni))

plot(ind$gni, ind$healthexp)   # X = independent, Y = dependent

# Multiple variables at once:
pairs(~ gni + healthexp + skul, data = ind)

On log transformation: Taking log10() of a right-skewed variable can normalize it and improve your model — but it makes results harder to explain to readers. Use it when the math requires it, not by default.

Once you've identified a strong correlation, linear regression quantifies the relationship — how much does Y change for every one-unit increase in X? In R, use lm() (linear model). The formula reads: dependent ~ independent.

options(scipen = 999)   # turns off scientific notation

ind_mod1 <- lm(formula = healthexp ~ gni, data = ind)
summary(ind_mod1)

ind_mod2 <- lm(formula = healthexp ~ gni + unemp, data = ind)
summary(ind_mod2)

The summary() output has a lot going on. Here's what to focus on:

R-squared: Tells you how much of the variation in Y is explained by X. An R² of 0.87 means 87% of the change in health expenditure is explained by GNI. Check the p-value to confirm it's statistically significant (below 0.05).

Coefficients: The coefficient for GNI was 0.107 — for every $1 increase in GNI, health expenditure is expected to rise by about 11 cents.

The coefficient always describes a one-unit change in the independent variable. Multiply by a round number to make the interpretation more intuitive for your audience.

A residual is the difference between the model's predicted value and the actual value. Large residuals are cases that behave unexpectedly — which is often where the most interesting stories are.

schools_pr <- schools %>%
  add_predictions(school_mod1) %>%
  add_residuals(school_mod1)

# Standardize to find true outliers (z-score > 2 or < -2)
schools_pr$z_score <- as.numeric(scale(schools_pr$resid))

plot(schools$poorpct, schools$p_passed,
     main = "Poverty vs. Test Pass Rate",
     xlab = "% in poverty", ylab = "% passed",
     pch = 19)
abline(school_mod1, col = "red")