# ST 453 Final Project - Output File
# Trace Brown
# 04/21/26

#
# This script loads all bootstrap results and produces plots/tables for the
# report.  Includes percentile, basic, normal, and BCa bootstrap confidence
# intervals, plus a side-by-side comparison with the parametric (Wald) CIs
# from the midterm.
#
# Usage: Rscript out_file.r <N>
#   N  -- number of bootstrap samples to load (seeds 1 through N)
#
# Example: Rscript out_file.r 1000

# No R packages are permitted for use in this assignment.



# -----------------------------------------------------------------------------
# Command line arguments
# -----------------------------------------------------------------------------

args = commandArgs(trailingOnly=TRUE)

if(length(args) < 1){
  stop("Usage: Rscript out_file.r <N>")
}

N = as.integer(args[1])

print(paste("Loading", N, "bootstrap results..."))



# -----------------------------------------------------------------------------
# Helper functions (needed for jackknife in BCa CI)
# -----------------------------------------------------------------------------

gaussian_elimination = function( M, b){
  
  p = dim(M)[1]
  aug = cbind( M, b)
  
  for(k in 1:(p-1)){
    
    max_idx = k
    for(i in (k+1):p){
      if(abs(aug[i,k]) > abs(aug[max_idx,k]))  max_idx = i
    }
    
    if(max_idx != k){
      temp = aug[k,]
      aug[k,] = aug[max_idx,]
      aug[max_idx,] = temp
    }
    
    if(abs(aug[k,k]) < 1e-12)  next
    
    for(i in (k+1):p){
      factor = aug[i,k] / aug[k,k]
      aug[i,] = aug[i,] - factor * aug[k,]
    }
  }
  
  return(list( M=aug[,1:p], b=aug[,p+1]))
}


back_substitution = function( M, b){
  
  p = dim(M)[1]
  x = rep( NA, p)
  
  x[p] = b[p] / M[p,p]
  for(i in (p-1):1){
    x[i] = b[i] / M[i,i]
    for(j in (i+1):p)  x[i] = x[i] - M[i,j] * x[j] / M[i,i]
  }
  
  return(x)
}


solve_system = function( A, b){
  row_echelon = gaussian_elimination( A, b)
  return( back_substitution( row_echelon$M, row_echelon$b))
}


logistic = function( z){
  return( 1 / (1 + exp(-z)))
}


newton_raphson_logistic = function( y, X, max_iter=200, tol=1e-8){
  
  n = dim(X)[1]
  p = dim(X)[2]
  beta = rep( 0, p)
  
  for(iter in 1:max_iter){
    
    eta = X %*% beta
    prob = logistic(eta)
    prob[prob < 1e-10] = 1e-10
    prob[prob > 1 - 1e-10] = 1 - 1e-10
    
    gradient = t(X) %*% (y - prob)
    W = as.vector(prob * (1 - prob))
    XtWX = t(X) %*% (W * X)
    
    delta = solve_system( XtWX, as.vector(gradient))
    beta_new = beta + delta
    
    if(sqrt(sum((beta_new - beta)^2)) < tol){
      
      prob_final = logistic(X %*% beta_new)
      prob_final[prob_final < 1e-10] = 1e-10
      prob_final[prob_final > 1 - 1e-10] = 1 - 1e-10
      W_final = as.vector(prob_final * (1 - prob_final))
      XtWX_final = t(X) %*% (W_final * X)
      
      cov_mat = matrix( 0, nrow=p, ncol=p)
      for(j in 1:p){
        e_j = rep( 0, p)
        e_j[j] = 1
        cov_mat[,j] = solve_system( XtWX_final, e_j)
      }
      
      se = sqrt(diag(cov_mat))
      return(list( beta=beta_new, cov_mat=cov_mat, se=se,
                   iter=iter, converged=TRUE))
    }
    
    beta = beta_new
  }
  
  prob_final = logistic(X %*% beta)
  prob_final[prob_final < 1e-10] = 1e-10
  prob_final[prob_final > 1 - 1e-10] = 1 - 1e-10
  W_final = as.vector(prob_final * (1 - prob_final))
  XtWX_final = t(X) %*% (W_final * X)
  
  cov_mat = matrix( 0, nrow=p, ncol=p)
  for(j in 1:p){
    e_j = rep( 0, p)
    e_j[j] = 1
    cov_mat[,j] = solve_system( XtWX_final, e_j)
  }
  
  se = sqrt(diag(cov_mat))
  return(list( beta=beta, cov_mat=cov_mat, se=se,
               iter=max_iter, converged=FALSE))
}



