#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 800


library(shiny)
library(latex2exp)

# --- Funzione di Plot ---
plot_hypothesis_comparison <- function(d_eff, n1, n2, s1, s2, alpha = 0.05, 
                                       tails = "two", dist_type = "t", 
                                       obs_stat = NULL, x_range = NULL, 
                                       is_std = FALSE, show_ci = FALSE,
                                       show_plots = c("H0", "H1")) {
    
    # Logica Cohen's d: se is_std è TRUE, forziamo medie e SD
    if(is_std) {
        curr_s1 <- 1
        curr_s2 <- 1
        curr_d_eff <- d_eff # Qui d_eff è interpretato come Cohen's d
    } else {
        curr_s1 <- s1
        curr_s2 <- s2
        curr_d_eff <- d_eff
    }
    
    s_pooled_sq <- ((n1 - 1) * curr_s1^2 + (n2 - 1) * curr_s2^2) / (n1 + n2 - 2)
    sp <- sqrt(s_pooled_sq)
    se <- sp * sqrt(1/n1 + 1/n2)
    df_val <- n1 + n2 - 2
    
    m0 <- 0
    m1 <- curr_d_eff 
    curr_se <- se
    curr_obs <- if(!is.null(obs_stat) && !is.na(obs_stat)) obs_stat else NULL
    
    dfunc <- if(dist_type == "z") dnorm else function(x, m, s) dt((x-m)/s, df_val)/s
    pfunc <- if(dist_type == "z") pnorm else function(q, m, s) pt((q-m)/s, df_val)
    qfunc <- if(dist_type == "z") qnorm else function(p, m, s) qt(p, df_val)*s + m
    
    # Soglie critiche su H0
    crit_low <- -Inf; crit_high <- Inf
    if (tails == "two") {
        crit_low  <- qfunc(alpha/2, m0, curr_se)
        crit_high <- qfunc(1 - alpha/2, m0, curr_se)
    } else if (tails == "right") {
        crit_high <- qfunc(1 - alpha, m0, curr_se)
    } else {
        crit_low  <- qfunc(alpha, m0, curr_se)
    }
    
    if (is.null(x_range)) {
        vals <- c(m0, m1, crit_low[is.finite(crit_low)], crit_high[is.finite(crit_high)])
        if(!is.null(curr_obs)) vals <- c(vals, curr_obs)
        x_range <- c(min(vals) - 4*curr_se, max(vals) + 4*curr_se)
    }
    
    n_active <- length(show_plots)
    if(n_active == 0) return(NULL)
    
    old_par <- par(mfrow = c(n_active, 1), mar = c(4, 4, 3, 2))
    on.exit(par(old_par))
    
    # --- 1. GRAFICO H0 ---
    if("H0" %in% show_plots) {
        y_max_h0 <- dfunc(m0, m0, curr_se)
        curve(dfunc(x, m0, curr_se), from = x_range[1], to = x_range[2], 
              ylab = "Densità", main = TeX(r"($H_0$)"), xaxt = "n", xlab = "")
        
        if(is.finite(crit_low)) {
            x_a <- seq(x_range[1], crit_low, length = 100)
            polygon(c(x_a, crit_low), c(dfunc(x_a, m0, curr_se), 0), density = 20, angle = 45, col = "black")
        }
        if(is.finite(crit_high)) {
            x_a <- seq(crit_high, x_range[2], length = 100)
            polygon(c(crit_high, x_a), c(0, dfunc(x_a, m0, curr_se)), density = 20, angle = 45, col = "black")
        }
        
        if (!is.null(curr_obs)) {
            y_obs <- dfunc(curr_obs, m0, curr_se)
            if (tails == "right" || (tails == "two" && curr_obs > m0)) {
                px <- seq(abs(curr_obs), x_range[2], length = 100)
                polygon(c(abs(curr_obs), px, x_range[2]), c(0, dfunc(px, m0, curr_se), 0), col = rgb(1,0,0,0.3), border = NA)
            }
            if (tails == "left" || (tails == "two" && curr_obs < m0)) {
                px <- seq(x_range[1], -abs(curr_obs), length = 100)
                polygon(c(x_range[1], px, -abs(curr_obs)), c(0, dfunc(px, m0, curr_se), 0), col = rgb(1,0,0,0.3), border = NA)
            }
            points(curr_obs, 0, col = "red", pch = 19, cex = 1.3)
            pv_val <- if(tails == "right") 1 - pfunc(curr_obs, m0, curr_se) else 
                if(tails == "left") pfunc(curr_obs, m0, curr_se) else 2 * (1 - pfunc(abs(curr_obs), m0, curr_se))
            text(curr_obs, y_obs, labels = paste("p =", round(pv_val, 4)), pos = 3, col = "red", font = 2)
            
            if (show_ci) {
                y_pos_h0 <- y_max_h0 * 0.20
                z_crit <- if(tails == "two") qfunc(1 - alpha/2, 0, curr_se) else qfunc(1 - alpha, 0, curr_se)
                if (tails == "two") {
                    ci_l <- curr_obs - z_crit; ci_h <- curr_obs + z_crit
                    lines(c(ci_l, ci_h), c(y_pos_h0, y_pos_h0), col = "darkred", lwd = 3)
                } else if (tails == "right") {
                    ci_l <- curr_obs - z_crit
                    lines(c(ci_l, x_range[2]-curr_se), c(y_pos_h0, y_pos_h0), col = "darkred", lwd = 3)
                    lines(c(x_range[2]-curr_se, x_range[2]), c(y_pos_h0, y_pos_h0), col = "darkred", lwd = 3, lty = 3)
                } else {
                    ci_h <- curr_obs + z_crit
                    lines(c(x_range[1]+curr_se, ci_h), c(y_pos_h0, y_pos_h0), col = "darkred", lwd = 3)
                    lines(c(x_range[1], x_range[1]+curr_se), c(y_pos_h0, y_pos_h0), col = "darkred", lwd = 3, lty = 3)
                }
            }
        }
        abline(v = c(m0, crit_low[is.finite(crit_low)], crit_high[is.finite(crit_high)]), lty = c(1, 3, 3), col = c("black", "blue", "blue"))
        axis(1)
    }
    
    # --- 2. GRAFICO H1 ---
    if("H1" %in% show_plots) {
        y_max_h1 <- dfunc(m1, m1, curr_se)
        curve(dfunc(x, m1, curr_se), from = x_range[1], to = x_range[2], 
              ylab = "Densità", xlab = if(is_std) "Cohen's d" else "Differenza tra Medie", main = TeX(r"($H_1$)"))
        
        xb_min <- max(x_range[1], if(is.finite(crit_low)) crit_low else x_range[1])
        xb_max <- min(x_range[2], if(is.finite(crit_high)) crit_high else x_range[2])
        if(xb_min < xb_max) polygon(c(xb_min, seq(xb_min, xb_max, length=100), xb_max), c(0, dfunc(seq(xb_min, xb_max, length=100), m1, curr_se), 0), col = "lightblue", border = NA)
        if(is.finite(crit_low)) polygon(c(x_range[1], seq(x_range[1], crit_low, length=100), crit_low), c(0, dfunc(seq(x_range[1], crit_low, length=100), m1, curr_se), 0), col = "lightgreen", border = NA)
        if(is.finite(crit_high)) polygon(c(crit_high, seq(crit_high, x_range[2], length=100), x_range[2]), c(0, dfunc(seq(crit_high, x_range[2], length=100), m1, curr_se), 0), col = "lightgreen", border = NA)
        
        if (show_ci && !is.null(curr_obs)) {
            y_pos_h1 <- y_max_h1 * 0.20
            z_crit <- if(tails == "two") qfunc(1 - alpha/2, 0, curr_se) else qfunc(1 - alpha, 0, curr_se)
            if (tails == "two") {
                ci_l <- curr_obs - z_crit; ci_h <- curr_obs + z_crit
                lines(c(ci_l, ci_h), c(y_pos_h1, y_pos_h1), col = "darkred", lwd = 3)
            } else if (tails == "right") {
                ci_l <- curr_obs - z_crit
                lines(c(ci_l, x_range[2]-curr_se), c(y_pos_h1, y_pos_h1), col = "darkred", lwd = 3)
                lines(c(x_range[2]-curr_se, x_range[2]), c(y_pos_h1, y_pos_h1), col = "darkred", lwd = 3, lty = 3)
            } else {
                ci_h <- curr_obs + z_crit
                lines(c(x_range[1]+curr_se, ci_h), c(y_pos_h1, y_pos_h1), col = "darkred", lwd = 3)
                lines(c(x_range[1], x_range[1]+curr_se), c(y_pos_h1, y_pos_h1), col = "darkred", lwd = 3, lty = 3)
            }
        }
        abline(v = c(m1, crit_low[is.finite(crit_low)], crit_high[is.finite(crit_high)]), lty = c(1, 3, 3), col = c("black", "blue", "blue"))
        legend("topright", legend = c("Beta (Err II)", "Power (1-Beta)"), fill = c("lightblue", "lightgreen"), bty = "n")
    }
}

