n      = sample_sizes,
target = seq(label_range[1], label_range[2], length.out = length(sample_sizes))
) |>
mutate(x = as.integer(round(target * n))) |>
inner_join(ci_wide, by = c("n", "x")) |>
transmute(n = factor(n, levels = n_levels), p_hat, Pi = Upper)
# 5.3 Label background without frame (argument renamed in ggplot2 4.0.0) ----
label_no_border <- if (packageVersion("ggplot2") >= "4.0.0") {
list(linewidth = 0)
} else {
list(label.size = 0)
}
# 5.4 Plot ----
conf_level_percent <- conf_level * 100
axis_breaks <- seq(0, 1, 0.1)
my_plot <- ggplot(ci_data, aes(x = p_hat, y = Pi,
group = interaction(n, Boundary))) +
geom_line(aes(color = n, linetype = Boundary)) +
do.call(geom_label, c(list(data = label_data,
mapping = aes(label = n, color = n, group = NULL),
fill = "white",
show.legend = FALSE),
label_no_border)) +
scale_x_continuous(name   = bquote("Observed Relative Frequency" ~ hat(p) ~ "in the Sample"),
limits = c(0, 1), breaks = axis_breaks, expand = c(0, 0)) +
scale_y_continuous(name   = bquote("Relative Frequency" ~ pi ~ "in the Population"),
limits = c(0, 1), breaks = axis_breaks, expand = c(0, 0)) +
coord_fixed(ratio = 1) +
labs(title    = bquote(.(conf_level_percent) * "% Confidence Intervals for" ~ pi ~
"based on Relative Frequencies" ~ hat(p) ~ "(Clopper-Pearson)"),
color    = "Sample Size n",
linetype = "Boundary") +
theme_minimal()
# 5.5 Show plot ----
print(my_plot)
# 6. Export to various formats with defined sizes ----
save_my_plot <- function(plot, dir, name, device) {
file <- file.path(dir, paste0(name, ".", device))
ggsave(filename = file,
plot     = plot,
device   = if (device == "svg") svglite::svglite else device,
width    = plot_width,
height   = plot_height,
units    = plot_units,
dpi      = plot_dpi)
message("Saved: ", file)
invisible(file)
}
saved_files <- vapply(plot_types, function(dev)
save_my_plot(my_plot, output_dir, plot_name, dev), character(1))
# 7. Clean workspace (Optional) ----
# Removes only the objects created by this script; the rest of the
# workspace stays untouched. Comment out to inspect the results afterwards.
rm(list = intersect(ls(), c(
"required_packages", "missing_packages", "conf_level", "sample_sizes",
"label_range", "plot_width", "plot_height", "plot_units", "plot_dpi",
"plot_name", "plot_types", "output_dir", "get_script_dir",
"cp_lower", "cp_upper", "ci_wide", "n_levels", "ci_data", "label_data",
"label_no_border", "conf_level_percent", "axis_breaks", "my_plot",
"save_my_plot", "saved_files"
)))
###############################################################################
# Author: Boris Kulig, Bertram Schäfer
# Base: R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
# Script-Name: A_Fundamental_Calc_CI_z_vs_t_Intervalls.R
# Script-Version: 2026-10-15
#
# Purpose:
# 1. Calculate and compare Confidence Interval (CI) widths between
#    Standard Normal (z) and t-distribution methods.
# 2. Visualize the convergence of the t-distribution to the normal distribution
#    for different sample sizes.
#
# Location independence:
# The script determines its own directory (Rscript, source(), RStudio "Source"
# or "Run") and saves all plots there. If the script location cannot be
# determined (e.g. code pasted into the console), the current working
# directory is used instead.
#
# Adjust the values in Section 2 to suit your own needs.
###############################################################################
# 0. Workspace ----
# Note: rm(list = ls()) at the start was removed on purpose. When the script
# is called via source(), it would delete the caller's workspace. All objects
# of this script are cleaned up in Section 6 instead (optional).
# 1. Check and load packages ----
required_packages <- c("ggplot2", "dplyr", "tidyr", "svglite")
missing_packages  <- required_packages[!vapply(required_packages,
requireNamespace,
logical(1),
quietly = TRUE)]
if (length(missing_packages) > 0) {
message("Installing missing packages: ", paste(missing_packages, collapse = ", "))
install.packages(missing_packages)
}
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(tidyr)
})
# ==============================================================================
# 2. User Input - ADJUST VALUES HERE ----
# ==============================================================================
# 2.1 Statistical Parameters ----
n_values    <- c(5, 10, 100)  # Sample sizes to analyze (each n >= 2)
se_value    <- 0.063          # Observed or assumed standard error
conf_level  <- 0.95           # Confidence level (0.95 = 95 %)
# 2.2 Plot Output Settings ----
plot_width  <- 36             # Plot width
plot_height <- 12             # Plot height
plot_units  <- "cm"           # "cm", "mm", or "in"
plot_dpi    <- 1200           # Resolution (raster formats only)
plot_name   <- "A_Distribution_Comparison"   # File name without extension
plot_types  <- c("jpeg", "pdf", "svg")       # Output formats
# 2.3 Output Directory ----
# NULL = directory of this script (fallback: working directory).
# Alternatively a fixed path, e.g. "C:/Projects/Output".
output_dir  <- NULL
# ==============================================================================
# 2.4 Input checks ----
stopifnot(
is.numeric(n_values), all(n_values >= 2), all(n_values == round(n_values)),
is.numeric(se_value), se_value > 0,
is.numeric(conf_level), conf_level > 0, conf_level < 1,
plot_units %in% c("cm", "mm", "in")
)
n_values <- sort(unique(n_values))
# 2.5 Determine script directory ----
get_script_dir <- function() {
# a) Started via source() (also RStudio "Source"); checked first, so that a
#    script sourced from another Rscript call resolves to its own folder
for (i in rev(seq_len(sys.nframe()))) {
ofile <- get0("ofile", envir = sys.frame(i), inherits = FALSE)
if (is.character(ofile) && length(ofile) == 1 && nzchar(ofile)) {
return(dirname(normalizePath(ofile)))
}
# b) Started via Rscript / R CMD BATCH
file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
if (length(file_arg) > 0) {
return(dirname(normalizePath(sub("^--file=", "", file_arg[1]))))
}
# c) Executed line by line in RStudio
if (requireNamespace("rstudioapi", quietly = TRUE) &&
rstudioapi::isAvailable()) {
path <- rstudioapi::getSourceEditorContext()$path
if (nzchar(path)) return(dirname(normalizePath(path)))
}
# d) Fallback
warning("Script location could not be determined; ",
"using working directory: ", getwd(), call. = FALSE)
getwd()
}
if (is.null(output_dir)) output_dir <- get_script_dir()
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
output_dir <- normalizePath(output_dir)
message("Output directory: ", output_dir)
# 3. Calculation for Statistics (Table) ----
# 3.1 Derived parameters ----
alpha       <- 1 - conf_level
ci_quantile <- 1 - alpha / 2          # e.g. 0.975 for 95 %
z_score     <- qnorm(ci_quantile)
# 3.2 Results (total CI width = 2 * quantile * SE) ----
results <- data.frame(
n       = n_values,
z_value = z_score,
t_value = qt(ci_quantile, df = n_values - 1)
) |>
mutate(
z_width         = 2 * z_value * se_value,
t_width         = 2 * t_value * se_value,
plus_in_percent = (t_width / z_width - 1) * 100
)
# 3.3 Formatting for display (decimal point, no thousands separator) ----
results_display <- results |>
transmute(
"n"          = sprintf("%d", as.integer(n)),
"z-Value"    = sprintf("%.2f", z_value),
"z-based-CI" = sprintf("%.3f", z_width),
"t-Value"    = sprintf("%.2f", t_value),
"t-based-CI" = sprintf("%.3f", t_width),
"Width Plus" = sprintf("%5.1f %%", plus_in_percent)
)
# 3.4 Print function for the table ----
print_clean_report <- function(data, se, conf) {
line_eq <- strrep("=", 56)
line_mi <- strrep("-", 56)
cat("\n\n", line_eq, "\n",
" Comparison of Total CI Widths (z vs. t)\n",
line_eq, "\n",
sprintf(" Standard Error SE = %.3f | Confidence Level = %.0f %%",
se, conf * 100), "\n",
line_mi, "\n", sep = "")
print(data, row.names = FALSE, right = TRUE)
cat("\n\n")
}
# 4. Calculation for Visualization (Plot) ----
# 4.1 Data preparation ----
x_vals <- seq(-4, 4, length.out = 1000)
make_label   <- function(n) sprintf("n = %d\n(df = %d)", as.integer(n), as.integer(n - 1))
facet_labels <- make_label(n_values)
plot_data_long <- expand.grid(x = x_vals, n = n_values) |>
mutate(
Normal = dnorm(x),
t_Dist = dt(x, df = n - 1),
Sample_Size_Label = factor(make_label(n), levels = facet_labels)
) |>
pivot_longer(c(Normal, t_Dist),
names_to = "Distribution_Type", values_to = "Density") |>
mutate(Distribution_Type = factor(Distribution_Type,
levels = c("Normal", "t_Dist"),
labels = c("Normal Distribution",
"t-Distribution")))
# 4.2 Create plot ----
plot_comparison <- ggplot(plot_data_long,
aes(x = x, y = Density,
color = Distribution_Type,
linetype = Distribution_Type)) +
geom_line(linewidth = 0.8) +
facet_wrap(~Sample_Size_Label, nrow = 1) +
scale_color_manual(values = c("Normal Distribution" = "black",
"t-Distribution"      = "red")) +
scale_linetype_manual(values = c("Normal Distribution" = "solid",
"t-Distribution"      = "dashed")) +
labs(
title    = "Comparison: Standard Normal vs. t-Distribution",
subtitle = "Visualization of convergence with increasing sample size (n)",
x        = "Standard Deviations (z / t)",
y        = "Density",
color    = "Distribution",
linetype = "Distribution"
) +
theme_minimal() +
theme(
plot.title       = element_text(hjust = 0.5, face = "bold", size = 14),
plot.subtitle    = element_text(hjust = 0.5, margin = margin(b = 15)),
legend.position  = "bottom",
strip.text       = element_text(face = "bold", size = 11),
panel.grid.minor = element_blank()
)
# 5. Create and Save Outputs ----
# 5.1 Display plot ----
print(plot_comparison)
# 5.2 Save function ----
save_my_plot <- function(plot, dir, name, device) {
file <- file.path(dir, paste0(name, ".", device))
ggsave(filename = file,
plot     = plot,
device   = if (device == "svg") svglite::svglite else device,
width    = plot_width,
height   = plot_height,
units    = plot_units,
dpi      = plot_dpi)
message("Saved: ", file)
invisible(file)
}
# 5.3 Save plots in all requested formats ----
saved_files <- vapply(plot_types, function(dev)
save_my_plot(plot_comparison, output_dir, plot_name, dev), character(1))
# 5.4 Print table ----
print_clean_report(results_display, se_value, conf_level)
# 6. Clean workspace (Optional) ----
# Removes only the objects created by this script; the rest of the
# workspace stays untouched. Comment out to inspect the results afterwards.
rm(list = intersect(ls(), c(
"required_packages", "missing_packages", "n_values", "se_value",
"conf_level", "plot_width", "plot_height", "plot_units", "plot_dpi",
"plot_name", "plot_types", "output_dir", "get_script_dir", "alpha",
"ci_quantile", "z_score", "results", "results_display",
"print_clean_report", "x_vals", "make_label", "facet_labels", "plot_data_long",
"plot_comparison", "save_my_plot", "saved_files"
)))
# Remove all graphics:
# graphics.off()
# Clear console (visual only):
# cat("\014")
################################################################################
# Autor: Boris Kulig, Bertram Schäfer ----
# Basis: R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life" ----
# Script-Name: D_Fundamental_Calc_CI_Clopper_and_Pearson_Nomogram.R ----
# Script-Version: 2026-10-15 ----
#
# Purpose ----
# The aim of the script is to generate a customized nomogram of the confidence
# intervals of the observed prevalence, taking into account the sample size and
# assuming a binomial distribution.
# The exact method according to Clopper and Pearson (Clopper, C.J., Pearson,
# E.S., 1934. The use of confidence or fiducial limits illustrated in the case
# of the binomial. Biometrika 26, 404–413.
# https://doi.org/10.1093/biomet/26.4.404) is used for calculation and
# presentation in order to deduce the parameter in the population from the
# observed prevalence in the sample (indirect inference). The confidence limits
# for π are calculated using ˆp with a confidence probability of 95%, i.e., the
# 95% CI for π. The numbers on the curves indicate the sample size. The
# confidence limits become narrower and more symmetrical as the sample size n
# increases, since the binomial distribution transitions to a normal
# distribution; for ˆp = 0.5, the confidence interval is symmetrical even for
# small n values. The figure can therefore be used to roughly estimate the n
# required to achieve a certain level of accuracy.
# For an exact calculation of the required sample size for a specific
# application, it is advisable to use the exact values from tables or better to
# calculate the values of the limits with R.
#
# Ziel des Skriptes ist es, ein für Sie angepasstes Nomogramm der
# Vertrauensgrenzen der beobachteten Prävalenz unter Berücksichtigung des
# Stichprobenumfangs unter Annahme der Binomialverteilung zu erzeugen. Zur
# Berechnung und Darstellung wird die exakte Methode nach Clopper und Pearson
# genutzt, um von der beobachteten Prävalenz in der Stichprobe auf den Parameter
# in der Grundgesamtheit zu schließen (indirekter Schluss). Es werden die
# Vertrauensgrenzen für π anhand von ˆp mit einer Vertrauenswahrscheinlichkeit
# von 95%, d. h. das 95%-KI für π berechnet. Die Zahlen an den Kurven bezeichnen
# den Stichprobenumfang. Die Vertrauensgrenzen werden mit zunehmendem
# Stichprobenumfang n enger und symmetrischer, da die Binomialverteilung in eine
# Normalverteilung übergeht, für ˆp = 0.5 ist der Vertrauensbereich auch bei
# kleinen n-Werten symmetrisch. Aus der Abbildung lässt sich somit grob das zum
# Erreichen einer bestimmten Genauigkeit notwendige n abschätzen.
# Für eine exakte Kalkulation des notwendigen Stichprobenumfangs für einen
# spezifischen Anwendungsfall sei angeraten, die exakten Werte aus Tabellen zu
# nutzen, oder besser die Werte der Grenzen mit R zu berechnen.
#
# Calculation note ----
# The Clopper-Pearson limits are computed directly from the beta distribution
#   Lower = qbeta(alpha/2,     x,     n - x + 1)   (0 for x = 0)
#   Upper = qbeta(1 - alpha/2, x + 1, n - x)       (1 for x = n)
# for all x = 0, ..., n. This is mathematically identical to
# binom.test(x, n)$conf.int, but vectorized and therefore much faster.
#
# Location independence ----
# The script determines its own directory (source(), RStudio "Source",
# Rscript, RStudio line-by-line) and saves all plots there. The working
# directory is NOT changed. If the location cannot be determined, the current
# working directory is used.
#
# Steps in the script ----
# 1. Check and load packages (missing ones are installed automatically)
# 2. Definition of parameters for calculation and graphing
#    Please change the values to your needs!
# 3. Input checks and output directory
# 4. Calculation of confidence intervals
# 5. Generating the graph with correct labeling
# 6. Export to various formats with defined sizes
# 7. Clean workspace (optional)
################################################################################
# 1. Check and load packages ----
required_packages <- c("ggplot2", "dplyr", "tidyr", "svglite")
missing_packages  <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages) > 0) {
message("Installing missing packages: ", paste(missing_packages, collapse = ", "))
install.packages(missing_packages)
}
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(tidyr)
})
# ==============================================================================
# 2. USER INPUT - ADJUST VALUES HERE ----
# ==============================================================================
# 2.1 Statistical Parameters ----
conf_level   <- 0.95                                 # Confidence level
sample_sizes <- c(10, 25, 50, 100, 250, 500, 1000)   # Sample sizes n
# 2.2 Label positions ----
# The n labels are placed on the upper limit curves at p_hat positions evenly
# spaced between these two values, so that they do not overlap.
label_range  <- c(0.2, 0.7)
# 2.3 Plot Output Settings ----
# The panel itself is square (coord_fixed), but because the title is long the
# whole plot gets a different size.
plot_width   <- 22
plot_height  <- 15
plot_units   <- "cm"                                  # "cm", "mm", or "in"
plot_dpi     <- 1200                                  # Raster formats only
plot_name    <- "Clopper_Pearson_Nomogram_CIs"        # File name w/o extension
plot_types   <- c("jpeg", "pdf", "svg")               # Output formats
# 2.4 Output Directory ----
# NULL = directory of this script (fallback: working directory).
# Alternatively a fixed path, e.g. "C:/Projects/Output".
output_dir   <- NULL
# ==============================================================================
# 3. Input checks and output directory ----
# 3.1 Input checks ----
stopifnot(
is.numeric(conf_level), conf_level > 0, conf_level < 1,
is.numeric(sample_sizes), all(sample_sizes >= 1),
all(sample_sizes == round(sample_sizes)),
length(label_range) == 2, all(label_range >= 0 & label_range <= 1),
plot_units %in% c("cm", "mm", "in")
)
sample_sizes <- sort(unique(as.integer(sample_sizes)))
# 3.2 Determine script directory ----
get_script_dir <- function() {
# a) Started via source() (also RStudio "Source"); checked first, so that a
#    script sourced from another Rscript call resolves to its own folder
for (i in rev(seq_len(sys.nframe()))) {
ofile <- get0("ofile", envir = sys.frame(i), inherits = FALSE)
if (is.character(ofile) && length(ofile) == 1 && nzchar(ofile)) {
return(dirname(normalizePath(ofile, winslash = "/")))
}
# b) Started via Rscript / R CMD BATCH
file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
if (length(file_arg) > 0) {
return(dirname(normalizePath(sub("^--file=", "", file_arg[1]), winslash = "/")))
}
# c) Executed line by line in RStudio
if (requireNamespace("rstudioapi", quietly = TRUE) &&
rstudioapi::isAvailable()) {
path <- rstudioapi::getSourceEditorContext()$path
if (nzchar(path)) return(dirname(normalizePath(path, winslash = "/")))
}
# d) Fallback
warning("Script location could not be determined; ",
"using working directory: ", getwd(), call. = FALSE)
getwd()
}
if (is.null(output_dir)) output_dir <- get_script_dir()
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
output_dir <- normalizePath(output_dir, winslash = "/")
message("Output directory: ", output_dir)
# 4. Calculation of confidence intervals ----
# 4.1 Exact Clopper-Pearson limits (vectorized) ----
# Shape parameters are kept > 0 (pmax) to avoid warnings; the edge cases
# x = 0 and x = n are then set to their exact values 0 and 1.
cp_lower <- function(x, n, conf) {
a <- (1 - conf) / 2
out <- qbeta(a, pmax(x, 1), n - x + 1)
out[x == 0] <- 0
out
}
cp_upper <- function(x, n, conf) {
a <- (1 - conf) / 2
out <- qbeta(1 - a, x + 1, pmax(n - x, 1))
out[x == n] <- 1
out
}
# 4.2 All possible outcomes x = 0, ..., n for every n ----
ci_wide <- data.frame(
n = rep(sample_sizes, times = sample_sizes + 1L),
x = sequence(sample_sizes + 1L, from = 0L)
) |>
mutate(p_hat = x / n,
Lower = cp_lower(x, n, conf_level),
Upper = cp_upper(x, n, conf_level))
# 5. Generating the graph with correct labeling ----
# 5.1 Long format for plotting ----
# Optional filter: drop p_hat = 0 and 1 so curves do not reach the edges.
# ci_wide <- ci_wide |> filter(p_hat > 0, p_hat < 1)
n_levels <- as.character(sample_sizes)
ci_data  <- ci_wide |>
pivot_longer(c(Lower, Upper), names_to = "Boundary", values_to = "Pi") |>
mutate(n = factor(n, levels = n_levels))
# 5.2 Label positions (closest grid point to the target on the upper curve) ----
label_data <- data.frame(
n      = sample_sizes,
target = seq(label_range[1], label_range[2], length.out = length(sample_sizes))
) |>
mutate(x = as.integer(round(target * n))) |>
inner_join(ci_wide, by = c("n", "x")) |>
transmute(n = factor(n, levels = n_levels), p_hat, Pi = Upper)
# 5.3 Label background without frame (argument renamed in ggplot2 4.0.0) ----
label_no_border <- if (packageVersion("ggplot2") >= "4.0.0") {
list(linewidth = 0)
} else {
list(label.size = 0)
}
# 5.4 Plot ----
conf_level_percent <- conf_level * 100
axis_breaks <- seq(0, 1, 0.1)
my_plot <- ggplot(ci_data, aes(x = p_hat, y = Pi,
group = interaction(n, Boundary))) +
geom_line(aes(color = n, linetype = Boundary)) +
do.call(geom_label, c(list(data = label_data,
mapping = aes(label = n, color = n, group = NULL),
fill = "white",
show.legend = FALSE),
label_no_border)) +
scale_x_continuous(name   = bquote("Observed Relative Frequency" ~ hat(p) ~ "in the Sample"),
limits = c(0, 1), breaks = axis_breaks, expand = c(0, 0)) +
scale_y_continuous(name   = bquote("Relative Frequency" ~ pi ~ "in the Population"),
limits = c(0, 1), breaks = axis_breaks, expand = c(0, 0)) +
coord_fixed(ratio = 1) +
labs(title    = bquote(.(conf_level_percent) * "% Confidence Intervals for" ~ pi ~
"based on Relative Frequencies" ~ hat(p) ~ "(Clopper-Pearson)"),
color    = "Sample Size n",
linetype = "Boundary") +
theme_minimal()
# 5.5 Show plot ----
print(my_plot)
# 6. Export to various formats with defined sizes ----
save_my_plot <- function(plot, dir, name, device) {
file <- file.path(dir, paste0(name, ".", device))
ggsave(filename = file,
plot     = plot,
device   = if (device == "svg") svglite::svglite else device,
width    = plot_width,
height   = plot_height,
units    = plot_units,
dpi      = plot_dpi)
message("Saved: ", file)
invisible(file)
}
saved_files <- vapply(plot_types, function(dev)
save_my_plot(my_plot, output_dir, plot_name, dev), character(1))
# 7. Clean workspace (Optional) ----
# Removes only the objects created by this script; the rest of the
# workspace stays untouched. Comment out to inspect the results afterwards.
rm(list = intersect(ls(), c(
"required_packages", "missing_packages", "conf_level", "sample_sizes",
"label_range", "plot_width", "plot_height", "plot_units", "plot_dpi",
"plot_name", "plot_types", "output_dir", "get_script_dir",
"cp_lower", "cp_upper", "ci_wide", "n_levels", "ci_data", "label_data",
"label_no_border", "conf_level_percent", "axis_breaks", "my_plot",
"save_my_plot", "saved_files"
)))
