# =============================================================================
# DAEB 2026/27 - Lecture 8: Dispersion and concentration measures in R;
#                           data cleaning
# Name: ____________________________      Date: 1 October 2026
#
# Before you start:
#   a. Open the RStudio project DAEB.Rproj (File > Open Project), or make sure
#      the working directory is the folder that contains the "data" and
#      "results" folders.
#   b. Check with getwd() that you are in the right place.
#   c. Files needed: data/ages.csv, data/earnings.csv, data/workers_raw.csv.
#      This script writes data/workers_clean.csv and four files in results/.
# =============================================================================


# --- 1. Packages and data ----------------------------------------------------

#install.packages("DescTools")   # run once: package for the Lorenz curve and Gini index
#install.packages("readr")       # run once, if library(readr) below fails

library(DescTools)  # inequality: Gini(), Lc()
library(dplyr)      # tables: count(), mutate(), summarise()
library(knitr)      # formatted tables: kable()
library(ggplot2)    # charts

ages <- read.csv("data/ages.csv")            # 50 students (Lecture 4)
earn <- read.csv("data/earnings.csv")        # 35 000 workers aged 25-30 (Lecture 5)


# --- 2. Dispersion measures with raw data ------------------------------------

# Range and interquartile range:  R = max - min,  IQR = Q3 - Q1
range(earn$earnings)
x <- earn$earnings
max(x) - min(x)                                       # range: raw data
IQR(x)                                                # interquartile range: raw data

# Mean deviation: 
mean(abs(x - mean(x)))                            # earnings
mean(abs(ages$age - mean(ages$age)))              # age

# Variance and standard deviation: R divides by N - 1, the course divides by N.
# var(x) and sd(x) give the sample versions. Write your own function:
sd_N <- function(x) {            # course formula: mean of squared deviations, divisor N
  v <- mean((x - mean(x))^2)     # variance with divisor N
  sqrt(v)                        # returned value
}
c(sd_N(ages$age), sd(ages$age))              # course vs R default (N - 1), N = 50: they differ
c(sd_N(x), sd(x))                            # with N = 35 000 the difference vanishes

# Relative dispersion:  RIQR = IQR / median,  CV = S / mean
riqr_earn <- IQR(x) / median(x)                  # RIQR: earnings
cv_earn   <- sd_N(x) / mean(x)                   # CV: earnings

riqr_age  <- IQR(ages$age) / median(ages$age)    # RIQR: age
cv_age    <- sd_N(ages$age) / mean(ages$age)     # CV: age
rel_disp <- data.frame(
  Variable = c("Earnings", "Age"),
  RIQR     = c(riqr_earn, riqr_age),
  CV       = c(cv_earn, cv_age)
)
kable(rel_disp, digits = 3)


# --- 3. Concentration measures -----------------------------------------------

# Lorenz curve

## Included: code for saving R based plots 
png("results/lorenz_earnings.png", width = 1600, height = 1200, res = 200)   # open the file (pixels; res = dpi)

lc <- Lc(x) 
plot(lc,
     main = "Lorenz curve - earnings of 35,000 workers",
     xlab = "Cumulative share of workers (p)",
     ylab = "Cumulative share of earnings (q)",
     col  = "darkred", lwd = 2)
dev.off()                                                                   # close the file: nothing is written before this

# Gini index
x_sorted <- sort(earn$earnings)
N_e      <- length(x_sorted)
q_i <- cumsum(x_sorted) / sum(x_sorted)   # cumulative share of earnings

c(bottom_half = q_i[N_e / 2], top_10pct = 1 - q_i[0.9 * N_e])   # shares held by the bottom 50% and the top 10%

Gini(x)           # Gini index (matches the course formula)


# --- 4. Data cleaning: import and inspect ------------------------------------

# The data set: 200 workers (hypothetical survey), file data/workers_raw.csv.
# Variables: id, wage (gross monthly wage, EUR), age, gender, education.
# Before running any code, open the file in RStudio (Files pane > View File)
# and keep a list of what looks wrong. The rest of the script goes through it.

# read.csv() is built into R; readr's read_csv() does the same job, faster and
# with a report of what it read (column types). 
library(readr)
workers <- read_csv("data/workers_raw.csv")

glimpse(workers)      # one line per variable: type and first values
summary(workers)      # numbers: min, quartiles, mean, max, NA's; text: only the length
# wage: minimum -999 and NA's; age: minimum 7, maximum 150; 203 rows instead of 200

# summary() says nothing useful about text. count() lists every distinct value:
workers |> count(gender)
workers |> count(education)
# two categories written six ways; three levels written six ways; NA in both


# --- 5. Problem 1: missing values --------------------------------------------

# The code -999 must become NA, otherwise it enters every mean.
workers <- workers |>
  mutate(wage = na_if(wage, -999))     # -999 -> NA; everything else unchanged
