############################################################
# Correlation and Regression Analysis of Songkorpus Features
# This script computes correlations, regression models,
# and diachronic trends for linguistic features in song lyrics
############################################################

# -----------------------------
# DATA PREPARATION (CORRELATION)
# -----------------------------

# Select relevant columns and exclude metadata
df_corr <- df_num[, !names(df_num) %in% c(
  "CO_SIGLE", "CO_TOKENCOUNT", "CO_TYPECOUNT", "CO_STTR",
  "CO_YEAR", "CO_SENTI", "CO_SENTI_TOKEN_RATIO",
  "CO_SENTI_POS", "CO_SENTI_NEG", "CO_PART", "CO_SENTIWORDS"
)]

# Filter out rows without modal particles or pronouns
df_corr <- df_corr[df_corr$CO_PART_TOKEN_RATIO != 0, ]
df_corr <- df_corr[df_corr$CO_PRON1ST != 0, ]
df_corr <- df_corr[df_corr$CO_PRON2ND != 0, ]
df_corr <- df_corr[df_corr$CO_PRON3RD != 0, ]

# Compute pronoun ratio (1st + 2nd person vs. 3rd person)
df_corr$PRON_DIFF <- (df_corr$CO_PRON1ST + df_corr$CO_PRON2ND) / df_corr$CO_PRON3RD

# Remove redundant columns
df_corr <- df_corr[, !(names(df_corr) %in% c("CO_ID", "CO_PRON1ST", "CO_PRON2ND", "CO_PRON3RD"))]

# Rename variables for readability
colnames(df_corr) <- c("MATTR", "MTLD", "SI", "SENTIMENT_POLARITY", "PTKM", "PRON")

# -----------------------------
# CORRELATION ANALYSIS
# -----------------------------

# Compute Spearman correlation matrix (robust to outliers)
cor_mat <- cor(df_corr, method = "spearman", use = "complete.obs")
print(cor_mat)

# Load visualization package
library(corrplot)

# Plot full correlation heatmap
corrplot(cor_mat,
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45)

# Plot selected subset of correlations
corrplot(cor_mat[c("MATTR", "MTLD", "SI", "PTKM", "PRON"),
                 c("MATTR", "MTLD", "SI", "PTKM", "PRON")],
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45)

# -----------------------------
# REGRESSION MODELS
# -----------------------------

# Fit multiple linear models to test relationships between variables

model <- lm(PTKM ~ SI + MATTR + PRON, data = df_corr)
summary(model)

model <- lm(SI ~ PTKM + MATTR + PRON, data = df_corr)
summary(model)

model <- lm(MATTR ~ PTKM + SI + PRON, data = df_corr)
summary(model)

model <- lm(PRON ~ PTKM + SI + MATTR, data = df_corr)
summary(model)

# -----------------------------
# DIACHRONIC ANALYSIS
# -----------------------------

library(ggplot2)
library(gridExtra)

# Prepare diachronic dataset (exclude metadata)
df_diachron <- df_num[, !names(df_num) %in% c(
  "CO_SIGLE", "CO_TOKENCOUNT", "CO_TYPECOUNT", "CO_STTR",
  "CO_SENTI", "CO_SENTI_TOKEN_RATIO",
  "CO_SENTI_POS", "CO_SENTI_NEG", "CO_PART", "CO_SENTIWORDS"
)]

# Apply same filtering as above
df_diachron <- df_diachron[df_diachron$CO_PART_TOKEN_RATIO != 0, ]
df_diachron <- df_diachron[df_diachron$CO_PRON1ST != 0, ]
df_diachron <- df_diachron[df_diachron$CO_PRON2ND != 0, ]
df_diachron <- df_diachron[df_diachron$CO_PRON3RD != 0, ]

# Restrict to songs from 1970 onwards
df_diachron <- df_diachron[df_diachron$CO_YEAR >= 1970, ]

# Compute pronoun ratio
df_diachron$PRON_DIFF <- (df_diachron$CO_PRON1ST + df_diachron$CO_PRON2ND) / df_diachron$CO_PRON3RD

# -----------------------------
# LINEAR MODELS (TIME TRENDS)
# -----------------------------

m1 <- lm(CO_MATTR ~ CO_YEAR, data = df_diachron)
m2 <- lm(CO_PART_TOKEN_RATIO ~ CO_YEAR, data = df_diachron)
m3 <- lm(CO_SENTI_SENTIWORDS_RATIO ~ CO_YEAR, data = df_diachron)
m4 <- lm(PRON_DIFF ~ CO_YEAR, data = df_diachron)

# Summaries of diachronic models
summary(m1)
summary(m2)
summary(m3)
summary(m4)

# -----------------------------
# VISUALIZATION
# -----------------------------

# MATTR over time
p1 <- ggplot(df_diachron, aes(x = CO_YEAR, y = CO_MATTR)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Moving Average Type-Token Ratio per Year",
       x = "",
       y = "MATTR") +
  theme_minimal()

# Modal particle density over time (filtered)
p2 <- ggplot(
  subset(df_diachron, CO_PART_TOKEN_RATIO <= 0.1),
  aes(x = CO_YEAR, y = CO_PART_TOKEN_RATIO)
) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Relative Modal Particle Rate by Year",
       x = "",
       y = "MP per Token") +
  theme_minimal()

# Sentiment over time (filtered)
p3 <- ggplot(
  subset(df_diachron, CO_SENTI_SENTIWORDS_RATIO <= 20000000000000),
  aes(x = CO_YEAR, y = CO_SENTI_SENTIWORDS_RATIO)
) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Sentiment per Year",
       x = "",
       y = "Combined Score") +
  theme_minimal()

# Pronoun ratio over time (filtered)
p4 <- ggplot(
  subset(df_diachron, PRON_DIFF <= 30),
  aes(x = CO_YEAR, y = PRON_DIFF)
) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Pronoun Ratio per Year",
       x = "Year",
       y = "1st/2nd vs. 3rd Person") +
  theme_minimal()

# Arrange plots in 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)