###############################################################################
# Basics of R: Exploring the Human Mortality Database
# GREAT LEAP Summer School – Monday, August 25, 2025
# Lecturer: Emre SARI
#
# Data Source:
#   Human Mortality Database (HMD) 
#   Max Planck Institute for Demographic Research (Germany),
#   University of California, Berkeley (USA),
#   French Institute for Demographic Studies (France).
#   Available at https://www.mortality.org (downloaded Aug 6, 2025).
#
# Goal of this script:
#   Learn the basics of R using real mortality data.
#   Practice importing, inspecting, and simple manipulations.
###############################################################################

### 0. BASIC FUNCTIONS (BEFORE IMPORTING DATA) ################################

# Logical comparisons (return TRUE or FALSE)
2 == 2        # equal to
3 != 2        # not equal to
5 < 10        # less than
5 > 10        # greater than
5 <= 5        # less than or equal to
5 >= 6        # greater than or equal to

# Assignment operator: save a value into an object
x <- 5
x            # returns 5
x + 10      # reuse x: add 10
y <- x * 2  # create a new object using x
y

# Combine values into a vector
ages <- c(0, 1, 5, 10)
ages
mean(ages)               # average age
ages * 2                 # multiply all values by 2

# Create sequences
1:5                 # 1, 2, 3, 4, 5
seq(0, 20, by = 5)  # 0, 5, 10, 15, 20
seq(2000, 2020, by = 5)  # years 2000, 2005, ...

# Repeat values
rep("hello", times = 3)     # "hello" "hello" "hello"
rep(1:3, times = 2)         # 1 2 3 1 2 3
rep(1:3, each = 2)          # 1 1 2 2 3 3

# Random sampling
sample(1:100, 5)         # pick 5 random numbers
sample(letters, 3)      # pick 3 random letters
sample(1:100, 5, replace = TRUE)  # allow repeats

  ### 1. SETTING UP THE ENVIRONMENT #############################################

# Working directory = the folder where your files are stored
# Adjust this path to match your own computer.
setwd("/Users/emre/Library/CloudStorage/OneDrive-NORCE/Projects/05. COST - GREAT LEAP/Summer School/Course materials")

# Start fresh: clear all objects from memory
rm(list = ls())

# List files in the working directory (to check if your file is there)
dir()

# Load tidyverse (includes readr, dplyr, ggplot2)
library(tidyverse)

### 2. IMPORTING DATA ##########################################################
# There are many ways to import data in R. Here are the most common ones:

# (A) CSV files (comma-separated values)
mortality_csv <- read.csv("NOR_d_short_idr.csv")       # Base R
mortality <- read_csv("NOR_d_short_idr.csv")           # tidyverse (faster, cleaner)

# (B) CSVs with different delimiters
# Not all "CSV" files use commas! Depending on the source, they may use tabs or semicolons.
# Always check the first few lines of the file in a text editor if something looks wrong.

# Tab-separated values (TSV) files (extension can be .csv or .tsv)
# Example:  col1\tcol2\tcol3
# mortality_tsv <- read_tsv("mortality.csv")   # readr guesses tab separator automatically

# Semicolon-separated values (very common in European CSV files)
# Example:  col1;col2;col3
# mortality_semicolon <- read_delim("mortality.csv", delim = ";")

# Tip: If read_csv() shows the entire row in one column, 
# try read_tsv() or read_delim() with the right delimiter.

# (C) Excel files
# install.packages("readxl")   # run once if not installed
# library(readxl)
# mortality_xlsx <- read_excel("mortality.xlsx", sheet = 1)

# (D) Stata, SPSS, or SAS data
# install.packages("haven")
# library(haven)
# mortality_dta <- read_dta("mortality.dta")           # Stata
# mortality_sav <- read_sav("mortality.sav")           # SPSS

# (E) R native formats
# load("mortality.RData")     # loads one or more R objects
# mortality_rds <- readRDS("mortality.rds")

# ---> Data removing
# Remove only specific objects
rm(mortality_csv)
# Or, if you want to be extra safe, remove both and reload just one
rm(mortality, mortality_csv)

