packages <- c(
  "readxl", "tidyverse", "janitor", "psych",
  "lavaan", "semTools", "naniar", "mice",
  "lme4", "lmerTest", "performance", "car",
  "broom.mixed", "tibble", "purrr", "tidyr", "dplyr"
)

installed_packages <- rownames(installed.packages())

for (pkg in packages) {
  if (!pkg %in% installed_packages) {
    install.packages(pkg)
  }
}

library(readxl)
library(tidyverse)
library(janitor)
library(psych)
library(lavaan)
library(semTools)
library(naniar)
library(mice)
library(lme4)
library(lmerTest)
library(performance)
library(car)
library(broom.mixed)
library(tibble)
library(purrr)
library(tidyr)
library(dplyr)



data <- read_excel("~/Desktop/Longitudinal associations/Missingbai_aztertzekoprest_long.xlsx")

data <- data %>% 
  clean_names()

names(data)
head(data)
str(data)
summary(data)

table(data$time, useNA = "ifany")
unique(data$time)

data <- data %>%
  mutate(
    student_id = as.factor(student_id),
    school = as.factor(school),
    grade = as.factor(grade),
    gender = as.factor(gender),
    repeat_student = as.factor(repeat_student)
  )


data <- data %>%
  mutate(
    time = as.numeric(time),
    time = factor(
      time,
      levels = c(1, 2, 3),
      labels = c("T1", "T2", "T3")
    )
  )


table(data$time, useNA = "ifany")
levels(data$time)

data <- data %>%
  mutate(
    time = relevel(time, ref = "T1"),
    grade = relevel(grade, ref = "G7"),
    school = relevel(school, ref = levels(school)[1])
  )

levels(data$time)
levels(data$grade)
levels(data$school)

missing_by_variable <- data %>%
  summarise(across(everything(), ~ mean(is.na(.)) * 100)) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "missing_percent"
  ) %>%
  arrange(desc(missing_percent))

missing_by_variable
print(missing_by_variable, n = Inf)

missing_by_time <- data %>%
  group_by(time) %>%
  summarise(
    n_rows = n(),
    n_students = n_distinct(student_id),
    missing_learning = mean(is.na(learning)) * 100,
    missing_participation = mean(is.na(participation)) * 100,
    missing_social_climate = mean(is.na(social_climate)) * 100,
    missing_self_efficacy = mean(is.na(self_efficacy_for_learning)) * 100,
    .groups = "drop"
  )

missing_by_time

student_missing_pattern <- data %>%
  group_by(student_id) %>%
  summarise(
    n_learning_available = sum(!is.na(learning)),
    n_participation_available = sum(!is.na(participation)),
    n_social_available = sum(!is.na(social_climate)),
    n_se_available = sum(!is.na(self_efficacy_for_learning)),
    .groups = "drop"
  )

table(student_missing_pattern$n_learning_available)
table(student_missing_pattern$n_participation_available)
table(student_missing_pattern$n_social_available)
table(student_missing_pattern$n_se_available)

se_items <- paste0("se_", 1:19)

se_items %in% names(data)
setdiff(se_items, names(data))

grep("^se_", names(data), value = TRUE)

alpha_by_time <- function(dataset, time_point, items) {
  
  temp <- dataset %>%
    filter(time == time_point) %>%
    select(all_of(items))
  
  psych::alpha(temp)
}

alpha_se_t1_19 <- alpha_by_time(data, "T1", se_items)
alpha_se_t2_19 <- alpha_by_time(data, "T2", se_items)
alpha_se_t3_19 <- alpha_by_time(data, "T3", se_items)

alpha_se_t1_19$total
alpha_se_t2_19$total
alpha_se_t3_19$total

alpha_se_t1_19$item.stats
alpha_se_t2_19$item.stats
alpha_se_t3_19$item.stats

se_items_18 <- setdiff(se_items, "se_1")

alpha_se_t1_18 <- alpha_by_time(data, "T1", se_items_18)
alpha_se_t2_18 <- alpha_by_time(data, "T2", se_items_18)
alpha_se_t3_18 <- alpha_by_time(data, "T3", se_items_18)

alpha_se_t1_18$total
alpha_se_t2_18$total
alpha_se_t3_18$total

set.seed(1234)

student_split <- data %>%
  distinct(student_id) %>%
  mutate(
    split_group = sample(c("EFA", "CFA"), size = n(), replace = TRUE)
  )

data_split <- data %>%
  left_join(student_split, by = "student_id")

table(data_split$split_group)

data_efa <- data_split %>%
  filter(split_group == "EFA")

