Nothing
# ==============================================================================
# stratifyR 2.0 — Professional Shiny Application
# Optimal Survey Stratification
#
# Required packages:
# install.packages(c("shiny", "bslib", "DT", "plotly", "readxl"))
# install.packages("~/Desktop/strataVerse/stratifyR 2.0/stratifyR 2.0-1",
# repos = NULL, type = "source")
# ==============================================================================
library(shiny)
library(bslib)
library(DT)
library(stratifyR)
has_plotly <- requireNamespace("plotly", quietly = TRUE)
has_readxl <- requireNamespace("readxl", quietly = TRUE)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
if (has_plotly) library(plotly)
if (has_ggplot2) library(ggplot2)
# ==============================================================================
# Static data
# ==============================================================================
builtin_choices <- c(
"Anaemia — Iron" = "anaemia|Iron",
"HIES — Income" = "hies|Income",
"Math — Final Marks" = "math|final_marks",
"Sugarcane — Production" = "sugarcane|Production"
)
distr_choices <- c(
"Weibull" = "weibull",
"Gamma" = "gamma",
"Log-Normal" = "lnorm",
"Normal" = "norm",
"Exponential" = "exp",
"Pareto" = "pareto",
"Cauchy" = "cauchy",
"Uniform" = "unif",
"Triangular" = "triangle",
"Right-Triangular" = "rtriangle"
)
distr_params <- list(
weibull = list(names = c("shape", "scale"), defaults = c(2.0, 1.0)),
gamma = list(names = c("shape", "rate"), defaults = c(2.0, 0.5)),
lnorm = list(names = c("meanlog", "sdlog"), defaults = c(0.0, 1.0)),
norm = list(names = c("mean", "sd"), defaults = c(0.0, 1.0)),
exp = list(names = c("rate"), defaults = c(1.0)),
pareto = list(names = c("shape", "scale"), defaults = c(1.5, 1.0)),
cauchy = list(names = c("location", "scale"), defaults = c(0.0, 1.0)),
unif = list(names = c("min", "max"), defaults = c(0.0, 1.0)),
triangle = list(names = c("min", "max", "mode"), defaults = c(0.0, 1.0, 0.5)),
rtriangle = list(names = c("min", "max", "mode"), defaults = c(0.0, 1.0, 1.0))
)
palettes <- list(
classic = c("#4E79A7","#F28E2B","#E15759","#76B7B2","#59A14F",
"#EDC948","#B07AA1","#FF9DA7","#9C755F","#BAB0AC"),
ocean = c("#005F73","#0A9396","#94D2BD","#E9D8A6","#EE9B00",
"#CA6702","#BB3E03","#AE2012","#9B2226","#001219"),
rose = c("#D64045","#E8A838","#9C2542","#F4C261","#B84A62",
"#E07B54","#FF7F5C","#FFA500","#C46E4B","#8B1A1A"),
forest = c("#2D6A4F","#40916C","#52B788","#74C69D","#95D5B2",
"#B7E4C7","#1B4332","#081C15","#D8F3DC","#6B9E7A"),
vivid = c("#E63946","#457B9D","#2DC653","#FF6B35","#7B2D8B",
"#F4A261","#264653","#2A9D8F","#E9C46A","#F77F00")
)
# ==============================================================================
# CSS
# ==============================================================================
app_css <- "
/* ── Header ──────────────────────────────────────────────────────────────── */
.navbar {
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 55%, #24527A 100%) !important;
border-bottom: 3px solid #E9A826 !important;
padding: 0.6rem 1.5rem !important;
box-shadow: 0 2px 12px rgba(0,0,0,0.35) !important;
}
.navbar-brand {
font-size: 1.85rem !important;
font-weight: 900 !important;
letter-spacing: -1px !important;
color: #FFFFFF !important;
display: flex !important;
align-items: center !important;
gap: 10px !important;
}
.brand-logo-wrap {
display: inline-flex;
align-items: center;
gap: 11px;
}
.brand-text-wrap {
display: inline-flex;
flex-direction: column;
line-height: 1;
}
.brand-name {
font-size: 1.75rem;
font-weight: 900;
letter-spacing: -1px;
color: #FFFFFF;
font-style: italic;
}
.brand-name .brand-badge { color: #0D1B2A; }
.brand-badge {
background: #E9A826;
color: #0D1B2A;
font-size: 0.58rem;
font-weight: 900;
padding: 2px 7px;
border-radius: 4px;
margin-left: 7px;
vertical-align: super;
letter-spacing: 0.8px;
}
.brand-sub {
font-size: 0.72rem;
color: rgba(255,255,255,0.5);
font-weight: 400;
letter-spacing: 0.5px;
margin-top: 2px;
}
/* ── Sidebar ──────────────────────────────────────────────────────────────── */
.bslib-sidebar-layout > .sidebar {
background: #FFFFFF !important;
border-right: 1px solid #E2E8F0 !important;
box-shadow: 3px 0 10px rgba(0,0,0,0.06) !important;
}
.sidebar-section-label {
font-size: 0.68rem;
font-weight: 700;
letter-spacing: 1.6px;
text-transform: uppercase;
color: #94A3B8;
margin: 1.1rem 0 0.45rem;
padding-bottom: 0.35rem;
border-bottom: 1px solid #F1F5F9;
}
/* ── Run button ───────────────────────────────────────────────────────────── */
#run_btn {
background: linear-gradient(135deg, #1B3A5C, #24527A) !important;
border: none !important;
color: #FFFFFF !important;
font-weight: 700 !important;
font-size: 0.9rem !important;
letter-spacing: 0.4px !important;
padding: 0.65rem 1rem !important;
border-radius: 8px !important;
width: 100%;
box-shadow: 0 4px 14px rgba(27,58,92,0.4) !important;
transition: all 0.18s ease !important;
}
#run_btn:hover {
transform: translateY(-1px) !important;
box-shadow: 0 6px 18px rgba(27,58,92,0.5) !important;
}
/* ── Status badge ─────────────────────────────────────────────────────────── */
.status-wrap { margin-top: 0.6rem; }
.status-pill {
display: inline-flex;
align-items: center;
gap: 6px;
font-size: 0.78rem;
font-weight: 600;
padding: 4px 12px;
border-radius: 20px;
width: 100%;
justify-content: center;
}
.s-idle { background:#F1F5F9; color:#64748B; }
.s-ready { background:#D1FAE5; color:#065F46; }
.s-error { background:#FEE2E2; color:#991B1B; }
.s-running { background:#DBEAFE; color:#1D4ED8; }
/* ── Metric cards ─────────────────────────────────────────────────────────── */
.metrics-row { display:flex; gap:12px; margin-bottom:1rem; }
.metric-card {
flex: 1;
background: #FFFFFF;
border-radius: 10px;
padding: 1rem 0.8rem 0.85rem;
text-align: center;
box-shadow: 0 1px 4px rgba(0,0,0,0.08), 0 0 0 1px rgba(0,0,0,0.04);
}
.metric-value {
font-size: 1.65rem;
font-weight: 800;
color: #1B3A5C;
line-height: 1.1;
}
.metric-label {
font-size: 0.84rem;
font-weight: 700;
letter-spacing: 0.6px;
text-transform: uppercase;
color: #94A3B8;
margin-top: 4px;
}
/* ── Nav tabs ─────────────────────────────────────────────────────────────── */
.nav-tabs {
border-bottom: none !important;
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 60%, #24527A 100%);
border-radius: 12px 12px 0 0;
padding: 0.55rem 0.6rem 0 0.6rem;
gap: 4px;
display: flex;
flex-wrap: wrap;
}
.nav-tabs .nav-link {
font-weight: 600 !important;
font-size: 0.92rem !important;
color: rgba(255,255,255,0.72) !important;
border: none !important;
padding: 0.65rem 1.15rem !important;
border-radius: 8px 8px 0 0 !important;
transition: background 0.18s, color 0.18s !important;
white-space: nowrap;
display: flex !important;
align-items: center !important;
gap: 0.38rem !important;
background: rgba(255,255,255,0.06) !important;
letter-spacing: 0.015em;
}
.nav-tabs .nav-link svg, .nav-tabs .nav-link i {
font-size: 0.92rem !important;
opacity: 0.8;
}
.nav-tabs .nav-link:hover:not(.active) {
background: rgba(233,168,38,0.18) !important;
color: #F6D280 !important;
}
.nav-tabs .nav-link:hover:not(.active) svg,
.nav-tabs .nav-link:hover:not(.active) i { opacity: 1; }
.nav-tabs .nav-link.active {
background: #E9A826 !important;
color: #0D1B2A !important;
font-weight: 700 !important;
box-shadow: 0 3px 10px rgba(233,168,38,0.35) !important;
}
.nav-tabs .nav-link.active svg,
.nav-tabs .nav-link.active i { opacity: 1 !important; color: #0D1B2A !important; }
.tab-content {
border: none !important;
border-radius: 0 0 12px 12px !important;
background: #FFFFFF !important;
box-shadow: 0 4px 16px rgba(13,27,42,0.10) !important;
}
/* ── Empty state ──────────────────────────────────────────────────────────── */
.empty-state {
display: flex;
flex-direction: column;
align-items: center;
justify-content: center;
min-height: 55vh;
color: #94A3B8;
text-align: center;
gap: 0.75rem;
}
.empty-icon { font-size: 3.5rem; opacity: 0.35; }
.empty-title { font-size: 1rem; font-weight: 700; color: #64748B; }
.empty-sub { font-size: 0.82rem; }
/* ── Accordion ────────────────────────────────────────────────────────────── */
.accordion-button {
background: linear-gradient(135deg, #134E4A 0%, #0F766E 100%) !important;
color: #FFFFFF !important;
font-weight: 700 !important;
font-size: 1.05rem !important;
letter-spacing: 0.2px !important;
padding: 0.90rem 1.15rem !important;
border-radius: 8px !important;
box-shadow: 0 2px 8px rgba(15,118,110,0.30) !important;
}
.accordion-button::after {
filter: invert(1) brightness(2) !important;
}
.accordion-button:not(.collapsed) {
background: linear-gradient(135deg, #042F2E 0%, #134E4A 100%) !important;
color: #E9A826 !important;
box-shadow: 0 2px 10px rgba(4,47,46,0.45) !important;
}
.accordion-button:not(.collapsed)::after {
filter: invert(0) sepia(1) saturate(4) hue-rotate(5deg) brightness(1.1) !important;
}
.accordion-item {
border: none !important;
margin-bottom: 6px !important;
border-radius: 10px !important;
overflow: visible !important;
box-shadow: 0 2px 8px rgba(27,58,92,0.12) !important;
}
.accordion-button {
border-radius: 10px !important;
}
.accordion-button:not(.collapsed) {
border-radius: 10px 10px 0 0 !important;
}
.accordion-collapse .accordion-body {
border-radius: 0 0 10px 10px !important;
}
.accordion-body {
padding: 0.85rem 1rem !important;
background: #FAFCFF !important;
color: #1E293B !important;
overflow: visible !important;
}
/* Ensure select dropdowns always render above everything */
.selectize-dropdown,
.dropdown-menu,
.selectize-dropdown-content {
z-index: 99999 !important;
position: absolute !important;
}
/* Ensure all form controls inside accordion are fully visible */
.accordion-body .form-select,
.accordion-body select,
.accordion-body .form-control,
.accordion-body input[type=text],
.accordion-body input[type=number] {
background-color: #FFFFFF !important;
color: #1E293B !important;
border: 1.5px solid #CBD5E1 !important;
border-radius: 6px !important;
font-size: 0.83rem !important;
}
.accordion-body label,
.accordion-body .form-check-label,
.accordion-body .shiny-input-container > label {
color: #374151 !important;
font-size: 0.8rem !important;
font-weight: 600 !important;
}
.accordion-body .form-check-input { border-color: #4A9EE8 !important; }
.accordion-body .radio label { color: #374151 !important; font-weight: 500 !important; }
/* ── Progress bar (sidebar shimmer) ──────────────────────────────────────── */
.progress-wrap {
margin-top: 10px;
border-radius: 8px;
overflow: hidden;
background: #E2E8F0;
height: 10px;
}
.progress-shimmer {
height: 100%;
width: 100%;
background: linear-gradient(90deg,
#134E4A 0%, #0F766E 30%, #E9A826 55%, #0F766E 75%, #134E4A 100%);
background-size: 300% 100%;
animation: shimmer 1.8s infinite linear;
border-radius: 6px;
}
@keyframes shimmer {
0% { background-position: 100% 0; }
100% { background-position: -100% 0; }
}
.progress-label {
font-size: 0.72rem;
color: #0F766E;
font-weight: 600;
margin-top: 5px;
text-align: center;
}
/* ── Shiny withProgress modal — custom theme ─────────────────────────────── */
#shiny-progress-container {
position: fixed !important;
bottom: 24px !important;
right: 24px !important;
left: auto !important;
width: 340px !important;
z-index: 9999;
}
.shiny-progress {
background: #0D1B2A !important;
border-radius: 14px !important;
padding: 18px 22px 16px !important;
box-shadow: 0 12px 40px rgba(0,0,0,0.55) !important;
border: 1.5px solid rgba(233,168,38,0.35) !important;
}
.shiny-progress .progress-message {
color: #E9A826 !important;
font-weight: 700 !important;
font-size: 0.92rem !important;
margin-bottom: 6px !important;
}
.shiny-progress .progress-detail {
color: rgba(255,255,255,0.65) !important;
font-size: 0.80rem !important;
margin-bottom: 10px !important;
}
.shiny-progress .progress {
height: 10px !important;
border-radius: 6px !important;
background: rgba(255,255,255,0.12) !important;
overflow: hidden !important;
margin-bottom: 0 !important;
}
.shiny-progress .progress-bar {
background: linear-gradient(90deg, #134E4A, #0F766E 50%, #E9A826) !important;
transition: width 0.4s cubic-bezier(0.4,0,0.2,1) !important;
border-radius: 6px !important;
}
/* ── Summary table ────────────────────────────────────────────────────────── */
.summary-meta {
font-size: 0.78rem;
color: #64748B;
margin-bottom: 0.75rem;
display: flex;
gap: 1.2rem;
flex-wrap: wrap;
}
.meta-item { display: flex; gap: 5px; align-items: center; }
.meta-badge {
background: #EFF6FF;
color: #1D4ED8;
font-size: 0.7rem;
font-weight: 700;
padding: 2px 8px;
border-radius: 10px;
}
.meta-badge.green { background:#D1FAE5; color:#065F46; }
.meta-badge.amber { background:#FEF3C7; color:#92400E; }
/* ── Design Comparison tab ────────────────────────────────────────────────── */
.dc-cards {
display: flex; flex-wrap: wrap; gap: 0.9rem; margin-bottom: 1.1rem;
}
.dc-card {
flex: 1 1 160px;
background: #FFFFFF;
border: 1.5px solid #E2E8F0;
border-radius: 10px;
padding: 0.85rem 1.1rem;
text-align: center;
}
.dc-card .dc-val {
font-size: 1.55rem; font-weight: 800; color: #1B3A5C;
line-height: 1.1; margin-bottom: 0.2rem;
}
.dc-card .dc-lbl {
font-size: 0.72rem; color: #64748B; letter-spacing: 0.3px;
}
.dc-card.dc-highlight { border-color: #059669; background: #F0FDF4; }
.dc-card.dc-highlight .dc-val { color: #065F46; }
.dc-card.dc-amber { border-color: #D97706; background: #FFFBEB; }
.dc-card.dc-amber .dc-val { color: #92400E; }
.dc-table {
width: 100%; border-collapse: collapse;
font-size: 0.875rem; margin: 0.4rem 0 1.1rem;
}
.dc-table th {
background: #1B3A5C; color: #E9A826;
padding: 8px 14px; text-align: left;
font-weight: 700; font-size: 0.80rem; letter-spacing: 0.3px;
}
.dc-table th:not(:first-child) { text-align: right; }
.dc-table td {
padding: 8px 14px; border-bottom: 1px solid #E2E8F0; color: #334155;
}
.dc-table td:not(:first-child) { text-align: right; font-family: 'Courier New', monospace; }
.dc-table tr.dc-row-opt td { background: #F0FDF4; font-weight: 700; color: #065F46; }
.dc-table tr.dc-row-prop td { background: #FAFAFA; }
.dc-table .dc-row-opt td:first-child::before { content: '★ '; }
.dc-saving-banner {
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 100%);
color: #FFFFFF; border-radius: 10px;
padding: 1rem 1.4rem; margin-bottom: 1.1rem;
display: flex; align-items: center; gap: 1.4rem; flex-wrap: wrap;
}
.dc-saving-banner .dsb-stat {
text-align: center; flex: 1; min-width: 100px;
}
.dc-saving-banner .dsb-val {
font-size: 1.9rem; font-weight: 800; color: #E9A826; line-height: 1;
}
.dc-saving-banner .dsb-lbl {
font-size: 0.73rem; color: rgba(255,255,255,0.72); margin-top: 2px;
}
.dc-saving-banner .dsb-div {
width: 1px; background: rgba(255,255,255,0.20); align-self: stretch;
}
/* ── Introduction page ────────────────────────────────────────────────────── */
.intro-hero {
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 60%, #24527A 100%);
color: #FFFFFF;
border-radius: 12px;
padding: 1.8rem 2.2rem 1.6rem;
margin-bottom: 1.4rem;
position: relative;
overflow: hidden;
}
.intro-hero::before {
content: '';
position: absolute; top: 0; right: 0; bottom: 0;
width: 38%;
background: radial-gradient(ellipse at top right, rgba(233,168,38,0.18) 0%, transparent 70%);
}
.intro-hero-title {
font-size: 1.6rem; font-weight: 800;
color: #E9A826; letter-spacing: -0.4px; margin-bottom: 0.35rem;
}
.intro-hero-sub {
font-size: 0.97rem; color: rgba(255,255,255,0.82);
max-width: 680px; line-height: 1.55;
}
.intro-section {
background: #FFFFFF;
border: 1.5px solid #E2E8F0;
border-radius: 10px;
padding: 1.35rem 1.6rem;
margin-bottom: 1.1rem;
}
.intro-section-title {
font-size: 1.02rem; font-weight: 700;
color: #1B3A5C; margin-bottom: 0.65rem;
display: flex; align-items: center; gap: 0.5rem;
}
.intro-section-title .num-badge {
background: #1B3A5C; color: #E9A826;
font-size: 0.78rem; font-weight: 800;
width: 22px; height: 22px; border-radius: 50%;
display: flex; align-items: center; justify-content: center;
flex-shrink: 0;
}
.intro-body { font-size: 0.91rem; color: #334155; line-height: 1.70; }
.intro-body p { margin: 0 0 0.7rem; }
.intro-body p:last-child { margin-bottom: 0; }
.intro-math {
background: #F8FAFC;
border-left: 3.5px solid #1B3A5C;
border-radius: 0 8px 8px 0;
padding: 0.7rem 1.2rem;
margin: 0.6rem 0;
font-size: 0.89rem;
overflow-x: auto;
}
.intro-key-result {
background: linear-gradient(135deg, #EFF6FF 0%, #F8FAFF 100%);
border: 1.5px solid #BFDBFE;
border-radius: 8px;
padding: 0.75rem 1.2rem;
margin: 0.7rem 0 0;
font-size: 0.88rem; color: #1E40AF;
font-weight: 600;
}
.intro-example-table {
width: 100%; border-collapse: collapse;
font-size: 0.85rem; margin: 0.6rem 0;
}
.intro-example-table th {
background: #1B3A5C; color: #E9A826;
padding: 7px 12px; text-align: center;
font-weight: 700; font-size: 0.80rem;
}
.intro-example-table td {
padding: 6px 12px; text-align: center;
border-bottom: 1px solid #E2E8F0; color: #334155;
}
.intro-example-table tr:last-child td {
border-bottom: none; font-weight: 700;
background: #F1F5F9; color: #1B3A5C;
}
.intro-example-table tr:nth-child(even) td { background: #F8FAFC; }
.intro-gain-box {
display: flex; gap: 1rem; flex-wrap: wrap; margin-top: 0.8rem;
}
.intro-gain-card {
flex: 1; min-width: 150px;
background: #FFFFFF; border: 1.5px solid #E2E8F0;
border-radius: 8px; padding: 0.7rem 1rem;
text-align: center;
}
.intro-gain-card .gval {
font-size: 1.4rem; font-weight: 800; color: #1B3A5C;
}
.intro-gain-card .glbl {
font-size: 0.73rem; color: #64748B; margin-top: 1px;
}
.intro-gain-card.highlight { border-color: #059669; background: #F0FDF4; }
.intro-gain-card.highlight .gval { color: #065F46; }
.intro-workflow {
display: grid; grid-template-columns: repeat(auto-fit, minmax(150px,1fr));
gap: 0.7rem; margin-top: 0.6rem;
}
.intro-step {
background: #F8FAFC; border: 1px solid #E2E8F0;
border-radius: 8px; padding: 0.75rem 0.9rem;
font-size: 0.84rem; color: #334155;
position: relative;
}
.intro-step .step-num {
font-size: 1.1rem; font-weight: 800; color: #E9A826;
display: block; margin-bottom: 0.2rem;
}
.intro-ref {
font-size: 0.78rem; color: #94A3B8;
border-top: 1px solid #E2E8F0; margin-top: 1rem;
padding-top: 0.6rem; line-height: 1.55;
}
/* ── H-Selector tab ───────────────────────────────────────────────────────── */
.hsel-banner {
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 100%);
color: #FFFFFF; border-radius: 10px;
padding: 1rem 1.5rem; margin-bottom: 1.1rem;
display: flex; align-items: center; justify-content: space-between;
flex-wrap: wrap; gap: 1rem;
}
.hsel-banner-title { font-size: 1.1rem; font-weight: 700; color: #E9A826; margin-bottom: 0.25rem; }
.hsel-banner-text { font-size: 0.87rem; color: rgba(255,255,255,0.80); max-width: 640px; }
.hsel-rec-box {
background: #F0FDF4; border: 2px solid #059669;
border-radius: 10px; padding: 0.85rem 1.4rem;
margin-bottom: 1rem; display: flex; align-items: center; gap: 1.4rem; flex-wrap: wrap;
}
.hsel-rec-h { font-size: 2.2rem; font-weight: 800; color: #065F46; line-height: 1; }
.hsel-rec-label { font-size: 0.82rem; color: #047857; font-weight: 600; margin-top: 2px; }
.hsel-rec-note { font-size: 0.82rem; color: #475569; }
.hsel-table {
width: 100%; border-collapse: collapse; font-size: 0.875rem; margin: 0.4rem 0 0.5rem;
}
.hsel-table th {
background: #1B3A5C; color: #E9A826;
padding: 8px 14px; font-weight: 700; font-size: 0.80rem; text-align: center;
}
.hsel-table td { padding: 7px 14px; border-bottom: 1px solid #E2E8F0; text-align: center; color: #334155; }
.hsel-table tr.hsel-cur-row td { background: #EFF6FF; font-weight: 700; color: #1D4ED8; }
.hsel-table tr.hsel-rec-row td { background: #F0FDF4; font-weight: 700; color: #065F46; }
.hsel-table tr.hsel-degen-row td { background: #FFFBEB; color: #92400E; font-style: italic; }
.hsel-table tr.hsel-stagnant-row td { background: #F8FAFC; color: #94A3B8; font-style: italic; }
.hsel-note { font-size: 0.76rem; color: #94A3B8; margin-top: 0.5rem; }
/* ── Boundary Methods Comparison ─────────────────────────────────────────── */
.bmc-banner {
background: linear-gradient(135deg, #0D1B2A 0%, #1B3A5C 100%);
color: #FFFFFF; border-radius: 10px;
padding: 1rem 1.5rem; margin-bottom: 1.1rem;
display: flex; align-items: center; justify-content: space-between;
flex-wrap: wrap; gap: 1rem;
}
.bmc-banner-title { font-size: 1.1rem; font-weight: 700; color: #E9A826; margin-bottom: 0.25rem; }
.bmc-banner-text { font-size: 0.87rem; color: rgba(255,255,255,0.80); max-width: 640px; }
.bmc-table {
width: 100%; border-collapse: collapse; font-size: 0.875rem; margin: 0.4rem 0 0.5rem;
}
.bmc-table th {
background: #1B3A5C; color: #E9A826;
padding: 8px 14px; font-weight: 700; font-size: 0.80rem; text-align: left;
}
.bmc-table th:not(:first-child) { text-align: right; }
.bmc-table td { padding: 8px 14px; border-bottom: 1px solid #E2E8F0; color: #334155; }
.bmc-table td:not(:first-child) { text-align: right; font-family: 'Courier New', monospace; }
.bmc-table tr.bmc-row-dp td { background: #F0FDF4; font-weight: 700; color: #065F46; }
.bmc-table tr.bmc-row-fail td { color: #94A3B8; font-style: italic; }
.bmc-badge-opt {
background: #D1FAE5; color: #065F46;
font-size: 0.68rem; font-weight: 700;
padding: 2px 7px; border-radius: 10px; margin-left: 6px;
}
.bmc-badge-fail {
background: #FEF2F2; color: #991B1B;
font-size: 0.68rem; font-weight: 700;
padding: 2px 7px; border-radius: 10px; margin-left: 6px;
}
/* ── Sample Size Calculator ──────────────────────────────────────────────── */
.ssc-result-box {
background: linear-gradient(135deg, #EFF6FF 0%, #DBEAFE 100%);
border: 2px solid #3B82F6; border-radius: 12px;
padding: 1.4rem 2rem; text-align: center; min-width: 180px;
}
.ssc-result-n { font-size: 3rem; font-weight: 900; color: #1B3A5C; line-height: 1; }
.ssc-result-lbl { font-size: 0.88rem; color: #475569; margin-top: 0.3rem; }
.ssc-results-row {
display: flex; gap: 1rem; flex-wrap: wrap; margin-bottom: 1rem; align-items: flex-start;
}
.ssc-input-panel {
background: #FFFFFF; border: 1.5px solid #E2E8F0; border-radius: 10px;
padding: 1rem 1.2rem; margin-bottom: 1rem;
display: flex; align-items: flex-end; gap: 1.4rem; flex-wrap: wrap;
}
/* ── Cost-Constrained Optimisation ──────────────────────────────────────── */
.cc-cost-grid {
display: grid; grid-template-columns: repeat(auto-fit, minmax(120px, 1fr));
gap: 0.6rem; margin: 0.8rem 0 1rem;
}
.cc-result-banner {
background: linear-gradient(135deg, #F0FDF4 0%, #DCFCE7 100%);
border: 2px solid #059669; border-radius: 10px;
padding: 1rem 1.4rem; margin-bottom: 1rem;
display: flex; gap: 2.5rem; flex-wrap: wrap; align-items: center;
}
.cc-result-stat { text-align: center; }
.cc-result-val { font-size: 1.8rem; font-weight: 800; color: #065F46; line-height: 1; }
.cc-result-lbl { font-size: 0.76rem; color: #047857; margin-top: 2px; }
.cc-result-banner { display: flex; gap: 1.6rem; flex-wrap: wrap; align-items: center; }
.cc-cost-card { background: #F8FAFC; border: 1px solid #E2E8F0; border-radius: 8px;
padding: 0.5rem 0.7rem; }
.cc-cost-card label { font-size: 0.72rem; font-weight: 600; color: #64748B;
display: block; margin-bottom: 0.25rem; }
/* ── Methodology panel ──────────────────────────────────────────────────────── */
.method-panel { background: #F8FAFC; border: 1.5px solid #CBD5E1;
border-radius: 10px; padding: 1rem 1.4rem; margin-bottom: 1rem; }
.method-panel summary { font-size: 0.82rem; font-weight: 700; color: #1B3A5C;
cursor: pointer; user-select: none; list-style: none; display: flex;
align-items: center; gap: 0.4rem; }
.method-panel summary::before { content: '▶'; font-size: 0.65rem; color: #64748B;
transition: transform 0.2s; display: inline-block; }
details[open] > summary::before { transform: rotate(90deg); }
.method-panel .method-body { margin-top: 0.8rem; font-size: 0.82rem;
color: #334155; line-height: 1.6; }
.method-panel .method-body p { margin: 0.4rem 0; }
.method-panel .method-body .method-formula { background: #EFF6FF;
border-left: 3px solid #3B82F6; border-radius: 0 6px 6px 0;
padding: 0.4rem 0.8rem; margin: 0.5rem 0; font-family: monospace;
font-size: 0.85rem; color: #1E3A5F; }
.method-panel .method-ref { font-size: 0.75rem; color: #64748B;
margin-top: 0.6rem; border-top: 1px solid #E2E8F0; padding-top: 0.5rem; }
/* ── Plot title input ─────────────────────────────────────────────────────── */
.title-input-wrap label { font-size: 0.78rem; font-weight: 600; color:#475569; }
.title-input-wrap input { font-size: 0.88rem; border-radius: 6px !important; }
/* ── Download button ──────────────────────────────────────────────────────── */
#dl_results {
width: 100%;
font-size: 0.8rem !important;
border-color: #CBD5E1 !important;
color: #475569 !important;
background: #F8FAFC !important;
}
#dl_results:hover { background: #F1F5F9 !important; }
/* ── Generate R Code button ───────────────────────────────────────────────── */
#gen_code_btn {
width: 100%;
font-size: 0.8rem !important;
border-color: #BFDBFE !important;
color: #1D4ED8 !important;
background: #EFF6FF !important;
margin-top: 6px;
}
#gen_code_btn:hover { background: #DBEAFE !important; }
/* ── R Code modal ─────────────────────────────────────────────────────────── */
.r-code-block {
background: #0F172A;
color: #E2E8F0;
font-family: 'Courier New', Consolas, monospace;
font-size: 0.82rem;
line-height: 1.65;
border-radius: 8px;
padding: 1.1rem 1.3rem;
white-space: pre;
overflow-x: auto;
max-height: 68vh;
overflow-y: auto;
}
.r-code-copy-btn {
font-size: 0.78rem !important;
padding: 4px 14px !important;
}
/* ── Main area background ─────────────────────────────────────────────────── */
.bslib-sidebar-layout > .main { background: #F8FAFD; padding: 1.2rem; }
/* ── Plot container ───────────────────────────────────────────────────────── */
.plot-container {
background: #FFFFFF;
border-radius: 10px;
padding: 0.5rem;
box-shadow: 0 1px 4px rgba(0,0,0,0.07);
}
/* ── Data preview ─────────────────────────────────────────────────────────── */
.data-preview {
background: #F8FAFD;
border: 1px solid #E2E8F0;
border-radius: 6px;
padding: 8px 10px;
font-size: 0.76rem;
color: #475569;
margin-top: 4px;
}
.data-preview strong { color: #1B3A5C; }
/* ── 2D Info panel below plot ─────────────────────────────────────────────── */
.info-dist-card {
background: linear-gradient(135deg, #EEF2FF 0%, #F8FAFF 100%);
border: 1.5px solid #B0C4DE;
border-radius: 10px;
padding: 0.9rem 1.2rem;
margin-top: 0.8rem;
}
.info-dist-name {
font-size: 1.30rem;
font-weight: 800;
color: #1a1a6e;
margin-bottom: 0.4rem;
letter-spacing: -0.2px;
}
.info-eqn {
background: #FFFFFF;
border: 1.5px solid #C8D8EE;
border-radius: 8px;
padding: 0.5rem 1.4rem 0.3rem;
margin: 0.45rem 0;
display: block;
text-align: center;
overflow-x: auto;
}
.info-eqn .MathJax_Display, .info-eqn mjx-container {
font-size: 1.25rem !important;
}
.info-params {
display: flex;
flex-wrap: wrap;
gap: 0.6rem;
margin-top: 0.6rem;
}
.param-badge {
background: #1B3A5C;
color: #E9A826;
font-size: 0.95rem;
font-weight: 700;
padding: 6px 16px;
border-radius: 20px;
font-family: 'Courier New', monospace;
letter-spacing: 0.2px;
white-space: nowrap;
}
.info-strata {
margin-top: 0.85rem;
display: flex;
flex-direction: column;
gap: 0.42rem;
}
.strata-row {
display: flex;
align-items: center;
gap: 0.65rem;
background: #FFFFFF;
border: 1px solid #E2E8F0;
border-radius: 9px;
padding: 0.6rem 1.1rem;
font-size: 1.05rem;
}
.strata-swatch {
width: 15px;
height: 15px;
border-radius: 4px;
flex-shrink: 0;
display: inline-block;
}
.strata-label {
font-weight: 800;
color: #1B3A5C;
min-width: 90px;
font-size: 1.05rem;
}
.strata-interval {
color: #C0392B;
font-weight: 600;
min-width: 160px;
font-family: 'Courier New', monospace;
font-size: 1.00rem;
white-space: nowrap;
}
.strata-stat {
color: #475569;
font-family: 'Courier New', monospace;
font-size: 1.05rem;
white-space: nowrap;
padding: 4px 13px;
background: #F1F5F9;
border-radius: 5px;
}
"
# ==============================================================================
# Theme
# ==============================================================================
app_theme <- bs_theme(
version = 5,
bg = "#F8FAFD",
fg = "#1E293B",
primary = "#1B3A5C",
secondary = "#4A6FA5",
success = "#0F9B6E",
warning = "#D97706",
danger = "#DC2626",
font_scale = 0.94
)
# ==============================================================================
# UI
# ==============================================================================
ui <- page_sidebar(
theme = app_theme,
window_title = "stratifyR 2.0",
tags$head(
tags$style(HTML(app_css)),
withMathJax()
),
title = div(class = "brand-logo-wrap",
# Logo: circle with density curve + 3 strata — clean, minimal, professional
HTML('
<svg width="48" height="48" viewBox="0 0 130 130"
xmlns="http://www.w3.org/2000/svg"
style="flex-shrink:0; filter:drop-shadow(0 2px 10px rgba(0,0,0,0.50));">
<defs>
<clipPath id="lgClip">
<polygon points="65,9 113.5,37 113.5,93 65,121 16.5,93 16.5,37"/>
</clipPath>
</defs>
<polygon points="65,3 118.7,34 118.7,96 65,127 11.3,96 11.3,34" fill="#1A2D50"/>
<g clip-path="url(#lgClip)">
<rect x="7" y="7" width="34" height="116" fill="#E9A826" fill-opacity="0.28"/>
<rect x="41" y="7" width="42" height="116" fill="#2DD4BF" fill-opacity="0.22"/>
<rect x="83" y="7" width="40" height="116" fill="#818CF8" fill-opacity="0.28"/>
<line x1="41" y1="7" x2="41" y2="123" stroke="rgba(255,255,255,0.50)" stroke-width="1.2" stroke-dasharray="3.5,3"/>
<line x1="83" y1="7" x2="83" y2="123" stroke="rgba(255,255,255,0.50)" stroke-width="1.2" stroke-dasharray="3.5,3"/>
<path d="M 8,96 C 12,96 16,88 22,74 C 27,62 32,44 39,31 C 44,22 49,19 55,25 C 62,34 68,52 76,66 C 83,78 90,87 98,92 C 106,95 114,96 122,96"
fill="none" stroke="rgba(255,255,255,0.95)" stroke-width="2.8" stroke-linecap="round"/>
</g>
<polygon points="65,5 117,35 117,95 65,125 13,95 13,35" fill="none" stroke="#E9A826" stroke-width="2.6"/>
</svg>
'),
div(class = "brand-text-wrap",
HTML('<div class="brand-name">stratifyR<span class="brand-badge">2.0-1</span></div>'),
div(class = "brand-sub", "Optimal Survey Stratification")
)
),
# ── Sidebar ─────────────────────────────────────────────────────────────────
sidebar = sidebar(
width = 310,
accordion(
open = c("acc_data", "acc_strat", "acc_plot"),
# ── 1. Data Source ───────────────────────────────────────────────────────
accordion_panel(
title = tagList(icon("database", style="margin-right:6px"), "Data Source"),
value = "acc_data",
radioButtons("data_source", NULL,
choices = c("Upload file" = "upload",
"Built-in dataset" = "builtin",
"Specify distribution" = "distr"),
selected = "builtin"),
# Upload
conditionalPanel("input.data_source == 'upload'",
fileInput("file_upload", NULL,
accept = c(".csv", ".xlsx", ".xls", ".txt"),
placeholder = "CSV, Excel or TXT…",
buttonLabel = "Browse"),
uiOutput("ui_col_select")
),
# Built-in
conditionalPanel("input.data_source == 'builtin'",
selectInput("builtin_ds", NULL, choices = builtin_choices),
uiOutput("ui_builtin_info")
),
# Distribution
conditionalPanel("input.data_source == 'distr'",
selectInput("distr_name", NULL, choices = distr_choices, selected = "weibull"),
uiOutput("ui_distr_params"),
fluidRow(
column(6, numericInput("initval", "Min value", value = 0, step = 0.1)),
column(6, numericInput("dist_range","Range (max−min)", value = 10, min = 0.001, step = 0.1))
),
numericInput("distr_N", "Population size (N)", value = 10000, min = 100, step = 500)
)
),
# ── 2. Stratification ────────────────────────────────────────────────────
accordion_panel(
title = tagList(icon("sliders", style="margin-right:6px"), "Stratification"),
value = "acc_strat",
fluidRow(
column(6, numericInput("h_strata", "Strata (h)", value = 3, min = 2, max = 10, step = 1)),
column(6, numericInput("n_sample", "Sample (n)", value = 300, min = 10, step = 10))
),
conditionalPanel("input.data_source != 'distr'",
checkboxInput("use_N", "Set population size (N)", value = FALSE),
conditionalPanel("input.use_N",
numericInput("pop_N", NULL, value = 10000, min = 100, step = 500)
)
),
selectInput("solver", "Solver",
choices = c("DP — Dynamic Programming (default)" = "dp",
"COBYLA — Fast Multi-start" = "cobyla",
"GLOBAL — DIRECT-L + COBYLA (best for skewed)" = "global"),
selected = "dp",
selectize = FALSE),
div(style = "font-size:0.75rem; color:#94A3B8; margin-top:2px;",
tags$strong("DP"), " — exact global optimum (recommended for empirical data). ",
tags$strong("COBYLA"), " — fast gradient-free search. ",
tags$strong("GLOBAL"), " — DIRECT-L global phase + COBYLA refinement; best for right-skewed or heavy-tailed distributions.")
),
# ── 3. Plot Options ──────────────────────────────────────────────────────
accordion_panel(
title = tagList(icon("chart-line", style="margin-right:6px"), "Plot Options"),
value = "acc_plot",
div(class = "title-input-wrap",
textInput("plot_title", "Plot title",
placeholder = "Leave blank for auto-generated title")
),
selectInput("palette", "Colour palette",
choices = c("Classic (Tableau)" = "classic",
"Ocean Blues" = "ocean",
"Rose & Amber" = "rose",
"Forest Greens" = "forest",
"Vivid" = "vivid"),
selected = "classic",
selectize = FALSE)
)
),
hr(style = "margin:0.9rem 0 0.75rem; border-color:#E2E8F0;"),
actionButton("run_btn", tagList(icon("play"), " Run Stratification"),
class = "btn btn-primary"),
div(class = "status-wrap", uiOutput("ui_status")),
hr(style = "margin:0.9rem 0 0.6rem; border-color:#E2E8F0;"),
downloadButton("dl_results", tagList(icon("download"), " Export Results (CSV)"),
class = "btn btn-outline-secondary btn-sm"),
actionButton("gen_code_btn", tagList(icon("code"), " Generate R Code"),
class = "btn btn-outline-primary btn-sm")
),
# ── Main content ─────────────────────────────────────────────────────────────
uiOutput("ui_main")
)
# ==============================================================================
# compare_designs — local implementation (does not require package reinstall)
# ==============================================================================
.S2_from_strata_app <- function(object) {
pop <- object$data_internal
if (!is.null(pop) && length(pop) >= 2L) {
N <- length(pop); mu <- mean(pop)
return(sum((pop - mu)^2) / N)
}
x_lo <- object$maxval * object$initval
x_hi <- object$maxval * object$finval
if (!is.finite(x_lo) || !is.finite(x_hi) || x_hi <= x_lo)
return(sum(object$Wh * object$Vh))
x_seq <- seq(x_lo, x_hi, length.out = 1024L)
params <- tryCatch({
fit <- object$fit
if (!is.null(fit$estimate)) fit$estimate
else if (!is.null(fit$fit$estimate)) fit$fit$estimate
else NULL
}, error = function(e) NULL)
if (is.null(params)) return(sum(object$Wh * object$Vh))
dens <- tryCatch(
switch(object$distr,
norm = stats::dnorm(x_seq, mean = params["mean"], sd = params["sd"]),
lnorm = stats::dlnorm(x_seq, meanlog = params["meanlog"], sdlog = params["sdlog"]),
gamma = stats::dgamma(x_seq, shape = params["shape"], rate = params["rate"]),
weibull = stats::dweibull(x_seq, shape = params["shape"], scale = params["scale"]),
exp = stats::dexp(x_seq, rate = params["rate"]),
cauchy = stats::dcauchy(x_seq, location = params["location"],scale = params["scale"]),
unif = stats::dunif(x_seq, min = params["min"], max = params["max"]),
rep(1 / (x_hi - x_lo), 1024L)
),
error = function(e) rep(1 / (x_hi - x_lo), 1024L)
)
dens[!is.finite(dens)] <- 0
dx <- diff(x_seq)
trap <- function(f) sum(dx * (head(f, -1L) + tail(f, -1L))) / 2
tot <- trap(dens)
if (tot < 1e-12) return(sum(object$Wh * object$Vh))
mu <- trap(x_seq * dens) / tot
mu2 <- trap(x_seq ^ 2L * dens) / tot
max(mu2 - mu^2, 0)
}
compare_designs <- function(object, n = NULL) {
if (!inherits(object, "strata"))
stop("'object' must be of class \"strata\".")
n_total <- if (!is.null(n)) as.integer(n) else object$nhTot
H <- nrow(object$h)
Wh <- object$Wh
Vh <- object$Vh
WSh <- object$WhShTot
V_within <- sum(Wh * Vh)
S2 <- .S2_from_strata_app(object)
V_srs <- S2 / n_total
V_prop <- V_within / n_total
V_opt <- WSh^2 / n_total
deff_prop <- if (V_srs > 0) V_prop / V_srs else NA_real_
deff_opt <- if (V_srs > 0) V_opt / V_srs else NA_real_
n_srs_equiv <- ceiling(S2 * n_total / WSh^2)
pct_saving <- if (n_srs_equiv > 0)
round(100 * (n_srs_equiv - n_total) / n_srs_equiv, 1) else NA_real_
list(
n = n_total,
H = H,
S2 = S2,
V_within = V_within,
V_srs = V_srs,
V_prop = V_prop,
V_opt = V_opt,
SE_srs = sqrt(max(V_srs, 0)),
SE_prop = sqrt(max(V_prop, 0)),
SE_opt = sqrt(max(V_opt, 0)),
deff_prop = deff_prop,
deff_opt = deff_opt,
n_srs_equiv = n_srs_equiv,
pct_saving = pct_saving,
gain_over_prop = if (V_opt > 0) V_prop / V_opt else NA_real_,
gain_over_srs = if (V_opt > 0) V_srs / V_opt else NA_real_,
WhShTot = WSh
)
}
# ==============================================================================
# Boundary Methods Comparison — helpers
# ==============================================================================
# ΣWₕSₕ from raw sorted data + full boundary vector (length H+1)
.whsh_from_data <- function(data, full_b) {
data <- sort(data)
N <- length(data)
H <- length(full_b) - 1L
total <- 0
for (h in seq_len(H)) {
lo <- full_b[h]; hi <- full_b[h + 1L]
sl <- if (h == 1L) data[data >= lo & data <= hi]
else data[data > lo & data <= hi]
Nh <- length(sl)
if (Nh < 2L) next
Wh <- Nh / N
Sh <- sqrt(sum((sl - mean(sl))^2) / Nh) # population SD
total <- total + Wh * Sh
}
total
}
# ΣWₕSₕ from density grid + full boundary vector (length H+1)
.whsh_from_dens <- function(x_seq, dens, full_b) {
H <- length(full_b) - 1L
dx <- diff(x_seq)
tot <- sum(dx * (head(dens, -1) + tail(dens, -1))) / 2
if (tot < 1e-12) return(NA_real_)
total <- 0
for (h in seq_len(H)) {
lo <- full_b[h]; hi <- full_b[h + 1L]
idx <- if (h < H) x_seq >= lo & x_seq < hi
else x_seq >= lo & x_seq <= hi
if (sum(idx) < 2L) next
xs <- x_seq[idx]; fs <- dens[idx]; dxh <- diff(xs)
Wh <- sum(dxh * (head(fs, -1) + tail(fs, -1))) / 2 / tot
if (Wh < 1e-12) next
mu <- sum(dxh * (head(xs * fs, -1) + tail(xs * fs, -1))) / 2 / (Wh * tot)
mu2 <- sum(dxh * (head(xs^2*fs, -1) + tail(xs^2*fs, -1))) / 2 / (Wh * tot)
Sh <- sqrt(max(mu2 - mu^2, 0))
total <- total + Wh * Sh
}
total
}
# Equal-width boundaries
.bmc_ew <- function(x_lo, x_hi, H) {
seq(x_lo, x_hi, length.out = H + 1L)
}
# Geometric (Gunning & Horgan 2004): b_i = x_lo × (x_hi/x_lo)^(i/H)
# Shifts data if x_lo ≤ 0 (log-equal-width in original space)
.bmc_geom <- function(x_lo, x_hi, H) {
shift <- 0
if (x_lo <= 0) { shift <- abs(x_lo) + 1; x_lo <- x_lo + shift; x_hi <- x_hi + shift }
if (x_lo <= 0 || x_hi <= x_lo) return(NULL)
(x_lo * (x_hi / x_lo)^(seq(0L, H) / H)) - shift
}
# Dalenius-Hodges cumulative √f: divide ∫√f(y)dy into H equal parts
.bmc_dh <- function(x_seq, dens, H) {
dx <- c(0, diff(x_seq))
sqrtf <- sqrt(pmax(dens, 0))
cum_sf <- cumsum(sqrtf * dx)
tot <- tail(cum_sf, 1)
if (tot < 1e-12) return(seq(x_seq[1], tail(x_seq,1), length.out = H + 1L))
targets <- seq(0, tot, length.out = H + 1L)
c(x_seq[1],
vapply(targets[seq(2L, H)], function(t) x_seq[which.min(abs(cum_sf - t))], numeric(1)),
tail(x_seq, 1))
}
# Equal-frequency (quantile) boundaries from raw data
.bmc_eq_data <- function(data, H) {
as.numeric(quantile(data, probs = seq(0, 1, 1/H), type = 2))
}
# Equal-frequency boundaries from density CDF (for distribution pathway)
.bmc_eq_dens <- function(x_seq, dens, H) {
dx <- c(diff(x_seq), 0)
cdf <- cumsum(dens * dx)
cdf <- cdf / max(cdf, 1e-12)
targets <- seq(0, 1, length.out = H + 1L)
c(x_seq[1],
vapply(targets[seq(2L, H)], function(t) x_seq[which.min(abs(cdf - t))], numeric(1)),
tail(x_seq, 1))
}
# ==============================================================================
# Server
# ==============================================================================
# ── LaTeX PDF equations for each stratifyR distribution ───────────────────────
STRAT_PDF_EQUATIONS <- list(
norm = "\\[ f(y) = \\frac{1}{\\sigma\\sqrt{2\\pi}}\\exp\\!\\left(-\\frac{(y-\\mu)^2}{2\\sigma^2}\\right), \\quad y\\in\\mathbb{R} \\]",
lnorm = "\\[ f(y) = \\frac{1}{y\\,\\sigma\\sqrt{2\\pi}}\\exp\\!\\left(-\\frac{(\\ln y-\\mu)^2}{2\\sigma^2}\\right), \\quad y>0 \\]",
gamma = "\\[ f(y) = \\frac{\\beta^{\\alpha}}{\\Gamma(\\alpha)}\\,y^{\\alpha-1}e^{-\\beta y}, \\quad y>0 \\]",
weibull = "\\[ f(y) = \\frac{k}{\\lambda}\\!\\left(\\frac{y}{\\lambda}\\right)^{\\!k-1}\\!\\exp\\!\\left(-\\!\\left(\\frac{y}{\\lambda}\\right)^{\\!k}\\right), \\quad y>0 \\]",
exp = "\\[ f(y) = \\lambda\\,e^{-\\lambda y}, \\quad y>0 \\]",
cauchy = "\\[ f(y) = \\frac{1}{\\pi\\gamma\\!\\left[1+\\left(\\dfrac{y-\\mu_0}{\\gamma}\\right)^{\\!2}\\right]}, \\quad y\\in\\mathbb{R} \\]",
unif = "\\[ f(y) = \\frac{1}{b-a}, \\quad a \\leq y \\leq b \\]",
pareto = "\\[ f(y) = \\frac{\\alpha\\,y_m^{\\alpha}}{y^{\\alpha+1}}, \\quad y>y_m \\]"
)
server <- function(input, output, session) {
# ── Palette helper ───────────────────────────────────────────────────────────
current_palette <- reactive({ palettes[[input$palette]] })
# ── Plot title helper ────────────────────────────────────────────────────────
plot_main <- reactive({
t <- trimws(input$plot_title)
if (nzchar(t)) t else NULL
})
# ── Uploaded file ────────────────────────────────────────────────────────────
uploaded_df <- reactive({
req(input$file_upload)
path <- input$file_upload$datapath
ext <- tolower(tools::file_ext(input$file_upload$name))
tryCatch({
if (ext %in% c("xlsx", "xls")) {
if (has_readxl) readxl::read_excel(path)
else stop("readxl not installed. Please install it to read Excel files.")
} else {
utils::read.csv(path, stringsAsFactors = FALSE)
}
}, error = function(e) { showNotification(e$message, type = "error"); NULL })
})
output$ui_col_select <- renderUI({
df <- uploaded_df(); req(df)
num_cols <- names(df)[vapply(df, is.numeric, logical(1))]
if (length(num_cols) == 0)
return(div(class="data-preview", "No numeric columns found."))
selectInput("data_col", "Variable column", choices = num_cols)
})
# ── Built-in dataset info ────────────────────────────────────────────────────
output$ui_builtin_info <- renderUI({
req(input$builtin_ds)
parts <- strsplit(input$builtin_ds, "\\|")[[1]]
ds_name <- parts[1]; col_name <- parts[2]
e <- new.env(parent = emptyenv())
tryCatch({
utils::data(list = ds_name, package = "stratifyR", envir = e)
y <- as.numeric(get(ds_name, envir = e)[[col_name]])
y <- y[is.finite(y)]
div(class = "data-preview",
tags$strong(sprintf("n = %d", length(y))),
sprintf(" | mean = %.2f | sd = %.2f | range = [%.2f, %.2f]",
mean(y), stats::sd(y), min(y), max(y))
)
}, error = function(e) NULL)
})
# ── Distribution parameter inputs ────────────────────────────────────────────
output$ui_distr_params <- renderUI({
req(input$distr_name)
info <- distr_params[[input$distr_name]]
pn <- info$names
pd <- info$defaults
np <- length(pn)
make_input <- function(i)
numericInput(paste0("dp_", pn[i]), pn[i], value = pd[i], step = 0.01)
if (np == 1) {
make_input(1)
} else if (np == 2) {
fluidRow(column(6, make_input(1)), column(6, make_input(2)))
} else {
tagList(
fluidRow(column(6, make_input(1)), column(6, make_input(2))),
make_input(3)
)
}
})
# ── Extract y vector from current source ─────────────────────────────────────
get_y <- reactive({
if (input$data_source == "upload") {
req(uploaded_df(), input$data_col)
y <- as.numeric(uploaded_df()[[input$data_col]])
} else {
parts <- strsplit(input$builtin_ds, "\\|")[[1]]
ds_name <- parts[1]; col_name <- parts[2]
e <- new.env(parent = emptyenv())
utils::data(list = ds_name, package = "stratifyR", envir = e)
y <- as.numeric(get(ds_name, envir = e)[[col_name]])
}
y[is.finite(y)]
})
# ── Core stratification ───────────────────────────────────────────────────────
strata_res <- reactiveVal(NULL)
run_err <- reactiveVal(NULL)
is_running <- reactiveVal(FALSE)
observeEvent(input$run_btn, {
strata_res(NULL); run_err(NULL); is_running(TRUE)
method_used <- input$solver
is_dp <- identical(method_used, "dp")
prog_msg <- switch(method_used,
"dp" = "▶ DP solver running…",
"cobyla" = "▶ COBYLA optimising…",
"global" = "▶ GLOBAL (DIRECT-L + COBYLA) running…")
shiny::withProgress(message = prog_msg, value = 0, {
tryCatch({
h <- as.integer(input$h_strata)
n_samp <- as.integer(input$n_sample)
method <- input$solver
shiny::incProgress(0.10, detail = "Validating inputs…")
if (input$data_source == "distr") {
# ── Distribution pathway ────────────────────────────────────────────
pn <- distr_params[[input$distr_name]]$names
params <- stats::setNames(
vapply(pn, function(p) {
v <- input[[paste0("dp_", p)]]
if (is.null(v)) 0 else as.numeric(v)
}, numeric(1)),
pn)
shiny::incProgress(0.20,
detail = switch(method_used,
"dp" = "Running DP algorithm (O(h·N²))…",
"cobyla" = "Running COBYLA (20 random starts)…",
"global" = "Running DIRECT-L global search + COBYLA refinement…"))
res <- strata.distr(
h = h,
initval = input$initval,
dist = input$dist_range,
distr = input$distr_name,
params = params,
n = n_samp,
N = as.integer(input$distr_N),
method = method
)
} else {
# ── Data pathway ────────────────────────────────────────────────────
y <- get_y()
if (n_samp > length(y))
stop(sprintf("Sample size n (%d) cannot exceed population size N (%d).", n_samp, length(y)))
shiny::incProgress(0.20,
detail = switch(method_used,
"dp" = "Running DP algorithm (O(h·N²))…",
"cobyla" = "Running COBYLA (20 random starts)…",
"global" = "Running DIRECT-L global search + COBYLA refinement…"))
res <- strata.data(
data = y,
h = h,
n = n_samp,
method = method
)
}
shiny::incProgress(0.60, detail = "Computing sample allocations…")
strata_res(res)
shiny::incProgress(0.10, detail = "Done!")
}, error = function(e) {
run_err(conditionMessage(e))
}, finally = {
is_running(FALSE)
})
})
})
# ── Status badge ──────────────────────────────────────────────────────────────
output$ui_status <- renderUI({
if (isTRUE(is_running())) {
tagList(
div(class = "status-pill s-running",
icon("spinner", class = "fa-spin"), "Computing stratification…"),
div(class = "progress-wrap",
div(class = "progress-shimmer")),
div(class = "progress-label",
switch(input$solver,
"dp" = "DP solver: O(h·N²) — larger strata counts take longer",
"cobyla" = "COBYLA: gradient-free boundary search (20 starts)…",
"global" = "GLOBAL: DIRECT-L phase 1 (500 evals) + COBYLA refinement…"))
)
} else if (!is.null(run_err())) {
div(class = "status-pill s-error", icon("circle-xmark"), run_err())
} else if (!is.null(strata_res())) {
div(class = "status-pill s-ready", icon("circle-check"), "Stratification complete")
} else {
div(class = "status-pill s-idle", icon("circle-info"), "Configure and press Run")
}
})
# ── Helper: stratum bounds table ─────────────────────────────────────────────
bounds_df <- reactive({
res <- strata_res(); req(res)
H <- nrow(res$h)
x_lo <- res$maxval * res$initval
x_hi <- max(res$maxval * res$finval, max(res$OSB))
int_osb <- res$OSB[res$OSB > x_lo & res$OSB < x_hi]
fb <- c(x_lo, int_osb, x_hi)
data.frame(
Stratum = paste0("S", seq_len(H)),
Lower = round(fb[seq_len(H)], 4),
Upper = round(fb[seq_len(H) + 1L], 4),
`Wh` = round(res$Wh, 4),
`Sh` = round(sqrt(pmax(res$Vh, 0)), 4),
`nh` = res$nh,
`Nh` = if (!is.null(res$Nh)) res$Nh else NA_integer_,
`fh` = round(res$fh, 4),
check.names = FALSE, stringsAsFactors = FALSE
)
})
# ── Introduction page content (static) ───────────────────────────────────────
intro_page <- div(
style = "max-width:900px; margin:0 auto; padding:0.4rem 0.2rem 1.5rem;",
# Hero
div(class = "intro-hero",
div(class = "intro-hero-title",
"Why Construct Optimal Strata?"),
div(class = "intro-hero-sub",
"Stratified random sampling divides a heterogeneous population into",
" homogeneous sub-groups (strata) before sampling.",
" Choosing boundaries optimally — rather than arbitrarily —",
" can reduce sampling variance by a factor of 2–5× for the same",
" total sample size, directly lowering survey costs.")
),
# ── Section 1: The Problem ────────────────────────────────────────────────
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "1"), "The Problem with Simple Random Sampling"),
div(class = "intro-body",
tags$p("In a simple random sample (SRS) of size ",
tags$em("n"), " from a population with variance ",
tags$em("S²"), ", the variance of the sample mean is:"),
div(class = "intro-math",
withMathJax("\\( V(\\bar{y}_{SRS}) = \\dfrac{S^2}{n} \\)")),
tags$p("When the population is ",
tags$strong("heterogeneous"), " — for example, a skewed income",
" distribution where most values are low but a few are very large —",
" a small sample may miss the rare high-value units entirely,",
" producing a biased or imprecise estimate.",
" Cochran (1977, Ch. 5) showed that stratification systematically",
" exploits known population structure to beat this limit.")
)
),
# ── Section 2: Stratified Sampling ───────────────────────────────────────
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "2"),
"Cochran's Framework — Stratified Random Sampling"),
div(class = "intro-body",
tags$p("Divide the population of size ", tags$em("N"),
" into ", tags$em("H"), " strata with sizes ",
tags$em("N₁, …, N_H"), ".",
" Draw independent samples of sizes ",
tags$em("n₁, …, n_H"), " from each stratum",
" (with ∑nₕ = n).",
" The stratified estimator of the population mean is:"),
div(class = "intro-math",
withMathJax(
"\\( \\bar{y}_{st} = \\sum_{h=1}^{H} W_h\\,\\bar{y}_h \\)")),
tags$p("where ", tags$em("Wₕ = Nₕ/N"),
" is the stratum weight and ",
tags$em("ȳₕ"), " is the stratum sample mean.",
" Its variance under ", tags$strong("Neyman (1934) optimal allocation"),
" (", tags$em("nₕ ∝ Wₕ Sₕ"), ") achieves:"),
div(class = "intro-math",
withMathJax(
"\\( V^*(\\bar{y}_{st}) = \\frac{\\left(\\sum_{h=1}^{H} W_h S_h\\right)^2}{n} \\)")),
div(class = "intro-key-result",
"Key insight: minimising V*(ȳ_st) for fixed n is equivalent to",
" minimising ΣWₕSₕ by choosing the strata boundaries. This is",
" precisely what stratifyR optimises.")
)
),
# ── Section 3: Optimal Boundaries ────────────────────────────────────────
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "3"),
"Constructing Optimal Strata Boundaries"),
div(class = "intro-body",
tags$p("Given the stratum weight–SD product ", tags$em("WₕSₕ"),
", the boundary construction problem is:"),
div(class = "intro-math",
withMathJax(
"\\[ \\min_{b_1,\\ldots,b_{H-1}} \\sum_{h=1}^{H} W_h(b)\\,S_h(b) \\]")),
tags$p(tags$strong("Dalenius (1950)"),
" proved the exact optimality conditions.",
" ", tags$strong("Dalenius & Hodges (1959)"),
" proposed the popular '√f rule' approximation.",
" ", tags$strong("Khan et al. (2008)"),
" reformulated the problem as a Dynamic Programming",
" problem that guarantees the",
tags$em(" globally optimal"), " strata boundaries.",
" stratifyR 2.0 implements all three approaches.",
" For continuous data the boundaries are found by integrating",
" over the fitted parametric density or the empirical distribution.")
)
),
# ── Section 4: Worked Example ─────────────────────────────────────────────
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "4"),
"Worked Example — Anaemia Survey (Iron, g/dL)"),
div(class = "intro-body",
tags$p("Population: N = 2,703 patients. Survey variable: serum iron.",
" Fitted distribution: Log-Normal. Sample size: n = 300.",
" Below are the H = 3 optimal strata from stratifyR:"),
tags$table(class = "intro-example-table",
tags$thead(tags$tr(
tags$th("Stratum"), tags$th("Interval (g/dL)"),
tags$th("Wₕ"), tags$th("Sₕ"), tags$th("WₕSₕ"), tags$th("nₕ (Neyman)")
)),
tags$tbody(
tags$tr(
tags$td("S1"), tags$td("[5.0, 12.3)"),
tags$td("0.342"), tags$td("2.14"), tags$td("0.732"), tags$td("67")),
tags$tr(
tags$td("S2"), tags$td("[12.3, 22.6)"),
tags$td("0.448"), tags$td("3.02"), tags$td("1.353"), tags$td("124")),
tags$tr(
tags$td("S3"), tags$td("[22.6, 68.0]"),
tags$td("0.210"), tags$td("7.86"), tags$td("1.651"), tags$td("109")),
tags$tr(
tags$td("Total"), tags$td("—"),
tags$td("1.000"), tags$td("—"),
tags$td("3.736"), tags$td("300"))
)
),
div(class = "intro-gain-box",
div(class = "intro-gain-card",
div(class = "gval", "0.0465"),
div(class = "glbl", "V*(ȳ) — stratified")),
div(class = "intro-gain-card",
div(class = "gval", "0.2180"),
div(class = "glbl", "V(ȳ) — SRS equivalent")),
div(class = "intro-gain-card highlight",
div(class = "gval", "4.7×"),
div(class = "glbl", "Efficiency gain over SRS"))
),
tags$p(style = "margin-top:0.8rem; font-size:0.83rem; color:#64748B;",
"Equivalently: to achieve the same precision as 300 stratified observations,",
" you would need ~1,410 observations under SRS — a 79% cost saving.")
)
),
# ── Section 5: How to Use ─────────────────────────────────────────────────
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "5"), "How to Use This App"),
div(class = "intro-body",
div(class = "intro-workflow",
div(class = "intro-step",
tags$span(class = "step-num", "① Data"),
"Choose a built-in dataset, upload your own CSV/Excel file,",
" or specify a parametric distribution directly."),
div(class = "intro-step",
tags$span(class = "step-num", "② Settings"),
"Set the number of strata H, sample size n, and solver.",
" DP gives the exact global optimum.",
" COBYLA is a fast gradient-free alternative.",
" GLOBAL (DIRECT-L + COBYLA) is the best choice for right-skewed",
" or heavy-tailed distributions when using a parametric pathway."),
div(class = "intro-step",
tags$span(class = "step-num", "③ Run"),
"Click ", tags$strong("Run Stratification"),
". The app finds optimal strata boundaries (OSB) and computes Neyman allocation."),
div(class = "intro-step",
tags$span(class = "step-num", "④ Explore"),
"Inspect results across eight analytical tabs:",
tags$ul(style = "margin:0.4rem 0 0 0.5rem; padding-left:1.2rem;",
tags$li(tags$strong("Optimal Boundaries"), " — summary table, strata statistics, and 2D density plot."),
tags$li(tags$strong("Design Comparison"), " — benchmark stratified vs SRS and proportional allocation."),
tags$li(tags$strong("H Selector"), " — plot ΣWₕSₕ and V*(n) against H to choose the number of strata."),
tags$li(tags$strong("Boundary Methods"), " — compare stratifyR against geometric, cumulative √f, and Lavallée-Hidiroglou boundary construction methods."),
tags$li(tags$strong("Sample Size"), " — back-calculate the n required to hit a target SE or CV under Neyman, proportional, and SRS."),
tags$li(tags$strong("Allocation"), " — compare Neyman, proportional, and equal allocation side-by-side for fixed boundaries."),
tags$li(tags$strong("Cost-Constrained"), " — find the variance-minimising allocation subject to a fieldwork budget and per-stratum costs."),
tags$li(tags$strong("Visualise"), " — 2D density, 3D surface, and interactive boundary explorer.")
)),
div(class = "intro-step",
tags$span(class = "step-num", "⑤ Export"),
"Download results as CSV or click",
tags$strong(" Generate R Code"),
" for a ready-to-run reproducible script.")
)
)
),
# References
div(class = "intro-ref",
tags$strong("References:"),
" Cochran, W.G. (1977). ",
tags$em("Sampling Techniques"), ", 3rd ed. Wiley. —",
" Neyman, J. (1934). On the two different aspects of the representative method.",
" ", tags$em("JRSS"), " 97, 558–625. —",
" Dalenius, T. & Hodges, J.L. (1959). Minimum variance stratification.",
" ", tags$em("JASA"), " 54, 88–101. —",
" Khan, M.G.M. et al. (2008). Determining optimum strata boundaries using",
" mathematical programming. ", tags$em("Survey Methodology"), " 34(2), 91–102. —",
" Gunning, P. & Horgan, J.M. (2004). A new algorithm for the construction of",
" stratum boundaries in skewed populations.",
" ", tags$em("Survey Methodology"), " 30(2), 159–166. —",
" Lavallée, P. & Hidiroglou, M.A. (1988). On the stratification of skewed",
" populations. ", tags$em("Survey Methodology"), " 14(1), 33–43. —",
" Baillargeon, S. & Rivest, L.-P. (2011). The construction of stratified designs",
" in R with the package stratification.",
" ", tags$em("Survey Methodology"), " 37(2), 109–129. —",
" Hansen, M.H., Hurwitz, W.N. & Madow, W.G. (1953). ",
tags$em("Sample Survey Methods and Theory"), ". Wiley."
)
)
# ── Main UI ───────────────────────────────────────────────────────────────────
output$ui_main <- renderUI({
res <- strata_res()
no_result_tab <- function(label)
div(class = "empty-state",
div(class = "empty-icon", "\U0001F4CA"),
div(class = "empty-title", "No results yet"),
div(class = "empty-sub",
"Run Stratification first — results will appear here."))
tagList(
# Metric strip (only when results exist)
if (!is.null(res)) {
gap_lbl <- if (!is.null(res$optimality_gap) && !is.na(res$optimality_gap))
formatC(res$optimality_gap, format = "e", digits = 2) else "N/A (DP)"
div(class = "metrics-row",
div(class = "metric-card",
div(class = "metric-value", nrow(res$h)),
div(class = "metric-label", "Strata")),
div(class = "metric-card",
div(class = "metric-value", res$nhTot),
div(class = "metric-label", "Total sample")),
div(class = "metric-card",
div(class = "metric-value", round(res$WhShTot, 4)),
div(class = "metric-label", "ΣWₕSₕ")),
div(class = "metric-card",
div(class = "metric-value", gap_lbl),
div(class = "metric-label", "Optimality gap"))
)
},
navset_card_tab(
# ── 1. Optimal Boundaries ────────────────────────────────────────────
nav_panel(tagList(icon("table-list"), " Optimal Boundaries"),
if (is.null(res)) no_result_tab()
else tagList(uiOutput("ui_summary"), uiOutput("ui_2d_info"))),
# ── 2. Design Comparison ─────────────────────────────────────────────
nav_panel(tagList(icon("scale-balanced"), " Design Comparison"),
if (is.null(res)) no_result_tab()
else uiOutput("ui_design_comp")),
# ── 3. H Selector ────────────────────────────────────────────────────
nav_panel(tagList(icon("layer-group"), " Strata Count (H)"),
uiOutput("ui_h_selector")),
# ── 4. Boundary Methods Comparison ───────────────────────────────────
nav_panel(tagList(icon("arrows-left-right"), " Boundary Methods"),
uiOutput("ui_bmc")),
# ── 5. Sample Size Calculator ─────────────────────────────────────────
nav_panel(tagList(icon("calculator"), " Sample Size"),
uiOutput("ui_ssc")),
# ── 6. Allocation Explorer ────────────────────────────────────────────
nav_panel(tagList(icon("chart-bar"), " Allocation"),
uiOutput("ui_alloc")),
# ── 7. Cost-Constrained ───────────────────────────────────────────────
nav_panel(tagList(icon("coins"), " Cost-Constrained"),
uiOutput("ui_cost")),
# ── 8. Visualise (dropdown) ──────────────────────────────────────────
nav_menu(
title = tagList(icon("chart-line"), " Visualise"),
align = "left",
nav_panel(tagList(icon("chart-area"), " 2D Density Plot"),
if (is.null(res)) no_result_tab()
else div(class = "plot-container",
plotlyOutput("plot_2d", height = "500px"))),
nav_panel(tagList(icon("cube"), " 3D Surface Plot"),
if (is.null(res)) no_result_tab()
else if (has_plotly)
div(class = "plot-container",
plotlyOutput("plot_3d", height = "720px"))
else div(class = "empty-state",
div("Install plotly:", tags$code("install.packages('plotly')")))),
nav_panel(tagList(icon("sliders"), " Boundary Explorer"),
uiOutput("ui_ib_v2"))
),
# ── 4. Notes ─────────────────────────────────────────────────────────
nav_panel(tagList(icon("circle-info"), " Notes"),
intro_page)
)
)
})
# ── H-Selector: reactive storage ─────────────────────────────────────────────
h_curve_res <- reactiveVal(NULL)
observeEvent(input$run_h_curve, {
res0 <- strata_res(); req(res0)
h_curve_res(NULL) # clear previous results immediately
n_val <- as.integer(input$n_sample)
max_H <- 12L
is_data <- !is.null(res0$data_internal) && length(res0$data_internal) >= 2L
# ── For the data pathway, extract the already-fitted distribution from res0
# so we can use strata.distr() in the loop — this avoids re-running the
# 9-distribution AIC selection on every H iteration (the main cause of
# slowness for heavy-tailed data like Cauchy).
use_cached_distr <- FALSE
if (is_data) {
cached_distr <- tryCatch(res0$distr, error = function(e) NULL)
cached_params <- tryCatch({
fit <- res0$fit
if (!is.null(fit$estimate)) fit$estimate
else if (!is.null(fit$fit$estimate)) fit$fit$estimate
else NULL
}, error = function(e) NULL)
cached_initval <- tryCatch(res0$initval, error = function(e) NULL)
cached_dist <- tryCatch(res0$dist, error = function(e) NULL)
cached_maxval <- tryCatch(res0$maxval, error = function(e) NULL)
cached_N <- tryCatch(res0$N, error = function(e) NULL)
use_cached_distr <- !is.null(cached_distr) && !is.null(cached_params) &&
!is.null(cached_initval) && !is.null(cached_dist)
}
withProgress(message = "H-Selector", value = 0, {
rows <- lapply(seq_len(max_H), function(h_i) {
incProgress(1 / max_H,
detail = sprintf("Computing H = %d of %d…", h_i, max_H))
tryCatch({
if (use_cached_distr) {
# Fast path: reuse the already-fitted distribution from the main run.
# strata.distr() skips the 9-distribution AIC fitting step entirely.
N_use <- if (!is.null(cached_N) && cached_N > 0) cached_N else
as.integer(input$distr_N)
r <- strata.distr(
h = h_i,
initval = cached_initval,
dist = cached_dist,
distr = cached_distr,
params = cached_params,
n = n_val,
N = N_use,
method = "cobyla")
} else if (is_data) {
# Fallback: no cached fit — run strata.data() as before
pop <- res0$data_internal
N_arg <- if (!is.null(res0$N) && res0$N > 0 && res0$N != length(pop))
res0$N else NULL
r <- if (is.null(N_arg))
strata.data(data = pop, h = h_i, n = n_val, method = "cobyla")
else
strata.data(data = pop, h = h_i, n = n_val, N = N_arg, method = "cobyla")
} else {
pnames <- distr_params[[input$distr_name]]$names
defs <- distr_params[[input$distr_name]]$defaults
params_v <- setNames(
sapply(seq_along(pnames), function(j) {
v <- input[[paste0("param_", pnames[j])]]
if (is.null(v)) defs[j] else as.numeric(v)
}),
pnames)
r <- strata.distr(
h = h_i,
initval = as.numeric(input$initval),
dist = as.numeric(input$dist_range),
distr = input$distr_name,
params = params_v,
n = n_val,
N = as.integer(input$distr_N),
method = "cobyla")
}
list(H = h_i, WhSh = r$WhShTot, V = r$WhShTot^2 / n_val,
converged = isTRUE(r$converged))
}, error = function(e) NULL)
})
})
h_curve_res(Filter(Negate(is.null), rows))
})
# ── Boundary Methods Comparison: reactive + observer ─────────────────────────
bmc_res <- reactiveVal(NULL)
observeEvent(input$run_bmc, {
res0 <- strata_res(); req(res0)
bmc_res(NULL)
# ── Lazy-install stratification (first time only) ──────────────────────
if (!requireNamespace("stratification", quietly = TRUE)) {
withProgress(message = "Installing stratification package…", value = 0.5, {
install.packages("stratification")
})
}
H <- nrow(res0$h)
n_val <- res0$nhTot
is_data <- !is.null(res0$data_internal) && length(res0$data_internal) >= 2L
dc <- .compute_density_2d(res0)
x_seq <- dc$x_seq; dens <- dc$dens
x_lo <- dc$x_lo; x_hi <- dc$x_hi
results <- withProgress(message = "Boundary Methods Comparison", value = 0, {
# ── 1. DP (stratifyR) — recompute with same formula as other methods ──
incProgress(0.1, detail = "DP (stratifyR)…")
int_osb <- res0$OSB[res0$OSB > x_lo & res0$OSB < x_hi]
dp_full_b <- c(x_lo, int_osb, x_hi)
# Recompute using the same .whsh_from_data / .whsh_from_dens formula so
# all methods are on an identical numerical footing.
dp_ws <- tryCatch(
if (is_data) .whsh_from_data(res0$data_internal, dp_full_b)
else .whsh_from_dens(x_seq, dens, dp_full_b),
error = function(e) res0$WhShTot) # fallback to stored value
solver_lbl <- toupper(if (!is.null(res0$method) && nzchar(res0$method)) res0$method else "DP")
dp_entry <- list(
method = sprintf("%s — stratifyR", solver_lbl),
WhSh = dp_ws,
V = dp_ws^2 / n_val,
is_optimal = TRUE,
converged = TRUE,
note = sprintf("stratifyR optimal boundaries (%s solver)", solver_lbl)
)
# ── 2. Geometric (Gunning & Horgan 2004) ─────────────────────────────
incProgress(0.2, detail = "Geometric (Gunning 2004)…")
geom_b <- tryCatch(.bmc_geom(x_lo, x_hi, H), error = function(e) NULL)
geom_ok <- !is.null(geom_b) && all(is.finite(geom_b)) && all(diff(geom_b) > 0)
geom_ws <- if (geom_ok) tryCatch(
if (is_data) .whsh_from_data(res0$data_internal, geom_b)
else .whsh_from_dens(x_seq, dens, geom_b),
error = function(e) NA_real_) else NA_real_
geom_entry <- list(
method = "Geometric (Gunning 2004)",
WhSh = geom_ws,
V = if (is.finite(geom_ws)) geom_ws^2 / n_val else NA_real_,
is_optimal = FALSE,
converged = geom_ok && is.finite(geom_ws),
note = "bᵢ = x_lo × (x_hi/x_lo)^(i/H) — log-equal-width spacing"
)
# ── 3. Cumulative √f (Dalenius & Hodges 1959) ─────────────────────────
incProgress(0.2, detail = "Cumulative √f (Dalenius-Hodges)…")
dh_b <- tryCatch(.bmc_dh(x_seq, dens, H), error = function(e) NULL)
dh_ok <- !is.null(dh_b) && all(is.finite(dh_b)) && all(diff(dh_b) >= 0)
dh_ws <- if (dh_ok) tryCatch(
if (is_data) .whsh_from_data(res0$data_internal, dh_b)
else .whsh_from_dens(x_seq, dens, dh_b),
error = function(e) NA_real_) else NA_real_
dh_entry <- list(
method = "Cum. √f (Dalenius-Hodges)",
WhSh = dh_ws,
V = if (is.finite(dh_ws)) dh_ws^2 / n_val else NA_real_,
is_optimal = FALSE,
converged = dh_ok && is.finite(dh_ws),
note = "Divide ∪√f(y)dy into H equal parts"
)
# ── 4. Lavallée-Hidiroglou (stratification package) ───────────────────
incProgress(0.3, detail = "Lavallée-Hidiroglou…")
lh_ws <- tryCatch({
# Get data vector: raw if available, otherwise simulate from fitted dist
y_lh <- if (is_data) {
res0$data_internal
} else {
set.seed(42L)
N_sim <- 5000L
params <- tryCatch({
fit <- res0$fit
if (!is.null(fit$estimate)) fit$estimate
else if (!is.null(fit$fit$estimate)) fit$fit$estimate
else NULL
}, error = function(e) NULL)
if (!is.null(params) && !is.null(res0$distr)) {
sim <- switch(res0$distr,
norm = stats::rnorm(N_sim, mean = params["mean"], sd = params["sd"]),
lnorm = stats::rlnorm(N_sim, meanlog = params["meanlog"], sdlog = params["sdlog"]),
gamma = stats::rgamma(N_sim, shape = params["shape"], rate = params["rate"]),
weibull = stats::rweibull(N_sim, shape = params["shape"], scale = params["scale"]),
exp = stats::rexp(N_sim, rate = params["rate"]),
cauchy = stats::rcauchy(N_sim, location = params["location"],scale = params["scale"]),
unif = stats::runif(N_sim, min = params["min"], max = params["max"]),
NULL)
if (!is.null(sim)) sim[sim >= x_lo & sim <= x_hi & is.finite(sim)] else NULL
} else NULL
}
if (is.null(y_lh) || length(y_lh) < H * 5L) stop("Insufficient data for L-H")
# strata.LH requires strictly positive data — shift if needed
lh_shift <- 0
if (min(y_lh, na.rm = TRUE) <= 0) {
lh_shift <- abs(min(y_lh, na.rm = TRUE)) + 1
y_lh <- y_lh + lh_shift
}
# strata.LH(x, Ls = H, n = ...) — Ls = strata count, n = sample size
lh_out <- suppressMessages(suppressWarnings(
stratification::strata.LH(x = y_lh, Ls = H, n = n_val)
))
# bh may be H-1 interior values OR H+1 including endpoints
bh_raw <- as.numeric(lh_out$bh) - lh_shift # shift back to original scale
lh_b <- if (length(bh_raw) == H - 1L) c(x_lo, bh_raw, x_hi)
else if (length(bh_raw) == H + 1L) bh_raw
else stop(sprintf("Unexpected boundary count from strata.LH: got %d, expected %d or %d",
length(bh_raw), H - 1L, H + 1L))
if (!all(is.finite(lh_b)) || any(diff(lh_b) <= 0))
stop("Non-monotone L-H boundaries")
if (is_data) .whsh_from_data(res0$data_internal, lh_b)
else .whsh_from_dens(x_seq, dens, lh_b)
}, error = function(e) {
message("L-H error: ", conditionMessage(e))
conditionMessage(e) # return the message string, not NA, so we can display it
})
lh_ok <- is.numeric(lh_ws) && length(lh_ws) == 1L && is.finite(lh_ws)
lh_err <- if (!lh_ok && is.character(lh_ws)) lh_ws else NULL
lh_entry <- list(
method = "Lavallée-Hidiroglou",
WhSh = if (lh_ok) lh_ws else NA_real_,
V = if (lh_ok) lh_ws^2 / n_val else NA_real_,
is_optimal = FALSE,
converged = lh_ok,
err_msg = lh_err,
note = "Iterative boundary-allocation algorithm (Lavallée & Hidiroglou 1988)"
)
incProgress(0.2, detail = "Done!")
list(dp = dp_entry, geom = geom_entry, dh = dh_entry, lh = lh_entry)
})
bmc_res(results)
})
# ── Design Comparison tab ─────────────────────────────────────────────────────
output$ui_design_comp <- renderUI({
res <- strata_res(); req(res)
cd <- tryCatch(
compare_designs(res),
error = function(e) e)
if (inherits(cd, "error")) {
return(div(class = "empty-state",
div(class = "empty-icon", "⚠"),
div(class = "empty-title", "Design Comparison failed"),
div(class = "empty-sub", conditionMessage(cd))))
}
fmt_v <- function(v) formatC(v, format = "e", digits = 4)
fmt_se <- function(v) formatC(v, format = "f", digits = 6)
fmt_deff <- function(v) if (is.finite(v)) sprintf("%.4f", v) else "—"
tagList(
# ── Headline banner ────────────────────────────────────────────────────
div(class = "dc-saving-banner",
div(class = "dsb-stat",
div(class = "dsb-val", sprintf("%.1f×", cd$gain_over_srs)),
div(class = "dsb-lbl", "More efficient than SRS")),
div(class = "dsb-div"),
div(class = "dsb-stat",
div(class = "dsb-val", sprintf("%.1f×", cd$gain_over_prop)),
div(class = "dsb-lbl", "More efficient than proportional")),
div(class = "dsb-div"),
div(class = "dsb-stat",
div(class = "dsb-val", format(cd$n_srs_equiv, big.mark = ",")),
div(class = "dsb-lbl", sprintf("SRS observations needed for same precision (vs n = %d)", cd$n))),
div(class = "dsb-div"),
div(class = "dsb-stat",
div(class = "dsb-val", sprintf("%.1f%%", cd$pct_saving)),
div(class = "dsb-lbl", "Sample-size saving over SRS"))
),
# ── Key metric cards ───────────────────────────────────────────────────
div(class = "dc-cards",
div(class = "dc-card",
div(class = "dc-val", formatC(cd$S2, format = "g", digits = 5)),
div(class = "dc-lbl", "S² population variance")),
div(class = "dc-card",
div(class = "dc-val", formatC(cd$V_within, format = "g", digits = 5)),
div(class = "dc-lbl", "ΣWₕSₕ² within-stratum component")),
div(class = "dc-card",
div(class = "dc-val", formatC(cd$WhShTot, format = "g", digits = 5)),
div(class = "dc-lbl", "ΣWₕSₕ Neyman objective")),
div(class = "dc-card dc-highlight",
div(class = "dc-val", sprintf("%.4f", cd$deff_opt)),
div(class = "dc-lbl", "DEFF Neyman vs SRS")),
div(class = "dc-card dc-amber",
div(class = "dc-val", sprintf("%.4f", cd$deff_prop)),
div(class = "dc-lbl", "DEFF Proportional vs SRS"))
),
# ── Comparison table ───────────────────────────────────────────────────
tags$table(class = "dc-table",
tags$thead(tags$tr(
tags$th("Design"),
tags$th("Allocation rule"),
tags$th("V(ȳ)"),
tags$th("SE(ȳ)"),
tags$th("DEFF"),
tags$th("Relative efficiency")
)),
tags$tbody(
tags$tr(
tags$td("Simple Random Sampling"),
tags$td("—"),
tags$td(fmt_v(cd$V_srs)),
tags$td(fmt_se(cd$SE_srs)),
tags$td("1.0000 (reference)"),
tags$td("1.00×")),
tags$tr(class = "dc-row-prop",
tags$td("Stratified — Proportional"),
tags$td(HTML("n<sub>h</sub> ∝ W<sub>h</sub>")),
tags$td(fmt_v(cd$V_prop)),
tags$td(fmt_se(cd$SE_prop)),
tags$td(fmt_deff(cd$deff_prop)),
tags$td(sprintf("%.2f×", cd$V_srs / cd$V_prop))),
tags$tr(class = "dc-row-opt",
tags$td("Stratified — Neyman (optimal)"),
tags$td(HTML("n<sub>h</sub> ∝ W<sub>h</sub>S<sub>h</sub>")),
tags$td(fmt_v(cd$V_opt)),
tags$td(fmt_se(cd$SE_opt)),
tags$td(fmt_deff(cd$deff_opt)),
tags$td(sprintf("%.2f×", cd$gain_over_srs)))
)
),
# ── Interpretation note ────────────────────────────────────────────────
div(style = paste0(
"background:#F8FAFC; border:1px solid #E2E8F0; border-radius:8px;",
" padding:0.75rem 1.1rem; font-size:0.82rem; color:#475569;",
" line-height:1.6;"),
tags$strong("How to read this: "),
sprintf(
paste0("With H = %d optimal strata and Neyman allocation (n = %d), ",
"the variance of the estimated mean is %.2f× smaller than under SRS ",
"and %.2f× smaller than under proportional allocation. ",
"To achieve the same precision with SRS alone you would need ",
"n = %s observations — a %.1f%% cost saving."),
cd$H, cd$n,
cd$gain_over_srs, cd$gain_over_prop,
format(cd$n_srs_equiv, big.mark = ","), cd$pct_saving),
tags$br(), tags$br(),
HTML("DEFF = V<sub>design</sub> / V<sub>SRS</sub> (Kish 1965). ",
"Values below 1 indicate stratification improves precision over SRS. ",
"★ marks the Neyman optimal design implemented by stratifyR.")
)
)
})
# ── Boundary Methods Comparison tab ──────────────────────────────────────────
output$ui_bmc <- renderUI({
res0 <- strata_res()
if (is.null(res0)) {
return(div(class = "empty-state",
div(class = "empty-icon", "\U0001F4CA"),
div(class = "empty-title", "No results yet"),
div(class = "empty-sub", "Run Stratification first.")))
}
rows <- bmc_res()
H <- nrow(res0$h)
n_val <- res0$nhTot
banner <- div(class = "bmc-banner",
div(
div(class = "bmc-banner-title", "Boundary Methods Comparison"),
div(class = "bmc-banner-text",
sprintf("Compares four boundary construction strategies for H = %d strata, n = %d.", H, n_val),
" ΣWₕSₕ and V*(n) are evaluated under Neyman allocation for each method.",
sprintf(" stratifyR (%s) is the reference — all gaps are relative to it.",
toupper(if (!is.null(res0$method) && nzchar(res0$method)) res0$method else "DP")))
),
actionButton("run_bmc",
tagList(icon(if (is.null(rows)) "play" else "rotate"),
if (is.null(rows)) " Run Comparison" else " Recompute"),
class = if (is.null(rows)) "btn btn-primary"
else "btn btn-outline-light btn-sm")
)
if (is.null(rows)) return(banner)
# ── Build data frame, sort by WhSh ascending ──────────────────────────
df <- data.frame(
Method = sapply(rows, `[[`, "method"),
WhSh = sapply(rows, `[[`, "WhSh"),
V = sapply(rows, `[[`, "V"),
is_opt = sapply(rows, `[[`, "is_optimal"),
ok = sapply(rows, `[[`, "converged"),
err_msg = sapply(rows, function(r) {
m <- r$err_msg; if (is.null(m)) "" else as.character(m) }),
stringsAsFactors = FALSE
)
dp_ws <- df$WhSh[df$is_opt]
if (length(dp_ws) != 1L || !is.finite(dp_ws)) dp_ws <- NA_real_
# gap_pct: % by which each method's ΣWₕSₕ exceeds DP (always ≥ 0 in theory)
df$gap_pct <- round(100 * (df$WhSh - dp_ws) / dp_ws, 2)
df$RE <- round(dp_ws^2 / df$WhSh^2, 4) # V_DP / V_method
df <- df[order(df$WhSh, na.last = TRUE), ]
tbl_rows <- lapply(seq_len(nrow(df)), function(i) {
is_dp <- df$is_opt[i]
ok <- df$ok[i]
cls <- if (is_dp) "bmc-row-dp" else if (!ok) "bmc-row-fail" else ""
ws_str <- if (is.finite(df$WhSh[i]))
formatC(df$WhSh[i], format = "f", digits = 5) else "—"
v_str <- if (is.finite(df$V[i]))
formatC(df$V[i], format = "e", digits = 4) else "—"
gap_str <- if (is_dp) "reference"
else if (is.finite(df$gap_pct[i])) sprintf("%+.2f%%", df$gap_pct[i])
else "—"
re_str <- if (is_dp) "1.000×"
else if (is.finite(df$RE[i])) sprintf("%.3f×", df$RE[i])
else "—"
err_msg <- df$err_msg[i]
has_err <- nzchar(err_msg)
badge <- if (is_dp) tags$span(class = "bmc-badge-opt", "★ optimal")
else if (!ok) tags$span(class = "bmc-badge-fail", "✗ failed")
else NULL
err_tip <- if (!ok && has_err)
tags$span(style = "font-size:0.72rem; color:#DC2626; display:block; margin-top:2px;",
err_msg)
else NULL
tags$tr(class = cls,
tags$td(tagList(df$Method[i], badge, err_tip)),
tags$td(ws_str),
tags$td(v_str),
tags$td(gap_str),
tags$td(re_str)
)
})
tagList(
banner,
div(class = "plot-container", style = "margin-bottom: 1rem;",
plotlyOutput("plot_bmc", height = "280px")),
tags$table(class = "bmc-table",
tags$thead(tags$tr(
tags$th("Method"),
tags$th("ΣWₕSₕ"),
tags$th("V*(n)"),
tags$th("Gap vs DP"),
tags$th("Rel. efficiency vs DP")
)),
tags$tbody(tbl_rows)
),
div(class = "hsel-note",
"All methods evaluated at the same H and n. Boundaries differ; allocation is Neyman-optimal for each set of boundaries."),
# ── Method explanations ────────────────────────────────────────────────
div(style = "margin-top: 1.4rem;",
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "1"), "Methods compared"),
div(class = "intro-body",
tags$p(tags$strong("DP / COBYLA / GLOBAL — stratifyR:"),
" Globally optimal boundaries via Dynamic Programming, COBYLA, or GLOBAL (Khan et al. 2008).",
" DP guarantees the global optimum on the empirical grid;",
" COBYLA is a fast gradient-free alternative;",
" GLOBAL (DIRECT-L + COBYLA) provides the best accuracy for skewed distributions.",
" This is the reference — all other methods are compared against it."),
tags$p(tags$strong("Geometric (Gunning & Horgan 2004):"),
" Places boundaries at geometrically spaced values:",
" bᵢ = x_lo × (x_hi / x_lo)^(i/H).",
" Equivalent to equal-width stratification in log space.",
" Near-optimal for Pareto-type and strongly right-skewed distributions."),
tags$p(tags$strong("Cumulative √f — Dalenius & Hodges (1959):"),
" Divides the cumulative square root of the density function, ∫√f(y) dy,",
" into H equal parts and places boundaries at the division points.",
" A classical and widely cited approximation that works well across many",
" distribution families."),
tags$p(tags$strong("Lavallée-Hidiroglou (1988):"),
" An iterative algorithm that alternates between optimising stratum boundaries",
" and recomputing Neyman allocation until the sample sizes stabilise.",
" Particularly effective for skewed business populations.",
" Implemented here via the ",
tags$em("stratification"), " R package (Baillargeon & Rivest 2011).",
" The package is installed on first use.")
)
),
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "2"), "How to read the results"),
div(class = "intro-body",
tags$p(
"All four methods produce H-1 interior boundary points.",
" For each set of boundaries, Neyman-optimal allocation is applied using",
" the same total sample size n, and ΣWₕSₕ is recomputed from scratch.",
" A smaller ΣWₕSₕ means lower variance for the same n."),
div(class = "intro-math",
withMathJax(
"\\[ V^*(\\bar{y}) = \\frac{(\\Sigma W_h S_h)^2}{n} \\]")),
tags$p(
tags$strong("Gap vs DP:"),
" percentage by which a method's ΣWₕSₕ exceeds the DP minimum.",
" Negative is impossible (DP is optimal); values near 0% mean the method",
" matches DP closely for this dataset."),
tags$p(
tags$strong("Rel. efficiency vs DP:"),
" V_DP / V_method. Values below 1 indicate the method is less efficient",
" than DP; a value of 0.95 means the method's variance is 1/0.95 ≈ 5%",
" larger than DP's.")
)
),
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "3"), "References"),
div(class = "intro-body",
div(class = "intro-ref",
tags$strong("DP method: "),
"Khan, M.G.M., Nand, N., & Ahmad, N. (2008). Determining optimum strata boundaries using mathematical programming. ",
tags$em("Survey Methodology"), " 34(2), 91–102.", tags$br(),
tags$strong("Geometric method: "),
"Gunning, P. & Horgan, J.M. (2004). A new algorithm for the construction of stratum boundaries in skewed populations. ",
tags$em("Survey Methodology"), " 30(2), 159–166.", tags$br(),
tags$strong("Cumulative √f: "),
"Dalenius, T. & Hodges, J.L. (1959). Minimum variance stratification. ",
tags$em("Journal of the American Statistical Association"), " 54, 88–101.", tags$br(),
tags$strong("Lavallée-Hidiroglou: "),
"Lavallée, P. & Hidiroglou, M. (1988). On the stratification of skewed populations. ",
tags$em("Survey Methodology"), " 14(1), 33–43.", tags$br(),
tags$strong("R package stratification: "),
"Baillargeon, S. & Rivest, L.-P. (2011). The construction of stratified designs in R with the package stratification. ",
tags$em("Survey Methodology"), " 37(1), 53–65."
)
)
)
)
)
})
output$plot_bmc <- renderPlotly({
req(has_plotly)
rows <- bmc_res(); req(!is.null(rows))
df <- data.frame(
Method = sapply(rows, `[[`, "method"),
WhSh = sapply(rows, `[[`, "WhSh"),
is_opt = sapply(rows, `[[`, "is_optimal"),
ok = sapply(rows, `[[`, "converged"),
stringsAsFactors = FALSE
)
df <- df[is.finite(df$WhSh), ]
if (nrow(df) == 0) return(plotly::plot_ly())
df <- df[order(df$WhSh), ]
bar_col <- ifelse(df$is_opt, "#059669",
ifelse(!df$ok, "#CBD5E1", "#4A6FA5"))
plotly::plot_ly(
x = df$WhSh,
y = factor(df$Method, levels = rev(df$Method)),
type = "bar",
orientation = "h",
marker = list(color = bar_col,
line = list(color = "white", width = 1.5)),
text = formatC(df$WhSh, format = "f", digits = 5),
textposition = "outside",
textfont = list(size = 11, color = "#334155"),
hovertemplate = "<b>%{y}</b><br>ΣWₕSₕ = %{x:.5f}<extra></extra>"
) |>
plotly::layout(
xaxis = list(
title = list(text = "ΣWₕSₕ (lower = better)", font = list(size = 12)),
zeroline = FALSE, gridcolor = "#F1F5F9",
showline = TRUE, linecolor = "#E2E8F0"),
yaxis = list(title = "", tickfont = list(size = 11)),
paper_bgcolor = "white", plot_bgcolor = "white",
margin = list(t = 20, b = 40, l = 240, r = 100),
showlegend = FALSE
)
})
# ── Sample Size Calculator ────────────────────────────────────────────────────
output$ui_ssc <- renderUI({
res <- strata_res()
if (is.null(res))
return(div(class="empty-state",
div(class="empty-icon","\U0001F4CA"),
div(class="empty-title","No results yet"),
div(class="empty-sub","Run Stratification first.")))
tagList(
div(class = "bmc-banner",
div(
div(class = "bmc-banner-title", "Precision / Sample Size Calculator"),
div(class = "bmc-banner-text",
"Enter a target precision and see how many observations are needed",
" under Neyman allocation, proportional allocation, and SRS.")
)
),
div(class = "ssc-input-panel",
div(radioButtons("ssc_metric", "Express target precision as:",
choices = c("Standard Error SE(ȳ)" = "se",
"Coefficient of Variation CV(%)" = "cv"),
selected = isolate(input$ssc_metric) %||% "se",
inline = TRUE)),
div(style = "flex:0 0 auto; min-width:160px;",
numericInput("ssc_target_val",
label = "Target value",
value = isolate(input$ssc_target_val) %||% 0.05,
min = 1e-6, step = 0.001))
),
uiOutput("ui_ssc_results"),
tags$details(class = "method-panel",
tags$summary("Methodology & References"),
div(class = "method-body",
tags$p(tags$strong("What this does:"),
" Given a target precision — expressed as a standard error SE(ȳ) or a coefficient",
" of variation CV(ȳ) = SE(ȳ)/|ȳ| × 100% — this tab back-calculates the minimum",
" sample size required under three designs."),
tags$p(tags$strong("Neyman allocation"), " is optimal: it distributes the sample",
" proportionally to stratum size and within-stratum variability (WₕSₕ), yielding",
" the smallest variance for a given n. The required sample size is:"),
div(class = "method-formula", "n_Neyman = ⌈ (ΣWₕSₕ)² / SE² ⌉"),
tags$p(tags$strong("Proportional allocation"), " (nₕ ∝ Wₕ) ignores within-stratum",
" variance; it is simpler to implement but requires more observations:"),
div(class = "method-formula", "n_prop = ⌈ ΣWₕVₕ / SE² ⌉"),
tags$p(tags$strong("Simple Random Sampling (SRS)"), " treats the population as",
" unstratified. Its variance is driven by the overall population variance S²:"),
div(class = "method-formula", "n_SRS = ⌈ S² / SE² ⌉"),
tags$p("The SE-vs-n curve shows how precision improves as n grows under each design,",
" with a dashed line marking the target n under Neyman allocation.",
" The savings percentages quantify how much smaller the Neyman sample is relative",
" to SRS and proportional allocation at the same precision."),
div(class = "method-ref",
tags$strong("References: "),
"Cochran, W.G. (1977). ",
tags$em("Sampling Techniques"), ", 3rd ed. Wiley. §5.4–5.5. — ",
"Lohr, S.L. (2010). ",
tags$em("Sampling: Design and Analysis"), ", 2nd ed. Brooks/Cole. §4.4. — ",
"Neyman, J. (1934). On the two different aspects of the representative method.",
" ", tags$em("Journal of the Royal Statistical Society"), ", 97(4), 558–625.")
)
)
)
})
output$ui_ssc_results <- renderUI({
res <- strata_res(); req(res)
metric <- input$ssc_metric; req(metric)
tval <- input$ssc_target_val; req(tval); req(is.finite(tval) && tval > 0)
WSh <- res$WhShTot
S2 <- .S2_from_strata_app(res)
V_within <- sum(res$Wh * res$Vh)
# Population mean (for CV → SE conversion)
mu <- if (!is.null(res$data_internal) && length(res$data_internal) >= 2L) {
mean(res$data_internal, na.rm = TRUE)
} else {
dc <- .compute_density_2d(res)
dx <- diff(dc$x_seq)
tot <- sum(dx * (head(dc$dens,-1) + tail(dc$dens,-1))) / 2
if (tot > 1e-12)
sum(dx * (head(dc$x_seq*dc$dens,-1) + tail(dc$x_seq*dc$dens,-1))) / 2 / tot
else NA_real_
}
# Convert CV → SE if needed
se_target <- if (metric == "se") {
tval
} else {
if (!is.finite(mu) || abs(mu) < 1e-12)
return(div(class="empty-state",
div(class="empty-icon","⚠"),
div(class="empty-title","Cannot compute SE"),
div(class="empty-sub","Population mean is near zero; please use SE mode.")))
(tval / 100) * abs(mu)
}
# Required n under each design
n_neyman <- ceiling(WSh^2 / se_target^2)
n_prop <- ceiling(V_within / se_target^2)
n_srs <- ceiling(S2 / se_target^2)
# SE at the current n for reference
n_cur <- res$nhTot
se_cur <- WSh / sqrt(n_cur)
# Savings
save_srs <- if (n_srs > 0) round(100*(n_srs - n_neyman)/n_srs, 1) else NA_real_
save_prop <- if (n_prop > 0) round(100*(n_prop - n_neyman)/n_prop, 1) else NA_real_
# Target label
tgt_lbl <- if (metric == "se") sprintf("SE = %.4g", tval)
else sprintf("CV = %.2f%%", tval)
tagList(
# ── Result boxes ─────────────────────────────────────────────────────
div(class = "ssc-results-row",
div(class = "ssc-result-box",
div(class = "ssc-result-n", format(n_neyman, big.mark=",")),
div(class = "ssc-result-lbl", sprintf("n required (Neyman) for %s", tgt_lbl))),
div(class = "dc-cards", style = "flex:1; align-items:flex-start;",
div(class = "dc-card",
div(class = "dc-val", format(n_prop, big.mark=",")),
div(class = "dc-lbl", "n — Proportional")),
div(class = "dc-card",
div(class = "dc-val", format(n_srs, big.mark=",")),
div(class = "dc-lbl", "n — SRS")),
div(class = "dc-card dc-highlight",
div(class = "dc-val", sprintf("%.1f%%", save_srs)),
div(class = "dc-lbl", "Saving vs SRS")),
div(class = "dc-card dc-amber",
div(class = "dc-val", sprintf("%.1f%%", save_prop)),
div(class = "dc-lbl", "Saving vs Proportional"))
)
),
# ── SE vs n curve ────────────────────────────────────────────────────
div(class = "plot-container", style = "margin-bottom:1rem;",
plotlyOutput("plot_ssc_curve", height = "300px")),
# ── Note about current run ────────────────────────────────────────────
div(class = "hsel-note",
sprintf("Current run: n = %d gives SE = %.5f (Neyman). Target: %s requires n = %s.",
n_cur, se_cur, tgt_lbl, format(n_neyman, big.mark=",")))
)
})
output$plot_ssc_curve <- renderPlotly({
req(has_plotly)
res <- strata_res(); req(res)
metric <- input$ssc_metric; req(metric)
tval <- input$ssc_target_val; req(tval); req(is.finite(tval) && tval > 0)
WSh <- res$WhShTot
S2 <- .S2_from_strata_app(res)
V_within <- sum(res$Wh * res$Vh)
n_cur <- res$nhTot
mu <- if (!is.null(res$data_internal) && length(res$data_internal) >= 2L)
mean(res$data_internal, na.rm = TRUE)
else NA_real_
se_target <- if (metric == "se") tval
else if (is.finite(mu) && abs(mu) > 1e-12) (tval/100)*abs(mu)
else return(plotly::plot_ly())
n_neyman <- ceiling(WSh^2 / se_target^2)
n_max <- max(n_neyman * 2L, n_cur * 2L, 50L)
n_seq <- unique(sort(c(seq(1L, n_max, length.out = 200L), n_neyman, n_cur)))
se_neyman <- WSh / sqrt(n_seq)
se_prop <- sqrt(V_within) / sqrt(n_seq)
se_srs <- sqrt(S2) / sqrt(n_seq)
y_label <- if (metric == "cv" && is.finite(mu) && abs(mu) > 1e-12)
"CV(ȳ) (%)" else "SE(ȳ)"
y_scale <- if (metric == "cv" && is.finite(mu) && abs(mu) > 1e-12)
100 / abs(mu) else 1
plotly::plot_ly() |>
plotly::add_trace(x=n_seq, y=se_srs *y_scale, name="SRS",
type="scatter", mode="lines",
line=list(color="#94A3B8", width=1.8, dash="dot")) |>
plotly::add_trace(x=n_seq, y=se_prop *y_scale, name="Proportional",
type="scatter", mode="lines",
line=list(color="#F28E2B", width=2)) |>
plotly::add_trace(x=n_seq, y=se_neyman *y_scale, name="Neyman (optimal)",
type="scatter", mode="lines",
line=list(color="#059669", width=2.5)) |>
plotly::add_trace(
x=c(n_neyman,n_neyman), y=c(0, se_target*y_scale*1.5),
type="scatter", mode="lines", showlegend=FALSE,
line=list(color="#C0392B", width=1.5, dash="dash")) |>
plotly::layout(
xaxis = list(title=list(text="Sample size n", font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9"),
yaxis = list(title=list(text=y_label, font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9"),
legend = list(orientation="h", x=0.5, xanchor="center",
y=-0.18, font=list(size=11)),
hovermode="x unified",
paper_bgcolor="white", plot_bgcolor="white",
margin=list(t=15, b=60, l=70, r=20),
annotations = list(list(
x=n_neyman, y=se_target*y_scale*1.55,
xanchor="center", yanchor="bottom",
text=sprintf("<b>n = %s</b>", format(n_neyman, big.mark=",")),
showarrow=FALSE, font=list(size=10, color="#C0392B"),
bgcolor="#FEF2F2", bordercolor="#C0392B", borderpad=3))
)
})
# ── Allocation Schemes Explorer ───────────────────────────────────────────────
output$ui_alloc <- renderUI({
res <- strata_res()
if (is.null(res))
return(div(class="empty-state",
div(class="empty-icon","\U0001F4CA"),
div(class="empty-title","No results yet"),
div(class="empty-sub","Run Stratification first.")))
H <- nrow(res$h)
n_val <- res$nhTot
Wh <- res$Wh
Vh <- res$Vh
Sh <- sqrt(pmax(Vh, 0))
WSh <- Wh * Sh
# ── Three allocations ─────────────────────────────────────────────────
# alloc_raw: falls back to equal allocation when sum(wts)=0 or NaN/NA
# (can happen when GLOBAL/COBYLA produces single-point strata with Sh=0)
alloc_raw <- function(wts) {
s <- sum(wts, na.rm = TRUE)
if (!is.finite(s) || s <= 0) {
e <- rep(floor(n_val / H), H); e[1L] <- e[1L] + n_val - sum(e)
return(pmax(as.integer(e), 1L))
}
r <- round(n_val * wts / s)
r[!is.finite(r) | r < 1L] <- 1L
r
}
neyman_undef <- {s <- sum(WSh, na.rm=TRUE); !is.finite(s) || s <= 0}
nh_neyman <- alloc_raw(WSh)
nh_prop <- alloc_raw(Wh)
nh_equal <- { e <- rep(floor(n_val / H), H); e[1] <- e[1] + n_val - sum(e); e }
# Adjust rounding so totals match n exactly
fix_total <- function(nh) {
nh[!is.finite(nh)] <- 1L # sanitise any leftover NaN/Inf
diff_n <- n_val - sum(nh, na.rm = TRUE)
if (isTRUE(diff_n > 0)) nh[order(nh)[seq_len( diff_n)]] <- nh[order(nh)[seq_len( diff_n)]] + 1L
if (isTRUE(diff_n < 0)) nh[order(-nh)[seq_len(-diff_n)]] <- nh[order(-nh)[seq_len(-diff_n)]] - 1L
pmax(nh, 1L)
}
nh_neyman <- fix_total(nh_neyman)
nh_prop <- fix_total(nh_prop)
# V = Σ Wh² Sh² / nh (ignoring FPC)
Wh2Sh2 <- Wh^2 * Vh
V_neyman <- sum(Wh2Sh2 / nh_neyman)
V_prop <- sum(Wh2Sh2 / nh_prop)
V_equal <- sum(Wh2Sh2 / nh_equal)
fmt_v <- function(v) formatC(v, format="e", digits=4)
fmt_se <- function(v) formatC(sqrt(max(v,0)), format="f", digits=5)
strata_lbl <- paste0("S", seq_len(H))
tagList(
if (neyman_undef)
div(class = "hsel-warn",
tags$strong("⚠ Neyman allocation undefined:"),
" all within-stratum variances are zero for this solution (e.g. single-point or",
" constant strata). The Neyman column shows equal allocation as a fallback."),
div(class = "bmc-banner",
div(
div(class = "bmc-banner-title", "Allocation Schemes Explorer"),
div(class = "bmc-banner-text",
sprintf("Compares Neyman, proportional, and equal allocation for H = %d strata, n = %d.", H, n_val),
" The strata boundaries are fixed at the stratifyR solution; only the within-stratum",
" sample sizes nₕ differ.")
)
),
# ── Summary cards ──────────────────────────────────────────────────
div(class = "dc-cards",
div(class = "dc-card dc-highlight",
div(class = "dc-val", fmt_v(V_neyman)),
div(class = "dc-lbl", "V*(ȳ) — Neyman")),
div(class = "dc-card",
div(class = "dc-val", fmt_v(V_prop)),
div(class = "dc-lbl", "V*(ȳ) — Proportional")),
div(class = "dc-card",
div(class = "dc-val", fmt_v(V_equal)),
div(class = "dc-lbl", "V*(ȳ) — Equal")),
div(class = "dc-card",
div(class = "dc-val", sprintf("%.3f×", V_prop / V_neyman)),
div(class = "dc-lbl", "Prop. / Neyman variance ratio")),
div(class = "dc-card",
div(class = "dc-val", sprintf("%.3f×", V_equal / V_neyman)),
div(class = "dc-lbl", "Equal / Neyman variance ratio"))
),
# ── Grouped bar chart ─────────────────────────────────────────────────
div(class = "plot-container", style = "margin-bottom:1rem;",
plotlyOutput("plot_alloc",
height = paste0(max(320L, min(600L,
ceiling(max(nh_neyman, nh_prop, nh_equal) * 1.5) + 120L)), "px"))),
# ── Table ─────────────────────────────────────────────────────────────
tags$table(class = "dc-table",
tags$thead(tags$tr(
tags$th("Stratum"),
tags$th("Wₕ"), tags$th("Sₕ"),
tags$th("nₕ Neyman"), tags$th("nₕ Proportional"), tags$th("nₕ Equal")
)),
tags$tbody(lapply(seq_len(H), function(h) {
tags$tr(
tags$td(tags$strong(strata_lbl[h])),
tags$td(round(Wh[h], 4)), tags$td(round(Sh[h], 4)),
tags$td(tags$strong(nh_neyman[h])),
tags$td(nh_prop[h]),
tags$td(nh_equal[h])
)
})),
tags$tfoot(tags$tr(style="font-weight:700; background:#F1F5F9;",
tags$td("Total"),
tags$td("1.0000"), tags$td("—"),
tags$td(sum(nh_neyman)), tags$td(sum(nh_prop)), tags$td(sum(nh_equal))
))
),
div(class = "hsel-note", style = "margin-top:0.6rem;",
sprintf("SE: Neyman = %s | Proportional = %s | Equal = %s",
fmt_se(V_neyman), fmt_se(V_prop), fmt_se(V_equal))),
tags$details(class = "method-panel",
tags$summary("Methodology & References"),
div(class = "method-body",
tags$p(tags$strong("What this does:"),
" For fixed strata boundaries (as determined by stratifyR), three allocation",
" rules are compared side-by-side. All use the same total sample size n; only",
" the within-stratum counts nₕ differ."),
tags$p(tags$strong("Neyman (optimal) allocation"), " minimises V*(ȳ) for a given n",
" by assigning more observations to strata that are larger and more variable:"),
div(class = "method-formula", "nₕ ∝ WₕSₕ → nₕ = n × WₕSₕ / ΣWₕSₕ"),
tags$p(tags$strong("Proportional allocation"), " assigns observations in proportion",
" to stratum size alone, ignoring within-stratum variance:"),
div(class = "method-formula", "nₕ ∝ Wₕ → nₕ = n × Wₕ"),
tags$p(tags$strong("Equal allocation"), " splits the sample evenly across strata",
" regardless of size or variability. Simple to administer, but often inefficient:"),
div(class = "method-formula", "nₕ = n / H (rounded)"),
tags$p("Variance under each allocation is computed as V*(ȳ) = Σ Wₕ²Vₕ/nₕ",
" (ignoring the finite population correction). The variance ratio (Proportional /",
" Neyman and Equal / Neyman) shows the cost in efficiency of the simpler designs.",
" A ratio of 1.20 means proportional allocation needs 20% more observations to",
" achieve the same precision as Neyman."),
div(class = "method-ref",
tags$strong("References: "),
"Neyman, J. (1934). On the two different aspects of the representative method.",
" ", tags$em("Journal of the Royal Statistical Society"), ", 97(4), 558–625. — ",
"Cochran, W.G. (1977). ",
tags$em("Sampling Techniques"), ", 3rd ed. Wiley. §5.5. — ",
"Lohr, S.L. (2010). ",
tags$em("Sampling: Design and Analysis"), ", 2nd ed. §4.3.")
)
)
)
})
output$plot_alloc <- renderPlotly({
req(has_plotly)
res <- strata_res(); req(res)
H <- nrow(res$h)
n_val <- res$nhTot
Wh <- res$Wh; Sh <- sqrt(pmax(res$Vh, 0)); WSh <- Wh * Sh
alloc_raw <- function(wts) {
s <- sum(wts, na.rm = TRUE)
if (!is.finite(s) || s <= 0) {
e <- rep(floor(n_val / H), H); e[1L] <- e[1L] + n_val - sum(e)
return(pmax(as.integer(e), 1L))
}
r <- round(n_val * wts / s)
r[!is.finite(r) | r < 1L] <- 1L
r
}
fix_total <- function(nh) {
nh[!is.finite(nh)] <- 1L
d <- n_val - sum(nh, na.rm = TRUE)
if (isTRUE(d > 0)) nh[order(nh)[seq_len( d)]] <- nh[order(nh)[seq_len( d)]] + 1L
if (isTRUE(d < 0)) nh[order(-nh)[seq_len(-d)]] <- nh[order(-nh)[seq_len(-d)]] - 1L
pmax(nh, 1L)
}
nh_n <- fix_total(alloc_raw(WSh))
nh_p <- fix_total(alloc_raw(Wh))
nh_e <- { e <- rep(floor(n_val/H), H); e[1] <- e[1] + n_val - sum(e); e }
s_lbl <- paste0("S", seq_len(H))
y_max <- max(nh_n, nh_p, nh_e)
y_range <- c(0, y_max * 1.28) # 28% headroom for outside labels
plotly::plot_ly(x = s_lbl, y = nh_n, type = "bar", name = "Neyman",
marker = list(color = "#059669"), textposition = "outside",
text = nh_n) |>
plotly::add_trace(x=s_lbl, y=nh_p, name="Proportional",
marker=list(color="#F28E2B"), text=nh_p, textposition="outside") |>
plotly::add_trace(x=s_lbl, y=nh_e, name="Equal",
marker=list(color="#94A3B8"), text=nh_e, textposition="outside") |>
plotly::layout(
barmode = "group",
xaxis = list(title=list(text="Stratum", font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9"),
yaxis = list(title=list(text="nₕ (sample size)", font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9",
range = y_range),
legend = list(orientation="h", x=0.5, xanchor="center",
y=-0.22, font=list(size=11)),
paper_bgcolor="white", plot_bgcolor="white",
margin=list(t=30, b=60, l=70, r=20)
)
})
# ── Cost-Constrained Optimisation ────────────────────────────────────────────
cost_res <- reactiveVal(NULL)
output$ui_cost <- renderUI({
res <- strata_res()
if (is.null(res))
return(div(class="empty-state",
div(class="empty-icon","\U0001F4CA"),
div(class="empty-title","No results yet"),
div(class="empty-sub","Run Stratification first.")))
H <- nrow(res$h)
n_val <- res$nhTot
Wh <- res$Wh; Sh <- sqrt(pmax(res$Vh, 0))
# Sensible default budget: cost-per-unit = 1 everywhere, budget = n
def_c <- 1
cost_inputs <- lapply(seq_len(H), function(h)
div(class = "cc-cost-card",
tags$label(sprintf("c%d (stratum %d)", h, h)),
numericInput(paste0("cc_cost_", h),
label = NULL,
value = isolate(input[[paste0("cc_cost_", h)]]) %||% def_c,
min = 0.001, step = 0.1)
)
)
tagList(
div(class = "bmc-banner",
div(
div(class = "bmc-banner-title", "Cost-Constrained Optimisation"),
div(class = "bmc-banner-text",
"Enter per-stratum fieldwork costs cₕ and a total budget C.",
" The cost-optimal allocation nₕ ∝ WₕSₕ/√cₕ minimises variance",
" subject to ∑cₕnₕ ≤ C.")
),
actionButton("run_cost",
tagList(icon("play"), " Compute"),
class = "btn btn-primary")
),
div(style = "background:#FFFFFF; border:1.5px solid #E2E8F0; border-radius:10px; padding:1rem 1.2rem; margin-bottom:1rem;",
div(style = "display:flex; gap:1.5rem; flex-wrap:wrap; align-items:flex-end;",
div(style = "flex:0 0 auto; min-width:160px;",
numericInput("cc_budget", "Total budget C",
value = isolate(input$cc_budget) %||% (def_c * n_val),
min = 1, step = 100)),
div(style = "flex:1;",
div(style = "font-size:0.78rem; font-weight:600; color:#1B3A5C; margin-bottom:0.4rem;",
"Per-stratum cost cₕ:"),
div(class = "cc-cost-grid", cost_inputs)
)
)
),
uiOutput("ui_cost_results"),
tags$details(class = "method-panel",
tags$summary("Methodology & References"),
div(class = "method-body",
tags$p(tags$strong("What this does:"),
" When different strata have different fieldwork costs per unit (e.g., rural strata",
" cost more to survey than urban ones), Neyman allocation is no longer optimal.",
" This tab finds the allocation that minimises variance for a given total budget C,",
" where the total cost is ΣcₕNₕ ≤ C."),
tags$p(tags$strong("Cost-optimal allocation"), " (Cochran 1977, §5.5) extends",
" Neyman by down-weighting expensive strata:"),
div(class = "method-formula", "nₕ ∝ WₕSₕ / √cₕ"),
tags$p("The total sample size under this allocation is:"),
div(class = "method-formula",
"n = C × Σ(WₕSₕ/√cₕ) / Σ(WₕSₕ√cₕ)"),
tags$p("The minimum achievable variance for budget C is given analytically by:"),
div(class = "method-formula", "V_min = (ΣWₕSₕ√cₕ)² / C"),
tags$p("For comparison, the tab also shows the ", tags$strong("Neyman allocation"),
" for the same budget — i.e., ignoring costs and distributing n ∝ WₕSₕ with",
" n chosen so that ΣcₕNₕ = C. When costs are unequal, cost-optimal allocation",
" always achieves a lower (or equal) variance than budget-constrained Neyman."),
tags$p(tags$strong("How to use: "),
"Set cₕ = 1 for all strata to recover standard Neyman allocation.",
" Set cₕ proportional to real fieldwork costs (e.g. travel time, interviewer",
" wages) to find the allocation that makes the most of your budget."),
div(class = "method-ref",
tags$strong("References: "),
"Cochran, W.G. (1977). ",
tags$em("Sampling Techniques"), ", 3rd ed. Wiley. §5.5 (Cost and variance in",
" stratified sampling). — ",
"Hansen, M.H., Hurwitz, W.N. & Madow, W.G. (1953). ",
tags$em("Sample Survey Methods and Theory"), ". Wiley. Vol. I, §6.9. — ",
"Lohr, S.L. (2010). ",
tags$em("Sampling: Design and Analysis"), ", 2nd ed. §4.5.")
)
)
)
})
observeEvent(input$run_cost, {
res <- strata_res(); req(res)
H <- nrow(res$h)
Wh <- res$Wh; Vh <- res$Vh; Sh <- sqrt(pmax(Vh, 0)); WSh <- Wh * Sh
C <- input$cc_budget; req(is.finite(C) && C > 0)
ch <- vapply(seq_len(H), function(h) {
v <- input[[paste0("cc_cost_", h)]]
if (is.null(v) || !is.finite(v) || v <= 0) 1 else v
}, numeric(1))
# Cost-optimal allocation: nₕ ∝ WₕSₕ / √cₕ (Cochran 1977, §5.5)
sqrt_ch <- sqrt(ch)
alloc_wts <- WSh / sqrt_ch
sum_wsh_sc <- sum(WSh * sqrt_ch) # Σ WₕSₕ√cₕ
if (sum_wsh_sc < 1e-12) { cost_res(NULL); return() }
nh_raw <- C * alloc_wts / sum_wsh_sc
nh_cost <- pmax(round(nh_raw), 1L)
n_total_cost <- sum(nh_cost)
# Neyman allocation for the same budget
# nₕ ∝ WₕSₕ → n × WₕSₕ/ΣWₕSₕ, where n = C / (ΣcₕWₕSₕ/ΣWₕSₕ) ...
sum_WSh <- sum(WSh)
sum_cWSh <- sum(ch * WSh)
n_neyman_budget <- if (sum_cWSh > 0) C * sum_WSh / sum_cWSh else NA_real_
nh_neyman_raw <- if (is.finite(n_neyman_budget)) n_neyman_budget * WSh / sum_WSh
else rep(NA_real_, H)
nh_neyman_b <- pmax(round(nh_neyman_raw), 1L)
# Variances: V = Σ Wh² Vh / nh
Wh2Vh <- Wh^2 * Vh
V_cost <- sum(Wh2Vh / nh_cost)
V_neyman_b <- if (all(is.finite(nh_neyman_b))) sum(Wh2Vh / nh_neyman_b) else NA_real_
# Theoretical minimum: (Σ WₕSₕ√cₕ)² / C
V_theory <- sum_wsh_sc^2 / C
cost_res(list(
H = H, C = C, ch = ch, Wh = Wh, Sh = Sh, WSh = WSh,
nh_cost = nh_cost, nh_neyman_b = nh_neyman_b,
n_total_cost = n_total_cost,
n_neyman_budget = n_neyman_budget,
V_cost = V_cost, V_neyman_b = V_neyman_b, V_theory = V_theory
))
})
output$ui_cost_results <- renderUI({
cr <- cost_res(); req(cr)
gain <- if (is.finite(cr$V_neyman_b) && cr$V_neyman_b > 0)
round(cr$V_neyman_b / cr$V_cost, 3) else NA_real_
tbl_rows <- lapply(seq_len(cr$H), function(h)
tags$tr(
tags$td(tags$strong(paste0("S", h))),
tags$td(round(cr$ch[h], 2)),
tags$td(round(cr$Wh[h], 4)),
tags$td(round(cr$Sh[h], 4)),
tags$td(tags$strong(cr$nh_cost[h])),
tags$td(if(is.finite(cr$nh_neyman_b[h])) cr$nh_neyman_b[h] else "—")
)
)
tagList(
div(class = "cc-result-banner",
div(class = "cc-result-stat",
div(class = "cc-result-val", format(cr$n_total_cost, big.mark=",")),
div(class = "cc-result-lbl", "Total n achieved")),
div(style="width:1px;background:rgba(5,150,105,0.25);align-self:stretch;"),
div(class = "cc-result-stat",
div(class = "cc-result-val", formatC(cr$V_cost, format="e", digits=3)),
div(class = "cc-result-lbl", "V*(ȳ) — cost-optimal")),
div(style="width:1px;background:rgba(5,150,105,0.25);align-self:stretch;"),
div(class = "cc-result-stat",
div(class = "cc-result-val",
if (is.finite(cr$V_neyman_b)) formatC(cr$V_neyman_b, format="e", digits=3) else "—"),
div(class = "cc-result-lbl", "V*(ȳ) — Neyman (same budget)")),
div(style="width:1px;background:rgba(5,150,105,0.25);align-self:stretch;"),
div(class = "cc-result-stat",
div(class = "cc-result-val",
if (is.finite(gain)) sprintf("%.3f×", gain) else "—"),
div(class = "cc-result-lbl", "Cost-opt. efficiency vs Neyman"))
),
div(class = "plot-container", style = "margin-bottom:1rem;",
plotlyOutput("plot_cost", height = "300px")),
tags$table(class = "dc-table",
tags$thead(tags$tr(
tags$th("Stratum"), tags$th("cₕ"), tags$th("Wₕ"), tags$th("Sₕ"),
tags$th("nₕ cost-optimal"), tags$th("nₕ Neyman (same C)")
)),
tags$tbody(tbl_rows),
tags$tfoot(tags$tr(style="font-weight:700;background:#F1F5F9;",
tags$td("Total"), tags$td("—"), tags$td("1.0000"), tags$td("—"),
tags$td(cr$n_total_cost),
tags$td(if(is.finite(sum(cr$nh_neyman_b))) sum(cr$nh_neyman_b) else "—")
))
),
div(class = "hsel-note", style = "margin-top:0.5rem;",
sprintf("Theoretical minimum V = (ΣWₕSₕ√cₕ)²/C = %.4e. Budget C = %.0f.",
cr$V_theory, cr$C))
)
})
output$plot_cost <- renderPlotly({
req(has_plotly)
cr <- cost_res(); req(cr)
s_lbl <- paste0("S", seq_len(cr$H))
ok_n <- is.finite(cr$nh_neyman_b)
y_max <- max(cr$nh_cost, cr$nh_neyman_b[ok_n], na.rm = TRUE)
y_range <- c(0, y_max * 1.28)
plotly::plot_ly(x=s_lbl, y=cr$nh_cost, type="bar", name="Cost-Optimal",
marker=list(color="#059669"), text=cr$nh_cost, textposition="outside") |>
plotly::add_trace(
x = s_lbl[ok_n], y = cr$nh_neyman_b[ok_n],
name = "Neyman (same budget)",
marker = list(color = "#F28E2B"),
text = cr$nh_neyman_b[ok_n], textposition = "outside") |>
plotly::layout(
barmode = "group",
xaxis = list(title=list(text="Stratum", font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9"),
yaxis = list(title=list(text="nₕ", font=list(size=12)),
zeroline=FALSE, gridcolor="#F1F5F9",
range = y_range),
legend = list(orientation="h", x=0.5, xanchor="center",
y=-0.22, font=list(size=11)),
paper_bgcolor="white", plot_bgcolor="white",
margin=list(t=30, b=60, l=70, r=20)
)
})
# ── H-Selector tab ───────────────────────────────────────────────────────────
output$ui_h_selector <- renderUI({
res0 <- strata_res()
if (is.null(res0)) {
return(div(class = "empty-state",
div(class = "empty-icon", "\U0001F4CA"),
div(class = "empty-title", "No results yet"),
div(class = "empty-sub", "Run Stratification first.")))
}
rows <- h_curve_res()
cur_H <- nrow(res0$h)
# ── Banner (always visible) ─────────────────────────────────────────────
banner <- div(class = "hsel-banner",
div(
div(class = "hsel-banner-title", "H-Selector Efficiency Curve"),
div(class = "hsel-banner-text",
"Runs stratifyR with COBYLA for H = 1 to 12 and plots ΣWₕSₕ vs H.",
" The elbow of the curve — detected automatically using the maximum",
" perpendicular-distance method — marks the recommended H.")
),
actionButton("run_h_curve",
tagList(icon(if (is.null(rows)) "play" else "rotate"),
if (is.null(rows)) " Compute H-Curve" else " Recompute"),
class = if (is.null(rows)) "btn btn-primary"
else "btn btn-outline-light btn-sm")
)
if (is.null(rows)) return(banner)
# ── Build summary data.frame ─────────────────────────────────────────────
df <- data.frame(
H = sapply(rows, `[[`, "H"),
WhSh = sapply(rows, `[[`, "WhSh"),
V = sapply(rows, `[[`, "V"),
stringsAsFactors = FALSE)
# ── Degenerate-result detection ──────────────────────────────────────────
# COBYLA can return WhSh ≈ 0 (boundary collapse) or WhSh that *increases*
# with H (the solver wandered to a worse region). Both are convergence
# failures. Flag them so they are excluded from elbow detection and
# annotated clearly in the table and plot.
first_valid_wsh <- df$WhSh[df$WhSh > 0][1]
degen_threshold <- if (!is.null(first_valid_wsh) && is.finite(first_valid_wsh))
first_valid_wsh * 0.02 else 1e-8
# A result is suspicious if WhSh is near-zero OR if it increased from the
# preceding reliable value (shouldn't happen for a correctly converged run).
df$ok <- df$WhSh > degen_threshold
# also flag increases vs the previous OK row
prev_ok_wsh <- NA_real_
for (i in seq_len(nrow(df))) {
if (df$ok[i]) {
if (!is.na(prev_ok_wsh) && df$WhSh[i] > prev_ok_wsh * 1.05)
df$ok[i] <- FALSE
else
prev_ok_wsh <- df$WhSh[i]
}
}
n_degen <- sum(!df$ok)
# % reduction vs previous H (among all rows for display)
df$gain_pct <- c(NA_real_,
round(100 * (df$WhSh[-nrow(df)] - df$WhSh[-1]) / df$WhSh[-nrow(df)], 1))
# ── Stagnation detection ──────────────────────────────────────────────────
# COBYLA often returns the SAME solution for consecutive H values in
# heavy-tailed settings (gain < 0.05% = rounds to "-0.0%").
# Flag runs of 2+ consecutive stagnant-reliable points so they render
# differently and the warning can mention them.
stagnant_threshold <- 0.05 # % gain below which we consider it stagnant
df$stagnant <- FALSE
for (i in seq_len(nrow(df))) {
if (df$ok[i] && !is.na(df$gain_pct[i]) &&
abs(df$gain_pct[i]) < stagnant_threshold)
df$stagnant[i] <- TRUE
}
# Only flag as stagnant if it's part of a run of >= 2
stag_run <- rle(df$stagnant & df$ok)
if (any(stag_run$values & stag_run$lengths >= 2)) {
idx <- 1L
for (k in seq_along(stag_run$lengths)) {
end_idx <- idx + stag_run$lengths[k] - 1L
if (!stag_run$values[k] || stag_run$lengths[k] < 2L)
df$stagnant[idx:end_idx] <- FALSE
idx <- end_idx + 1L
}
} else {
df$stagnant[] <- FALSE
}
n_stagnant <- sum(df$stagnant)
# % cumulative reduction vs H = 1
df$cum_pct <- round(100 * (df$WhSh[1] - df$WhSh) / df$WhSh[1], 1)
# ── Elbow detection — only on reliable rows ───────────────────────────────
df_ok <- df[df$ok, , drop = FALSE]
rec_H <- if (nrow(df_ok) >= 3L) {
h_sc <- (df_ok$H - df_ok$H[1]) / max(df_ok$H[nrow(df_ok)] - df_ok$H[1], 1e-12)
w_sc <- (df_ok$WhSh - df_ok$WhSh[nrow(df_ok)]) /
max(df_ok$WhSh[1] - df_ok$WhSh[nrow(df_ok)], 1e-12)
df_ok$H[which.min(h_sc + w_sc)]
} else if (nrow(df_ok) >= 1L) max(df_ok$H) else max(df$H)
# ── Heavy-tailed distribution warning ────────────────────────────────────
fitted_distr <- tryCatch(res0$distr, error = function(e) "")
is_heavy_tailed <- fitted_distr %in% c("cauchy", "pareto")
heavy_warn <- if (is_heavy_tailed || n_degen > 0 || n_stagnant > 0) {
div(style = paste0(
"background:#FFF3CD; border:1px solid #FFC107; border-radius:8px;",
" padding:0.75rem 1rem; margin-bottom:0.8rem; font-size:0.85rem;",
" color:#856404; display:flex; gap:0.6rem; align-items:flex-start;"),
icon("triangle-exclamation", style = "margin-top:2px; flex-shrink:0;"),
div(
if (is_heavy_tailed)
tags$p(style = "margin:0 0 0.3rem 0;",
tags$strong("Heavy-tailed distribution detected (", fitted_distr, ")."),
" COBYLA is documented to perform poorly on Cauchy and Pareto data —",
" it converges to poor local optima and can produce degenerate boundary solutions.",
" The H-Selector curve may be unreliable for these distributions.",
" For more accurate results, use DP in the main run and treat this curve as indicative only.")
else NULL,
if (n_degen > 0)
tags$p(style = "margin:0 0 0.3rem 0;",
tags$strong(n_degen, " degenerate result(s) detected"),
sprintf(" (H = %s) and excluded from the elbow calculation.",
paste(df$H[!df$ok], collapse = ", ")),
" These rows are marked ⚠ in the table. Degenerate solutions occur when COBYLA",
" collapses all stratum boundaries to nearly the same point (ΣWₕSₕ ≈ 0),",
" or when the objective worsened rather than improved with additional strata.")
else NULL,
if (n_stagnant > 0)
tags$p(style = "margin:0;",
tags$strong(n_stagnant, " stagnant result(s) detected"),
sprintf(" (H = %s).", paste(df$H[df$stagnant], collapse = ", ")),
" COBYLA returned essentially the same solution for these consecutive H values",
" (< 0.05% improvement), shown as grey ◆ on the plot.",
" This is a known COBYLA limitation for heavy-tailed distributions —",
" the solver gets stuck in the same local optimum regardless of how many strata are requested.")
else NULL
)
)
} else NULL
# ── Recommendation box ───────────────────────────────────────────────────
rec_note <- if (rec_H == cur_H)
"Your current H sits at the elbow — a good choice."
else
sprintf("The curve bends most sharply at H = %d. Beyond this, each extra stratum returns diminishing variance reduction.", rec_H)
rec_box <- div(class = "hsel-rec-box",
div(
div(class = "hsel-rec-h", paste0("H = ", rec_H)),
div(class = "hsel-rec-label",
if (n_degen > 0 || n_stagnant > 0) "Recommended strata (reliable points only)"
else "Recommended strata")
),
div(class = "hsel-rec-note", rec_note)
)
# ── Table rows ───────────────────────────────────────────────────────────
tbl_rows <- lapply(seq_len(nrow(df)), function(i) {
is_cur <- df$H[i] == cur_H
is_rec <- df$H[i] == rec_H
is_degen <- !df$ok[i]
is_stagnant <- df$stagnant[i]
row_cls <- if (is_degen) "hsel-degen-row"
else if (is_stagnant) "hsel-stagnant-row"
else if (is_cur && is_rec) "hsel-cur-row hsel-rec-row"
else if (is_cur) "hsel-cur-row"
else if (is_rec) "hsel-rec-row"
else ""
tag_lbl <- if (is_degen) " ⚠ degenerate"
else if (is_stagnant) " ◆ stagnant"
else if (is_cur && is_rec) " ◄ current ★ rec"
else if (is_cur) " ◄ current"
else if (is_rec) " ★ rec"
else ""
tags$tr(class = row_cls,
tags$td(tagList(as.character(df$H[i]),
tags$small(
style = if (is_degen) "color:#DC6C00; margin-left:4px;"
else if (is_stagnant) "color:#94A3B8; margin-left:4px;"
else "color:#94A3B8; margin-left:4px;",
tag_lbl))),
tags$td(if (is_degen) tags$em(style="color:#DC6C00;",
formatC(df$WhSh[i], format = "f", digits = 5))
else formatC(df$WhSh[i], format = "f", digits = 5)),
tags$td(formatC(df$V[i], format = "e", digits = 4)),
tags$td(if (is_degen || is.na(df$gain_pct[i])) "—"
else sprintf("%.1f%%", df$gain_pct[i])),
tags$td(if (is_degen) "—"
else sprintf("%.1f%%", df$cum_pct[i]))
)
})
tagList(
banner,
heavy_warn,
rec_box,
div(class = "plot-container", style = "margin-bottom: 1rem;",
plotlyOutput("plot_h_curve", height = "360px")),
tags$table(class = "hsel-table",
tags$thead(tags$tr(
tags$th("H"),
tags$th("ΣWₕSₕ"),
tags$th("V*(n)"),
tags$th("Marginal gain"),
tags$th("Cumulative reduction vs H=1")
)),
tags$tbody(tbl_rows)
),
div(class = "hsel-note",
"COBYLA solver used throughout (analytical pathway, reusing the fitted distribution from the main run).",
" Results may differ slightly from a full DP run.",
if (n_degen > 0) " Degenerate rows (⚠) are excluded from the elbow calculation." else NULL,
if (n_stagnant > 0) " Stagnant rows (◆) indicate COBYLA convergence to the same local optimum." else NULL),
# ── Method explanation ─────────────────────────────────────────────────
div(style = "margin-top: 1.4rem;",
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "1"), "Why plot ΣWₕSₕ against H?"),
div(class = "intro-body",
tags$p(
"Under Neyman (optimal) allocation the variance of the stratified mean is",
" ", tags$em("V*(ȳ) = (ΣWₕSₕ)² / n"), ".",
" Minimising ΣWₕSₕ over the strata boundaries is therefore the core",
" objective of stratifyR. As H increases, more strata are available to",
" adapt to the shape of the population distribution, and ΣWₕSₕ decreases."),
tags$p(
"However, the improvement is not linear — early strata deliver large gains",
" while each additional stratum beyond the elbow contributes only marginally.",
" The H-Selector curve makes this trade-off explicit: pick the H that sits",
" at the bend, and you capture most of the attainable variance reduction",
" at minimal survey complexity.")
)
),
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "2"), "The Elbow Method — formalised via Kneedle"),
div(class = "intro-body",
tags$p(
"The ", tags$strong("elbow method"), " is a general visual heuristic:",
" inspect the curve and pick the H at which it bends most sharply,",
" after which each additional stratum contributes diminishing returns.",
" Because 'by eye' detection is subjective, Satopaa et al. (2011) proposed",
" the ", tags$strong("Kneedle algorithm"), " as a formal, automated way to",
" locate that same elbow point mathematically.",
" stratifyR uses Kneedle to remove ambiguity from the selection."),
tags$p(
"Given the sequence of points ", tags$em("(H₁, f₁), …, (Hₘ, fₘ)"),
" where ", tags$em("fₕ = ΣWₕSₕ"), ", the algorithm proceeds in three steps:"),
tags$p(tags$strong("Step 1 — Normalise both axes to [0, 1]:")),
div(class = "intro-math",
withMathJax(
"\\[ \\tilde{H}_i = \\frac{H_i - H_1}{H_m - H_1}, \\qquad
\\tilde{f}_i = \\frac{f_i - f_m}{f_1 - f_m} \\]")),
tags$p(
"After normalisation the first point is at (0, 1) and the last at (1, 0)."),
tags$p(tags$strong("Step 2 — Reference line:")),
tags$p(
"Draw the straight line connecting",
" (0, 1) to (1, 0). Its equation is ",
tags$em("x + y = 1"), "."),
tags$p(tags$strong("Step 3 — Find the point of maximum perpendicular distance:")),
div(class = "intro-math",
withMathJax(
"\\[ d_i = \\frac{\\left|\\,\\tilde{H}_i + \\tilde{f}_i - 1\\,\\right|}{\\sqrt{2}} \\]")),
tags$p(
"For a convex-decreasing curve (ΣWₕSₕ always decreases, steeply at first),",
" every interior point lies ", tags$em("below"), " the diagonal, so",
" ", tags$em("d"), " simplifies to"),
div(class = "intro-math",
withMathJax(
"\\[ d_i = \\frac{1 - \\tilde{H}_i - \\tilde{f}_i}{\\sqrt{2}} \\]")),
div(class = "intro-key-result",
"The recommended H is ",
tags$em("argmin (\\(\\tilde{H}_i + \\tilde{f}_i\\))"),
" over the interior points i = 2, …, m−1.",
" This is the Kneedle formalisation of the elbow: the H at which the curve",
" departs furthest from a straight-line interpolation between the two",
" extremes. Both the concept (elbow) and its detection (Kneedle) point to",
" the same answer — Kneedle simply removes the subjectivity of visual inspection.")
)
),
div(class = "intro-section",
div(class = "intro-section-title",
div(class = "num-badge", "3"), "References"),
div(class = "intro-body",
div(class = "intro-ref",
tags$strong("Kneedle algorithm: "),
"Satopaa, V., Albrecht, J., Irwin, D., & Raghavan, B. (2011).",
" Finding a 'kneedle' in a haystack: Detecting knee points in system behaviour.",
" ", tags$em("31st IEEE International Conference on Distributed Computing Systems Workshops"),
", pp. 166–171.", tags$br(),
tags$strong("Optimal stratification: "),
"Cochran, W.G. (1977). ",
tags$em("Sampling Techniques"), ", 3rd ed. Wiley. — Ch. 5.",
tags$br(),
tags$strong("Neyman allocation: "),
"Neyman, J. (1934). On the two different aspects of the representative method. ",
tags$em("Journal of the Royal Statistical Society"), " 97, 558–625.",
tags$br(),
tags$strong("Optimum boundaries: "),
"Khan, M.G.M., Nand, N., & Ahmad, N. (2008). Determining optimum strata",
" boundaries using mathematical programming. ",
tags$em("Survey Methodology"), " 34(2), 91–102."
)
)
)
)
)
})
output$plot_h_curve <- renderPlotly({
req(has_plotly)
rows <- h_curve_res(); req(!is.null(rows) && length(rows) >= 2L)
res0 <- strata_res(); req(res0)
cur_H <- nrow(res0$h)
df <- data.frame(
H = sapply(rows, `[[`, "H"),
WhSh = sapply(rows, `[[`, "WhSh"),
V = sapply(rows, `[[`, "V"),
stringsAsFactors = FALSE)
df$gain_pct <- c(NA_real_,
round(100 * (df$WhSh[-nrow(df)] - df$WhSh[-1]) / df$WhSh[-nrow(df)], 1))
# ── Degenerate + stagnant detection (mirror renderUI logic) ─────────────
first_valid <- df$WhSh[df$WhSh > 0][1]
degen_thr <- if (!is.null(first_valid) && is.finite(first_valid))
first_valid * 0.02 else 1e-8
df$ok <- df$WhSh > degen_thr
prev_ok <- NA_real_
for (i in seq_len(nrow(df))) {
if (df$ok[i]) {
if (!is.na(prev_ok) && df$WhSh[i] > prev_ok * 1.05) df$ok[i] <- FALSE
else prev_ok <- df$WhSh[i]
}
}
# Stagnation: gain < 0.05% AND part of a run of >= 2
df$stagnant <- FALSE
for (i in seq_len(nrow(df))) {
if (df$ok[i] && !is.na(df$gain_pct[i]) && abs(df$gain_pct[i]) < 0.05)
df$stagnant[i] <- TRUE
}
stag_run <- rle(df$stagnant & df$ok)
idx <- 1L
for (k in seq_along(stag_run$lengths)) {
end_idx <- idx + stag_run$lengths[k] - 1L
if (!stag_run$values[k] || stag_run$lengths[k] < 2L)
df$stagnant[idx:end_idx] <- FALSE
idx <- end_idx + 1L
}
# "Good" = ok and not stagnant — used for the solid line and elbow
df$good <- df$ok & !df$stagnant
df_good <- df[df$good, , drop = FALSE]
df_ok <- df[df$ok, , drop = FALSE] # for annotation lookup
# Elbow on good points only
rec_H <- if (nrow(df_good) >= 3L) {
h_sc <- (df_good$H - df_good$H[1]) / max(df_good$H[nrow(df_good)] - df_good$H[1], 1e-12)
w_sc <- (df_good$WhSh - df_good$WhSh[nrow(df_good)]) /
max(df_good$WhSh[1] - df_good$WhSh[nrow(df_good)], 1e-12)
df_good$H[which.min(h_sc + w_sc)]
} else if (nrow(df_good) >= 1L) max(df_good$H)
else if (nrow(df_ok) >= 1L) max(df_ok$H)
else max(df$H)
is_special <- df$H == cur_H | df$H == rec_H
# Point colours: degenerate=orange×, stagnant=grey◆, current=red●, elbow=green●, normal=blue●
pt_col <- ifelse(!df$ok, "#F59E0B",
ifelse(df$stagnant, "#94A3B8",
ifelse(df$H == cur_H, "#C0392B",
ifelse(df$H == rec_H, "#059669", "#4A6FA5"))))
pt_sz <- ifelse(is_special & df$ok & !df$stagnant, 14, 8)
pt_sym <- ifelse(!df$ok, "x",
ifelse(df$stagnant, "diamond", "circle"))
hover_txt <- sprintf(
paste0("<b>H = %d</b>%s<br>",
"ΣWₕSₕ = %.5f<br>",
"V*(n) = %.6f%s<extra></extra>"),
df$H,
ifelse(!df$ok, " ⚠ degenerate",
ifelse(df$stagnant, " ◆ stagnant (COBYLA stuck)", "")),
df$WhSh, df$V,
ifelse(is.na(df$gain_pct) | !df$ok | df$stagnant, "",
sprintf("<br>Reduction = −%.1f%%", df$gain_pct)))
# Gain labels only on good (non-special, non-stagnant, non-degen) points
gain_lbl <- ifelse(!df$ok | df$stagnant | is_special | is.na(df$gain_pct), "",
sprintf("−%.1f%%", df$gain_pct))
txt_pos <- ifelse(df$H %% 2 == 0, "top center", "bottom center")
fig <- plotly::plot_ly()
# Solid line through good (non-stagnant) points only
if (nrow(df_good) >= 2L) {
fig <- fig |> plotly::add_trace(
x = df_good$H, y = df_good$WhSh,
type = "scatter", mode = "lines",
line = list(color = "#CBD5E1", width = 2.2),
showlegend = FALSE, hoverinfo = "none")
}
# Dashed grey line through stagnant points (if any) — connects them visually
# but signals they're unreliable
df_stag <- df[df$stagnant, , drop = FALSE]
# include the last good point before the stagnant run as anchor
if (nrow(df_stag) >= 1L && nrow(df_good) >= 1L) {
anchor_H <- max(df_good$H[df_good$H < min(df_stag$H)], na.rm = TRUE)
anchor_row <- df_good[df_good$H == anchor_H, , drop = FALSE]
if (nrow(anchor_row) == 1L)
df_stag_line <- rbind(anchor_row[, c("H","WhSh")], df_stag[, c("H","WhSh")])
else
df_stag_line <- df_stag[, c("H","WhSh")]
if (nrow(df_stag_line) >= 2L) {
fig <- fig |> plotly::add_trace(
x = df_stag_line$H, y = df_stag_line$WhSh,
type = "scatter", mode = "lines",
line = list(color = "#CBD5E1", width = 1.5, dash = "dot"),
showlegend = FALSE, hoverinfo = "none")
}
}
# All points
fig <- fig |> plotly::add_trace(
x = df$H, y = df$WhSh,
type = "scatter", mode = "markers+text",
marker = list(color = pt_col, size = pt_sz, symbol = pt_sym,
line = list(color = "white", width = 2)),
text = gain_lbl,
textposition = txt_pos,
textfont = list(size = 10, color = "#64748B"),
hovertemplate = hover_txt,
showlegend = FALSE)
# Annotations for special points — only on good (non-stagnant, non-degen) rows
annots <- list()
if (cur_H %in% df_good$H) {
y_cur <- df_good$WhSh[df_good$H == cur_H]
if (cur_H == rec_H) {
annots <- c(annots, list(list(
x = cur_H, y = y_cur, xanchor = "center", yanchor = "top",
text = sprintf("<b>H = %d (current = elbow)</b>", cur_H),
showarrow = TRUE, arrowhead = 0,
arrowcolor = "#059669", arrowwidth = 1.5, ax = 0, ay = 40,
font = list(size = 11, color = "#065F46"),
bgcolor = "#F0FDF4", bordercolor = "#059669", borderwidth = 1.5, borderpad = 5)))
} else {
annots <- c(annots, list(list(
x = cur_H, y = y_cur, xanchor = "center", yanchor = "top",
text = sprintf("<b>H = %d (current)</b>", cur_H),
showarrow = TRUE, arrowhead = 0,
arrowcolor = "#C0392B", arrowwidth = 1.5, ax = 0, ay = 40,
font = list(size = 11, color = "#C0392B"),
bgcolor = "#FEF2F2", bordercolor = "#C0392B", borderwidth = 1.5, borderpad = 5)))
}
}
if (rec_H != cur_H && rec_H %in% df_good$H) {
annots <- c(annots, list(list(
x = rec_H, y = df_good$WhSh[df_good$H == rec_H],
xanchor = "center", yanchor = "top",
text = sprintf("<b>H = %d (elbow)</b>", rec_H),
showarrow = TRUE, arrowhead = 0,
arrowcolor = "#059669", arrowwidth = 1.5, ax = 0, ay = 40,
font = list(size = 11, color = "#065F46"),
bgcolor = "#F0FDF4", bordercolor = "#059669", borderwidth = 1.5, borderpad = 5)))
}
fig |> plotly::layout(
xaxis = list(
title = list(text = "Number of Strata (H)", font = list(size = 13)),
tickmode = "linear", tick0 = 1, dtick = 1,
zeroline = FALSE, gridcolor = "#F1F5F9",
showline = TRUE, linecolor = "#E2E8F0"),
yaxis = list(
title = list(text = "ΣWₕSₕ (Neyman objective)", font = list(size = 13)),
zeroline = FALSE, gridcolor = "#F1F5F9",
showline = TRUE, linecolor = "#E2E8F0"),
hovermode = "closest",
paper_bgcolor = "white", plot_bgcolor = "white",
margin = list(t = 20, b = 60, l = 80, r = 30),
annotations = annots
)
})
# ── Summary tab ───────────────────────────────────────────────────────────────
output$ui_summary <- renderUI({
res <- strata_res(); req(res)
solver_lbl <- toupper(res$method)
distr_lbl <- if (!is.null(res$distr) && nzchar(res$distr)) res$distr else "data-fitted"
conv_cls <- if (isTRUE(res$converged)) "green" else "amber"
conv_lbl <- if (isTRUE(res$converged)) "Converged" else "Not converged"
tagList(
div(class = "summary-meta",
div(class = "meta-item",
span("Solver:"),
span(class = "meta-badge", solver_lbl)),
div(class = "meta-item",
span("Distribution:"),
span(class = "meta-badge", distr_lbl)),
div(class = "meta-item",
span("Status:"),
span(class = paste("meta-badge", conv_cls), conv_lbl)),
div(class = "meta-item",
span("Total variance V*:"),
span(class = "meta-badge",
formatC(res$WhShTot^2 / res$nhTot, format = "e", digits = 3)))
),
DT::datatable(
bounds_df(),
rownames = FALSE,
class = "compact stripe hover",
options = list(
dom = "t",
pageLength = 15,
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) |>
DT::formatStyle("Stratum", fontWeight = "bold", color = "#1B3A5C") |>
DT::formatStyle(c("Lower","Upper"), color = "#C0392B", fontWeight = "600")
)
})
# ── Shared density computation helper ────────────────────────────────────────
.compute_density_2d <- function(res, n_pts = 512L) {
pop_data <- res$data_internal
distr <- res$distr
maxval <- res$maxval
OSB <- res$OSB
x_lo <- maxval * res$initval
x_hi <- max(maxval * res$finval, max(OSB))
full_B <- c(x_lo, OSB[OSB > x_lo & OSB < x_hi], x_hi)
x_seq <- seq(x_lo, x_hi, length.out = n_pts)
params <- tryCatch({
fit <- res$fit
if (!is.null(fit$estimate)) fit$estimate
else if (!is.null(fit$fit$estimate)) fit$fit$estimate
else NULL
}, error = function(e) NULL)
dens <- tryCatch({
p2 <- params
switch(distr,
norm = stats::dnorm(x_seq, mean = p2["mean"], sd = p2["sd"]),
lnorm = stats::dlnorm(x_seq, meanlog = p2["meanlog"], sdlog = p2["sdlog"]),
gamma = stats::dgamma(x_seq, shape = p2["shape"], rate = p2["rate"]),
weibull = stats::dweibull(x_seq, shape = p2["shape"], scale = p2["scale"]),
exp = stats::dexp(x_seq, rate = p2["rate"]),
cauchy = stats::dcauchy(x_seq, location = p2["location"], scale = p2["scale"]),
unif = stats::dunif(x_seq, min = p2["min"], max = p2["max"]),
{
if (!is.null(pop_data)) {
d2 <- density(pop_data, n = n_pts, from = x_lo, to = x_hi)
stats::approx(d2$x, d2$y, xout = x_seq, rule = 2)$y
} else rep(1 / (x_hi - x_lo), n_pts)
}
)
}, error = function(e) rep(1 / (x_hi - x_lo), n_pts))
dens[!is.finite(dens)] <- 0
list(x_seq = x_seq, dens = dens, full_B = full_B, x_lo = x_lo, x_hi = x_hi,
pop_data = pop_data)
}
# ── Build shared plotly 2D chart ──────────────────────────────────────────────
.build_plotly_2d <- function(res, pal, ttl = NULL, show_stats_box = FALSE) {
req(has_plotly)
H <- nrow(res$h)
OSB <- res$OSB
cols <- rep_len(pal, H)
dc <- .compute_density_2d(res)
x_seq <- dc$x_seq; dens <- dc$dens; full_B <- dc$full_B
pop_data <- dc$pop_data
d_max <- max(dens, na.rm = TRUE)
x_lo <- dc$x_lo; x_hi <- dc$x_hi
x_rng <- x_hi - x_lo
fig <- plotly::plot_ly()
# ── Histogram ──────────────────────────────────────────────────────────────
if (!is.null(pop_data) && length(pop_data) >= 10L) {
hh <- hist(pop_data, breaks = "Sturges", plot = FALSE)
fig <- fig |>
plotly::add_bars(
x = hh$mids, y = hh$density,
name = "Histogram",
marker = list(color = "rgba(91,158,214,0.40)",
line = list(color = "rgba(91,158,214,0.75)", width = 0.7)),
width = diff(hh$breaks)[1],
hovertemplate = "Midpoint: %{x:.2f}<br>Density: %{y:.5f}<extra></extra>"
)
}
# ── Stratum fills ──────────────────────────────────────────────────────────
for (h in seq_len(H)) {
lo_h <- full_B[h]; hi_h <- full_B[h + 1L]
idx <- x_seq >= lo_h & x_seq <= hi_h
if (sum(idx) < 2L) next
wh <- res$Wh[h]
sh <- round(sqrt(max(res$Vh[h], 0)), 4)
nh <- res$nh[h]
fig <- fig |>
plotly::add_trace(
x = c(x_seq[idx], rev(x_seq[idx])),
y = c(dens[idx], rep(0, sum(idx))),
type = "scatter", mode = "lines", fill = "toself",
fillcolor = grDevices::adjustcolor(cols[h], alpha.f = 0.38),
line = list(color = "transparent"),
name = sprintf("Stratum %d", h),
hovertemplate = sprintf(
"<b>Stratum %d</b><br>Wₕ = %.4f<br>Sₕ = %.4f<br>nₕ = %d<extra></extra>",
h, wh, sh, nh)
)
}
# ── Density curve ──────────────────────────────────────────────────────────
distr_lbl <- if (!is.null(res$distr) && nzchar(res$distr)) res$distr else "KDE"
fig <- fig |>
plotly::add_trace(
x = x_seq, y = dens, type = "scatter", mode = "lines",
line = list(color = "#1a1a2e", width = 2.6),
name = sprintf("f(y) — %s", distr_lbl),
hovertemplate = "y = %{x:.3f}<br>f(y) = %{y:.5f}<extra></extra>"
)
# ── Boundary lines ─────────────────────────────────────────────────────────
for (i in seq_along(OSB)) {
fig <- fig |>
plotly::add_trace(
x = rep(OSB[i], 2), y = c(0, d_max * 1.08),
type = "scatter", mode = "lines",
line = list(color = "#C0392B", width = 2, dash = "dash"),
showlegend = FALSE, hoverinfo = "none"
)
}
# ── Boundary annotations (horizontal labels at top, matching screenshot) ───
sub_h <- c("₁","₂","₃","₄","₅",
"₆","₇","₈","₉")
osb_annots <- lapply(seq_along(OSB), function(i) {
list(x = OSB[i], y = d_max * 0.52,
xref = "x", yref = "y",
xanchor = "center", yanchor = "middle",
text = sprintf("b%s = %.3f", sub_h[min(i, 9L)], OSB[i]),
textangle = -90,
showarrow = FALSE,
font = list(size = 10, color = "#C0392B"),
bgcolor = "rgba(255,255,255,0.60)",
borderpad = 2)
})
# ── Optional stats box (for Interactive Explorer) ─────────────────────────
stats_annot <- if (show_stats_box) {
list(list(
x = 0.99, y = 0.99, xref = "paper", yref = "paper",
xanchor = "right", yanchor = "top",
text = sprintf(
paste0("<b>Optimal solution</b><br>",
"ΣWₕSₕ = %.5f<br>",
"V(n) = %.6f<br>",
"n<sub>total</sub> = %d"),
res$WhShTot, res$WhShTot^2 / res$nhTot, res$nhTot),
showarrow = FALSE,
font = list(size = 12, color = "#C0392B"),
bgcolor = "rgba(255,240,240,0.92)",
bordercolor = "#C0392B", borderwidth = 1.2, borderpad = 8
))
} else list()
all_annots <- c(osb_annots, stats_annot)
title_txt <- if (!is.null(ttl) && nzchar(ttl)) ttl else
sprintf("%d-Strata Optimal Solution", H)
fig |> plotly::layout(
title = list(text = title_txt,
font = list(size = 16, color = "#1A2B3C", family = "sans-serif"),
x = 0.5, xanchor = "center"),
xaxis = list(title = list(text = "Y (variable value)", font = list(size = 14)),
tickfont = list(size = 13),
zeroline = FALSE, gridcolor = "#EEEEEE"),
yaxis = list(title = list(text = "Density f(y)", font = list(size = 14)),
tickfont = list(size = 13),
zeroline = FALSE, gridcolor = "#EEEEEE",
range = c(0, d_max * 1.15)),
paper_bgcolor = "white",
plot_bgcolor = "white",
hovermode = "x unified",
legend = list(orientation = "h", x = 0.01, y = -0.20,
xanchor = "left", font = list(size = 12)),
margin = list(t = 60, b = 70, l = 65, r = 25),
annotations = all_annots
)
}
# ── 2D plot (plotly, clean — matches screenshot) ──────────────────────────────
output$plot_2d <- renderPlotly({
req(has_plotly)
res <- strata_res(); req(res)
.build_plotly_2d(res, current_palette(), plot_main(), show_stats_box = TRUE)
})
# ── 2D info panel (distribution + strata summary) ─────────────────────────────
output$ui_2d_info <- renderUI({
res <- strata_res(); req(res)
distr <- if (!is.null(res$distr) && nzchar(res$distr)) res$distr else "data"
dist_disp <- switch(distr,
norm = "Normal Distribution",
lnorm = "Log-Normal Distribution",
gamma = "Gamma Distribution",
weibull = "Weibull Distribution",
exp = "Exponential Distribution",
cauchy = "Cauchy Distribution",
unif = "Uniform Distribution",
pareto = "Pareto Distribution",
data = "Data-fitted (non-parametric)",
paste0(tools::toTitleCase(distr), " Distribution")
)
params <- tryCatch({
fit <- res$fit
if (!is.null(fit$estimate)) fit$estimate
else if (!is.null(fit$fit$estimate)) fit$fit$estimate
else NULL
}, error = function(e) NULL)
greek <- switch(distr,
norm = c(mean = "μ", sd = "σ"),
lnorm = c(meanlog = "μ", sdlog = "σ"),
gamma = c(shape = "α", rate = "β"),
weibull = c(shape = "k", scale = "λ"),
exp = c(rate = "λ"),
cauchy = c(location = "μ₀", scale = "γ"),
unif = c(min = "a", max = "b"),
NULL
)
# LaTeX equation via MathJax
eq_latex <- STRAT_PDF_EQUATIONS[[distr]]
# Parameters: "param_name (symbol) = value"
param_badges <- if (!is.null(params) && length(params) > 0L) {
lapply(seq_along(params), function(i) {
nm <- names(params)[i]
val <- as.numeric(params[i])
sym <- if (!is.null(greek) && nm %in% names(greek)) greek[[nm]] else nm
tags$span(class = "param-badge",
sprintf("%s (%s) = %.4f", nm, sym, val))
})
} else NULL
bd <- bounds_df()
pal <- current_palette()
cols <- rep_len(pal, nrow(bd))
strata_rows <- lapply(seq_len(nrow(bd)), function(i) {
row <- bd[i, ]
div(class = "strata-row",
tags$span(class = "strata-swatch",
style = sprintf("background:%s;", cols[i])),
tags$span(class = "strata-label", sprintf("Stratum %d", i)),
tags$span(class = "strata-interval",
sprintf("[%s – %s]",
formatC(row$Lower, format = "f", digits = 2, big.mark = ","),
formatC(row$Upper, format = "f", digits = 2, big.mark = ","))),
tags$span(class = "strata-stat", sprintf("Wh = %.4f", row$Wh)),
tags$span(class = "strata-stat", sprintf("Sh = %.4f", row$Sh)),
tags$span(class = "strata-stat", sprintf("nh = %d", row$nh))
)
})
tagList(
div(class = "info-dist-card",
div(class = "info-dist-name", dist_disp),
if (!is.null(eq_latex))
withMathJax(div(class = "info-eqn", HTML(eq_latex))),
if (!is.null(param_badges))
div(class = "info-params", param_badges)
),
div(class = "info-strata", strata_rows)
)
})
# ── 3D plot ────────────────────────────────────────────────────────────────────
output$plot_3d <- renderPlotly({
req(has_plotly)
res <- strata_res(); req(res)
stratifyR:::plot.strata(res, type = "3d",
palette = current_palette(),
main = plot_main())
})
# ── Interactive Explorer v2 — Shiny-native sliders ───────────────────────────
#
# Layout: plotly 2D density plot (reacts to sliders)
# stats strip (current ΣWhSh | optimal | gap)
# sliderInput controls panel (one per internal boundary)
#
# Completely avoids plotly's built-in slider coordinate system.
# ---- helper: compute ΣWhSh for any set of boundary values ------------------
.ib_whsh <- function(res, cur_b) {
dc <- .compute_density_2d(res)
x_lo <- dc$x_lo; x_hi <- dc$x_hi
x_seq <- dc$x_seq; dens <- dc$dens
fB <- c(x_lo, sort(cur_b), x_hi)
Hh <- length(cur_b) + 1L
pop <- dc$pop_data
if (!is.null(pop) && length(pop) >= 2L) {
ds <- sort(pop)
N_pop <- if (!is.null(res$N) && res$N > 0) res$N else length(pop)
obj <- 0
for (h in seq_len(Hh)) {
sl <- if (h == 1L) ds[ds >= fB[h] & ds <= fB[h + 1L]]
else ds[ds > fB[h] & ds <= fB[h + 1L]]
if (length(sl) < 2L) next
obj <- obj + (length(sl) / N_pop) * sqrt(stats::var(sl))
}
return(obj)
}
# density-integral fallback
tot <- sum(diff(x_seq) * (head(dens, -1) + tail(dens, -1))) / 2
Wh_v <- numeric(Hh); Sh_v <- numeric(Hh)
for (h in seq_len(Hh)) {
idx <- x_seq >= fB[h] & x_seq < fB[h + 1L]
if (sum(idx) < 2L) next
xs <- x_seq[idx]; fs <- dens[idx]; dx <- diff(xs)
Wh_v[h] <- sum(dx * (head(fs, -1) + tail(fs, -1))) / 2 / max(tot, 1e-12)
if (Wh_v[h] < 1e-12) next
mu <- sum(dx * (head(xs * fs, -1) + tail(xs * fs, -1))) / 2 /
Wh_v[h] / max(tot, 1e-12)
mu2 <- sum(dx * (head(xs ^ 2 * fs, -1) + tail(xs ^ 2 * fs, -1))) / 2 /
Wh_v[h] / max(tot, 1e-12)
Sh_v[h] <- sqrt(max(mu2 - mu ^ 2, 0))
}
sum(Wh_v * Sh_v)
}
# ---- reactive: current slider boundary values --------------------------------
ib_cur_b <- reactive({
res <- strata_res(); req(res)
H <- nrow(res$h)
b <- sapply(seq_len(H - 1L), function(i) {
v <- input[[paste0("ib_b_", i)]]
if (is.null(v)) res$OSB[i] else v
})
sort(b)
})
# ---- tab container: renders plot + stats strip + slider widgets ---------------
output$ui_ib_v2 <- renderUI({
res <- strata_res()
if (is.null(res)) {
return(div(class = "empty-state",
div(class = "empty-icon", "\U0001F4CA"),
div(class = "empty-title", "No results yet"),
div(class = "empty-sub", "Run stratification first.")))
}
if (!has_plotly) {
return(div(class = "empty-state",
div("Install plotly:", tags$code("install.packages('plotly')"))))
}
H <- nrow(res$h)
n_sl <- H - 1L
x_lo <- res$maxval * res$initval
x_hi <- max(res$maxval * res$finval, max(res$OSB))
x_rng <- x_hi - x_lo
# full_B has H+1 entries: lo, b1…b(H-1), hi
opt_B <- res$OSB[res$OSB > x_lo & res$OSB < x_hi]
full_B <- c(x_lo, opt_B, x_hi)
sliders_ui <- lapply(seq_len(n_sl), function(i) {
lo_i <- round(full_B[i] + 0.02 * x_rng, 4)
hi_i <- round(full_B[i + 2L] - 0.02 * x_rng, 4)
if (lo_i >= hi_i) lo_i <- round(full_B[i] + x_rng * 0.001, 4)
stp <- max(0.0001, round(x_rng / 400, 4))
div(style = "flex: 1 1 160px; padding: 0 10px;",
sliderInput(
inputId = paste0("ib_b_", i),
label = HTML(sprintf(
"<b style='color:#1B3A5C'>b<sub>%d</sub></b>", i)),
min = lo_i,
max = hi_i,
value = round(res$OSB[i], 4),
step = stp,
width = "100%"
)
)
})
tagList(
div(class = "plot-container",
plotlyOutput("plot_ib_v2", height = "460px")),
div(style = "margin-top: 0.55rem;",
uiOutput("ui_ib_stats")),
div(style = paste0(
"margin-top: 0.6rem; background: #F3F4F6;",
" border-radius: 10px; padding: 0.85rem 1.1rem;",
" border: 1px solid #E2E8F0;"),
div(style = paste0(
"font-size: 0.68rem; font-weight: 700;",
" letter-spacing: 1.5px; text-transform: uppercase;",
" color: #94A3B8; margin-bottom: 0.55rem;"),
"■ Boundary Controls — drag to explore"),
div(style = "display: flex; flex-wrap: wrap; gap: 4px;",
sliders_ui)
)
)
})
# ---- reactive plot -----------------------------------------------------------
output$plot_ib_v2 <- renderPlotly({
req(has_plotly)
res <- strata_res(); req(res)
cur_b <- ib_cur_b()
H <- nrow(res$h)
cols <- rep_len(current_palette(), H)
dc <- .compute_density_2d(res)
x_seq <- dc$x_seq; dens <- dc$dens
pop_data <- dc$pop_data
d_max <- max(dens, na.rm = TRUE)
x_lo <- dc$x_lo; x_hi <- dc$x_hi
full_B_cur <- c(x_lo, cur_b, x_hi)
opt_OSB <- res$OSB[res$OSB > x_lo & res$OSB < x_hi]
fig <- plotly::plot_ly()
# Histogram (if raw data available)
if (!is.null(pop_data) && length(pop_data) >= 10L) {
hh <- hist(pop_data, breaks = "Sturges", plot = FALSE)
fig <- fig |>
plotly::add_bars(
x = hh$mids, y = hh$density, name = "Histogram",
marker = list(color = "rgba(91,158,214,0.38)",
line = list(color = "rgba(91,158,214,0.72)", width = 0.6)),
width = diff(hh$breaks)[1],
hovertemplate = "Mid: %{x:.2f} Density: %{y:.5f}<extra></extra>")
}
# Stratum fills using CURRENT slider boundaries
for (h in seq_len(H)) {
lo_h <- full_B_cur[h]; hi_h <- full_B_cur[h + 1L]
idx <- x_seq >= lo_h & x_seq <= hi_h
if (sum(idx) < 2L) next
fig <- fig |>
plotly::add_trace(
x = c(x_seq[idx], rev(x_seq[idx])),
y = c(dens[idx], rep(0, sum(idx))),
type = "scatter", mode = "lines", fill = "toself",
fillcolor = grDevices::adjustcolor(cols[h], alpha.f = 0.38),
line = list(color = "transparent"),
name = sprintf("Stratum %d", h),
hovertemplate = sprintf("<b>Stratum %d</b><extra></extra>", h))
}
# Density curve
distr_lbl <- if (!is.null(res$distr) && nzchar(res$distr)) res$distr else "KDE"
fig <- fig |>
plotly::add_trace(
x = x_seq, y = dens, type = "scatter", mode = "lines",
line = list(color = "#1a1a2e", width = 2.5),
name = sprintf("f(y) — %s", distr_lbl),
hovertemplate = "y=%{x:.3f} f(y)=%{y:.5f}<extra></extra>")
# Optimal boundary reference lines (green dotted)
for (osb in opt_OSB) {
fig <- fig |>
plotly::add_trace(
x = rep(osb, 2), y = c(0, d_max * 1.08),
type = "scatter", mode = "lines",
line = list(color = "rgba(16,120,80,0.48)", width = 1.6, dash = "dot"),
showlegend = FALSE, hoverinfo = "none")
}
# Current slider boundaries (solid red)
for (bv in cur_b) {
fig <- fig |>
plotly::add_trace(
x = rep(bv, 2), y = c(0, d_max * 1.08),
type = "scatter", mode = "lines",
line = list(color = "#C0392B", width = 2.2),
showlegend = FALSE, hoverinfo = "none")
}
# Rotated boundary value labels
sub_h <- c("₁","₂","₃","₄","₅",
"₆","₇","₈","₉")
bnd_annots <- lapply(seq_along(cur_b), function(i) {
list(x = cur_b[i], y = d_max * 0.52,
xref = "x", yref = "y",
xanchor = "center", yanchor = "middle",
text = sprintf("b%s=%.3f", sub_h[min(i, 9L)], cur_b[i]),
textangle = -90, showarrow = FALSE,
font = list(size = 10, color = "#C0392B"),
bgcolor = "rgba(255,255,255,0.65)", borderpad = 2)
})
# Stats box (top-right corner)
cur_obj <- tryCatch(.ib_whsh(res, cur_b), error = function(e) NA_real_)
opt_obj <- res$WhShTot
pct_gap <- if (is.finite(cur_obj) && opt_obj > 0)
100 * (cur_obj - opt_obj) / opt_obj else NA_real_
gap_col <- if (!is.na(pct_gap) && abs(pct_gap) <= 0.5) "#065F46" else "#C0392B"
gap_str <- if (is.finite(pct_gap)) sprintf("%+.2f%%", pct_gap) else "N/A"
stats_annot <- list(
x = 0.99, y = 0.99, xref = "paper", yref = "paper",
xanchor = "right", yanchor = "top",
text = sprintf(
paste0("<b>Current:</b> ΣW<sub>h</sub>S<sub>h</sub> = %.5f<br>",
"<b>Optimal:</b> ΣW<sub>h</sub>S<sub>h</sub> = %.5f<br>",
"<span style='color:%s'><b>Gap = %s</b></span>"),
cur_obj, opt_obj, gap_col, gap_str),
showarrow = FALSE,
font = list(size = 11.5, color = "#1A2B3C"),
bgcolor = "rgba(255,248,240,0.94)",
bordercolor = "#C0392B", borderwidth = 1.3, borderpad = 9)
ttl <- { t <- plot_main()
if (!is.null(t) && nzchar(t)) t else
sprintf("%d-Strata Boundary Explorer", H) }
fig |> plotly::layout(
title = list(text = ttl,
font = list(size = 15, color = "#1A2B3C"),
x = 0.5, xanchor = "center"),
xaxis = list(title = list(text = "Y (variable value)", font = list(size = 13)),
tickfont = list(size = 12),
zeroline = FALSE, gridcolor = "#EEEEEE"),
yaxis = list(title = list(text = "Density f(y)", font = list(size = 13)),
tickfont = list(size = 12),
zeroline = FALSE, gridcolor = "#EEEEEE",
range = c(0, d_max * 1.18)),
paper_bgcolor = "white",
plot_bgcolor = "white",
hovermode = "x unified",
legend = list(orientation = "h", x = 0.01, y = -0.18,
xanchor = "left", font = list(size = 11)),
margin = list(t = 55, b = 70, l = 65, r = 25),
annotations = c(bnd_annots, list(stats_annot))
)
})
# ---- stats strip below plot --------------------------------------------------
output$ui_ib_stats <- renderUI({
res <- strata_res(); req(res)
cur_b <- ib_cur_b()
cur_obj <- tryCatch(.ib_whsh(res, cur_b), error = function(e) NA_real_)
opt_obj <- res$WhShTot
n_tot <- res$nhTot
pct_gap <- if (is.finite(cur_obj) && opt_obj > 0)
round(100 * (cur_obj - opt_obj) / opt_obj, 2) else NA_real_
gap_cls <- if (!is.na(pct_gap) && abs(pct_gap) <= 0.5) "green" else "amber"
gap_lbl <- if (is.finite(pct_gap)) sprintf("%+.2f%%", pct_gap) else "N/A"
div(class = "summary-meta",
div(class = "meta-item",
span("Current ΣWₕSₕ:"),
span(class = "meta-badge",
if (is.finite(cur_obj)) formatC(cur_obj, digits = 5, format = "f") else "—")),
div(class = "meta-item",
span("Optimal ΣWₕSₕ:"),
span(class = "meta-badge green",
formatC(opt_obj, digits = 5, format = "f"))),
div(class = "meta-item",
span("Gap vs optimal:"),
span(class = paste("meta-badge", gap_cls), gap_lbl)),
div(class = "meta-item",
span("V*(n) current:"),
span(class = "meta-badge",
if (is.finite(cur_obj))
formatC(cur_obj ^ 2 / n_tot, digits = 6, format = "f") else "—"))
)
})
if (FALSE) {
# ── placeholder – never executed ──────────────────────────────────────────
n_sliders <- 0L
# ── Inline objective functions (replicate package internals exactly) ─────────
data_sorted <- if (!is.null(pop_data) && length(pop_data) >= 2L)
sort(pop_data) else NULL
obj_data <- function(b_test) {
fB <- c(x_lo, sort(b_test), x_hi)
Hh <- length(b_test) + 1L
obj <- 0
for (h in seq_len(Hh)) {
sl <- if (h == 1L)
data_sorted[data_sorted >= fB[h] & data_sorted <= fB[h + 1L]]
else
data_sorted[data_sorted > fB[h] & data_sorted <= fB[h + 1L]]
Nh <- length(sl)
if (Nh < 2L) next
obj <- obj + (Nh / N_pop) * sqrt(stats::var(sl))
}
obj
}
obj_dens <- function(b_test) {
fB <- c(x_lo, sort(b_test), x_hi)
Hh <- length(b_test) + 1L
tot <- sum(diff(x_seq) * (head(dens, -1) + tail(dens, -1))) / 2
Wh <- numeric(Hh); Sh <- numeric(Hh)
for (h in seq_len(Hh)) {
idx <- x_seq >= fB[h] & x_seq < fB[h + 1L]
if (sum(idx) < 2L) next
xs <- x_seq[idx]; fs <- dens[idx]; dx <- diff(xs)
Wh[h] <- sum(dx * (head(fs, -1) + tail(fs, -1))) / 2 / max(tot, 1e-12)
if (Wh[h] < 1e-12) next
mu <- sum(dx * (head(xs * fs, -1) + tail(xs * fs, -1))) /
2 / Wh[h] / max(tot, 1e-12)
mu2 <- sum(dx * (head(xs^2 * fs, -1) + tail(xs^2 * fs, -1))) /
2 / Wh[h] / max(tot, 1e-12)
Sh[h] <- sqrt(max(mu2 - mu^2, 0))
}
sum(Wh * Sh)
}
obj_fn <- if (!is.null(data_sorted)) obj_data else obj_dens
# ── Precompute slider grids + step labels ────────────────────────────────────
b_grids <- vector("list", n_sliders)
steps_list <- vector("list", n_sliders)
for (bi in seq_len(n_sliders)) {
b_lo_i <- full_B[bi] + 0.02 * x_rng
b_hi_i <- full_B[bi + 2L] - 0.02 * x_rng
if (b_lo_i >= b_hi_i) b_lo_i <- full_B[bi] + 1e-6
b_grid <- seq(b_lo_i, b_hi_i, length.out = n_steps)
b_grids[[bi]] <- b_grid
steps_list[[bi]] <- lapply(seq_len(n_steps), function(si) {
b_test <- sort(replace(internal_OSB, bi, b_grid[si]))
obj_val <- tryCatch(obj_fn(b_test), error = function(e) NA_real_)
if (!is.finite(obj_val)) obj_val <- 0
# annotations[1] = dynamic stats box (index 0=header, 1=stats)
list(
label = "", # hide rail tick labels; currentvalue prefix shows the value
method = "update",
args = list(
list(),
list(
`annotations[1].text` = sprintf(
"<b>b<sub>%d</sub> = %.4f</b><br>ΣW<sub>h</sub>S<sub>h</sub> = %.5f<br>V(n) = %.6f",
bi, b_grid[si], obj_val, obj_val^2 / n_total)
)
)
)
})
}
# ── Slider layout sizing (bottom grey-box, full width) ───────────────────────
px_per_slider <- 110L # px allocated per slider
margin_b_px <- 100L + n_sliders * px_per_slider # total bottom margin (100px base for x-axis title clearance)
h_px <- max(580L, 440L + n_sliders * px_per_slider) # 440 = 355 plot + 85 top, keeps plot area fixed despite larger margin base
step_paper <- px_per_slider / h_px # paper-unit height per slider
bottom_y <- 82.0 / h_px # gap at very bottom (clear x-axis label)
sub_h <- c("₁","₂","₃","₄","₅",
"₆","₇","₈","₉")
slider_specs <- lapply(seq_len(n_sliders), function(bi) {
opt_step <- which.min(abs(b_grids[[bi]] - OSB[bi])) - 1L
y_pos <- bottom_y + (n_sliders - bi) * step_paper
list(
active = opt_step,
currentvalue = list(
prefix = paste0("b", sub_h[min(bi, 9L)], " = "),
font = list(size = 11, color = "#C0392B"),
visible = TRUE,
offset = 8
),
x = 0.10,
y = y_pos,
len = 0.90,
xanchor = "left",
yanchor = "top",
pad = list(t = 24, b = 4),
ticklen = 0,
steps = steps_list[[bi]]
)
})
# ── Build base figure ────────────────────────────────────────────────────────
fig <- plotly::plot_ly()
if (!is.null(pop_data) && length(pop_data) >= 10L) {
hh <- hist(pop_data, breaks = "Sturges", plot = FALSE)
fig <- fig |>
plotly::add_bars(
x = hh$mids, y = hh$density, name = "Histogram",
marker = list(color = "rgba(91,158,214,0.40)",
line = list(color = "rgba(91,158,214,0.75)", width = 0.7)),
width = diff(hh$breaks)[1],
hovertemplate = "Midpoint: %{x:.2f}<br>Density: %{y:.5f}<extra></extra>"
)
}
for (h in seq_len(H)) {
lo_h <- full_B[h]; hi_h <- full_B[h + 1L]
idx <- x_seq >= lo_h & x_seq <= hi_h
if (sum(idx) < 2L) next
wh <- res$Wh[h]
sh <- round(sqrt(max(res$Vh[h], 0)), 4)
nh <- res$nh[h]
fig <- fig |>
plotly::add_trace(
x = c(x_seq[idx], rev(x_seq[idx])),
y = c(dens[idx], rep(0, sum(idx))),
type = "scatter", mode = "lines", fill = "toself",
fillcolor = grDevices::adjustcolor(cols[h], alpha.f = 0.38),
line = list(color = "transparent"),
name = sprintf("Stratum %d", h),
hovertemplate = sprintf(
"<b>Stratum %d</b><br>Wₕ=%.4f Sₕ=%.4f nₕ=%d<extra></extra>",
h, wh, sh, nh)
)
}
distr_lbl <- if (!is.null(res$distr) && nzchar(res$distr)) res$distr else "KDE"
fig <- fig |>
plotly::add_trace(
x = x_seq, y = dens, type = "scatter", mode = "lines",
line = list(color = "#1a1a2e", width = 2.6),
name = sprintf("f(y) — %s", distr_lbl),
hovertemplate = "y=%{x:.3f} f(y)=%{y:.5f}<extra></extra>"
)
for (i in seq_along(OSB)) {
fig <- fig |>
plotly::add_trace(
x = rep(OSB[i], 2), y = c(0, d_max * 1.08),
type = "scatter", mode = "lines",
line = list(color = "#C0392B", width = 2, dash = "dash"),
showlegend = FALSE, hoverinfo = "none"
)
}
# ── Annotations — ORDER MATTERS: index 0 = header, 1 = stats (slider target)
ibe_header <- list(
x = 0.655, y = 0.99, xref = "paper", yref = "paper",
xanchor = "left", yanchor = "top",
text = sprintf(
"<b style='font-size:12px'>⚙ Boundary Explorer</b><br><span style='color:#555'>%d strata │ n = %d</span>",
H, n_total),
showarrow = FALSE,
font = list(size = 11.5, color = "#1A2B3C"),
bgcolor = "#DDE8FA",
bordercolor = "#3B5FAE",
borderwidth = 2,
borderpad = 9
)
# Stats box goes BELOW the legend — estimate legend height from item count
n_legend_items <- H + 2L # Histogram + H strata + density line
legend_h_paper <- (n_legend_items * 26L + 20L) / h_px # realistic row height
stats_y <- max(0.05, 0.78 - legend_h_paper - 0.10) # clear gap below legend
init_txt <- sprintf(
"<b>Optimal boundaries</b><br>ΣW<sub>h</sub>S<sub>h</sub> = %.5f<br>V(n) = %.6f",
res$WhShTot, res$WhShTot^2 / n_total)
dyn_stats <- list(
x = 0.655, y = stats_y, xref = "paper", yref = "paper",
xanchor = "left", yanchor = "top",
text = init_txt, showarrow = FALSE,
font = list(size = 11, color = "#C0392B"),
bgcolor = "rgba(255,240,240,0.92)",
bordercolor = "#C0392B", borderwidth = 1.2, borderpad = 8
)
sub_h2 <- c("₁","₂","₃","₄","₅",
"₆","₇","₈","₉")
osb_annots <- lapply(seq_along(OSB), function(i) {
list(x = OSB[i], y = d_max * 0.52,
xref = "x", yref = "y",
xanchor = "center", yanchor = "middle",
text = sprintf("b%s = %.3f", sub_h2[min(i, 9L)], OSB[i]),
textangle = -90,
showarrow = FALSE,
font = list(size = 10, color = "#C0392B"),
bgcolor = "rgba(255,255,255,0.60)",
borderpad = 2)
})
# annotation[0]=header, [1]=dyn_stats (slider updates this), [2+]=osb labels
all_annots <- c(list(ibe_header, dyn_stats), osb_annots)
ttl <- { t <- plot_main(); if (!is.null(t) && nzchar(t)) t else
sprintf("%d-Strata Boundary Explorer", H) }
# Grey box behind slider area
margin_b_paper <- margin_b_px / h_px
slider_box <- list(
type = "rect", layer = "below",
xref = "paper", yref = "paper",
x0 = -0.05, x1 = 1.05,
y0 = -0.005, y1 = margin_b_paper + 0.005,
fillcolor = "#F3F4F6",
line = list(color = "#D1D5DB", width = 1)
)
fig |> plotly::layout(
height = h_px,
title = list(text = ttl,
font = list(size = 15, color = "#1A2B3C"),
x = 0.31, xanchor = "center"),
xaxis = list(domain = c(0, 0.62),
title = list(text = "Y (variable value)", font = list(size = 14)),
tickfont = list(size = 12),
zeroline = FALSE, gridcolor = "#EEEEEE"),
yaxis = list(title = list(text = "Density f(y)", font = list(size = 14)),
tickfont = list(size = 12),
zeroline = FALSE, gridcolor = "#EEEEEE",
range = c(0, d_max * 1.15)),
sliders = slider_specs,
shapes = list(slider_box),
legend = list(orientation = "v",
x = 0.655, y = 0.78,
xanchor = "left", yanchor = "top",
bgcolor = "rgba(255,255,255,0.85)",
bordercolor = "#DDDDDD", borderwidth = 1,
font = list(size = 11)),
paper_bgcolor = "white",
plot_bgcolor = "white",
hovermode = "x unified",
margin = list(t = 65, b = margin_b_px, l = 95, r = 20),
annotations = all_annots
)
} # end if(FALSE) placeholder
# ── Generate R Code ───────────────────────────────────────────────────────────
.build_r_code <- function(input, res) {
h <- as.integer(input$h_strata)
n <- as.integer(input$n_sample)
mtd <- input$solver
header <- paste0(
"# ================================================================\n",
"# stratifyR 2.0 — Reproducible R Code\n",
"# Generated: ", format(Sys.time(), "%Y-%m-%d %H:%M"), "\n",
"# ================================================================\n\n",
"library(stratifyR)\n"
)
if (input$data_source == "distr") {
# ---- distribution pathway -----------------------------------------------
pn <- distr_params[[input$distr_name]]$names
pv <- vapply(pn, function(p) {
v <- input[[paste0("dp_", p)]]
if (is.null(v)) 0 else as.numeric(v)
}, numeric(1))
p_lines <- paste(
sprintf(" %-10s = %s", pn, formatC(pv, format = "g", digits = 6)),
collapse = ",\n")
data_call <- sprintf(
"\nres <- strata.distr(\n h = %dL,\n initval = %s,\n dist = %s,\n distr = \"%s\",\n params = c(\n%s\n ),\n n = %dL,\n N = %dL,\n method = \"%s\"\n)",
h,
formatC(as.numeric(input$initval), format = "g", digits = 6),
formatC(as.numeric(input$dist_range), format = "g", digits = 6),
input$distr_name,
p_lines,
n,
as.integer(input$distr_N),
mtd)
} else if (input$data_source == "builtin") {
# ---- built-in dataset pathway -------------------------------------------
parts <- strsplit(input$builtin_ds, "\\|")[[1]]
ds_name <- parts[1]; col_name <- parts[2]
N_arg <- if (!is.null(res) && !is.null(res$N) && res$N > 0 && res$N != res$nhTot)
sprintf(",\n N = %dL", as.integer(res$N)) else ""
data_call <- sprintf(
"\n# Load built-in data\ndata(\"%s\", package = \"stratifyR\")\ny <- %s$%s\n\nres <- strata.data(\n data = y,\n h = %dL,\n n = %dL%s,\n method = \"%s\"\n)",
ds_name, ds_name, col_name,
h, n, N_arg, mtd)
} else {
# ---- uploaded file pathway -----------------------------------------------
col_name <- if (!is.null(input$data_col)) input$data_col else "your_column"
N_arg <- if (isTRUE(input$use_N) && !is.null(input$pop_N))
sprintf(",\n N = %dL", as.integer(input$pop_N)) else ""
data_call <- sprintf(
"\n# Replace the path with your actual file\ny_df <- read.csv(\"your_data.csv\", stringsAsFactors = FALSE)\ny <- y_df$%s\n\nres <- strata.data(\n data = y,\n h = %dL,\n n = %dL%s,\n method = \"%s\"\n)",
col_name, h, n, N_arg, mtd)
}
# ---- result annotation (if a run has been completed) --------------------
result_block <- ""
if (!is.null(res)) {
x_lo <- res$maxval * res$initval
x_hi <- max(res$maxval * res$finval, max(res$OSB))
int_osb <- res$OSB[res$OSB > x_lo & res$OSB < x_hi]
osb_lines <- paste(
sprintf("# b%d = %.4f", seq_along(int_osb), int_osb),
collapse = "\n")
Nh_lines <- paste(
sprintf("# Stratum %d: nh = %d", seq_along(res$nh), res$nh),
collapse = "\n")
conv_txt <- if (isTRUE(res$converged)) "YES" else
if (is.null(res$converged)) "N/A (DP)" else "NO"
result_block <- paste0(
"\n\n",
"# ── Result summary (from app run) ────────────────────────\n",
sprintf("# Method : %s\n", toupper(mtd)),
sprintf("# H (strata) : %d\n", nrow(res$h)),
sprintf("# Total n : %d\n", res$nhTot),
sprintf("# ΣWhSh : %.6f\n", res$WhShTot),
sprintf("# V*(n) : %.8f\n", res$WhShTot^2 / res$nhTot),
"# Converged : ", conv_txt, "\n",
if (!is.na(res$optimality_gap) && is.finite(res$optimality_gap))
sprintf("# Optimality gap : %.4f%%\n", 100 * res$optimality_gap) else "",
"#\n",
"# Optimal strata boundaries (OSB):\n",
osb_lines, "\n",
"#\n",
"# Neyman sample allocation:\n",
Nh_lines, "\n",
"# ─────────────────────────────────────────────────────────"
)
}
# ---- visualisation calls ------------------------------------------------
viz_block <- paste0(
"\n\n",
"summary(res)\n\n",
"plot(res) # 2D density + strata\n",
"plot(res, type = \"3d\") # 3D surface\n",
"plot(res, type = \"interactive\") # boundary explorer\n"
)
paste0(header, data_call, result_block, viz_block)
}
observeEvent(input$gen_code_btn, {
res <- strata_res()
code <- .build_r_code(input, res)
showModal(modalDialog(
title = tagList(
icon("code", style = "color:#1D4ED8; margin-right:6px;"),
"Generated R Code"
),
size = "l",
easyClose = TRUE,
footer = tagList(
tags$button(
id = "copy_r_code",
class = "btn btn-primary r-code-copy-btn",
onclick = paste0(
"var el=document.getElementById('r_code_pre');",
"navigator.clipboard.writeText(el.innerText).then(function(){",
" document.getElementById('copy_r_code').innerHTML='✓ Copied!';",
" setTimeout(function(){",
" document.getElementById('copy_r_code').innerHTML='📋 Copy to clipboard';",
" },2000);",
"});"
),
HTML("📋 Copy to clipboard")
),
modalButton("Close")
),
tags$pre(
id = "r_code_pre",
class = "r-code-block",
code
)
))
})
# ── CSV download ───────────────────────────────────────────────────────────────
output$dl_results <- downloadHandler(
filename = function()
paste0("stratifyR_", Sys.Date(), "_h", input$h_strata, "_n", input$n_sample, ".csv"),
content = function(file) {
df <- bounds_df()
utils::write.csv(df, file, row.names = FALSE)
}
)
}
# ==============================================================================
shinyApp(ui, server, options = list(launch.browser = TRUE))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.