# -----------------------------------------------------------------------------
# Load bootstrap results
# -----------------------------------------------------------------------------

first_result = readRDS("bootstrap_results/boot_seed_1.rds")
p             = first_result$p
n             = first_result$n
beta_hat_real = first_result$beta_hat_real
se_real       = first_result$se_real
X_means       = first_result$X_means
X_sds         = first_result$X_sds

coef_names = c("Intercept", "Age", "Ejection Fraction", "Serum Creatinine", 
               "Time")

# Allocate storage matrices
beta_boot_mat = matrix( NA, nrow=N, ncol=p)
se_boot_mat   = matrix( NA, nrow=N, ncol=p)
converged_vec = rep( NA, N)
iter_vec      = rep( NA, N)
events_vec    = rep( NA, N)

# Load all bootstrap results
num_loaded = 0
for(k in 1:N){
  
  fname = paste0("bootstrap_results/boot_seed_", k, ".rds")
  
  if(!file.exists(fname)){
    print(paste("WARNING: missing file for seed", k))
    next
  }
  
  res = readRDS(fname)
  beta_boot_mat[k,] = res$beta_boot
  se_boot_mat[k,]   = res$se_boot
  converged_vec[k]  = res$converged
  iter_vec[k]       = res$iter
  events_vec[k]     = res$boot_events
  num_loaded        = num_loaded + 1
}

print(paste("Successfully loaded", num_loaded, "of", N, "bootstrap results"))



# -----------------------------------------------------------------------------
# Bootstrap summary statistics
# -----------------------------------------------------------------------------

beta_boot_mean = colMeans(beta_boot_mat, na.rm=T)
beta_boot_med  = apply(beta_boot_mat, 2, median, na.rm=T)
beta_boot_sd   = apply(beta_boot_mat, 2, sd, na.rm=T)

# Bootstrap-estimated bias of the MLE
boot_bias = beta_boot_mean - beta_hat_real

# Mean of the per-bootstrap Hessian-based SEs (diagnostic, not used for CIs)
se_boot_meanH = colMeans(se_boot_mat, na.rm=T)



# -----------------------------------------------------------------------------
# Wald (parametric) 95% CI from the original fit, for comparison
# -----------------------------------------------------------------------------

alpha    = 0.05
z_crit   = qnorm(1 - alpha/2)
wald_lower = beta_hat_real - z_crit * se_real
wald_upper = beta_hat_real + z_crit * se_real
wald_width = wald_upper - wald_lower



# -----------------------------------------------------------------------------
# Bootstrap confidence intervals
# -----------------------------------------------------------------------------

# (1) Percentile bootstrap: alpha/2 and 1-alpha/2 quantiles of bootstrap dist
perc_lower = apply(beta_boot_mat, 2, quantile, probs=alpha/2,   na.rm=T)
perc_upper = apply(beta_boot_mat, 2, quantile, probs=1-alpha/2, na.rm=T)
perc_width = perc_upper - perc_lower


# (2) Basic (pivotal/reflection) bootstrap CI
#     [2*theta_hat - q_(1-alpha/2), 2*theta_hat - q_(alpha/2)]
basic_lower = 2 * beta_hat_real - perc_upper
basic_upper = 2 * beta_hat_real - perc_lower
basic_width = basic_upper - basic_lower


# (3) Normal-approximation bootstrap CI: theta_hat +/- z * se_boot
norm_lower = beta_hat_real - z_crit * beta_boot_sd
norm_upper = beta_hat_real + z_crit * beta_boot_sd
norm_width = norm_upper - norm_lower


# (4) BCa (bias-corrected and accelerated) bootstrap CI
#     Requires the jackknife estimates to compute the acceleration constant.

print("Computing jackknife replicates for BCa CI...")

dat = read.csv("heart_failure_clinical_records_dataset.csv")

X_real = cbind( 1,
                dat$age,
                dat$ejection_fraction,
                dat$serum_creatinine,
                dat$time)
y_real = dat$DEATH_EVENT