# --- UI ---
ui <- fluidPage(
    titlePanel("Test per due campioni indipendenti"),
    sidebarLayout(
        sidebarPanel(
            checkboxGroupInput("show_plots", "Visualizza Grafici:",
                               choices = list("Ipotesi Nulla (H0)" = "H0", 
                                              "Ipotesi Alternativa (H1)" = "H1"),
                               selected = c("H0", "H1")),
            hr(),
            checkboxInput("is_std", "Modalità Cohen's d (Mean=0, SD=1)", FALSE),
            checkboxInput("show_ci", "Visualizza IC nel Grafico", TRUE),
            hr(),
            radioButtons("dist", "Distribuzione:", choices = list("Z" = "z", "T" = "t"), selected = "t"),
            radioButtons("tails", "Direzione Test:", choices = list("Bilaterale (≠)" = "two", "Destro (>)" = "right", "Sinistro (<)" = "left")),
            numericInput("alpha", "Alpha:", value = 0.05, step = 0.01),
            hr(),
            fluidRow(column(6, numericInput("n1", "n1:", 30)), column(6, numericInput("n2", "n2:", 30))),
            
            # Pannello condizionale per SD: scompare se Cohen's d è attivo
            conditionalPanel(
                condition = "input.is_std == false",
                fluidRow(column(6, numericInput("s1", "SD 1:", 10)), column(6, numericInput("s2", "SD 2:", 10)))
            ),
            
            numericInput("d_eff", "Differenza Attesa / d:", value = 5),
            numericInput("obs_stat", "Differenza Osservata / d_obs:", value = 4.5),
            hr(),
            checkboxInput("fix_axes", "Blocca Assi X", FALSE),
            conditionalPanel("input.fix_axes", 
                             fluidRow(column(6, numericInput("xmin", "Min:", -30)), column(6, numericInput("xmax", "Max:", 30))))
        ),
        mainPanel(
            plotOutput("hypPlot", height = "580px"),
            hr(),
            tableOutput("summaryTable")
        )
    )
)