efa_se <- function(dataset, time_point, items) {
  
  temp <- dataset %>%
    filter(time == time_point) %>%
    select(all_of(items))
  
  poly_cor <- psych::polychoric(temp)$rho
  
  psych::fa(
    poly_cor,
    nfactors = 1,
    n.obs = nrow(temp),
    fm = "minres"
  )
}

efa_t1 <- efa_se(data_efa, "T1", se_items_18)
efa_t2 <- efa_se(data_efa, "T2", se_items_18)
efa_t3 <- efa_se(data_efa, "T3", se_items_18)

print(efa_t1$loadings, cutoff = .30)
print(efa_t2$loadings, cutoff = .30)
print(efa_t3$loadings, cutoff = .30)

data_cfa <- data_split %>%
  filter(split_group == "CFA")

se_model <- '
SE =~ se_2 + se_3 + se_4 + se_5 + se_6 + se_7 + se_8 + se_9 +
      se_10 + se_11 + se_12 + se_13 + se_14 + se_15 +
      se_16 + se_17 + se_18 + se_19
'

cfa_by_time <- function(dataset, time_point) {
  
  temp <- dataset %>%
    filter(time == time_point)
  
  lavaan::cfa(
    se_model,
    data = temp,
    estimator = "WLSMV",
    ordered = se_items_18
  )
}

fit_t1_cfa <- cfa_by_time(data_cfa, "T1")
fit_t2_cfa <- cfa_by_time(data_cfa, "T2")
fit_t3_cfa <- cfa_by_time(data_cfa, "T3")

fitMeasures(fit_t1_cfa, c("cfi", "tli", "rmsea", "srmr"))
fitMeasures(fit_t2_cfa, c("cfi", "tli", "rmsea", "srmr"))
fitMeasures(fit_t3_cfa, c("cfi", "tli", "rmsea", "srmr"))

min_se_items <- ceiling(length(se_items_18) * 0.80)

data <- data %>%
  rowwise() %>%
  mutate(
    se_valid_items = sum(!is.na(c_across(all_of(se_items_18)))),
    se_mean_final = if_else(
      se_valid_items >= min_se_items,
      mean(c_across(all_of(se_items_18)), na.rm = TRUE),
      NA_real_
    )
  ) %>%
  ungroup()

summary(data$se_mean_final)

learning_items <- paste0("l_", 1:7)
participation_items <- paste0("p_", 1:4)
social_items <- paste0("sc_", 1:5)

setdiff(learning_items, names(data))
setdiff(participation_items, names(data))
setdiff(social_items, names(data))
setdiff(se_items_18, names(data))

alpha_summary <- tibble(
  variable = rep(c("Learning", "Participation", "Social climate", "Self-efficacy"), each = 3),
  time = rep(c("T1", "T2", "T3"), times = 4),
  alpha = c(
    alpha_by_time(data, "T1", learning_items)$total$raw_alpha,
    alpha_by_time(data, "T2", learning_items)$total$raw_alpha,
    alpha_by_time(data, "T3", learning_items)$total$raw_alpha,
    
    alpha_by_time(data, "T1", participation_items)$total$raw_alpha,
    alpha_by_time(data, "T2", participation_items)$total$raw_alpha,
    alpha_by_time(data, "T3", participation_items)$total$raw_alpha,
    
    alpha_by_time(data, "T1", social_items)$total$raw_alpha,
    alpha_by_time(data, "T2", social_items)$total$raw_alpha,
    alpha_by_time(data, "T3", social_items)$total$raw_alpha,
    
    alpha_by_time(data, "T1", se_items_18)$total$raw_alpha,
    alpha_by_time(data, "T2", se_items_18)$total$raw_alpha,
    alpha_by_time(data, "T3", se_items_18)$total$raw_alpha
  )
)

alpha_summary

student_level <- data %>%
  distinct(student_id, school, grade, gender, repeat_student)

n_distinct(student_level$student_id)

student_level %>%
  count(school)

student_level %>%
  count(grade)

student_level %>%
  count(gender)

student_level %>%
  count(school, grade, gender) %>%
  arrange(school, grade, gender)

descriptives_main <- data %>%
  group_by(time) %>%
  summarise(
    learning_mean = mean(learning, na.rm = TRUE),
    learning_sd = sd(learning, na.rm = TRUE),
    
    participation_mean = mean(participation, na.rm = TRUE),
    participation_sd = sd(participation, na.rm = TRUE),
    
    social_climate_mean = mean(social_climate, na.rm = TRUE),
    social_climate_sd = sd(social_climate, na.rm = TRUE),
    
    se_mean = mean(se_mean_final, na.rm = TRUE),
    se_sd = sd(se_mean_final, na.rm = TRUE),
    
    n_learning = sum(!is.na(learning)),
    n_participation = sum(!is.na(participation)),
    n_social = sum(!is.na(social_climate)),
    n_se = sum(!is.na(se_mean_final)),
    .groups = "drop"
  )

