Data Fundamentals

* Load Stata dataset
use "data.dta", clear

* Load CSV file
import delimited "data.csv", clear

* Load Excel file
import excel "data.xlsx", sheet("Sheet1") firstrow clear

pacman::p_load(haven, readr, readxl)

# Load Stata dataset
data <- read_dta("data.dta")

# Load CSV file
data <- read_csv("data.csv")

# Load Excel file
data <- read_excel("data.xlsx", sheet = "Sheet1")

import pandas as pd

# Load Stata dataset
data = pd.read_stata("data.dta")

# Load CSV file
data = pd.read_csv("data.csv")

# Load Excel file
data = pd.read_excel("data.xlsx", sheet_name="Sheet1")

* Load only specific variables
use income age education using "data.dta", clear

* Load only observations meeting a condition
use "data.dta" if year == 2020, clear

* Combine both
use income age using "data.dta" if state == "CA", clear

# Load specific columns from CSV
data <- read_csv("data.csv", col_select = c(income, age, education))

# Filter while reading (with data.table for large files)
pacman::p_load(data.table)
data <- fread("data.csv", select = c("income", "age"))[year == 2020]

# Load specific columns
data = pd.read_csv("data.csv", usecols=["income", "age", "education"])

# Load and filter (for very large files, use chunks)
data = pd.read_stata("data.dta", columns=["income", "age"])
data = data[data["year"] == 2020]

* Save as Stata dataset
save "cleaned_data.dta", replace

* Export to CSV
export delimited "cleaned_data.csv", replace

# Save as Stata dataset
write_dta(data, "cleaned_data.dta")

# Export to CSV
write_csv(data, "cleaned_data.csv")

# Save as Stata dataset
data.to_stata("cleaned_data.dta", write_index=False)

# Export to CSV
data.to_csv("cleaned_data.csv", index=False)

* Example: Download FARS data (fatal crash data from NHTSA)
* Each year's data is a separate zip file

local years "2000 2001 2002 2003 2004 2005"

foreach year of local years {
    * Download the zip file
    copy "https://static.nhtsa.gov/nhtsa/downloads/FARS/`year'/National/FARS`year'NationalCSV.zip" ///
         "data/raw/FARS`year'.zip", replace

    * Unzip it
    unzipfile "data/raw/FARS`year'.zip", replace
}

* Now you have reproducible data acquisition!

# Example: Download FARS data (fatal crash data from NHTSA)
pacman::p_load(tidyverse)

years <- 2000:2005

for (year in years) {
  url <- glue::glue(
    "https://static.nhtsa.gov/nhtsa/downloads/FARS/{year}/National/FARS{year}NationalCSV.zip"
  )
  destfile <- glue::glue("data/raw/FARS{year}.zip")

  # Download
  download.file(url, destfile, mode = "wb")

  # Unzip
  unzip(destfile, exdir = glue::glue("data/raw/FARS{year}"))
}

# Now you have reproducible data acquisition!

# Example: Download FARS data (fatal crash data from NHTSA)
import requests
import zipfile
from pathlib import Path

years = range(2000, 2006)

for year in years:
    url = f"https://static.nhtsa.gov/nhtsa/downloads/FARS/{year}/National/FARS{year}NationalCSV.zip"
    zip_path = Path(f"data/raw/FARS{year}.zip")

    # Download
    response = requests.get(url)
    zip_path.write_bytes(response.content)

    # Unzip
    with zipfile.ZipFile(zip_path, 'r') as z:
        z.extractall(f"data/raw/FARS{year}")

# Now you have reproducible data acquisition!

* Overview of dataset
describe

* Summary statistics
summarize

* First few rows
list in 1/10

* Check for missing values
misstable summarize

* Detailed info on specific variables
codebook income education

# Overview of dataset
str(data)
glimpse(data)

# Summary statistics
summary(data)

# First few rows
head(data, 10)

# Check for missing values
colSums(is.na(data))

# Overview of dataset
data.info()
data.dtypes

# Summary statistics
data.describe()

# First few rows
data.head(10)

# Check for missing values
data.isna().sum()

* Create new variable
gen log_income = log(income)

* Conditional assignment
gen high_earner = (income > 100000)

* Modify variable
replace age = . if age < 0