colSums(is.na(workers))                # missing values per column
summary(workers)      # numbers: min, quartiles, mean, max, NA's; text: only the length

# With NA present, R refuses to guess:
mean(workers$wage)                     # NA: R will not compute a mean it cannot trust
mean(workers$wage, na.rm = TRUE)       # mean of the observed wages
# na.rm = TRUE drops the NA's for that computation only. Report how many were
# dropped. The cleaner route is to declare the codes when reading:
#   read_csv("data/workers_raw.csv", na = c("", "NA", "-999"))


# --- 6. Problem 2: duplicated rows -------------------------------------------

nrow(workers)                                   # rows in the file
n_distinct(workers$id)                          # distinct workers
workers |> count(id) |> filter(n > 1)           # which ids appear more than once

workers <- workers |> distinct()                # keep one copy of each identical row
nrow(workers)
# distinct() removes rows identical in EVERY column. The same id with two
# different wages is a conflict, not a duplicate. Check count(id) after distinct().


# --- 7. Problem 3: impossible values -----------------------------------------

workers |> filter(age < 16 | age > 70)          # | means "or"
workers |> arrange(desc(wage)) |> head(3)       # the three highest wages

# We cannot know the true age of the 150-year-old: NA. The wage of 25 000 is
# probably 2 500 with an extra zero, but "probably" is a guess: NA too, or
# correct it and write down that you did.
workers <- workers |>
  mutate(age  = if_else(age < 16 | age > 70, NA, age),    # if_else(condition, value if TRUE, value if FALSE)
         wage = if_else(wage > 10000, NA, wage))
colSums(is.na(workers))


# --- 8. Problem 4: inconsistent text -----------------------------------------

# trimws() removes spaces at the ends (safety net: read_csv() already did it),
# tolower() puts everything in lower case, case_when() maps each spelling to one label.
workers <- workers |>
  mutate(gender = tolower(trimws(gender)),                     # "Female" -> "female", "M" -> "m"
         gender = case_when(gender %in% c("female", "f") ~ "Female",
                            gender %in% c("male", "m")   ~ "Male"))
workers |> count(gender)
# %in% asks "is the value one of these?". Anything not listed in case_when()
# becomes NA, including the NA's already there: extra NA's mean a missed spelling.

workers <- workers |>
  mutate(education = tolower(trimws(education)),
         education = case_when(education == "basic"     ~ "Basic",
                               education == "secondary" ~ "Secondary",
                               education == "higher"    ~ "Higher"))
workers |> count(education)


# --- 9. Factors: telling R a variable is categorical -------------------------

# gender and education are still text (chr). A factor has a fixed set of
# levels, in an order you choose: tables and charts follow that order, a value
# outside the levels becomes NA, and summary() counts each level.
workers <- workers |>
  mutate(gender    = factor(gender, levels = c("Female", "Male")),                       # nominal
         education = factor(education, levels = c("Basic", "Secondary", "Higher"),
                            ordered = TRUE))                                             # ordinal: levels are ranked


# Check: the clean data set
str(workers)
summary(workers)
colSums(is.na(workers))
# 200 workers, two numeric variables with a plausible range, two factors with
# the right levels, and the missing values counted.

# --- 10. Save the clean file -------------------------------------------------

write_csv(workers, "data/workers_clean.csv")     # a new file; the raw file stays untouched
# A CSV stores text, not factors: after read_csv("data/workers_clean.csv") the
# two factors come back as chr. Redefine them with factor() in the analysis script.

# --- 12. From cleaning to analysis: measures by group ------------------------

wage_by_educ <- workers |>                         # store the table in an object: saved below
  group_by(education) |>                           # one row of results per level
  summarise(n        = n(),                        # workers in the group
            n_wage   = sum(!is.na(wage)),          # ... with an observed wage
            mean     = mean(wage, na.rm = TRUE),   # Lecture 6, now with the NA's dropped explicitly
            median   = median(wage, na.rm = TRUE),
            sd       = sd_N(na.omit(wage)),        # our function from section 2; na.omit() drops the NA's first
            cv       = sd / mean)                  # coefficient of variation by group
wage_by_educ

ggplot(filter(workers, !is.na(education), !is.na(wage)),
       aes(x = education, y = wage)) +
  geom_boxplot(coef = Inf, width = 0.5, fill = "grey90") +   # whiskers to min and max, as in the course
  labs(x = NULL, y = "Gross monthly wage (EUR)")


# --- 13. Save a table of descriptive statistics ------------------------------

# A table built with summarise() is a data frame: save it by code, into results/,
# never by copying numbers from the console. Rerunning overwrites the files, so
# the tables in the report always match the data.
write_csv(wage_by_educ, "results/wage_by_education.csv")        # 1. CSV: opens in any spreadsheet or word processor