# --- Server ---
server <- function(input, output) {
    calc_data <- reactive({
        # Gestione SD per il calcolo SE
        curr_s1 <- if(input$is_std) 1 else input$s1
        curr_s2 <- if(input$is_std) 1 else input$s2
        
        s_pooled_sq <- ((input$n1 - 1) * curr_s1^2 + (input$n2 - 1) * curr_s2^2) / (input$n1 + input$n2 - 2)
        sp <- sqrt(s_pooled_sq)
        df_val <- input$n1 + input$n2 - 2
        se <- sp * sqrt(1/input$n1 + 1/input$n2)
        
        qfunc <- if(input$dist == "z") qnorm else function(p) qt(p, df_val)
        pfunc <- if(input$dist == "z") pnorm else function(q) pt(q, df_val)
        
        if(input$tails == "two") {
            c_low <- qfunc(input$alpha/2) * se
            c_high <- qfunc(1 - input$alpha/2) * se
            beta <- pfunc(c_high/se - input$d_eff/se) - pfunc(c_low/se - input$d_eff/se)
        } else if(input$tails == "right") {
            c_high <- qfunc(1 - input$alpha) * se
            beta <- pfunc(c_high/se - input$d_eff/se)
        } else {
            c_low <- qfunc(input$alpha) * se
            beta <- 1 - pfunc(c_low/se - input$d_eff/se)
        }
        
        pv <- NA
        ci_txt <- "-"
        if(!is.na(input$obs_stat)) {
            t_stat <- input$obs_stat / se
            pv <- if(input$tails == "right") 1 - pfunc(t_stat) else 
                if(input$tails == "left") pfunc(t_stat) else 2 * (1 - pfunc(abs(t_stat)))
            
            z_crit_ci <- if(input$tails == "two") qfunc(1 - input$alpha/2) else qfunc(1 - input$alpha)
            if(input$tails == "two") {
                ci_txt <- paste0("[", round(input$obs_stat - z_crit_ci*se, 3), " ; ", round(input$obs_stat + z_crit_ci*se, 3), "]")
            } else if(input$tails == "right") {
                ci_txt <- paste0("[", round(input$obs_stat - z_crit_ci*se, 3), " ; +∞)")
            } else {
                ci_txt <- paste0("(-∞ ; ", round(input$obs_stat + z_crit_ci*se, 3), "]")
            }
        }
        
        list(sp=sp, se=se, df=df_val, beta=beta, pwr=1-beta, pv=pv, ci_txt=ci_txt,
             d_theo=input$d_eff/sp, d_obs=input$obs_stat/sp)
    })
    
    output$hypPlot <- renderPlot({
        req(input$show_plots)
        plot_hypothesis_comparison(input$d_eff, input$n1, input$n2, input$s1, input$s2,
                                   input$alpha, input$tails, input$dist, input$obs_stat,
                                   if(input$fix_axes) c(input$xmin, input$xmax) else NULL,
                                   is_std = input$is_std, show_ci = input$show_ci,
                                   show_plots = input$show_plots)
    })
    
    output$summaryTable <- renderTable({
        res <- calc_data()
        df_res <- data.frame(
            Descrizione = c("Errore Standard (SE)", "Power (1-Beta)", "Errore II tipo (Beta)", 
                            "Cohen's d [vero]", "Cohen's d [osservato]"),
            Valore = c(round(res$se, 4), round(res$pwr, 4), round(res$beta, 4),
                       round(res$d_theo, 4), if(!is.na(res$d_obs)) round(res$d_obs, 4) else "-")
        )
        
        if(!is.na(input$obs_stat)) {
            df_res <- rbind(df_res, 
                            data.frame(Descrizione = paste0("IC ", (1-input$alpha)*100, "%"), Valore = res$ci_txt),
                            data.frame(Descrizione = "P-Value", Valore = as.character(round(res$pv, 5))))
        }
        df_res
    }, striped = TRUE, bordered = TRUE)
}

shinyApp(ui, server)