library(tidyverse)
library(gt)
library(plotly)
library(htmltools)

set.seed(1234)

# -----------------------------------
# 1. CREEM UNA POBLACIÓ "REAL"
# -----------------------------------

N <- 20000

poblacio <- data.frame(
  id = 1:N,
  sexe = sample(c("Noi", "Noia"), N, replace = TRUE, prob = c(0.51, 0.49)),
  etapa = sample(c("ESO", "Batx", "FP"), N, replace = TRUE, prob = c(0.60, 0.25, 0.15)),
  centre = sample(paste0("Centre_", 1:20), N, replace = TRUE)
)

# Variable latent relacionada amb l'estat emocional
poblacio$malestar <- with(
  poblacio,
  50 +
    ifelse(sexe == "Noi", 6, 0) +
    ifelse(etapa == "FP", 3, ifelse(etapa == "Batx", 1, 0)) +
    rnorm(N, mean = 0, sd = 10)
)

# Limitem per evitar valors massa extrems
poblacio$malestar <- pmin(pmax(poblacio$malestar, 0), 100)

# -----------------------------------
# 2. PARÀMETRES REALS DE LA POBLACIÓ
# -----------------------------------

mitjana_real <- mean(poblacio$malestar)

cat("Mitjana real de malestar a la població:", round(mitjana_real, 2), "\n\n")

prop_sexe <- prop.table(table(poblacio$sexe))
prop_etapa <- prop.table(table(poblacio$etapa))

print(prop_sexe)
print(prop_etapa)
cat("\n")

# -----------------------------------
# 3. PREPAREM ELS PARÀMETRES DE RESPOSTA
# -----------------------------------

# Efecte fix dels centres (el mantenim constant a totes les iteracions)
efecte_centre <- rnorm(20, 0, 0.7)
names(efecte_centre) <- paste0("Centre_", 1:20)

# Model logístic de resposta
linpred <- -4 +
  0.06 * poblacio$malestar +
  efecte_centre[poblacio$centre]

prob_resposta <- 1 / (1 + exp(-linpred))

# Estrats poblacionals per a la postestratificació
poblacio$estrat <- interaction(poblacio$sexe, poblacio$etapa, drop = TRUE)
dist_real <- prop.table(table(poblacio$estrat))
estrats_poblacio <- names(dist_real)

# -----------------------------------
# 4. MONTE CARLO
# -----------------------------------

B <- 1000
n_sub <- 1000

resultats_mc <- vector("list", B)

for (b in 1:B) {
  
  # Simulem qui respon voluntàriament en aquesta iteració
  respon <- rbinom(N, 1, prob_resposta)
  enquesta_voluntaria <- poblacio[respon == 1, ]
  
  # Per seguretat, si hi ha massa poques respostes, saltem la iteració
  if (nrow(enquesta_voluntaria) < n_sub) {
    resultats_mc[[b]] <- data.frame(
      iteracio = b,
      metode = c("Enquesta voluntària",
                 "Submostra aleatòria de l'enquesta voluntària",
                 "Submostra corregida (pesos) de l'enquesta voluntària"),
      mitjana = NA_real_
    )
    next
  }
  
  # Mitjana de l'enquesta voluntària
  mitjana_ev <- mean(enquesta_voluntaria$malestar)
  
  # Submostra aleatòria dins l'enquesta voluntària
  idx_sub <- sample(seq_len(nrow(enquesta_voluntaria)), n_sub, replace = FALSE)
  submostra_aleatoria <- enquesta_voluntaria[idx_sub, ]
  
  mitjana_sub <- mean(submostra_aleatoria$malestar)
  
  # Postestratificació per sexe x etapa
  submostra_aleatoria$estrat <- interaction(
    submostra_aleatoria$sexe,
    submostra_aleatoria$etapa,
    drop = TRUE
  )
  
  dist_sub <- prop.table(table(submostra_aleatoria$estrat))
  
  # Alineem estrats
  estrats <- union(estrats_poblacio, names(dist_sub))
  dr <- dist_real[match(estrats, names(dist_real))]
  ds <- dist_sub[match(estrats, names(dist_sub))]
  
  dr[is.na(dr)] <- 0
  ds[is.na(ds)] <- 0
  
  pesos_estrat <- dr / ds
  pesos_estrat[!is.finite(pesos_estrat)] <- NA
  
  submostra_aleatoria$pes <- pesos_estrat[match(submostra_aleatoria$estrat, estrats)]
  
  mitjana_ponderada <- with(
    submostra_aleatoria,
    weighted.mean(malestar, w = pes, na.rm = TRUE)
  )
  
  resultats_mc[[b]] <- data.frame(
    iteracio = b,
    metode = c("Enquesta voluntària",
               "Submostra aleatòria de l'enquesta voluntària",
               "Submostra corregida (pesos) de l'enquesta voluntària"),
    mitjana = c(mitjana_ev, mitjana_sub, mitjana_ponderada)
  )
}

