Time to fatigue fracture of glass fibre composites using survival analysis

# Contoh kode R: Kaplan-Meier Survival and Cox Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(survival)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: waktu hingga kejadian dengan sensor kanan ----
n <- 200
d <- data.frame(stress_amplitude = rnorm(n, 50, 10),
                test_temperature = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$stress_amplitude))
z2 <- as.numeric(scale(d$test_temperature))
b1 <- runif(1, 0.4, 0.75) * sample(c(-1, 1), 1, prob = c(0.35, 0.65))   # log-HR per 1 SD
b2 <- runif(1, 0.1, 0.35) * sample(c(-1, 1), 1)
k <- runif(1, 1.1, 1.8)                            # bentuk Weibull (hazard naik)
med0 <- max(100, 1)                             # median waktu pada profil rata-rata
tk <- med0 * (-log(runif(n)) / (log(2) * exp(b1 * z1 + b2 * z2)))^(1 / k)
sensor <- pmin(runif(n, 0.3, 3.5) * med0, 3 * med0)  # putus pantau + akhir studi
d$time_to_fatigue_fracture <- round(pmax(pmin(tk, sensor), 0.05 * med0), 2)
d$status <- as.integer(tk <= sensor)                # 1 = kejadian, 0 = tersensor
d$kelompok <- ifelse(d$stress_amplitude > median(d$stress_amplitude), "High", "Low")
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: Kaplan-Meier, uji log-rank, regresi Cox ----
d$kelompok <- factor(d$kelompok, levels = c("Low", "High"))
km <- survfit(Surv(time_to_fatigue_fracture, status) ~ kelompok, data = d)
print(km)
km_all <- survfit(Surv(time_to_fatigue_fracture, status) ~ 1, data = d)
lr <- survdiff(Surv(time_to_fatigue_fracture, status) ~ kelompok, data = d)   # uji log-rank
print(lr)
p_lr <- pchisq(lr$chisq, df = length(lr$n) - 1, lower.tail = FALSE)
cox <- coxph(Surv(time_to_fatigue_fracture, status) ~ stress_amplitude + test_temperature, data = d)
sc <- summary(cox)
print(sc)
zph <- cox.zph(cox)                                 # uji asumsi proportional hazards
print(zph)
cidx <- unname(sc$concordance[1])                   # C-index (Harrell)
# HR per 1 simpangan baku supaya mudah dibandingkan antarvariabel
sdx <- c(sd(d$stress_amplitude), sd(d$test_temperature))
hr_sd <- exp(coef(cox) * sdx)
ci_sd <- exp(confint(cox) * sdx)
print(round(cbind(HR_per_SD = hr_sd, ci_sd), 3))
median_all <- unname(summary(km_all)$table["median"])
if (!is.finite(median_all)) median_all <- unname(summary(km_all)$table["rmean"])  # cadangan: rerata terbatas
med_kel <- summary(km)$table[, "median"]
p_x1 <- sc$coefficients["stress_amplitude", "Pr(>|z|)"]


# ---- Grafik ----
# 1) kurva Kaplan-Meier bertangga dengan pita CI 95% dan tanda sensor
kel <- rep(levels(d$kelompok), km$strata)
dk <- data.frame(kel = kel, t = km$time, s = km$surv, lo = km$lower, hi = km$upper, cen = km$n.censor)
dk$lo[is.na(dk$lo)] <- 0
dk$hi[is.na(dk$hi)] <- dk$s[is.na(dk$hi)]
tangga <- function(g) {
  a <- rbind(data.frame(kel = g$kel[1], t = 0, s = 1, lo = 1, hi = 1, cen = 0), g)
  m <- nrow(a)
  data.frame(kel = a$kel[1], t = rep(a$t, each = 2)[-1], s = rep(a$s, each = 2)[-(2 * m)],
             lo = rep(a$lo, each = 2)[-(2 * m)], hi = rep(a$hi, each = 2)[-(2 * m)])
}
ds <- do.call(rbind, lapply(split(dk, dk$kel), tangga))
lab_kel <- c(Low = "Low stress_amplitude", High = "High stress_amplitude")
ds$kel <- factor(lab_kel[ds$kel], levels = lab_kel)
dc <- subset(dk, cen > 0)
dc$kel <- factor(lab_kel[dc$kel], levels = lab_kel)
p1 <- ggplot(ds, aes(t, s, colour = kel, fill = kel)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, colour = NA) +
  geom_line(linewidth = 1) +
  geom_point(data = dc, aes(t, s), shape = 3, size = 1.8, show.legend = FALSE) +
  geom_hline(yintercept = 0.5, linetype = "dotted", colour = warna$redup) +
  annotate("text", x = max(ds$t) * 0.98, y = 0.97, hjust = 1, size = 3.2, colour = warna$teks,
           label = sprintf("Log-rank p %s", if (p_lr < 0.001) "< 0.001" else sprintf("= %.3f", p_lr))) +
  scale_colour_manual(values = pal[c(1, 2)], name = NULL) +
  scale_fill_manual(values = pal[c(1, 2)], name = NULL) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(title = "Kaplan-Meier survival (+ = censored)", x = "time_to_fatigue_fracture", y = "Survival probability")
# 2) forest plot hazard ratio per 1 SD
dh <- data.frame(var = factor(make.unique(c("stress_amplitude", "test_temperature"), sep = " "),
                              levels = rev(make.unique(c("stress_amplitude", "test_temperature"), sep = " "))),
                 hr = hr_sd, lo = ci_sd[, 1], hi = ci_sd[, 2])
dh$teks <- sprintf("%.2f (%.2f-%.2f)", dh$hr, dh$lo, dh$hi)
p2 <- ggplot(dh, aes(hr, var)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = warna$redup) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.18, colour = pal[3], linewidth = 0.9) +
  geom_point(size = 3.2, colour = pal[3]) +
  geom_text(aes(label = teks), vjust = -1.3, size = 3, colour = warna$teks) +
  scale_x_log10() +
  scale_y_discrete(expand = expansion(add = 0.7)) +
  labs(title = "Cox hazard ratios per 1 SD (95% CI)", x = "Hazard ratio (log scale)", y = NULL)
# 3) residual Schoenfeld terskala untuk stress_amplitude (cek proportional hazards)
dz <- data.frame(t = zph$time, r = zph$y[, 1])
p3 <- ggplot(dz, aes(t, r)) +
  geom_hline(yintercept = coef(cox)[1], linetype = "dashed", colour = warna$redup) +
  geom_point(colour = pal[4], alpha = 0.5, size = 1.3) +
  geom_smooth(method = "loess", formula = y ~ x, colour = pal[2], fill = pal[2], alpha = 0.15, linewidth = 0.9) +
  labs(title = sprintf("Schoenfeld: stress_amplitude (p = %.3f)", zph$table[1, "p"]),
       x = "time_to_fatigue_fracture", y = "Scaled residual")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3) + plot_layout(heights = c(1.3, 1))

print(p_hasil)