descriptives_main

cor_vars <- c(
  "learning",
  "participation",
  "social_climate",
  "se_mean_final"
)

cor_by_time <- function(dataset, time_point) {
  
  temp <- dataset %>%
    filter(time == time_point) %>%
    select(all_of(cor_vars))
  
  psych::corr.test(
    temp,
    use = "pairwise",
    method = "pearson"
  )
}

cor_t1 <- cor_by_time(data, "T1")
cor_t2 <- cor_by_time(data, "T2")
cor_t3 <- cor_by_time(data, "T3")

cor_t1$r
cor_t2$r
cor_t3$r

create_correlation_table <- function(dataset, time_point, vars) {
  
  temp <- dataset %>%
    filter(time == time_point) %>%
    select(all_of(vars))
  
  variable_pairs <- combn(vars, 2, simplify = FALSE)
  
  correlation_table <- map_dfr(variable_pairs, function(pair) {
    
    var1 <- pair[1]
    var2 <- pair[2]
    
    pair_data <- temp %>%
      select(all_of(c(var1, var2))) %>%
      drop_na()
    
    test <- cor.test(
      pair_data[[var1]],
      pair_data[[var2]],
      method = "pearson"
    )
    
    tibble(
      Time = time_point,
      Variable_1 = var1,
      Variable_2 = var2,
      n = nrow(pair_data),
      df = unname(test$parameter),
      Pearson_r = unname(test$estimate),
      CI_lower = test$conf.int[1],
      CI_upper = test$conf.int[2],
      p_value = test$p.value
    )
  })
  
  return(correlation_table)
}

table_3_t1 <- create_correlation_table(data, "T1", cor_vars)
table_4_t2 <- create_correlation_table(data, "T2", cor_vars)
table_5_t3 <- create_correlation_table(data, "T3", cor_vars)

table_3_t1
table_4_t2
table_5_t3

table_3_t1_round <- table_3_t1 %>%
  mutate(
    Pearson_r = round(Pearson_r, 3),
    CI_lower = round(CI_lower, 3),
    CI_upper = round(CI_upper, 3),
    p_value = ifelse(p_value < .001, "< .001", sprintf("%.3f", p_value))
  )

table_4_t2_round <- table_4_t2 %>%
  mutate(
    Pearson_r = round(Pearson_r, 3),
    CI_lower = round(CI_lower, 3),
    CI_upper = round(CI_upper, 3),
    p_value = ifelse(p_value < .001, "< .001", sprintf("%.3f", p_value))
  )

table_5_t3_round <- table_5_t3 %>%
  mutate(
    Pearson_r = round(Pearson_r, 3),
    CI_lower = round(CI_lower, 3),
    CI_upper = round(CI_upper, 3),
    p_value = ifelse(p_value < .001, "< .001", sprintf("%.3f", p_value))
  )

table_3_t1_round
table_4_t2_round
table_5_t3_round

data_model <- data %>%
  group_by(student_id) %>%
  mutate(
    participation_between = mean(participation, na.rm = TRUE),
    social_between = mean(social_climate, na.rm = TRUE),
    se_between = mean(se_mean_final, na.rm = TRUE),
    
    participation_within = participation - participation_between,
    social_within = social_climate - social_between,
    se_within = se_mean_final - se_between
  ) %>%
  ungroup()

summary(data_model$participation_within)
summary(data_model$participation_between)

summary(data_model$social_within)
summary(data_model$social_between)

summary(data_model$se_within)
summary(data_model$se_between)

required_vars <- c(
  "learning", "time", "grade", "school", "student_id",
  "participation_within", "participation_between",
  "social_within", "social_between",
  "se_within", "se_between"
)

setdiff(required_vars, names(data_model))

model_main <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 | student_id),
  data = data_model,
  REML = TRUE
)

summary(model_main)

confint(model_main, method = "Wald")

r2(model_main)
icc(model_main)

model_school_random <- lmer(
  learning ~ time + grade +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 | school) +
    (1 | student_id),
  data = data_model,
  REML = TRUE
)

summary(model_school_random)

isSingular(model_school_random)

VarCorr(model_school_random)