X_scaled = X_real
for(j in 2:p){
  X_scaled[,j] = (X_real[,j] - X_means[j-1]) / X_sds[j-1]
}

# Leave-one-out jackknife
jack_beta = matrix( NA, nrow=n, ncol=p)
for(i in 1:n){
  fit_jack = newton_raphson_logistic( y_real[-i], X_scaled[-i, , drop=FALSE])
  jack_beta[i,] = fit_jack$beta
  if(i %% 50 == 0)  print(paste("  jackknife", i, "of", n))
}

jack_mean = colMeans(jack_beta)


# BCa parameters
# Bias-correction:    z0 = Phi^{-1}(P*(beta* < beta_hat))
# Acceleration:       a  = sum((mean - x_jack)^3) / (6 * sum((mean - x_jack)^2)^(3/2))
z0 = qnorm(colMeans(beta_boot_mat < matrix(beta_hat_real, nrow=N, ncol=p, byrow=T),
                    na.rm=T))

a_const = rep( NA, p)
for(j in 1:p){
  num = sum((jack_mean[j] - jack_beta[,j])^3)
  den = 6 * (sum((jack_mean[j] - jack_beta[,j])^2))^(3/2)
  a_const[j] = num / den
}

# BCa-adjusted percentile boundaries
bca_lower = rep( NA, p)
bca_upper = rep( NA, p)
for(j in 1:p){
  
  z_lo = qnorm(alpha/2)
  z_hi = qnorm(1 - alpha/2)
  
  alpha_lo = pnorm(z0[j] + (z0[j] + z_lo) / (1 - a_const[j] * (z0[j] + z_lo)))
  alpha_hi = pnorm(z0[j] + (z0[j] + z_hi) / (1 - a_const[j] * (z0[j] + z_hi)))
  
  bca_lower[j] = quantile(beta_boot_mat[,j], probs=alpha_lo, na.rm=T)
  bca_upper[j] = quantile(beta_boot_mat[,j], probs=alpha_hi, na.rm=T)
}
bca_width = bca_upper - bca_lower



# -----------------------------------------------------------------------------
# Bootstrap hypothesis tests (H0: beta_j = 0)
# Two-sided percentile-based p-value derived by inverting the percentile CI:
#   p = 2 * min( P*(beta_b* <= 0), P*(beta_b* >= 0) )
# -----------------------------------------------------------------------------

boot_pvals = rep( NA, p)
for(j in 1:p){
  prop_below = mean(beta_boot_mat[,j] <= 0, na.rm=T)
  prop_above = mean(beta_boot_mat[,j] >= 0, na.rm=T)
  boot_pvals[j] = 2 * min(prop_below, prop_above)
  
  # Floor at 1/N (a bootstrap p-value cannot be smaller than the resolution
  # of the bootstrap distribution)
  if(boot_pvals[j] == 0)  boot_pvals[j] = 1 / N
}

# Wald p-values from the original fit (for comparison)
wald_pvals = rep( NA, p)
for(j in 1:p){
  z_stat = beta_hat_real[j] / se_real[j]
  wald_pvals[j] = 2 * (1 - pnorm(abs(z_stat)))
}



# -----------------------------------------------------------------------------
# Bootstrap sampling distribution plots
# -----------------------------------------------------------------------------

print("Generating bootstrap distribution plots...")

pdf("bootstrap_distributions.pdf", width=10, height=8)

par(mfrow=c(2,3))

for(j in 1:p){
  
  # Use empirical range padded for the normal overlay
  xrange = range(c(beta_boot_mat[,j], 
                   beta_hat_real[j] - 4*beta_boot_sd[j],
                   beta_hat_real[j] + 4*beta_boot_sd[j]), na.rm=T)
  grid = seq( xrange[1], xrange[2], length.out=200)
  
  hist(beta_boot_mat[,j], freq=F, breaks=floor(sqrt(N)),
       main=coef_names[j],
       xlab=expression(hat(beta)^"*"),
       xlim=xrange)
  
  # Original estimate (green solid line)
  abline( v=beta_hat_real[j], col="green", lwd=3)
  
  # Bootstrap mean (blue dashed line)
  abline( v=beta_boot_mean[j], col="blue", lwd=2, lty=2)
  
  # Percentile CI bounds (purple dashed lines)
  abline( v=perc_lower[j], col="purple", lwd=2, lty=2)
  abline( v=perc_upper[j], col="purple", lwd=2, lty=2)
  
  # Wald CI bounds (red dotted lines)
  abline( v=wald_lower[j], col="red", lwd=2, lty=3)
  abline( v=wald_upper[j], col="red", lwd=2, lty=3)
  
  # Normal density centered at the original estimate, with bootstrap SD
  lines( grid, dnorm(grid, mean=beta_hat_real[j], sd=beta_boot_sd[j]), 
         lwd=3)
  
  if(j == 1){
    legend("topright",
           legend=c("Estimate", "Boot Mean", "Perc CI", "Wald CI", "Normal"),
           col=c("green", "blue", "purple", "red", "black"),
           lty=c(1, 2, 2, 3, 1),
           lwd=c(3, 2, 2, 2, 3),
           cex=0.6)
  }
}