* Drop observations
drop if missing(income)

* Keep only certain variables
keep id income age

pacman::p_load(dplyr)

data <- data %>%
  mutate(
    # Create new variable
    log_income = log(income),
    # Conditional assignment
    high_earner = income > 100000,
    # Modify variable
    age = if_else(age < 0, NA_real_, age)
  ) %>%
  # Drop observations
  filter(!is.na(income)) %>%
  # Keep only certain variables
  select(id, income, age)

import numpy as np

# Create new variable
data["log_income"] = np.log(data["income"])

# Conditional assignment
data["high_earner"] = data["income"] > 100000

# Modify variable
data.loc[data["age"] < 0, "age"] = np.nan

# Drop observations
data = data.dropna(subset=["income"])

# Keep only certain variables
data = data[["id", "income", "age"]]

* Recode numeric values into categories
recode education (1/8 = 1 "Less than HS") ///
                 (9/11 = 2 "Some HS") ///
                 (12 = 3 "HS Graduate") ///
                 (13/15 = 4 "Some College") ///
                 (16/20 = 5 "College+"), gen(educ_cat)

* Simple recode
recode age (18/29 = 1) (30/49 = 2) (50/64 = 3) (65/100 = 4), gen(age_group)

data <- data %>%
  mutate(
    educ_cat = case_when(
      education <= 8 ~ "Less than HS",
      education <= 11 ~ "Some HS",
      education == 12 ~ "HS Graduate",
      education <= 15 ~ "Some College",
      education >= 16 ~ "College+"
    ),
    age_group = cut(age, breaks = c(18, 30, 50, 65, 100),
                    labels = c("18-29", "30-49", "50-64", "65+"),
                    right = FALSE)
  )

# Using pd.cut for numeric ranges
data["age_group"] = pd.cut(data["age"],
                           bins=[18, 30, 50, 65, 100],
                           labels=["18-29", "30-49", "50-64", "65+"])

# Using np.select for complex recoding
conditions = [
    data["education"] <= 8,
    data["education"] <= 11,
    data["education"] == 12,
    data["education"] <= 15,
    data["education"] >= 16
]
choices = ["Less than HS", "Some HS", "HS Graduate", "Some College", "College+"]
# Note: default="" required for numpy 1.25+ when choices are strings
data["educ_cat"] = np.select(conditions, choices, default="")

* i() = unit identifier, j() = what becomes the new column
reshape long income_, i(id) j(year)
rename income_ income

pacman::p_load(tidyr)

data_long <- data %>%
  pivot_longer(
    cols = starts_with("income_"),
    names_to = "year",
    names_prefix = "income_",
    values_to = "income"
  ) %>%
  mutate(year = as.numeric(year))

# Melt wide to long
data_long = pd.melt(
    data,
    id_vars=["id"],
    value_vars=["income_2020", "income_2021"],
    var_name="year",
    value_name="income"
)
# Clean up year column
data_long["year"] = data_long["year"].str.replace("income_", "").astype(int)

* i() = unit identifier, j() = column whose values become new variable names
reshape wide income, i(id) j(year)

* This creates: income2020, income2021, etc.

pacman::p_load(tidyr)

data_wide <- data_long %>%
  pivot_wider(
    id_cols = id,
    names_from = year,
    names_prefix = "income_",
    values_from = income
  )

# Pivot long to wide
data_wide = data_long.pivot(
    index="id",
    columns="year",
    values="income"
).reset_index()

# Rename columns
data_wide.columns = ["id"] + [f"income_{y}" for y in data_wide.columns[1:]]

* Start with crash-level data
use "fars_crashes.dta", clear

* Collapse to state-year level
collapse (count) total_crashes=crash_id ///
         (sum) drunk_crashes=drunk_driver ///
         (sum) total_fatalities=fatalities, ///
         by(state year)

* Now you have one row per state-year
* Ready to merge with policy data

pacman::p_load(dplyr)

# Start with crash-level data
crashes <- read_dta("fars_crashes.dta")

# Collapse to state-year level
state_year <- crashes %>%
  group_by(state, year) %>%
  summarize(
    total_crashes = n(),
    drunk_crashes = sum(drunk_driver, na.rm = TRUE),
    total_fatalities = sum(fatalities, na.rm = TRUE),
    .groups = "drop"
  )