# For this lecture, we will continue with the tidyverse import:
mortality <- read_csv("NOR_d_short_idr.csv")

### 3. INSPECTING THE DATA #####################################################

# 3A. Base tidyverse approaches
glimpse(mortality)     # column names, types, and first values
head(mortality, 5)     # first 5 rows
tail(mortality, 5)     # last 5 rows
summary(mortality)     # quick summary for numeric columns
names(mortality)       # list all column names
nrow(mortality)        # number of rows
ncol(mortality)        # number of columns

# Count how many rows contain missing values in total deaths
length(which(is.na(mortality$cause)))

# 3B. Using psych package for more detailed descriptive stats
# install.packages("psych") # run this once if not installed
library(psych)
describe(mortality)

# 3C. Calculate basic statistics (average, min, max) for numeric variables
# Goal: What is the average number of deaths per sex?

# Step 1: Filter data by sex = 1 (usually male)
male_data <- mortality[mortality$sex == 1, ]

# Bonus: By using tidyverse
mortality %>%
  filter(sex == 1) %>%
  summarise(mean_deaths = mean(total, na.rm = TRUE))
# %>% is called the pipe operator. It comes from the tidyverse, and it helps you 
# read code from left to right, like a sentence.
# You can read the above as:
#    “Take the mortality data, then filter it for males, then calculate 
#     the mean number of deaths.”

# Step 2: Calculate statistics for total deaths
mean(male_data$total, na.rm = TRUE)
sd(male_data$total, na.rm = TRUE)
min(male_data$total, na.rm = TRUE)
max(male_data$total, na.rm = TRUE)

# Step 3: Now do the same for sex = 2 (usually female)
female_data <- mortality[mortality$sex == 2, ]
mean(female_data$total, na.rm = TRUE)

# 3D. Frequency tables (categorical variable)
#Goal: How many times does each cause appear?

# Option 1: Frequency of a single variable r Copy Edit
# Count how many times each cause appears
table(mortality$cause)

# Option 2: Two-way table (cross-tabulation)
# How many times each sex and cause combination occurs
table(mortality$sex, mortality$cause)

# Option 3: Three-way table
# Sex x Year x Cause — might be large, but useful to show structure
table(mortality$sex, mortality$year, mortality$cause)

# Tip: Use only a subset for display:
# For example: only year 1996
table(mortality$sex[mortality$year == 1996], 
      mortality$cause[mortality$year == 1996])

# Option 4: Include missing values
# Add useNA = "ifany" to include NA counts
table(mortality$cause, useNA = "ifany")

# Option 5: Show proportions / percentages
# Get proportions (relative frequencies)
prop.table(table(mortality$cause))

# Two-way proportion table (e.g., within sex)
prop.table(table(mortality$sex, mortality$cause), margin = 1)# Row-wise (by sex)

# Option 6: Save table as an object and inspect
# Save as object
my_table <- table(mortality$sex, mortality$cause)

# Inspect structure
str(my_table)

# Convert to a data frame if needed
as.data.frame(my_table)

# 3E. Cross-tabulation (combine two variables)
#Goal: How many deaths were recorded for each sex in 1996? 

# Step 1: Filter data for year 1996
mort_1996 <- mortality[mortality$year == 1996, ]

# Step 2: Separate by sex
male_1996 <- mort_1996[mort_1996$sex == 1, ]
female_1996 <- mort_1996[mort_1996$sex == 2, ]

# Step 3: Sum deaths
sum(male_1996$total, na.rm = TRUE)
sum(female_1996$total, na.rm = TRUE)

#Bonus:
# Select specific columns by position
mortality[ , c(1, 3, 6)]

# Select columns by name
mortality[ , c("year", "sex", "total")]

# Or with dplyr (recommended)
mortality %>%
  filter(year == 1996) %>%
  select(year, sex, total)

# Another example: filter by multiple values
mortality %>%
  filter(year == 1996, sex == 3) %>%
  select(year, total)

### 4. DERIVED VARIABLES ######################################################

# Create a logical variable: is total above 1000?
mortality <- mortality %>%
  mutate(high_mortality = total > 1000)