dev.off()
print("Saved: bootstrap_distributions.pdf")



# -----------------------------------------------------------------------------
# Diagnostic plots
# -----------------------------------------------------------------------------

print("Generating diagnostic plots...")

pdf("bootstrap_diagnostics.pdf", width=10, height=8)

par(mfrow=c(2,2))

# Plot 1: Bootstrap-estimated bias
barplot(boot_bias, names.arg=c("Int","Age","EF","SC","Time"),
        main="Bootstrap-Estimated Bias",
        ylab="Mean(boot) - beta_hat",
        col="steelblue")
abline( h=0, col="red", lty=2, lwd=2)

# Plot 2: SE comparison: Wald vs Bootstrap
barplot( rbind(se_real, beta_boot_sd), beside=T,
         names.arg=c("Int","Age","EF","SC","Time"),
         main="Standard Errors: Wald vs Bootstrap",
         ylab="Standard Error",
         col=c("steelblue", "salmon"),
         legend.text=c("Wald (Hessian)", "Bootstrap SD"))

# Plot 3: 95% CI width comparison across all four bootstrap CI types + Wald
ci_widths = rbind(wald_width, perc_width, basic_width, norm_width, bca_width)
barplot(ci_widths, beside=T,
        names.arg=c("Int","Age","EF","SC","Time"),
        main="95% CI Width Comparison",
        ylab="CI Width",
        col=c("steelblue","salmon","goldenrod","darkgreen","purple"),
        legend.text=c("Wald","Percentile","Basic","Normal","BCa"),
        args.legend=list(cex=0.7))

# Plot 4: P-values comparison
pval_mat = rbind(wald_pvals, boot_pvals)
ymax = max(0.06, max(pval_mat) * 1.2)
barplot( pval_mat, beside=T,
         names.arg=c("Int","Age","EF","SC","Time"),
         main="P-values: Wald vs Bootstrap (H0: beta = 0)",
         ylab="P-value",
         col=c("steelblue","salmon"),
         legend.text=c("Wald","Bootstrap"),
         ylim=c(0, ymax))
abline( h=0.05, col="red", lty=2, lwd=2)

dev.off()
print("Saved: bootstrap_diagnostics.pdf")



# -----------------------------------------------------------------------------
# Justification of N
# 
# For percentile-based bootstrap CIs, the relevant precision is the MC variance
# of the empirical alpha/2 and 1-alpha/2 quantiles.  With N = 1000 and 
# alpha = 0.05, the 2.5th percentile is estimated from order statistics near
# rank 25, and the 97.5th near rank 975, which gives stable tail estimates.
#
# More formally, the MC SE of an empirical quantile is approximately
#   sqrt(p*(1-p) / N) / f(q_p)
# For p = 0.025 and N = 1000, sqrt(p*(1-p)/N) = 0.0049, so quantile rank
# uncertainty is small relative to the distribution.
#
# For the bootstrap mean, MC SE = SD / sqrt(N) is roughly 3% of the
# bootstrap SD with N = 1000.
# -----------------------------------------------------------------------------

mc_se_quantile_rank = sqrt(alpha/2 * (1 - alpha/2) / N)
mc_se_mean = beta_boot_sd / sqrt(N)



# -----------------------------------------------------------------------------
# Save bootstrap summary
# -----------------------------------------------------------------------------

print("Writing summary text files...")

sink("bootstrap_summary.txt")