# Now you have one row per state-year

import pandas as pd

# Start with crash-level data
crashes = pd.read_stata("fars_crashes.dta")

# Collapse to state-year level
state_year = crashes.groupby(["state", "year"]).agg(
    total_crashes=("crash_id", "count"),
    drunk_crashes=("drunk_driver", "sum"),
    total_fatalities=("fatalities", "sum")
).reset_index()

# Now you have one row per state-year

* Mean across all observations
egen income_mean = mean(income)

* Mean within groups (by state)
bysort state: egen income_mean_state = mean(income)

* Other useful egen functions:
bysort state: egen income_max = max(income)
bysort state: egen income_sd = sd(income)
bysort state: egen income_count = count(income)
bysort state: egen income_p50 = pctile(income), p(50)

pacman::p_load(dplyr)

# Add group-level statistics while keeping all rows
data <- data %>%
  group_by(state) %>%
  mutate(
    income_mean_state = mean(income, na.rm = TRUE),
    income_max = max(income, na.rm = TRUE),
    income_sd = sd(income, na.rm = TRUE),
    income_count = n(),
    income_p50 = median(income, na.rm = TRUE)
  ) %>%
  ungroup()

# Add group-level statistics while keeping all rows
data["income_mean_state"] = data.groupby("state")["income"].transform("mean")
data["income_max"] = data.groupby("state")["income"].transform("max")
data["income_sd"] = data.groupby("state")["income"].transform("std")
data["income_count"] = data.groupby("state")["income"].transform("count")
data["income_p50"] = data.groupby("state")["income"].transform("median")

* Collapse to state-level means
collapse (mean) income age, by(state)

* Multiple statistics
collapse (mean) income_mean=income ///
         (sd) income_sd=income ///
         (count) n=income, by(state)

* Different statistics for different variables
collapse (mean) income (max) tax_rate (sum) population, by(state year)

# Collapse to state level
state_data <- data %>%
  group_by(state) %>%
  summarize(
    income_mean = mean(income, na.rm = TRUE),
    income_sd = sd(income, na.rm = TRUE),
    n = n()
  )

# Different statistics for different variables
state_year_data <- data %>%
  group_by(state, year) %>%
  summarize(
    income = mean(income, na.rm = TRUE),
    tax_rate = max(tax_rate, na.rm = TRUE),
    population = sum(population, na.rm = TRUE),
    .groups = "drop"
  )

# Collapse to state level
state_data = data.groupby("state").agg(
    income_mean=("income", "mean"),
    income_sd=("income", "std"),
    n=("income", "count")
).reset_index()

# Different statistics for different variables
state_year_data = data.groupby(["state", "year"]).agg(
    income=("income", "mean"),
    tax_rate=("tax_rate", "max"),
    population=("population", "sum")
).reset_index()

* BEFORE merging: Check for unexpected duplicates
duplicates report person_id      // Should say "all unique" for 1:1
duplicates report state          // OK to have duplicates for m:1

* Count rows before and after
count
local n_before = r(N)

merge m:1 state using "state_data.dta"

count
local n_after = r(N)

* If n_after > n_before in a m:1 merge, something is wrong!
di "Rows before: `n_before'"
di "Rows after: `n_after'"
if `n_after' > `n_before' {
    di as error "WARNING: Row count increased - check for duplicate keys!"
}

# BEFORE merging: Check for duplicates
individuals %>% count(person_id) %>% filter(n > 1)  # Should be empty for 1:1
state_data %>% count(state) %>% filter(n > 1)       # Should be empty

# Count rows before and after
n_before <- nrow(individuals)

merged <- merge(individuals, state_data, by = "state", all.x = TRUE)

n_after <- nrow(merged)

# Check if row count changed unexpectedly
if (n_after > n_before) {
  warning("Row count increased! Check for duplicate keys in state_data")
  # Find the problem
  state_data %>% count(state) %>% filter(n > 1)
}

# BEFORE merging: Check for duplicates
print(individuals.groupby("person_id").size().loc[lambda x: x > 1])  # Should be empty
print(state_data.groupby("state").size().loc[lambda x: x > 1])       # Should be empty

# Count rows before and after
n_before = len(individuals)