model_rs_participation <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 + participation_within | student_id),
  data = data_model,
  REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(model_rs_participation)
isSingular(model_rs_participation)

model_rs_social <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 + social_within | student_id),
  data = data_model,
  REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(model_rs_social)
isSingular(model_rs_social)

model_rs_se <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 + se_within | student_id),
  data = data_model,
  REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(model_rs_se)
isSingular(model_rs_se)

model_main_ml <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 | student_id),
  data = data_model,
  REML = FALSE
)

model_rs_all_ml <- lmer(
  learning ~ time + grade + school +
    participation_within + participation_between +
    social_within + social_between +
    se_within + se_between +
    (1 + participation_within + social_within + se_within || student_id),
  data = data_model,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(model_rs_all_ml)
isSingular(model_rs_all_ml)

anova(model_main_ml, model_rs_all_ml)

model_rs_all_reml <- update(model_rs_all_ml, REML = TRUE)

summary(model_rs_all_reml)
r2(model_rs_all_reml)
icc(model_rs_all_reml)

model_rs_all_no_int_ml <- model_rs_all_ml

model_rs_all_int_ml <- lmer(
  learning ~ time * participation_within +
    time * social_within +
    time * se_within +
    grade + school +
    participation_between +
    social_between +
    se_between +
    (1 + participation_within + social_within + se_within || student_id),
  data = data_model,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(model_rs_all_int_ml)
isSingular(model_rs_all_int_ml)

anova(model_rs_all_no_int_ml, model_rs_all_int_ml)

model_final <- update(model_rs_all_int_ml, REML = TRUE)

summary(model_final)
r2(model_final)
icc(model_final)
confint(model_final, method = "Wald")

table6 <- tidy(
  model_final,
  effects = "fixed",
  conf.int = TRUE,
  conf.method = "Wald"
)

table6_clean <- table6 %>%
  select(
    term,
    estimate,
    std.error,
    conf.low,
    conf.high,
    statistic,
    df,
    p.value
  ) %>%
  rename(
    Predictor = term,
    B = estimate,
    SE = std.error,
    `95% CI Lower` = conf.low,
    `95% CI Upper` = conf.high,
    t = statistic,
    df = df,
    p = p.value
  )

table6_clean <- table6_clean %>%
  mutate(
    B = round(B, 3),
    SE = round(SE, 3),
    `95% CI Lower` = round(`95% CI Lower`, 3),
    `95% CI Upper` = round(`95% CI Upper`, 3),
    t = round(t, 2),
    df = round(df, 0),
    p = case_when(
      p < .001 ~ "< .001",
      TRUE ~ sprintf("%.3f", p)
    )
  )

table6_clean <- table6_clean %>%
  mutate(
    Predictor = case_when(
      Predictor == "(Intercept)" ~ "Intercept",
      Predictor == "timeT2" ~ "Time (T2 vs T1)",
      Predictor == "timeT3" ~ "Time (T3 vs T1)",
      Predictor == "gradeG8" ~ "Grade 8 (vs 7)",
      Predictor == "gradeG9" ~ "Grade 9 (vs 7)",
      Predictor == "gradeG10" ~ "Grade 10 (vs 7)",
      Predictor == "schoolG" ~ "School B (vs A)",
      Predictor == "schoolH" ~ "School C (vs A)",
      Predictor == "schoolL" ~ "School D (vs A)",
      Predictor == "participation_within" ~ "Participation (within)",
      Predictor == "participation_between" ~ "Participation (between)",
      Predictor == "social_within" ~ "Social climate (within)",
      Predictor == "social_between" ~ "Social climate (between)",
      Predictor == "se_within" ~ "Self-efficacy (within)",
      Predictor == "se_between" ~ "Self-efficacy (between)",
      Predictor == "timeT2:participation_within" ~ "T2 × Participation (within)",
      Predictor == "timeT3:participation_within" ~ "T3 × Participation (within)",
      Predictor == "timeT2:social_within" ~ "T2 × Social climate (within)",
      Predictor == "timeT3:social_within" ~ "T3 × Social climate (within)",
      Predictor == "timeT2:se_within" ~ "T2 × Self-efficacy (within)",
      Predictor == "timeT3:se_within" ~ "T3 × Self-efficacy (within)",
      TRUE ~ Predictor
    )
  )

table6_clean
print(table6_clean, n = Inf)

lm_aux <- lm(
  learning ~ time * participation_within +
    time * social_within +
    time * se_within +
    grade + school +
    participation_between +
    social_between +
    se_between,
  data = data_model
)

vif_results <- car::vif(lm_aux)

vif_results