cat("======================================================================\n")
cat("ST 453 Final Project - Bootstrap Study Summary\n")
cat("Trace Brown\n")
cat("======================================================================\n\n")

cat("BOOTSTRAP PARAMETERS\n")
cat("----------------------------------------------------------------------\n")
cat("Bootstrap method: Pairs/case resampling\n")
cat("                  (resample (X, y) rows with replacement, n at a time)\n")
cat(paste("Number of bootstrap samples: N =", N, "\n"))
cat(paste("Sample size per bootstrap: n =", n, "(matches original)\n"))
cat(paste("Number of parameters: p =", p, "\n"))
cat(paste("Convergence rate:", round(mean(converged_vec, na.rm=T)*100, 2), 
          "%\n"))
cat(paste("Mean events per bootstrap sample:", 
          round(mean(events_vec, na.rm=T), 2),
          "(real data:", sum(y_real), ")\n"))
cat(paste("Random seeds: 1 through", N, "\n\n"))

cat("JUSTIFICATION OF N\n")
cat("----------------------------------------------------------------------\n")
cat(paste("With N =", N, "bootstrap samples, the 2.5%/97.5% quantiles\n"))
cat(paste("are estimated from order statistics near ranks",
          round(N * alpha/2), "and", round(N * (1-alpha/2)), ",\n"))
cat("giving stable percentile-CI boundaries.\n\n")
cat(paste("MC SE of quantile rank position:",
          round(mc_se_quantile_rank, 6), "\n"))
cat("MC SE of bootstrap mean estimate (~SD / sqrt(N)):\n")
for(j in 1:p){
  cat(sprintf("  %-20s: %.6f\n", coef_names[j], mc_se_mean[j]))
}
cat("\n")

cat("BOOTSTRAP COEFFICIENT SUMMARY\n")
cat("----------------------------------------------------------------------\n")
cat(sprintf("%-20s %10s %10s %10s %10s %10s\n",
            "Coefficient", "Estimate", "BootMean", "BootMed", "BootSD", "Bias"))
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%-20s %10.4f %10.4f %10.4f %10.4f %10.4f\n",
              coef_names[j], beta_hat_real[j], beta_boot_mean[j],
              beta_boot_med[j], beta_boot_sd[j], boot_bias[j]))
}
cat("\n")

cat("STANDARD ERRORS: WALD (HESSIAN) vs BOOTSTRAP\n")
cat("----------------------------------------------------------------------\n")
cat(sprintf("%-20s %12s %12s %10s\n", 
            "Coefficient", "Wald SE", "Boot SD", "Ratio"))
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%-20s %12.4f %12.4f %10.4f\n",
              coef_names[j], se_real[j], beta_boot_sd[j],
              beta_boot_sd[j] / se_real[j]))
}
cat("\n")

cat("95% CONFIDENCE INTERVALS (5 METHODS)\n")
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%s\n", coef_names[j]))
  cat(sprintf("  %-12s [%9.4f, %9.4f]  width = %.4f\n",
              "Wald:", wald_lower[j], wald_upper[j], wald_width[j]))
  cat(sprintf("  %-12s [%9.4f, %9.4f]  width = %.4f\n",
              "Percentile:", perc_lower[j], perc_upper[j], perc_width[j]))
  cat(sprintf("  %-12s [%9.4f, %9.4f]  width = %.4f\n",
              "Basic:", basic_lower[j], basic_upper[j], basic_width[j]))
  cat(sprintf("  %-12s [%9.4f, %9.4f]  width = %.4f\n",
              "Normal:", norm_lower[j], norm_upper[j], norm_width[j]))
  cat(sprintf("  %-12s [%9.4f, %9.4f]  width = %.4f\n",
              "BCa:", bca_lower[j], bca_upper[j], bca_width[j]))
}
cat("\n")

cat("BCa BIAS-CORRECTION AND ACCELERATION CONSTANTS\n")
cat("----------------------------------------------------------------------\n")
cat(sprintf("%-20s %10s %10s\n", "Coefficient", "z0", "a"))
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%-20s %10.4f %10.4f\n", coef_names[j], z0[j], a_const[j]))
}
cat("\n")

