library(shiny)
library(MASS)
library(ggplot2)
library(gridExtra)
# ============================================================
# FUNÇÃO: BUSCA DE RAÍZES
# ============================================================
find_roots <- function(f, x_min, x_max, n_points = 1000, tol = 1e-6) {
x_grid <- seq(x_min, x_max, length.out = n_points)
y_grid <- f(x_grid)
roots <- numeric(0)
for (i in 1:(length(x_grid) - 1)) {
if (y_grid[i] * y_grid[i + 1] < 0) {
root <- tryCatch(
uniroot(f, c(x_grid[i], x_grid[i + 1]), tol = tol)$root,
error = function(e) NA
)
if (!is.na(root)) roots <- c(roots, root)
}
}
if (length(roots) > 1) {
roots <- roots[order(roots)]
roots <- roots[c(TRUE, diff(roots) > tol * 10)]
}
return(roots)
}
# ============================================================
# FUNÇÃO: WALD POST-HOC
# ============================================================
wald_posthoc <- function(res1, res2 = NULL, adjust = "bonferroni") {
dias_list <- list()
labels <- character(0)
if (!is.null(res1$result$all_dias)) {
for (i in seq_along(res1$result$all_dias)) {
dias_list[[length(dias_list) + 1]] <- res1$result$all_dias[[i]]
labels <- c(labels, paste0("G1-DIAS", i))
}
} else if (!is.null(res1$result$dias_x)) {
dias_list[[length(dias_list) + 1]] <- list(
dias_x = res1$result$dias_x,
dias_y_original = res1$result$dias_y_original,
cov_xy = res1$result$cov_xy
)
labels <- c(labels, "G1")
}
if (!is.null(res2)) {
if (!is.null(res2$result$all_dias)) {
for (i in seq_along(res2$result$all_dias)) {
dias_list[[length(dias_list) + 1]] <- res2$result$all_dias[[i]]
labels <- c(labels, paste0("G2-DIAS", i))
}
} else if (!is.null(res2$result$dias_x)) {
dias_list[[length(dias_list) + 1]] <- list(
dias_x = res2$result$dias_x,
dias_y_original = res2$result$dias_y_original,
cov_xy = res2$result$cov_xy
)
labels <- c(labels, "G2")
}
}
n <- length(dias_list)
if (n < 2) return(NULL)
pairs <- combn(n, 2)
ncomp <- ncol(pairs)
results <- data.frame(
Comparison = character(ncomp),
W = numeric(ncomp),
df = integer(ncomp),
p_raw = numeric(ncomp),
D2 = numeric(ncomp),
stringsAsFactors = FALSE
)
for (i in seq_len(ncomp)) {
idx1 <- pairs[1, i]
idx2 <- pairs[2, i]
d1 <- dias_list[[idx1]]
d2 <- dias_list[[idx2]]
mu1 <- c(d1$dias_x, d1$dias_y_original)
mu2 <- c(d2$dias_x, d2$dias_y_original)
cov1 <- d1$cov_xy
cov2 <- d2$cov_xy
delta <- mu1 - mu2
cov_combined <- cov1 + cov2
W <- as.numeric(t(delta) %*% solve(cov_combined) %*% delta)
p <- 1 - pchisq(W, 2)
D2 <- sqrt(W)
results$Comparison[i] <- paste0(labels[idx1], " vs ", labels[idx2])
results$W[i] <- W
results$df[i] <- 2
results$p_raw[i] <- p
results$D2[i] <- D2
}
if (adjust == "bonferroni") {
results$p_adj <- pmin(results$p_raw * ncomp, 1)
} else if (adjust %in% c("holm", "hochberg", "hommel", "BH", "BY")) {
results$p_adj <- p.adjust(results$p_raw, method = adjust)
} else {
results$p_adj <- results$p_raw
}
results$Significance <- ifelse(results$p_adj < 0.05, "SIGNIFICANT", "NOT SIGNIFICANT")
return(results)
}
wald_posthoc_text <- function(res1, res2 = NULL, adjust = "bonferroni") {
results <- wald_posthoc(res1, res2, adjust)
if (is.null(results)) return("No comparisons available (need at least 2 DIAS).")
out <- paste0(
"========================================\n",
" POST-HOC WALD TESTS (Pairwise Comparisons)\n",
" Adjustment: ", toupper(adjust), "\n",
"========================================\n\n"
)
for (i in 1:nrow(results)) {
out <- paste0(out,
sprintf("--- %s ---\n", results$Comparison[i]),
sprintf("W = %.3f, df = %d\n", results$W[i], results$df[i]),
sprintf("p (raw) = %.6f\n", results$p_raw[i]),
sprintf("p (adjusted) = %.6f\n", results$p_adj[i]),
sprintf("D² = %.3f (%s)\n", results$D2[i], ifelse(results$D2[i] < 0.5, "Small",
ifelse(results$D2[i] < 1, "Medium",
ifelse(results$D2[i] < 2, "Large", "Very Large")))),
sprintf("Result: %s\n\n", results$Significance[i])
)
}
out <- paste0(out, "========================================\n")
return(out)
}
# ============================================================
# MODELO LINEAR
# ============================================================
line <- function(x, intercept, slope) {
intercept + slope * x
}
dias_individual_linear <- function(a1, b1, a2, b2) {
if (abs(b1 - b2) < 1e-12) return(NA)
return((a2 - a1) / (b1 - b2))
}
dias_linear_text <- function(x, y1, y2, invert_curve1 = FALSE) {
if (invert_curve1) y1 <- 100 - y1
m1 <- lm(y1 ~ x)
m2 <- lm(y2 ~ x)
p1 <- coef(m1); s1 <- summary(m1)$coefficients[, 2]; c1 <- vcov(m1)
p2 <- coef(m2); s2 <- summary(m2)$coefficients[, 2]; c2 <- vcov(m2)
names(p1) <- c("intercept", "slope")
names(p2) <- c("intercept", "slope")
Cov_all <- matrix(0, 4, 4)
Cov_all[1:2, 1:2] <- c1; Cov_all[3:4, 3:4] <- c2
px <- dias_individual_linear(p1[1], p1[2], p2[1], p2[2])
if (is.na(px)) {
return(list(text = "Error: Lines are parallel. No intersection point.", result = NULL))
}
h <- 1e-6; J <- matrix(0, 2, 4); mt <- c(p1, p2)
for (i in 1:4) {
mp <- mt; mp[i] <- mp[i] + h
J[1, i] <- (dias_individual_linear(mp[1], mp[2], mp[3], mp[4]) - px) / h
}
py_current <- line(px, p1[1], p1[2])
for (i in 1:4) {
mp <- mt; mp[i] <- mp[i] + h
J[2, i] <- (line(px, mp[1], mp[2]) - py_current) / h
}
py_original <- if (invert_curve1) 100 - py_current else py_current
Cov_XY <- J %*% Cov_all %*% t(J)
sx <- sqrt(Cov_XY[1, 1]); sy_dias <- sqrt(Cov_XY[2, 2])
out <- paste0(
"===== FITTED PARAMETERS (Linear Regression) =====\n",
sprintf("Curve 1: Intercept = %.2f ± %.2f | Slope = %.4f ± %.4f\n", p1[1], s1[1], p1[2], s1[2]),
sprintf("Curve 2: Intercept = %.2f ± %.2f | Slope = %.4f ± %.4f\n\n", p2[1], s2[1], p2[2], s2[2]),
"===== DIAS (Linear) =====\n",
sprintf("DIAS-X: %.3f ± %.4f\n", px, sx),
sprintf("95%% CI (X): [%.3f, %.3f]\n", px - 1.96 * sx, px + 1.96 * sx),
sprintf("DIAS-Y: %.1f ± %.2f%%\n", py_original, sy_dias),
sprintf("95%% CI (Y): [%.1f, %.1f]%%\n\n", py_original - 1.96 * sy_dias, py_original + 1.96 * sy_dias),
"----- DIAS Covariance Matrix (X, Y) -----\n",
sprintf(" Var(DIAS-X) = %.6f Cov(DIAS-X, DIAS-Y) = %.4f\n", Cov_XY[1, 1], Cov_XY[1, 2]),
sprintf(" Cov(DIAS-X, DIAS-Y) = %.4f Var(DIAS-Y) = %.4f\n", Cov_XY[2, 1], Cov_XY[2, 2])
)
return(list(
text = out,
result = list(
dias_x = px, se_x = sx,
dias_y_original = py_original, se_y = sy_dias,
cov_xy = Cov_XY,
model1 = m1, model2 = m2,
x = x, y1 = y1, y2 = y2
)
))
}
# ------------------------------------------------------------
# MODELO SIGMOIDE 3P (find_roots)
# ------------------------------------------------------------
sigmoide3 <- function(x, top, bottom, lec50) {
bottom + (top - bottom) / (1 + 10^(lec50 - x))
}
dias_individual_3p <- function(top1, bot1, lec501, top2, bot2, lec502) {
diff_func <- function(x) {
sigmoide3(x, top1, bot1, lec501) - sigmoide3(x, top2, bot2, lec502)
}
x_center <- mean(c(lec501, lec502))
x_range <- 4
x_min <- min(c(lec501, lec502)) - x_range
x_max <- max(c(lec501, lec502)) + x_range
roots <- find_roots(diff_func, x_min, x_max)
if (length(roots) == 0) {
return(NA)
}
return(roots[1])
}
dias_3p_text <- function(x, y1, y2, invert_curve1 = FALSE,
chute_top1 = 100, chute_bottom1 = 0, chute_lec501 = -7,
chute_top2 = 100, chute_bottom2 = 0, chute_lec502 = -7) {
if (invert_curve1) y1 <- 100 - y1
m1 <- try(nls(y1 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
start = list(top = chute_top1, bottom = chute_bottom1, lec50 = chute_lec501),
control = nls.control(maxiter = 500, minFactor = 1e-10)), silent = TRUE)
m2 <- try(nls(y2 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
start = list(top = chute_top2, bottom = chute_bottom2, lec50 = chute_lec502),
control = nls.control(maxiter = 500, minFactor = 1e-10)), silent = TRUE)
if (inherits(m1, "try-error") || inherits(m2, "try-error")) {
return(list(text = "Error: Unable to fit the sigmoid curve.", result = NULL))
}
p1 <- coef(m1); s1 <- summary(m1)$coefficients[, 2]; c1 <- vcov(m1)
p2 <- coef(m2); s2 <- summary(m2)$coefficients[, 2]; c2 <- vcov(m2)
Cov_all <- matrix(0, 6, 6); Cov_all[1:3, 1:3] <- c1; Cov_all[4:6, 4:6] <- c2
px <- dias_individual_3p(p1[1], p1[2], p1[3], p2[1], p2[2], p2[3])
if (is.na(px)) {
return(list(text = "Error: Curves do not intersect within the search range. No DIAS calculated.", result = NULL))
}
h <- 1e-6; J <- matrix(0, 2, 6); mt <- c(p1, p2)
for (i in 1:6) { mp <- mt; mp[i] <- mp[i] + h; J[1, i] <- (dias_individual_3p(mp[1], mp[2], mp[3], mp[4], mp[5], mp[6]) - px) / h }
py_current <- sigmoide3(px, p1[1], p1[2], p1[3])
for (i in 1:6) { mp <- mt; mp[i] <- mp[i] + h; J[2, i] <- (sigmoide3(px, mp[1], mp[2], mp[3]) - py_current) / h }
py_original <- if (invert_curve1) 100 - py_current else py_current
Cov_XY <- J %*% Cov_all %*% t(J)
sx <- sqrt(Cov_XY[1, 1]); sy_dias <- sqrt(Cov_XY[2, 2])
out <- paste0(
"===== FITTED PARAMETERS (3-Parameter Sigmoid) =====\n",
sprintf("Curve 1: Top = %.2f ± %.2f | Bottom = %.2f ± %.2f | LogEC50 = %.3f ± %.4f\n", p1[1], s1[1], p1[2], s1[2], p1[3], s1[3]),
sprintf("Curve 2: Top = %.2f ± %.2f | Bottom = %.2f ± %.2f | LogEC50 = %.3f ± %.4f\n\n", p2[1], s2[1], p2[2], s2[2], p2[3], s2[3]),
"===== DIAS (3 Parameters) =====\n",
sprintf("DIAS-X: %.3f ± %.4f\n", px, sx),
sprintf("95%% CI (X): [%.3f, %.3f]\n", px - 1.96 * sx, px + 1.96 * sx),
sprintf("DIAS-Y: %.1f ± %.2f%%\n", py_original, sy_dias),
sprintf("95%% CI (Y): [%.1f, %.1f]%%\n\n", py_original - 1.96 * sy_dias, py_original + 1.96 * sy_dias),
"----- DIAS Covariance Matrix (X, Y) -----\n",
sprintf(" Var(DIAS-X) = %.6f Cov(DIAS-X, DIAS-Y) = %.4f\n", Cov_XY[1, 1], Cov_XY[1, 2]),
sprintf(" Cov(DIAS-X, DIAS-Y) = %.4f Var(DIAS-Y) = %.4f\n", Cov_XY[2, 1], Cov_XY[2, 2])
)
return(list(text = out, result = list(dias_x = px, se_x = sx, dias_y_original = py_original, se_y = sy_dias, cov_xy = Cov_XY,
model1 = m1, model2 = m2, x = x, y1 = y1, y2 = y2)))
}
# ============================================================
# MODELO SIGMOIDE 4P
# ============================================================
sigmoide4 <- function(x, top, bottom, lec50, hs) {
bottom + (top - bottom) / (1 + 10^((lec50 - x) * hs))
}
dias_individual_4p <- function(top1, bot1, lec501, hs1, top2, bot2, lec502, hs2) {
diff_func <- function(x) {
sigmoide4(x, top1, bot1, lec501, hs1) - sigmoide4(x, top2, bot2, lec502, hs2)
}
x_center <- mean(c(lec501, lec502))
x_range <- 4
x_min <- min(c(lec501, lec502)) - x_range
x_max <- max(c(lec501, lec502)) + x_range
roots <- find_roots(diff_func, x_min, x_max)
if (length(roots) == 0) {
return(NA)
}
return(roots[1])
}
dias_4p_text <- function(x, y1, y2, invert_curve1 = FALSE,
chute_top1 = 100, chute_bottom1 = 0, chute_lec501 = -7, chute_hs1 = 1,
chute_top2 = 100, chute_bottom2 = 0, chute_lec502 = -7, chute_hs2 = 1) {
if (invert_curve1) y1 <- 100 - y1
m1 <- try(nls(y1 ~ bottom + (top - bottom) / (1 + 10^((lec50 - x) * hs)),
start = list(top = chute_top1, bottom = chute_bottom1,
lec50 = chute_lec501, hs = chute_hs1),
control = nls.control(maxiter = 500, minFactor = 1e-10)), silent = TRUE)
m2 <- try(nls(y2 ~ bottom + (top - bottom) / (1 + 10^((lec50 - x) * hs)),
start = list(top = chute_top2, bottom = chute_bottom2,
lec50 = chute_lec502, hs = chute_hs2),
control = nls.control(maxiter = 500, minFactor = 1e-10)), silent = TRUE)
if (inherits(m1, "try-error") || inherits(m2, "try-error")) {
return(list(text = "Error: Unable to fit the sigmoid curve.", result = NULL))
}
p1 <- coef(m1); s1 <- summary(m1)$coefficients[, 2]; c1 <- vcov(m1)
p2 <- coef(m2); s2 <- summary(m2)$coefficients[, 2]; c2 <- vcov(m2)
Cov_all <- matrix(0, 8, 8); Cov_all[1:4, 1:4] <- c1; Cov_all[5:8, 5:8] <- c2
px <- dias_individual_4p(p1[1], p1[2], p1[3], p1[4], p2[1], p2[2], p2[3], p2[4])
if (is.na(px)) {
return(list(text = "Error: Curves do not intersect within the search range. No DIAS calculated.", result = NULL))
}
h <- 1e-6; J <- matrix(0, 2, 8); mt <- c(p1, p2)
for (i in 1:8) { mp <- mt; mp[i] <- mp[i] + h; J[1, i] <- (dias_individual_4p(mp[1], mp[2], mp[3], mp[4], mp[5], mp[6], mp[7], mp[8]) - px) / h }
py_current <- sigmoide4(px, p1[1], p1[2], p1[3], p1[4])
for (i in 1:8) { mp <- mt; mp[i] <- mp[i] + h; J[2, i] <- (sigmoide4(px, mp[1], mp[2], mp[3], mp[4]) - py_current) / h }
py_original <- if (invert_curve1) 100 - py_current else py_current
Cov_XY <- J %*% Cov_all %*% t(J)
sx <- sqrt(Cov_XY[1, 1]); sy_dias <- sqrt(Cov_XY[2, 2])
out <- paste0(
"===== FITTED PARAMETERS (4-Parameter Sigmoid) =====\n",
sprintf("Curve 1: Top = %.2f ± %.2f | Bottom = %.2f ± %.2f | LogEC50 = %.3f ± %.4f | HillSlope = %.4f ± %.4f\n", p1[1], s1[1], p1[2], s1[2], p1[3], s1[3], p1[4], s1[4]),
sprintf("Curve 2: Top = %.2f ± %.2f | Bottom = %.2f ± %.2f | LogEC50 = %.3f ± %.4f | HillSlope = %.4f ± %.4f\n\n", p2[1], s2[1], p2[2], s2[2], p2[3], s2[3], p2[4], s2[4]),
"===== DIAS (4 Parameters) =====\n",
sprintf("DIAS-X: %.3f ± %.4f\n", px, sx),
sprintf("95%% CI (X): [%.3f, %.3f]\n", px - 1.96 * sx, px + 1.96 * sx),
sprintf("DIAS-Y: %.1f ± %.2f%%\n", py_original, sy_dias),
sprintf("95%% CI (Y): [%.1f, %.1f]%%\n\n", py_original - 1.96 * sy_dias, py_original + 1.96 * sy_dias),
"----- DIAS Covariance Matrix (X, Y) -----\n",
sprintf(" Var(DIAS-X) = %.6f Cov(DIAS-X, DIAS-Y) = %.4f\n", Cov_XY[1, 1], Cov_XY[1, 2]),
sprintf(" Cov(DIAS-X, DIAS-Y) = %.4f Var(DIAS-Y) = %.4f\n", Cov_XY[2, 1], Cov_XY[2, 2])
)
return(list(text = out, result = list(dias_x = px, se_x = sx, dias_y_original = py_original, se_y = sy_dias, cov_xy = Cov_XY,
model1 = m1, model2 = m2, x = x, y1 = y1, y2 = y2)))
}
# ============================================================
# MODELO COSSENO
# ============================================================
cos_model <- function(x, mesor, amplitude, acrophase, period = 24) {
mesor + amplitude * cos(2 * pi * (x - acrophase) / period)
}
dias_cos_text <- function(x, y1, y2, invert_curve1 = FALSE,
x_min = 0, x_max = 24, period = 24,
expected_intersections = NULL,
chute_mesor1 = NULL, chute_amplitude1 = NULL, chute_acrophase1 = NULL,
chute_mesor2 = NULL, chute_amplitude2 = NULL, chute_acrophase2 = NULL,
chute_period = NULL) {
if (invert_curve1) y1 <- 100 - y1
if (is.null(chute_mesor1)) chute_mesor1 <- mean(y1)
if (is.null(chute_amplitude1)) chute_amplitude1 <- (max(y1) - min(y1)) / 2
if (is.null(chute_acrophase1)) chute_acrophase1 <- 0
if (is.null(chute_mesor2)) chute_mesor2 <- mean(y2)
if (is.null(chute_amplitude2)) chute_amplitude2 <- (max(y2) - min(y2)) / 2
if (is.null(chute_acrophase2)) chute_acrophase2 <- 0
if (is.null(chute_period)) chute_period <- period
start1 <- list(mesor = chute_mesor1, amplitude = chute_amplitude1, acrophase = chute_acrophase1)
start2 <- list(mesor = chute_mesor2, amplitude = chute_amplitude2, acrophase = chute_acrophase2)
m1 <- try(nls(y1 ~ mesor + amplitude * cos(2 * pi * (x - acrophase) / period),
start = start1, control = nls.control(maxiter = 500)), silent = TRUE)
m2 <- try(nls(y2 ~ mesor + amplitude * cos(2 * pi * (x - acrophase) / period),
start = start2, control = nls.control(maxiter = 500)), silent = TRUE)
if (inherits(m1, "try-error") || inherits(m2, "try-error")) {
return(list(text = "Error: Unable to fit the cosine curve.", result = NULL))
}
p1 <- coef(m1); c1 <- vcov(m1)
p2 <- coef(m2); c2 <- vcov(m2)
diff_func <- function(x) {
cos_model(x, p1[1], p1[2], p1[3], period) - cos_model(x, p2[1], p2[2], p2[3], period)
}
roots <- find_roots(diff_func, x_min, x_max)
warning_msg <- NULL
if (!is.null(expected_intersections) && length(roots) > expected_intersections) {
roots <- roots[1:expected_intersections]
warning_msg <- paste0("⚠️ More intersections found than expected. Keeping the first ",
expected_intersections, ".")
}
results <- list()
for (i in seq_along(roots)) {
px <- roots[i]
py_current <- cos_model(px, p1[1], p1[2], p1[3], period)
py_original <- if (invert_curve1) 100 - py_current else py_current
Cov_all <- matrix(0, 6, 6)
Cov_all[1:3, 1:3] <- c1
Cov_all[4:6, 4:6] <- c2
h <- 1e-6
J <- matrix(0, 2, 6)
mt <- c(p1, p2)
for (j in 1:6) {
mp <- mt; mp[j] <- mp[j] + h
diff_func_pert <- function(x) {
cos_model(x, mp[1], mp[2], mp[3], period) - cos_model(x, mp[4], mp[5], mp[6], period)
}
root_pert <- tryCatch(
uniroot(diff_func_pert, c(px - 0.5, px + 0.5), tol = 1e-6)$root,
error = function(e) px
)
J[1, j] <- (root_pert - px) / h
}
for (j in 1:3) {
mp <- mt; mp[j] <- mp[j] + h
y_pert <- cos_model(px, mp[1], mp[2], mp[3], period)
J[2, j] <- (y_pert - py_current) / h
}
for (j in 4:6) {
J[2, j] <- 0
}
Cov_XY <- J %*% Cov_all %*% t(J)
sx <- sqrt(Cov_XY[1, 1])
sy_dias <- sqrt(Cov_XY[2, 2])
results[[i]] <- list(
dias_x = px,
se_x = sx,
dias_y_original = py_original,
se_y = sy_dias,
cov_xy = Cov_XY
)
}
out <- paste0(
"===== FITTED PARAMETERS (Cosine) =====\n",
sprintf("Curve 1: Mesor = %.2f ± %.2f | Amplitude = %.2f ± %.2f | Acrophase = %.2f ± %.2f\n",
p1[1], sqrt(diag(c1)[1]), p1[2], sqrt(diag(c1)[2]), p1[3], sqrt(diag(c1)[3])),
sprintf("Curve 2: Mesor = %.2f ± %.2f | Amplitude = %.2f ± %.2f | Acrophase = %.2f ± %.2f\n\n",
p2[1], sqrt(diag(c2)[1]), p2[2], sqrt(diag(c2)[2]), p2[3], sqrt(diag(c2)[3])),
if (!is.null(warning_msg)) paste0(warning_msg, "\n\n") else "",
sprintf("Number of intersections found: %d\n\n", length(results))
)
for (i in seq_along(results)) {
r <- results[[i]]
out <- paste0(
out,
sprintf("===== DIAS %d (Cosine) =====\n", i),
sprintf("DIAS-X: %.3f ± %.4f\n", r$dias_x, r$se_x),
sprintf("95%% CI (X): [%.3f, %.3f]\n", r$dias_x - 1.96 * r$se_x, r$dias_x + 1.96 * r$se_x),
sprintf("DIAS-Y: %.1f ± %.2f%%\n", r$dias_y_original, r$se_y),
sprintf("95%% CI (Y): [%.1f, %.1f]%%\n\n", r$dias_y_original - 1.96 * r$se_y, r$dias_y_original + 1.96 * r$se_y),
"----- DIAS Covariance Matrix (X, Y) -----\n",
sprintf(" Var(DIAS-X) = %.6f Cov(DIAS-X, DIAS-Y) = %.4f\n", r$cov_xy[1, 1], r$cov_xy[1, 2]),
sprintf(" Cov(DIAS-X, DIAS-Y) = %.4f Var(DIAS-Y) = %.4f\n", r$cov_xy[2, 1], r$cov_xy[2, 2])
)
}
return(list(
text = out,
result = list(
all_dias = results,
model1 = m1, model2 = m2,
x = x, y1 = y1, y2 = y2,
warning = warning_msg,
roots = roots
)
))
}
# ------------------------------------------------------------
# FUNÇÃO: WALD
# ------------------------------------------------------------
dias_wald_text <- function(res1, res2) {
r1 <- res1$result; r2 <- res2$result
if (is.null(r1) || is.null(r2)) return("Error: One of the DIAS values was not calculated.")
if (!is.null(r1$all_dias) && !is.null(r2$all_dias)) {
d1 <- r1$all_dias[[1]]
d2 <- r2$all_dias[[1]]
mu1 <- c(d1$dias_x, d1$dias_y_original)
mu2 <- c(d2$dias_x, d2$dias_y_original)
cov1 <- d1$cov_xy; cov2 <- d2$cov_xy
} else {
mu1 <- c(r1$dias_x, r1$dias_y_original)
mu2 <- c(r2$dias_x, r2$dias_y_original)
cov1 <- r1$cov_xy; cov2 <- r2$cov_xy
}
delta <- mu1 - mu2; cov_combined <- cov1 + cov2
W <- as.numeric(t(delta) %*% solve(cov_combined) %*% delta)
p_value <- 1 - pchisq(W, 2)
D2 <- sqrt(W)
effect <- if (D2 < 0.5) "Small" else if (D2 < 1) "Medium" else if (D2 < 2) "Large" else "Very Large"
sig <- if (p_value < 0.05) "SIGNIFICANT (p < 0.05)" else "NOT SIGNIFICANT (p ≥ 0.05)"
paste0(
"========================================\n BIDIMENSIONAL WALD TEST (DIAS)\n (Analogous to Hotelling's T²)\n========================================\n\n",
sprintf("Group 1: DIAS-X = %.3f ± %.4f, DIAS-Y = %.1f ± %.2f%%\n", mu1[1], sqrt(cov1[1, 1]), mu1[2], sqrt(cov1[2, 2])),
sprintf("Group 2: DIAS-X = %.3f ± %.4f, DIAS-Y = %.1f ± %.2f%%\n\n", mu2[1], sqrt(cov2[1, 1]), mu2[2], sqrt(cov2[2, 2])),
sprintf("W (T²) statistic: %.3f\nDegrees of freedom: 2\np-value: %.6f\nSignificance: %s\n\n", W, p_value, sig),
sprintf("D² (Mahalanobis) distance: %.3f\nClassification: %s\n========================================\n", D2, effect)
)
}
# RESPIDUOS
get_fit_stats <- function(model, y) {
if (inherits(model, "lm")) {
r2 <- summary(model)$r.squared
adj_r2 <- summary(model)$adj.r.squared
aic <- AIC(model)
bic <- BIC(model)
res_std <- summary(model)$sigma
} else if (inherits(model, "nls")) {
ss_res <- sum(residuals(model)^2)
ss_tot <- sum((y - mean(y))^2)
r2 <- 1 - (ss_res / ss_tot)
adj_r2 <- NA
aic <- AIC(model)
bic <- BIC(model)
res_std <- sqrt(ss_res / (length(y) - length(coef(model))))
} else {
return(NULL)
}
list(r2 = r2, adj_r2 = adj_r2, aic = aic, bic = bic, res_std = res_std)
}
get_reliability <- function(r2) {
if (r2 > 0.95) {
return("Very reliable")
} else if (r2 > 0.85) {
return("Reliable")
} else {
return("Poorly reliable")
}
}
# ============================================================
# INTERFACE
# ============================================================
ui <- fluidPage(
titlePanel("DIAS - Dias' Intercept Analysis System"),
tabsetPanel(
id = "main_tabs",
tabPanel(
title = "Model",
h2("Welcome to DIAS"),
h4("Select the curve type for analysis:"),
br(),
fluidRow(
column(4,
div(style = "border: 2px solid #2c7fb8; border-radius: 10px; padding: 15px; text-align: center;",
h4("Sigmoid 3"), p("HillSlope = 1", style = "font-size: 12px;"),
actionButton("btn_sig3", "Select", style = "background-color: #2c7fb8; color: white; width: 100%;")
)
),
column(4,
div(style = "border: 2px solid #2c7fb8; border-radius: 10px; padding: 15px; text-align: center;",
h4("Sigmoid 4"), p("Free HillSlope", style = "font-size: 12px;"),
actionButton("btn_sig4", "Select", style = "background-color: #2c7fb8; color: white; width: 100%;")
)
),
column(4,
div(style = "border: 2px solid #2c7fb8; border-radius: 10px; padding: 15px; text-align: center;",
h4("Linear"), p("Simple Linear Regression", style = "font-size: 12px;"),
actionButton("btn_lin", "Select", style = "background-color: #2c7fb8; color: white; width: 100%;")
)
)
),
br(),
fluidRow(
column(4,
div(style = "border: 2px solid #2c7fb8; border-radius: 10px; padding: 15px; text-align: center;",
h4("Cosine"), p("Circadian Rhythm", style = "font-size: 12px;"),
p("Multiple intersections", style = "font-size: 10px; color: #888; font-style: italic;"),
actionButton("btn_cos", "Select", style = "background-color: #2c7fb8; color: white; width: 100%;")
)
)
)
),
tabPanel(
title = "Analysis", value = "tab_analysis",
sidebarLayout(
sidebarPanel(
h4(uiOutput("selected_model")), hr(),
h4("Group 1"),
textInput("name1", "Group 1 Name:", value = "Group 1"),
radioButtons("invert1", "Invert Curve 1 (100% → 0%)?", choices = c("Yes" = TRUE, "No" = FALSE), selected = TRUE, inline = TRUE),
textAreaInput("x1", "X:", placeholder = "Ex: -9, -8.5, -8, -7.5, -7", rows = 2),
textAreaInput("y1_1", "Y1:", placeholder = "Ex: 12.04, 0, 10.47, 16.75", rows = 2),
textAreaInput("y2_1", "Y2:", placeholder = "Ex: 0, 17.45, 45.89, 52.14", rows = 2),
hr(),
h4("Group 2"),
textInput("name2", "Group 2 Name:", value = "Group 2"),
radioButtons("invert2", "Invert Curve 1 (100% → 0%)?", choices = c("Yes" = TRUE, "No" = FALSE), selected = TRUE, inline = TRUE),
textAreaInput("x2", "X:", placeholder = "Ex: -9, -8.5, -8, -7.5, -7", rows = 2),
textAreaInput("y1_2", "Y1:", placeholder = "Ex: 10.0, 2.0, 12.0, 18.0", rows = 2),
textAreaInput("y2_2", "Y2:", placeholder = "Ex: 2.0, 19.0, 47.0, 54.0", rows = 2),
conditionalPanel(
condition = "output.model_is_cosine == true",
hr(),
h4("Cosine Parameters"),
numericInput("cos_x_min", "X min:", value = 0, step = 0.5),
numericInput("cos_x_max", "X max:", value = 24, step = 0.5),
numericInput("cos_period", "Period (T):", value = 24, step = 1, min = 1),
numericInput("cos_n_int", "Expected intersections:", value = NA, step = 1, min = 1, max = 10)
),
hr(),
tags$details(
tags$summary(style = "font-weight: bold; cursor: pointer; color: #2c7fb8;",
"⚙️ Advanced Settings (Initial Guesses)"),
div(style = "padding: 10px; border: 1px solid #ddd; border-radius: 5px; margin-top: 10px;",
uiOutput("chutes_ui")
)
),
hr(),
radioButtons("wald_type", "Wald Test Type:",
choices = c(
"Automatic (recommended)" = "auto",
"Standard (1 DIAS)" = "standard"
),
selected = "auto"),
checkboxInput("show_posthoc", "Show post-hoc comparisons (all pairs)", FALSE),
conditionalPanel(
condition = "input.show_posthoc == true",
selectInput("posthoc_adjust", "Adjustment method:",
choices = c("bonferroni", "holm", "hochberg", "hommel", "BH", "BY", "none"),
selected = "bonferroni")
),
hr(),
uiOutput("ellipse_selector"),
actionButton("calculate", "Calculate DIAS", width = "100%", style = "color: #fff; background-color: #2c7fb8; font-size: 16px;"),
hr(),
h4("Bootstrap Paramétrico"),
numericInput("n_boot", "Número de reamostragens:",
value = 500, min = 100, max = 2000, step = 100),
actionButton("run_boot", "Run Bootstrap", width = "100%",
style = "color: #fff; background-color: #2c7fb8; font-size: 14px;"),
verbatimTextOutput("boot_result"),
downloadButton("downloadBoot", "Bootstrap Result (.txt)"),
hr(),
h4("Export Results"),
uiOutput("downloadDIAS1_btn"), uiOutput("downloadDIAS2_btn"),
downloadButton("downloadWald", "Wald Test (.txt)"),
downloadButton("downloadEllipse", "Ellipse Plot (.pdf)")
),
mainPanel(
tabsetPanel(
tabPanel("Numerical Results",
h3(textOutput("title_dias1")), verbatimTextOutput("dias1"),
h3(textOutput("title_dias2")), verbatimTextOutput("dias2"),
h3("Wald Test"), verbatimTextOutput("wald")
),
tabPanel("Ellipse Plot",
h3(textOutput("ellipse_title")),
plotOutput("ellipse_plot", height = "500px")
),
tabPanel("Fit & Residuals",
h3("Diagnostic Plots for Group 1"),
fluidRow(
column(6, plotOutput("plot_fit1", height = "250px")),
column(6, plotOutput("plot_resid1", height = "250px"))
),
fluidRow(
column(6, plotOutput("plot_qq1", height = "250px")),
column(6, verbatimTextOutput("stats1"))
),
hr(),
h3("Diagnostic Plots for Group 2"),
fluidRow(
column(6, plotOutput("plot_fit2", height = "250px")),
column(6, plotOutput("plot_resid2", height = "250px"))
),
fluidRow(
column(6, plotOutput("plot_qq2", height = "250px")),
column(6, verbatimTextOutput("stats2"))
)
)
)
)
)
)
)
)
# ------------------------------------------------------
# BOOTSTRAP PARAMÉTRICO
# ------------------------------------------------------
simulate_from_model <- function(x, coefs, model_id, sigma) {
if (model_id == 3) {
mu <- coefs[2] + (coefs[1] - coefs[2]) / (1 + 10^(coefs[3] - x))
} else if (model_id == 4) {
mu <- coefs[2] + (coefs[1] - coefs[2]) / (1 + 10^((coefs[3] - x) * coefs[4]))
} else if (model_id == 7) {
mu <- coefs[1] + coefs[2] * x
} else if (model_id == 8) {
mu <- coefs[1] + coefs[2] * cos(2 * pi * (x - coefs[3]) / 12)
} else {
mu <- rep(NA, length(x))
}
mu + rnorm(length(x), 0, sigma)
}
fit_model <- function(x, y1, y2, model_id) {
if (model_id == 3) {
res <- dias_3p_text(x, y1, y2, invert_curve1 = FALSE,
chute_top1 = 100, chute_bottom1 = 0, chute_lec501 = -7,
chute_top2 = 100, chute_bottom2 = 0, chute_lec502 = -7)
} else if (model_id == 4) {
res <- dias_4p_text(x, y1, y2, invert_curve1 = FALSE,
chute_top1 = 100, chute_bottom1 = 0, chute_lec501 = -7, chute_hs1 = 1,
chute_top2 = 100, chute_bottom2 = 0, chute_lec502 = -7, chute_hs2 = 1)
} else if (model_id == 7) {
res <- dias_linear_text(x, y1, y2, invert_curve1 = FALSE)
} else if (model_id == 8) {
res <- dias_cos_text(x, y1, y2, invert_curve1 = FALSE,
x_min = 0, x_max = 12, period = 12,
chute_mesor1 = mean(y1),
chute_amplitude1 = (max(y1) - min(y1)) / 2,
chute_acrophase1 = 0,
chute_mesor2 = mean(y2),
chute_amplitude2 = (max(y2) - min(y2)) / 2,
chute_acrophase2 = 0)
} else {
return(NULL)
}
return(res$result)
}
# SERVER
server <- function(input, output, session) {
selected_model <- reactiveVal(NULL)
observeEvent(input$btn_sig3, { selected_model(3); updateTabsetPanel(session, "main_tabs", selected = "tab_analysis") })
observeEvent(input$btn_sig4, { selected_model(4); updateTabsetPanel(session, "main_tabs", selected = "tab_analysis") })
observeEvent(input$btn_lin, { selected_model(7); updateTabsetPanel(session, "main_tabs", selected = "tab_analysis") })
observeEvent(input$btn_cos, { selected_model(8); updateTabsetPanel(session, "main_tabs", selected = "tab_analysis") })
output$selected_model <- renderUI({
req(selected_model())
if (selected_model() == 3) HTML("📊 Model: 3-Parameter Sigmoid (HillSlope = 1)")
else if (selected_model() == 4) HTML("📊 Model: 4-Parameter Sigmoid (Free HillSlope)")
else if (selected_model() == 7) HTML("📊 Model: Linear (Simple Linear Regression)")
else if (selected_model() == 8) HTML("📊 Model: Cosine (Multiple intersections)")
})
output$model_is_cosine <- reactive({ selected_model() == 8 })
outputOptions(output, "model_is_cosine", suspendWhenHidden = FALSE)
has_multiple_dias <- reactive({
req(res1(), res2())
r1 <- res1()$result; r2 <- res2()$result
n1 <- if (!is.null(r1$all_dias)) length(r1$all_dias) else 1
n2 <- if (!is.null(r2$all_dias)) length(r2$all_dias) else 1
(n1 > 1 || n2 > 1)
})
output$ellipse_selector <- renderUI({
if (has_multiple_dias()) {
req(res1(), res2())
r1 <- res1()$result; r2 <- res2()$result
n1 <- if (!is.null(r1$all_dias)) length(r1$all_dias) else 1
n2 <- if (!is.null(r2$all_dias)) length(r2$all_dias) else 1
max_dias <- max(n1, n2)
choices <- setNames(1:max_dias, paste("DIAS", 1:max_dias))
tagList(
hr(),
h4("Ellipse Plot Options"),
selectInput("ellipse_dias_select", "Select DIAS pair to plot:", choices = choices, selected = 1)
)
}
})
output$chutes_ui <- renderUI({
modelo <- selected_model()
if (is.null(modelo)) return(NULL)
if (modelo == 3) {
tagList(
h5("Group 1"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_top1_3", "Top:", value = 100), numericInput("chute_bottom1_3", "Bottom:", value = 0), numericInput("chute_lec501_3", "LogEC50:", value = -7)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_top2_3", "Top:", value = 100), numericInput("chute_bottom2_3", "Bottom:", value = 0), numericInput("chute_lec502_3", "LogEC50:", value = -7))
),
h5("Group 2"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_top1_3_2", "Top:", value = 100), numericInput("chute_bottom1_3_2", "Bottom:", value = 0), numericInput("chute_lec501_3_2", "LogEC50:", value = -7)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_top2_3_2", "Top:", value = 100), numericInput("chute_bottom2_3_2", "Bottom:", value = 0), numericInput("chute_lec502_3_2", "LogEC50:", value = -7))
)
)
} else if (modelo == 4) {
tagList(
h5("Group 1"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_top1_4", "Top:", value = 100), numericInput("chute_bottom1_4", "Bottom:", value = 0), numericInput("chute_lec501_4", "LogEC50:", value = -7), numericInput("chute_hs1_4", "HillSlope:", value = 1)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_top2_4", "Top:", value = 100), numericInput("chute_bottom2_4", "Bottom:", value = 0), numericInput("chute_lec502_4", "LogEC50:", value = -7), numericInput("chute_hs2_4", "HillSlope:", value = 1))
),
h5("Group 2"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_top1_4_2", "Top:", value = 100), numericInput("chute_bottom1_4_2", "Bottom:", value = 0), numericInput("chute_lec501_4_2", "LogEC50:", value = -7), numericInput("chute_hs1_4_2", "HillSlope:", value = 1)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_top2_4_2", "Top:", value = 100), numericInput("chute_bottom2_4_2", "Bottom:", value = 0), numericInput("chute_lec502_4_2", "LogEC50:", value = -7), numericInput("chute_hs2_4_2", "HillSlope:", value = 1))
)
)
} else if (modelo == 7) {
p("Linear regression does not require initial guesses.")
} else if (modelo == 8) {
tagList(
h5("Group 1"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_mesor1_cos", "Mesor (mean):", value = 50), numericInput("chute_amplitude1_cos", "Amplitude:", value = 25), numericInput("chute_acrophase1_cos", "Acrophase (peak time):", value = 6)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_mesor2_cos", "Mesor (mean):", value = 50), numericInput("chute_amplitude2_cos", "Amplitude:", value = 25), numericInput("chute_acrophase2_cos", "Acrophase (peak time):", value = 6))
),
h5("Group 2"),
fluidRow(
column(6, h6("Curve 1 (Y1)"), numericInput("chute_mesor1_cos_2", "Mesor (mean):", value = 50), numericInput("chute_amplitude1_cos_2", "Amplitude:", value = 25), numericInput("chute_acrophase1_cos_2", "Acrophase (peak time):", value = 6)),
column(6, h6("Curve 2 (Y2)"), numericInput("chute_mesor2_cos_2", "Mesor (mean):", value = 50), numericInput("chute_amplitude2_cos_2", "Amplitude:", value = 25), numericInput("chute_acrophase2_cos_2", "Acrophase (peak time):", value = 6))
),
h6("Period (T):"),
numericInput("chute_period_cos", "Period (T):", value = 24)
)
}
})
parse_data <- function(x_str, y1_str, y2_str) {
list(
x = as.numeric(unlist(strsplit(gsub(" ", "", x_str), ","))),
y1 = as.numeric(unlist(strsplit(gsub(" ", "", y1_str), ","))),
y2 = as.numeric(unlist(strsplit(gsub(" ", "", y2_str), ",")))
)
}
res1 <- eventReactive(input$calculate, {
d <- parse_data(input$x1, input$y1_1, input$y2_1)
inv <- as.logical(input$invert1)
if (selected_model() == 3) {
dias_3p_text(d$x, d$y1, d$y2, inv, chute_top1 = input$chute_top1_3, chute_bottom1 = input$chute_bottom1_3, chute_lec501 = input$chute_lec501_3, chute_top2 = input$chute_top2_3, chute_bottom2 = input$chute_bottom2_3, chute_lec502 = input$chute_lec502_3)
} else if (selected_model() == 4) {
dias_4p_text(d$x, d$y1, d$y2, inv, chute_top1 = input$chute_top1_4, chute_bottom1 = input$chute_bottom1_4, chute_lec501 = input$chute_lec501_4, chute_hs1 = input$chute_hs1_4, chute_top2 = input$chute_top2_4, chute_bottom2 = input$chute_bottom2_4, chute_lec502 = input$chute_lec502_4, chute_hs2 = input$chute_hs2_4)
} else if (selected_model() == 7) {
dias_linear_text(d$x, d$y1, d$y2, inv)
} else if (selected_model() == 8) {
dias_cos_text(d$x, d$y1, d$y2, inv, x_min = input$cos_x_min, x_max = input$cos_x_max, period = input$cos_period, expected_intersections = input$cos_n_int, chute_mesor1 = input$chute_mesor1_cos, chute_amplitude1 = input$chute_amplitude1_cos, chute_acrophase1 = input$chute_acrophase1_cos, chute_mesor2 = input$chute_mesor2_cos, chute_amplitude2 = input$chute_amplitude2_cos, chute_acrophase2 = input$chute_acrophase2_cos, chute_period = input$chute_period_cos)
}
})
res2 <- eventReactive(input$calculate, {
d <- parse_data(input$x2, input$y1_2, input$y2_2)
inv <- as.logical(input$invert2)
if (selected_model() == 3) {
dias_3p_text(d$x, d$y1, d$y2, inv, chute_top1 = input$chute_top1_3_2, chute_bottom1 = input$chute_bottom1_3_2, chute_lec501 = input$chute_lec501_3_2, chute_top2 = input$chute_top2_3_2, chute_bottom2 = input$chute_bottom2_3_2, chute_lec502 = input$chute_lec502_3_2)
} else if (selected_model() == 4) {
dias_4p_text(d$x, d$y1, d$y2, inv, chute_top1 = input$chute_top1_4_2, chute_bottom1 = input$chute_bottom1_4_2, chute_lec501 = input$chute_lec501_4_2, chute_hs1 = input$chute_hs1_4_2, chute_top2 = input$chute_top2_4_2, chute_bottom2 = input$chute_bottom2_4_2, chute_lec502 = input$chute_lec502_4_2, chute_hs2 = input$chute_hs2_4_2)
} else if (selected_model() == 7) {
dias_linear_text(d$x, d$y1, d$y2, inv)
} else if (selected_model() == 8) {
dias_cos_text(d$x, d$y1, d$y2, inv, x_min = input$cos_x_min, x_max = input$cos_x_max, period = input$cos_period, expected_intersections = input$cos_n_int, chute_mesor1 = input$chute_mesor1_cos_2, chute_amplitude1 = input$chute_amplitude1_cos_2, chute_acrophase1 = input$chute_acrophase1_cos_2, chute_mesor2 = input$chute_mesor2_cos_2, chute_amplitude2 = input$chute_amplitude2_cos_2, chute_acrophase2 = input$chute_acrophase2_cos_2, chute_period = input$chute_period_cos)
}
})
output$title_dias1 <- renderText({ paste0("DIAS - ", input$name1) })
output$title_dias2 <- renderText({ paste0("DIAS - ", input$name2) })
output$dias1 <- renderPrint({ req(res1()); cat(res1()$text) })
output$dias2 <- renderPrint({ req(res2()); cat(res2()$text) })
output$wald <- renderPrint({
req(res1(), res2())
wald_type <- input$wald_type
r1 <- res1()$result; r2 <- res2()$result
n1 <- if (!is.null(r1$all_dias)) length(r1$all_dias) else 1
n2 <- if (!is.null(r2$all_dias)) length(r2$all_dias) else 1
has_multiple <- (n1 > 1 || n2 > 1)
if (wald_type == "auto") {
if (has_multiple) {
cat("⚠️ Multiple DIAS detected. Using bivariate Wald tests with post-hoc adjustment.\n\n")
}
cat(dias_wald_text(res1(), res2()))
} else if (wald_type == "standard") {
if (has_multiple) cat("⚠️ Multiple DIAS detected. Standard Wald compares only the first DIAS.\n\n")
cat(dias_wald_text(res1(), res2()))
}
if (input$show_posthoc) {
cat("\n", wald_posthoc_text(res1(), res2(), adjust = input$posthoc_adjust))
}
})
output$downloadDIAS1_btn <- renderUI({ req(res1()); downloadButton("downloadDIAS1", paste0("DIAS - ", input$name1, " (.txt)")) })
output$downloadDIAS1 <- downloadHandler(
filename = function() { paste0("DIAS_", gsub(" ", "_", input$name1), ".txt") },
content = function(file) { writeLines(res1()$text, file) }
)
output$downloadDIAS2_btn <- renderUI({ req(res2()); downloadButton("downloadDIAS2", paste0("DIAS - ", input$name2, " (.txt)")) })
output$downloadDIAS2 <- downloadHandler(
filename = function() { paste0("DIAS_", gsub(" ", "_", input$name2), ".txt") },
content = function(file) { writeLines(res2()$text, file) }
)
output$downloadWald <- downloadHandler(
filename = function() { "Wald_Test.txt" },
content = function(file) {
wald_type <- input$wald_type
r1 <- res1()$result; r2 <- res2()$result
n1 <- if (!is.null(r1$all_dias)) length(r1$all_dias) else 1
n2 <- if (!is.null(r2$all_dias)) length(r2$all_dias) else 1
has_multiple <- (n1 > 1 || n2 > 1)
out <- ""
if (wald_type == "auto") {
if (has_multiple) {
out <- paste0("⚠️ Multiple DIAS detected. Using bivariate Wald tests with post-hoc adjustment.\n\n")
}
out <- paste0(out, dias_wald_text(res1(), res2()))
} else if (wald_type == "standard") {
if (has_multiple) {
out <- paste0("⚠️ Multiple DIAS detected. Standard Wald compares only the first DIAS.\n\n")
}
out <- paste0(out, dias_wald_text(res1(), res2()))
}
if (input$show_posthoc) out <- paste0(out, "\n", wald_posthoc_text(res1(), res2(), adjust = input$posthoc_adjust))
writeLines(out, file)
}
)
boot_text <- reactiveVal(NULL)
observeEvent(input$run_boot, {
req(res1(), res2(), selected_model())
r1 <- res1()$result
r2 <- res2()$result
if (is.null(r1) || is.null(r2)) return(NULL)
model_id <- selected_model()
B <- input$n_boot
extract_dias <- function(r) {
if (!is.null(r$all_dias)) {
d <- r$all_dias[[1]]
return(list(mu = c(d$dias_x, d$dias_y_original),
cov = d$cov_xy))
} else {
return(list(mu = c(r$dias_x, r$dias_y_original),
cov = r$cov_xy))
}
}
d1_obs <- extract_dias(r1)
d2_obs <- extract_dias(r2)
delta_obs <- d1_obs$mu - d2_obs$mu
cov_obs <- d1_obs$cov + d2_obs$cov
W_obs <- as.numeric(t(delta_obs) %*% solve(cov_obs) %*% delta_obs)
fit1 <- r1$model1
fit2 <- r1$model2
fit3 <- r2$model1
fit4 <- r2$model2
beta1 <- (coef(fit1) + coef(fit3)) / 2
beta2 <- (coef(fit2) + coef(fit4)) / 2
n1 <- length(r1$x)
n2 <- length(r2$x)
p1 <- length(coef(fit1))
p2 <- length(coef(fit2))
sig1 <- sqrt((sum(residuals(fit1)^2) + sum(residuals(fit3)^2)) / (n1 + n2 - 2 * p1))
sig2 <- sqrt((sum(residuals(fit2)^2) + sum(residuals(fit4)^2)) / (n1 + n2 - 2 * p2))
W_boot <- numeric(B)
withProgress(message = 'Running bootstrap...', value = 0, {
for (b in 1:B) {
incProgress(1/B, detail = paste("Iteration", b, "of", B))
sim_y1_g1 <- simulate_from_model(r1$x, beta1, model_id, sig1)
sim_y2_g1 <- simulate_from_model(r1$x, beta2, model_id, sig2)
sim_y1_g2 <- simulate_from_model(r2$x, beta1, model_id, sig1)
sim_y2_g2 <- simulate_from_model(r2$x, beta2, model_id, sig2)
sim_res1 <- fit_model(r1$x, sim_y1_g1, sim_y2_g1, model_id)
sim_res2 <- fit_model(r2$x, sim_y1_g2, sim_y2_g2, model_id)
if (is.null(sim_res1) || is.null(sim_res2)) {
W_boot[b] <- NA
next
}
sim_d1 <- extract_dias(sim_res1)
sim_d2 <- extract_dias(sim_res2)
delta_sim <- sim_d1$mu - sim_d2$mu
cov_sim <- sim_d1$cov + sim_d2$cov
W_boot[b] <- as.numeric(t(delta_sim) %*% solve(cov_sim) %*% delta_sim)
}
})
W_boot <- W_boot[!is.na(W_boot)]
p_boot <- mean(W_boot >= W_obs)
txt <- paste0(
"===== BOOTSTRAP PARAMÉTRICO =====\n",
sprintf("Reamostragens válidas: %d\n", length(W_boot)),
sprintf("W observado: %.3f\n", W_obs),
sprintf("p-valor bootstrap: %.4f\n", p_boot)
)
boot_text(txt)
})
output$boot_result <- renderPrint({
cat(boot_text())
})
output$downloadBoot <- downloadHandler(
filename = function() { "bootstrap_result.txt" },
content = function(file) {
writeLines(boot_text(), file)
}
)
# ELIPSES
draw_ellipse <- function(mean, cov, level = 0.95, col = "black", lwd = 2) {
radius <- sqrt(qchisq(level, 2))
eig <- eigen(cov)
angles <- seq(0, 2 * pi, length.out = 200)
lines(t(mean + radius * eig$vectors %*% diag(sqrt(eig$values)) %*% rbind(cos(angles), sin(angles))), col = col, lwd = lwd)
}
output$ellipse_title <- renderText({
if (has_multiple_dias()) {
paste0("Bidimensional Comparison – DIAS ", input$ellipse_dias_select, " (95% Ellipses)")
} else {
"Bidimensional Comparison (95% Ellipses)"
}
})
plot_ellipse <- function() {
r1 <- res1()$result
r2 <- res2()$result
if (is.null(r1) && is.null(r2)) return(NULL)
has_multi <- FALSE
if (!is.null(r1$all_dias) || !is.null(r2$all_dias)) has_multi <- TRUE
sel <- 1
if (has_multi) {
if (!is.null(input$ellipse_dias_select)) {
sel <- as.numeric(input$ellipse_dias_select)
if (is.na(sel) || sel < 1) sel <- 1
}
}
extract_dias <- function(r, s) {
if (is.null(r)) return(NULL)
if (!is.null(r$all_dias)) {
if (length(r$all_dias) >= s) {
d <- r$all_dias[[s]]
return(list(mu = c(d$dias_x, d$dias_y_original), cov = d$cov_xy))
} else return(NULL)
} else {
return(list(mu = c(r$dias_x, r$dias_y_original), cov = r$cov_xy))
}
}
d1 <- extract_dias(r1, sel)
d2 <- extract_dias(r2, sel)
mu_list <- list(); cov_list <- list(); cols <- c(); labs <- c()
if (!is.null(d1)) {
mu_list[[length(mu_list)+1]] <- d1$mu
cov_list[[length(cov_list)+1]] <- d1$cov
cols <- c(cols, "blue"); labs <- c(labs, input$name1)
}
if (!is.null(d2)) {
mu_list[[length(mu_list)+1]] <- d2$mu
cov_list[[length(cov_list)+1]] <- d2$cov
cols <- c(cols, "red"); labs <- c(labs, input$name2)
}
if (length(mu_list) == 0) return(NULL)
x_min <- Inf; x_max <- -Inf; y_min <- Inf; y_max <- -Inf
for (i in seq_along(mu_list)) {
x_min <- min(x_min, mu_list[[i]][1] - 3*sqrt(cov_list[[i]][1,1]))
x_max <- max(x_max, mu_list[[i]][1] + 3*sqrt(cov_list[[i]][1,1]))
y_min <- min(y_min, mu_list[[i]][2] - 3*sqrt(cov_list[[i]][2,2]))
y_max <- max(y_max, mu_list[[i]][2] + 3*sqrt(cov_list[[i]][2,2]))
}
main_title <- ifelse(has_multi,
paste0("Bidimensional Comparison – DIAS ", sel, " (95% Ellipses)"),
"Bidimensional Comparison (95% Ellipses)")
plot(0, type="n", xlim=c(x_min, x_max), ylim=c(y_min, y_max),
xlab="DIAS-X", ylab="DIAS-Y", main=main_title)
for (i in seq_along(mu_list)) {
draw_ellipse(mu_list[[i]], cov_list[[i]], col=cols[i], lwd=2.5)
points(mu_list[[i]][1], mu_list[[i]][2], col=cols[i],
pch=ifelse(i==1,16,1), cex=2)
arrows(mu_list[[i]][1] - sqrt(cov_list[[i]][1,1]), mu_list[[i]][2],
mu_list[[i]][1] + sqrt(cov_list[[i]][1,1]), mu_list[[i]][2],
col=cols[i], angle=90, code=3, length=0.1)
arrows(mu_list[[i]][1], mu_list[[i]][2] - sqrt(cov_list[[i]][2,2]),
mu_list[[i]][1], mu_list[[i]][2] + sqrt(cov_list[[i]][2,2]),
col=cols[i], angle=90, code=3, length=0.1)
}
if (length(labs) > 0) {
legend("topright", legend=labs, col=cols,
pch=ifelse(seq_along(labs)==1,16,1), lwd=2.5, bty="n")
}
}
output$ellipse_plot <- renderPlot({ plot_ellipse() })
output$downloadEllipse <- downloadHandler(
filename = function() { "DIAS_Ellipses.pdf" },
content = function(file) {
pdf(file, width = 8, height = 6)
plot_ellipse()
dev.off()
}
)
# ============================================================
# GRÁFICOS DE RESÍDUOS
# ============================================================
valid_model1 <- reactive({
req(res1())
r <- tryCatch(res1()$result, error=function(e) NULL)
if (is.null(r) || !is.list(r)) return(NULL)
if (!all(c("model1","model2","x","y1","y2") %in% names(r))) return(NULL)
if (!inherits(r$model1,"nls") && !inherits(r$model1,"lm")) return(NULL)
return(list(model1=r$model1, model2=r$model2, x=r$x, y1=r$y1, y2=r$y2))
})
valid_model2 <- reactive({
req(res2())
r <- tryCatch(res2()$result, error=function(e) NULL)
if (is.null(r) || !is.list(r)) return(NULL)
if (!all(c("model1","model2","x","y1","y2") %in% names(r))) return(NULL)
if (!inherits(r$model1,"nls") && !inherits(r$model1,"lm")) return(NULL)
return(list(model1=r$model1, model2=r$model2, x=r$x, y1=r$y1, y2=r$y2))
})
output$plot_fit1 <- renderPlot({
tryCatch({
req(valid_model1())
vm <- valid_model1()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
d <- data.frame(x=vm$x, y=vm$y1)
m <- vm$model1
if (inherits(m,"nls")) {
new_x <- seq(min(d$x), max(d$x), length.out=200)
pred <- tryCatch(predict(m, newdata=list(x=new_x), se.fit=TRUE),
error=function(e) list(fit=predict(m,newdata=list(x=new_x)), se.fit=rep(NA,length(new_x))))
if (!is.list(pred) || !("fit" %in% names(pred))) pred <- list(fit=predict(m,newdata=list(x=new_x)), se.fit=rep(NA,length(new_x)))
if (all(is.na(pred$se.fit))) {
pred_df <- data.frame(x=new_x, fit=pred$fit)
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_line(data=pred_df,aes(x=x,y=fit),color="red",size=1)+
labs(x="X",y="Y",title="Curve 1: Data + Fit (no CI available)")+theme_minimal()
} else {
ci_u <- pred$fit+1.96*pred$se.fit; ci_l <- pred$fit-1.96*pred$se.fit
pred_df <- data.frame(x=new_x, fit=pred$fit, lower=ci_l, upper=ci_u)
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_ribbon(data=pred_df,aes(x=x,ymin=lower,ymax=upper),fill="red",alpha=0.15,inherit.aes=FALSE)+
geom_line(data=pred_df,aes(x=x,y=fit),color="red",size=1)+
labs(x="X",y="Y",title="Curve 1: Data + Fit (95% CI)")+theme_minimal()
}
} else if (inherits(m,"lm")) {
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_smooth(method="lm",se=TRUE,color="red",fill="red",alpha=0.15)+
labs(x="X",y="Y",title="Curve 1: Data + Linear Fit (95% CI)")+theme_minimal()
} else {
ggplot()+annotate("text",x=0.5,y=0.5,label="Model type not recognized.")+theme_void()
}
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$plot_resid1 <- renderPlot({
tryCatch({
req(valid_model1()); vm <- valid_model1()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
m <- vm$model1
if (!inherits(m,"nls") && !inherits(m,"lm")) return(ggplot()+annotate("text",x=0.5,y=0.5,label="Model not available.")+theme_void())
d <- data.frame(fitted=fitted(m), resid=residuals(m))
ggplot(d,aes(fitted,resid))+geom_point(size=2.5,color="darkgreen")+
geom_hline(yintercept=0,color="red",linetype="dashed",size=1)+
geom_smooth(method="loess",se=FALSE,color="blue",size=1)+
labs(x="Fitted",y="Residuals",title="Curve 1: Residuals vs Fitted")+theme_minimal()
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$plot_qq1 <- renderPlot({
tryCatch({
req(valid_model1()); vm <- valid_model1()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
m <- vm$model1
if (!inherits(m,"nls") && !inherits(m,"lm")) return(ggplot()+annotate("text",x=0.5,y=0.5,label="Model not available.")+theme_void())
resid_vals <- residuals(m)
ggplot(data.frame(resid=resid_vals),aes(sample=resid))+
stat_qq(size=2.5,color="darkblue")+stat_qq_line(color="red",size=1)+
labs(x="Theoretical Quantiles",y="Sample Quantiles",title="Curve 1: Q-Q Plot")+theme_minimal()
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$stats1 <- renderPrint({
tryCatch({
req(valid_model1()); vm <- valid_model1()
if (is.null(vm)) { cat("No valid fit available."); return() }
m <- vm$model1; y <- vm$y1
if (!inherits(m,"nls") && !inherits(m,"lm")) { cat("Model not available."); return() }
stats <- get_fit_stats(m,y)
if (is.null(stats)) { cat("Could not compute statistics."); return() }
reliability <- get_reliability(stats$r2)
cat(sprintf("R² = %.4f\n", stats$r2))
cat(sprintf("AIC = %.2f\n", stats$aic))
cat(sprintf("BIC = %.2f\n", stats$bic))
cat(sprintf("Residual SD = %.4f\n", stats$res_std))
cat(sprintf("\nReliability: %s\n", reliability))
}, error=function(e) cat("Error computing statistics:", e$message))
})
output$plot_fit2 <- renderPlot({
tryCatch({
req(valid_model2()); vm <- valid_model2()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
d <- data.frame(x=vm$x, y=vm$y1)
m <- vm$model1
if (inherits(m,"nls")) {
new_x <- seq(min(d$x), max(d$x), length.out=200)
pred <- tryCatch(predict(m, newdata=list(x=new_x), se.fit=TRUE),
error=function(e) list(fit=predict(m,newdata=list(x=new_x)), se.fit=rep(NA,length(new_x))))
if (!is.list(pred) || !("fit" %in% names(pred))) pred <- list(fit=predict(m,newdata=list(x=new_x)), se.fit=rep(NA,length(new_x)))
if (all(is.na(pred$se.fit))) {
pred_df <- data.frame(x=new_x, fit=pred$fit)
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_line(data=pred_df,aes(x=x,y=fit),color="red",size=1)+
labs(x="X",y="Y",title="Curve 2: Data + Fit (no CI available)")+theme_minimal()
} else {
ci_u <- pred$fit+1.96*pred$se.fit; ci_l <- pred$fit-1.96*pred$se.fit
pred_df <- data.frame(x=new_x, fit=pred$fit, lower=ci_l, upper=ci_u)
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_ribbon(data=pred_df,aes(x=x,ymin=lower,ymax=upper),fill="red",alpha=0.15,inherit.aes=FALSE)+
geom_line(data=pred_df,aes(x=x,y=fit),color="red",size=1)+
labs(x="X",y="Y",title="Curve 2: Data + Fit (95% CI)")+theme_minimal()
}
} else if (inherits(m,"lm")) {
ggplot(d,aes(x,y))+geom_point(size=2.5,color="blue")+
geom_smooth(method="lm",se=TRUE,color="red",fill="red",alpha=0.15)+
labs(x="X",y="Y",title="Curve 2: Data + Linear Fit (95% CI)")+theme_minimal()
} else {
ggplot()+annotate("text",x=0.5,y=0.5,label="Model type not recognized.")+theme_void()
}
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$plot_resid2 <- renderPlot({
tryCatch({
req(valid_model2()); vm <- valid_model2()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
m <- vm$model1
if (!inherits(m,"nls") && !inherits(m,"lm")) return(ggplot()+annotate("text",x=0.5,y=0.5,label="Model not available.")+theme_void())
d <- data.frame(fitted=fitted(m), resid=residuals(m))
ggplot(d,aes(fitted,resid))+geom_point(size=2.5,color="darkgreen")+
geom_hline(yintercept=0,color="red",linetype="dashed",size=1)+
geom_smooth(method="loess",se=FALSE,color="blue",size=1)+
labs(x="Fitted",y="Residuals",title="Curve 2: Residuals vs Fitted")+theme_minimal()
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$plot_qq2 <- renderPlot({
tryCatch({
req(valid_model2()); vm <- valid_model2()
if (is.null(vm)) return(ggplot()+annotate("text",x=0.5,y=0.5,label="No valid fit available.")+theme_void())
m <- vm$model1
if (!inherits(m,"nls") && !inherits(m,"lm")) return(ggplot()+annotate("text",x=0.5,y=0.5,label="Model not available.")+theme_void())
resid_vals <- residuals(m)
ggplot(data.frame(resid=resid_vals),aes(sample=resid))+
stat_qq(size=2.5,color="darkblue")+stat_qq_line(color="red",size=1)+
labs(x="Theoretical Quantiles",y="Sample Quantiles",title="Curve 2: Q-Q Plot")+theme_minimal()
}, error=function(e) ggplot()+annotate("text",x=0.5,y=0.5,label=paste("Error:",e$message))+theme_void())
})
output$stats2 <- renderPrint({
tryCatch({
req(valid_model2()); vm <- valid_model2()
if (is.null(vm)) { cat("No valid fit available."); return() }
m <- vm$model1; y <- vm$y1
if (!inherits(m,"nls") && !inherits(m,"lm")) { cat("Model not available."); return() }
stats <- get_fit_stats(m,y)
if (is.null(stats)) { cat("Could not compute statistics."); return() }
reliability <- get_reliability(stats$r2)
cat(sprintf("R² = %.4f\n", stats$r2))
cat(sprintf("AIC = %.2f\n", stats$aic))
cat(sprintf("BIC = %.2f\n", stats$bic))
cat(sprintf("Residual SD = %.4f\n", stats$res_std))
cat(sprintf("\nReliability: %s\n", reliability))
}, error=function(e) cat("Error computing statistics:", e$message))
})
}
shinyApp(ui = ui, server = server)