# Create age group labels based on age-specific columns
mortality <- mortality %>%
  mutate(age_group_0_14 = d0 + d1 + d5 + d10,       # Ages 0–14
         age_group_15_64 = d15 + d20 + d25 + d30 + d35 + d40 + d45 + d50 + d55 + d60,  # Ages 15–64
         age_group_65plus = d65 + d70 + d75 + d80 + d85p)                              # 65+

# Use case_when() for age group tags (based on deaths in specific columns)
mortality <- mortality %>%
  mutate(main_age_group = case_when(
    d0 + d1 + d5 + d10 > 0 ~ "0-14",
    d15 + d20 + d25 + d30 + d35 + d40 + d45 + d50 + d55 + d60 > 0 ~ "15-64",
    d65 + d70 + d75 + d80 + d85p > 0 ~ "65+",
    TRUE ~ "unknown"
  ))

# Check result
glimpse(mortality)

### 5. GROUPING AND SUMMARISING ###############################################

# Mean number of deaths by sex
mortality %>%
  group_by(sex) %>%
  summarise(mean_deaths = mean(total, na.rm = TRUE))

# Sum of deaths by year and sex
mortality %>%
  group_by(year, sex) %>%
  summarise(total_deaths = sum(total, na.rm = TRUE)) %>%
  arrange(year)

### 4. VERY BASIC PLOTS ########################################################
###############################################################################
# Plotting Deaths in R – Step by Step (Basic Learners)
# Based on Human Mortality Database (HMD) data
###############################################################################

### STEP 1: Basic Plot of One Variable
# This plots the total deaths for each row in the dataset.
# It doesn’t mean much yet, but it shows what a simple plot looks like.

plot(mortality$total)

# Extra: Other simple plot types (same data)
hist(mortality$total)
boxplot(mortality$total)

### STEP 2: Add Labels and a Title
# Always label your axes and give your plot a title.

plot(mortality$total,
     main = "Total Deaths (Raw Order)",
     xlab = "Observation Number",
     ylab = "Total Deaths")

### STEP 3: Plot Total Deaths by Year
# Plotting deaths against actual years gives more meaning.
# But there are multiple entries per year, so it’s still messy.

plot(mortality$year, mortality$total,
     main = "Total Deaths Over Time (Raw)",
     xlab = "Year",
     ylab = "Total Deaths")

### STEP 4: Clean It Up – Summarise Deaths by Year (Base R)
# First, create a summary using tapply()

deaths_by_year <- tapply(mortality$total, mortality$year, sum, na.rm = TRUE)

# Now plot just one value per year — much cleaner.

plot(names(deaths_by_year), deaths_by_year,
     type = "l",                      # "l" = line
     main = "Total Deaths Per Year",
     xlab = "Year",
     ylab = "Total Deaths")

### STEP 5: Make It Nicer – Add Points, Color, Thickness

plot(names(deaths_by_year), deaths_by_year,
     type = "b",                     # "b" = both line and points
     col = "darkgreen",              # line color
     pch = 19,                        # point symbol
     lwd = 2,                         # line width
     main = "Total Deaths Per Year",
     xlab = "Year",
     ylab = "Total Deaths")

### 6. FOCUSED EXERCISE: WHO DIES MORE AND WHEN? ##############################

# Filter for years after 1990
mortality_90 <- mortality[mortality$year > 1990, ]

# Separate by sex
male_90 <- mortality_90[mortality_90$sex == 1, ]
female_90 <- mortality_90[mortality_90$sex == 2, ]

# Calculate average deaths
mean(male_90$total, na.rm = TRUE)
mean(female_90$total, na.rm = TRUE)

# Find min and max year (should all be >1990)
min(male_90$year)
max(male_90$year)

# Create basic plots
plot(male_90$year, male_90$total,
     type = "l",
     col = "blue",
     main = "Male Deaths After 1990",
     xlab = "Year",
     ylab = "Total Deaths")

plot(female_90$year, female_90$total,
     type = "l",
     col = "red",
     main = "Female Deaths After 1990",
     xlab = "Year",
     ylab = "Total Deaths")