resultats_mc <- bind_rows(resultats_mc)

# Eliminem files NA si n'hi hagués
resultats_mc <- resultats_mc %>%
  filter(!is.na(mitjana))

# -----------------------------------
# 5. TAULA RESUM DEL MONTE CARLO
# -----------------------------------

taula_ic_presentacio <- resultats_mc %>%
  group_by(metode) %>%
  summarise(
    n_iteracions = n(),
    mitjana_mc = mean(mitjana, na.rm = TRUE),
    li_2.5 = quantile(mitjana, 0.025, na.rm = TRUE),
    p50 = quantile(mitjana, 0.50, na.rm = TRUE),
    ls_97.5 = quantile(mitjana, 0.975, na.rm = TRUE),
    sd = sd(mitjana, na.rm = TRUE),
    biaix = mean(mitjana, na.rm = TRUE) - mitjana_real,
    error_absolut = abs(mean(mitjana, na.rm = TRUE) - mitjana_real),
    .groups = "drop"
  ) %>%
  mutate(
    mitjana_real = mitjana_real,
    IC_95 = paste0("[", round(li_2.5, 2), ", ", round(ls_97.5, 2), "]")
  ) %>%
  select(
    metode,
    n_iteracions,
    mitjana_real,
    mitjana_mc,
    biaix,
    error_absolut,
    sd,
    li_2.5,
    p50,
    ls_97.5,
    IC_95
  )

# Afegim una fila de referència per a la població real
fila_poblacio_real <- tibble(
  metode = "Població real",
  n_iteracions = NA_integer_,
  mitjana_real = mitjana_real,
  mitjana_mc = mitjana_real,
  biaix = 0,
  error_absolut = 0,
  sd = NA_real_,
  li_2.5 = NA_real_,
  p50 = mitjana_real,
  ls_97.5 = NA_real_,
  IC_95 = "-"
)

taula_ic_presentacio <- bind_rows(fila_poblacio_real, taula_ic_presentacio)