merged = individuals.merge(state_data, on="state", how="left")

n_after = len(merged)

# Check if row count changed unexpectedly
if n_after > n_before:
    print("WARNING: Row count increased! Check for duplicate keys")
    print(state_data.groupby("state").size().loc[lambda x: x > 1])

* Load master data
use "individuals.dta", clear

* Merge: many individuals per state → m:1
merge m:1 state using "state_data.dta"

* ALWAYS check the results
tab _merge

/*
   _merge |      Freq.
----------+------------
        1 |          1    ← Only in master (NY - no state data)
        2 |          1    ← Only in using (TX - no individuals)
        3 |          3    ← Matched (MA, CA)
*/

* Investigate unmatched
list if _merge == 1   // NY had no state data
list if _merge == 2   // TX had no individuals

* Keep matched only (after understanding why others didn't match)
keep if _merge == 3
drop _merge

pacman::p_load(dplyr, haven)

individuals <- read_dta("individuals.dta")
state_data <- read_dta("state_data.dta")

# Left merge: keep all individuals, add state data where available
merged <- merge(individuals, state_data, by = "state", all.x = TRUE)

# Check for unmatched
merged %>% filter(is.na(min_wage))  # NY had no state data

# Inner merge: keep only matched
merged <- merge(individuals, state_data, by = "state")

individuals = pd.read_stata("individuals.dta")
state_data = pd.read_stata("state_data.dta")

# Left join: keep all individuals, add state data where available
merged = individuals.merge(state_data, on="state", how="left", indicator=True)

# Check merge results
print(merged["_merge"].value_counts())

# Check for unmatched
print(merged[merged["min_wage"].isna()])  # NY had no state data

# Inner join: keep only matched
merged = individuals.merge(state_data, on="state", how="inner")

* Check if your variable is actually numeric or string
browse, nolabel   // Shows actual values, not labels

* Common string cleaning operations:

* 1. Remove periods: "D.C." → "DC"
replace state = subinstr(state, ".", "", .)

* 2. Remove leading/trailing spaces: " MA " → "MA"
replace state = trim(state)

* 3. Extract substring: "Washington DC" → "DC"
replace state = substr(state, -2, 2)  // Last 2 characters

* 4. Convert case: "ma" → "MA"
replace state = upper(state)

* 5. Rename to match other dataset
rename statefip state

* 6. Convert string to numeric (or vice versa)
destring state_code, replace
tostring state_fips, replace

pacman::p_load(stringr)

# Common string cleaning operations:

# 1. Remove periods
data$state <- str_replace_all(data$state, "\\.", "")

# 2. Remove leading/trailing spaces
data$state <- str_trim(data$state)

# 3. Extract substring
data$state <- str_sub(data$state, -2, -1)  # Last 2 characters

# 4. Convert case
data$state <- str_to_upper(data$state)

# 5. Rename to match other dataset
data <- data %>% rename(state = statefip)

# 6. Convert types
data$state_code <- as.numeric(data$state_code)
data$state_fips <- as.character(data$state_fips)

# Common string cleaning operations:

# 1. Remove periods
data["state"] = data["state"].str.replace(".", "", regex=False)

# 2. Remove leading/trailing spaces
data["state"] = data["state"].str.strip()

# 3. Extract substring
data["state"] = data["state"].str[-2:]  # Last 2 characters

# 4. Convert case
data["state"] = data["state"].str.upper()

# 5. Rename to match other dataset
data = data.rename(columns={"statefip": "state"})

# 6. Convert types
data["state_code"] = pd.to_numeric(data["state_code"])
data["state_fips"] = data["state_fips"].astype(str)

* Append multiple years of data
use "data_2020.dta", clear
append using "data_2021.dta"
append using "data_2022.dta"

# Bind rows from multiple dataframes
data_all <- bind_rows(
  read_dta("data_2020.dta"),
  read_dta("data_2021.dta"),
  read_dta("data_2022.dta")
)

# Concatenate multiple dataframes
data_all = pd.concat([
    pd.read_stata("data_2020.dta"),
    pd.read_stata("data_2021.dta"),
    pd.read_stata("data_2022.dta")
], ignore_index=True)

* This INCLUDES missing values (dangerous!)
gen high_income = (income > 100000)

* This is what you actually want
gen high_income = (income > 100000) if !missing(income)

* Or equivalently
gen high_income = (income > 100000 & income < .)

* The ordering is: all numbers < . < .a < .b < ... < .z
* So . > 999999999 evaluates to TRUE

# R handles this more intuitively - NA propagates
high_income <- income > 100000
# Returns NA where income is NA, not TRUE

# But still be explicit about NA handling
high_income <- income > 100000 & !is.na(income)

# Python/pandas handles this intuitively - NaN propagates
high_income = data["income"] > 100000
# Returns NaN where income is NaN, not True

# But still be explicit about NA handling
high_income = (data["income"] > 100000) & data["income"].notna()

* DANGEROUS: If you have >32,767 unique beneficiaries,
* later IDs will all be assigned the same number!
egen bene_id = group(beneficiary_name)

* SAFE: Explicitly use a long integer
egen long bene_id = group(beneficiary_name)

* Or for very large datasets (>2 billion unique values)
egen double bene_id = group(beneficiary_name)

* Always verify you got the right number of groups
distinct beneficiary_name
assert r(ndistinct) == r(N)  // This will fail if top-coded

# R doesn't have this problem - factors work fine
data$bene_id <- as.numeric(factor(data$beneficiary_name))

# Or use dplyr
data <- data %>%
  mutate(bene_id = as.numeric(factor(beneficiary_name)))

# Python doesn't have this problem - uses 64-bit integers by default
data["bene_id"] = pd.factorize(data["beneficiary_name"])[0] + 1

# Verify
assert data["bene_id"].nunique() == data["beneficiary_name"].nunique()

* float has ~7 digits of precision
* This can cause merge failures on IDs!
gen float id = 1234567890
* id is now 1234567936 (corrupted!)

* Always use long or double for IDs
gen long id = 1234567890      // integers up to ~2 billion
gen double id = 1234567890    // integers up to ~9 quadrillion

* Check if a variable might have precision issues
summarize id
* If max is > 16777216 and storage type is float, you have a problem

* Fix: compress optimizes storage types automatically
compress

# R uses double precision by default, so this is less of an issue
# But be careful with very large integers
id <- 1234567890  # Fine in R

# For extremely large integers, use bit64 package
pacman::p_load(bit64)
big_id <- as.integer64("9223372036854775807")

# Python uses 64-bit integers by default, so this is less of an issue
id = 1234567890  # Fine in Python

# pandas also handles large integers well
data["id"] = pd.to_numeric(data["id"], downcast=None)  # Keep full precision

# Check dtypes
print(data.dtypes)

* DANGEROUS: Order of tied observations is random
sort state
gen observation_number = _n  // Different every time you run!

* SAFE: Break ties with a unique identifier
sort state id
gen observation_number = _n  // Now reproducible

* ALTERNATIVE: Use stable sort (preserves existing order for ties)
sort state, stable

* Best practice: Always sort by enough variables to eliminate ties
* or use stable when you don't care about tie-breaking

# R's order() is stable by default (preserves original order for ties)
data <- data[order(data$state), ]

# To be explicit about tie-breaking:
data <- data[order(data$state, data$id), ]

# dplyr::arrange() is also stable
pacman::p_load(dplyr)
data <- data %>% arrange(state, id)

# pandas sort_values() has kind='stable' option
data = data.sort_values('state', kind='stable')

# Better: explicitly break ties
data = data.sort_values(['state', 'id'])

# Note: default sort is 'quicksort' which is unstable

* Suppose you have: income, age, education
summarize inc  // Works! Stata abbreviates to "income"

* Later, you add: income_spouse
summarize inc  // ERROR: "inc" is ambiguous (income vs income_spouse)
* Or worse: no error, wrong variable used silently

* SOLUTION: Disable variable abbreviation at the start of every do-file
set varabbrev off

* Now Stata requires exact variable names
summarize income  // Must use full name

# R does NOT allow variable abbreviation by default
# You must use exact column names
summary(data$income)  # Must be exact

# Partial matching only happens with $ on lists, not data frames
# This is generally not a problem in R

# Python does NOT allow column abbreviation
# You must use exact column names
data['income'].describe()  # Must be exact

# Autocomplete in IDEs helps with long names
# No risk of silent abbreviation errors