
###################################################.
# OLM and RMA regression 
# Cicla tra tutte le colonne di AE (x) e tutte le colonne di TE (y), 
# calcola sia la regressione OLS che la RMA (Reduced Major Axis), e salva in un data frame tutte le statistiche rilevanti.
# L’obiettivo sarà:
# OLS: slope, intercept, R², p-value, RMSE, ecc.
# RMA: slope, intercept, R² (e magari anche deviazione dalla slope=1).
# Usare la colonna TIMESTAMP_START per riconoscere data e ora.
# Fare le regressioni sui dati originali (semiorari).
# Aggregare giornalmente (medie o somme, a seconda del tipo di variabile), e rifare le regressioni sui dati giornalieri.
# Salvare in un unico data frame sia i risultati semiorari che giornalieri, aggiungendo una colonna t_scale_temporale che indica il livello di aggregazione.
# LOOP: 
# Ogni dataset contiene le colonne TIMESTAMP_START, variabili AE e variabili TE.
# Il loop:
# Converte il timestamp
# Separa AE e TE
# Fa regressioni semiorarie e giornaliere con la funzione EBC_calc()
# Salva i risultati in un unico risultati_finali con colonna stazione_id.

# per ogni stazione:
# Fare le regressioni OLS + RMA non solo sull’intero dataset, ma anche all’interno di ogni livello di ciascun fattore presente in proxies.
# Ogni fattore ha 2+ livelli (stagioni, classi di RH, stabilità atmosferica, ecc.).
# L’output deve includere anche il nome del fattore e il livello su cui è stata fatta la regressione.


# FUNZIONE  ----------------------------

library(lmodel2)
library(broom)
library(dplyr)
library(lubridate)

# ---- Funzione per regressione ----
calc_reg <- function(x, y) { # EBC regressions function
  
  # OLS
  model <- lm(y ~ x)
  ols_stats <- summary(model)
  ols_slope <- ols_stats$coefficients[2,1]
  ols_intercept <- ols_stats$coefficients[1,1]
  ols_r2 <- ols_stats$r.squared
  ols_diff <- abs(ols_slope - 1)
  # TODO?
  # ols_slope_SE = ols_stats$coefficients[2,2],
  # ols_slope_p = ols_stats$coefficients[2,4],
  # ols_rmse = sqrt(mean(residuals(ols_fit)^2)),
  
  # OLS, intercept == 0
  model_i0 <- lm(y ~ x + 0)
  ols_stats_i0 <- summary(model_i0)
  ols_slope_i0 <- ols_stats_i0$coefficients[1]
  ols_r2_i0 <- ols_stats_i0$r.squared
  ols_diff_i0 <- abs(ols_slope_i0 - 1)
  
  # RMA
  rma <- lmodel2(y ~ x, range.x = "interval", range.y = "interval")
  rma_slope <- rma$'regression.results'[rma$'regression.results'$'Method' == "RMA", "Slope"]
  rma_r2 <- rma$rsquare
  rma_diff <- abs(rma_slope - 1)
  
  return(data.frame(
    OLS_slope = ols_slope,
    OLS_r2 = ols_r2,
    OLS_intercept = ols_intercept,
    OLS_slope_d1 = ols_diff,
    OLS_slope_i0 = ols_slope_i0,
    OLS_r2_i0 = ols_r2_i0,
    OLS_slope_d1_i0 = ols_diff_i0,
    RMA_slope = rma_slope,
    RMA_r2 = rma_r2,
    RMA_slope_d1 = rma_diff
  ))
}