cat("HYPOTHESIS TESTS (H0: beta_j = 0)\n")
cat("----------------------------------------------------------------------\n")
cat(sprintf("%-20s %10s %14s %14s %6s %6s\n",
            "Coefficient", "z-stat", "Wald p", "Boot p", "Wsig", "Bsig"))
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  z_stat = beta_hat_real[j] / se_real[j]
  
  wsig = ""
  if(wald_pvals[j] < 0.001)  wsig = "***"
  else if(wald_pvals[j] < 0.01)  wsig = "**"
  else if(wald_pvals[j] < 0.05)  wsig = "*"
  
  bsig = ""
  if(boot_pvals[j] < 0.001)  bsig = "***"
  else if(boot_pvals[j] < 0.01)  bsig = "**"
  else if(boot_pvals[j] < 0.05)  bsig = "*"
  
  cat(sprintf("%-20s %10.4f %14s %14s %6s %6s\n",
              coef_names[j], z_stat,
              format(wald_pvals[j], scientific=T, digits=3),
              format(boot_pvals[j], scientific=T, digits=3),
              wsig, bsig))
}
cat("\n")
cat("Significance codes: *** p<.001, ** p<.01, * p<.05\n")
cat("Note: Bootstrap p-values are floored at 1/N (resolution limit).\n")

sink()
print("Saved: bootstrap_summary.txt")



# -----------------------------------------------------------------------------
# Save real-data fit summary with all CI types
# -----------------------------------------------------------------------------

sink("bootstrap_real_data_fit.txt")

cat("======================================================================\n")
cat("Logistic Regression on Heart Failure Data with Bootstrap Inference\n")
cat("======================================================================\n\n")

cat("DATA DESCRIPTION\n")
cat("----------------------------------------------------------------------\n")
cat("Source: UCI Machine Learning Repository\n")
cat("URL: https://archive.ics.uci.edu/dataset/519/heart+failure+clinical+records\n")
cat("Reference: Ahmad et al. (2017), PLoS ONE 12(7):e0181001\n")
cat("Analysis: Chicco and Jurman (2020), BMC Med Inform Decis Mak 20:16\n")
cat(paste("Sample size: n =", n, "\n"))
cat(paste("Number of deaths:", sum(y_real),
          "(", round(mean(y_real)*100, 1), "%)\n\n"))

cat("MODEL\n")
cat("----------------------------------------------------------------------\n")
cat("Response: DEATH_EVENT (0 = survived, 1 = died)\n")
cat("Predictors: Age, Ejection Fraction, Serum Creatinine, Follow-up Time\n\n")
cat("Model: log(p/(1-p)) = b0 + b1*Age + b2*EF + b3*SC + b4*Time\n")
cat("(Predictors standardized to mean 0, SD 1.)\n\n")

cat("BOOTSTRAP STRATEGY\n")
cat("----------------------------------------------------------------------\n")
cat("Method:    Pairs/case bootstrap (resample (X, y) rows with replacement)\n")
cat(paste("N:         ", N, "\n"))
cat(paste("Estimator: Newton-Raphson MLE on each bootstrap sample\n\n"))

cat("COEFFICIENT ESTIMATES\n")
cat("----------------------------------------------------------------------\n")
cat(sprintf("%-20s %10s %12s %12s\n",
            "Coefficient", "Estimate", "Wald SE", "Boot SD"))
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%-20s %10.4f %12.4f %12.4f\n",
              coef_names[j], beta_hat_real[j], se_real[j], beta_boot_sd[j]))
}
cat("\n")

cat("95% PERCENTILE BOOTSTRAP CIs (PRIMARY)\n")
cat("----------------------------------------------------------------------\n")
for(j in 1:p){
  cat(sprintf("%-20s %8.4f  (95%% CI: %8.4f, %8.4f)\n",
              coef_names[j], beta_hat_real[j],
              perc_lower[j], perc_upper[j]))
}
cat("\n")

cat("INTERPRETATION (on standardized scale)\n")
cat("----------------------------------------------------------------------\n")
cat("Each coefficient gives the change in log-odds of death per 1-SD\n")
cat("change in the predictor (with the other three held at their means).\n")
cat("Negative coefficients (Ejection Fraction, Time) are protective.\n")
cat("Positive coefficients (Age, Serum Creatinine) increase mortality risk.\n")

sink()
print("Saved: bootstrap_real_data_fit.txt")


print("")
print("======================================================================")
print("All output files generated successfully!")
print("======================================================================")