# Taula GT
taula_gt <- taula_ic_presentacio %>%
  gt() %>%
  tab_header(
    title = md("**Monte Carlo de les estimacions de la mitjana de malestar**"),
    subtitle = paste0(
      "1000 iteracions. La fila 'Població real' és la referència de comparació."
    )
  ) %>%
  cols_label(
    metode = "Mètode",
    n_iteracions = "N iter.",
    mitjana_real = "Mitjana real",
    mitjana_mc = "Mitjana estimada",
    biaix = "Biaix",
    error_absolut = "Error absolut",
    sd = "Desv. estàndard",
    li_2.5 = "P2,5",
    p50 = "P50",
    ls_97.5 = "P97,5",
    IC_95 = "Interval 95%"
  ) %>%
  fmt_number(
    columns = c(mitjana_real, mitjana_mc, biaix, error_absolut, sd, li_2.5, p50, ls_97.5),
    decimals = 2
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  # Estil per a la fila de la Població Real (tota la fila)
  tab_style(
    style = list(
      cell_fill(color = "#EAF2F8"),
      cell_text(weight = "bold")
    ),
    locations = cells_body(
      rows = metode == "Població real"
    )
  ) %>%
  # Estil per a la columna Mitjana Real (tota la columna)
  tab_style(
    style = list(
      cell_fill(color = "#EAF2F8"),
      cell_text(weight = "bold")
    ),
    locations = cells_body(
      columns = mitjana_real
    )
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels(everything())
  )

taula_gt

# ---------------------------------------------------------
# 1. EL MIRALL DEFORMAT: DENSITATS (REAL VS VOLUNTÀRIA)
# ---------------------------------------------------------

# Creem un dataframe per a les densitats
densitat_real <- density(poblacio$malestar)
# Agafem una iteració qualsevol de l'enquesta voluntària (la darrera de l'experiment)
enquesta_v <- poblacio[rbinom(N, 1, prob_resposta) == 1, ]
densitat_vol <- density(enquesta_v$malestar)

fig1 <- plot_ly() %>%
  add_lines(x = densitat_real$x, y = densitat_real$y, name = "Població real (la veritat)", 
            line = list(color = '#2E86C1', width = 3), fill = 'tozeroy', alpha = 0.2) %>%
  add_lines(x = densitat_vol$x, y = densitat_vol$y, name = "Enquesta voluntària (el biaix)", 
            line = list(color = '#CB4335', width = 3, dash = 'dot'), fill = 'tozeroy', alpha = 0.2) %>%
  layout(
    title = "El mirall deformat: com l'autoselecció desplaça la realitat",
    xaxis = list(title = "Nivell de malestar (0-100)"),
    yaxis = list(title = "Densitat"),
    annotations = list(
      list(x = mitjana_real, y = max(densitat_real$y), 
           text = "Centre real", showarrow = T, arrowhead = 2,
           ax = -40, ay = -40), # Fletxa cap a l'esquerra
      
      list(x = mean(enquesta_v$malestar), y = max(densitat_vol$y), 
           text = "Centre enquesta voluntària", showarrow = T, arrowhead = 2,
           ax = 40, ay = -40)  # Fletxa cap a la dreta
    )
  )


# ---------------------------------------------------------
# 2. LA PONDERACIÓ MÀGICA: SCATTER PLOT AMB "FLETXA D'INTENT"
# ---------------------------------------------------------

# Agafem una mostra aleatòria pura de la població i la comparem amb la voluntària
mostra_pura <- poblacio[sample(1:N, 500), ]
mostra_vol  <- enquesta_v[sample(1:nrow(enquesta_v), 500), ]

# Calculem la mitjana ponderada de la voluntària (com fas al teu loop)
# Per simplicitat visual, usem el valor de la teva taula final (56.88 aprox)
mitjana_pond_v <- taula_ic_presentacio$mitjana_mc[taula_ic_presentacio$metode == "Submostra corregida (pesos) de l'enquesta voluntària"]

fig2 <- plot_ly() %>%
  # Punts població real
  add_markers(data = mostra_pura, x = ~id, y = ~malestar, name = "Mostra aleatòria real",
              marker = list(color = '#2E86C1', opacity = 0.3, size = 4)) %>%
  # Punts enquesta voluntària
  add_markers(data = mostra_vol, x = ~id, y = ~malestar, name = "Enquesta voluntària",
              marker = list(color = '#CB4335', opacity = 0.6, size = 5)) %>%
  # Línia real
  add_segments(x = 0, xend = N, y = mitjana_real, yend = mitjana_real, 
               name = "Mitjana real", line = list(color = 'blue', width = 2)) %>%
  # Línia ponderada (l'intent de rescat)
  add_segments(x = 0, xend = N, y = mitjana_pond_v, yend = mitjana_pond_v, 
               name = "Intent de ponderació", line = list(color = 'black', width = 2, dash = 'dash')) %>%
  layout(
    title = "La ponderació màgica: intentant moure una pared",
    xaxis = list(title = "Identificador individu sintètic"),
    yaxis = list(title = "Malestar"),
    annotations = list(
      list(
        x = N/2, y = mitjana_real,
        ax = N/2, ay = mean(enquesta_v$malestar),
        xref = "x", yref = "y", axref = "x", ayref = "y",
        text = "L'esforç inútil de ponderar",
        showarrow = T, arrowhead = 3, color = "black"
      )
    )
  )

fig1
fig2