Windowing Strategy for the Additive p-Control Chart

knitr::opts_chunk$set(fig.width=16, fig.height=8, fig.align="center", out.width="100%", dpi=350)
required_pkgs <- c("readxl", "dplyr", "tidyr", "ggplot2", "knitr")
optional_pkgs <- c("magick", "writexl")
install_if_missing <- function(pkgs){for(pkg in pkgs){if(!requireNamespace(pkg, quietly=TRUE)){install.packages(pkg, dependencies=TRUE)}; suppressPackageStartupMessages(library(pkg, character.only=TRUE))}}
install_if_missing(required_pkgs)
for(pkg in optional_pkgs){if(requireNamespace(pkg, quietly=TRUE)){suppressPackageStartupMessages(library(pkg, character.only=TRUE))}}
XLSX_PATH <- "GW_Comparative.xlsx"
OUT_DIR <- "outputs_GW_RHTML"
dir.create(OUT_DIR, showWarnings=FALSE, recursive=TRUE)
IMG_WIDTH_PX <- 5200; IMG_HEIGHT_PX <- 1900; IMG_DPI <- 360
IMG_FACET_W <- 18; IMG_FACET_H <- 12; IMG_BAR_W <- 16; IMG_BAR_H <- 7
N_VENT <- 12
COL_PI <- "black"; COL_CONTROL <- "darkblue"; COL_CENTER <- "gray40"; COL_SPEC <- "darkgreen"; COL_OOC_CTRL <- "red"; COL_OOC_SPEC <- "purple4"
if(!file.exists(XLSX_PATH)){stop(paste0("The file '", XLSX_PATH, "' was not found. Place GW_Comparative.xlsx in the same directory as this RHTML file or modify XLSX_PATH."))}
sheets_available <- readxl::excel_sheets(XLSX_PATH)
knitr::kable(data.frame(Available_sheets=sheets_available), caption="Sheets found in the Excel file")

validate_columns <- function(df, required_cols, sheet_name){
  missing_cols <- setdiff(required_cols, names(df))
  if(length(missing_cols)>0){stop(paste0("In sheet '", sheet_name, "' the following columns are missing: ", paste(missing_cols, collapse=", "), "\nAvailable columns: ", paste(names(df), collapse=", ")))}
  invisible(TRUE)
}
safe_numeric <- function(x){suppressWarnings(as.numeric(x))}

plot_control_from_cols <- function(x, pi, ucl, cl, lcl, out_png=NULL, width_px=IMG_WIDTH_PX, height_px=IMG_HEIGHT_PX, res_dpi=IMG_DPI){
  ok <- !(is.na(pi) & is.na(ucl) & is.na(cl) & is.na(lcl)); n_residuals <- sum(!ok)
  x <- x[ok]; pi <- pi[ok]; ucl <- ucl[ok]; cl <- cl[ok]; lcl <- lcl[ok]
  n <- length(x); if(n<2) stop("There are not enough data points to plot the general chart.")
  y_all <- c(pi,ucl,cl,lcl); y_min <- min(y_all,na.rm=TRUE); y_max <- max(y_all,na.rm=TRUE); pad <- .02*(y_max-y_min); if(!is.finite(pad)||pad==0) pad <- 1e-6; ylim <- c(y_min-pad,y_max+pad)
  ooc_up <- which(pi>ucl); ooc_down <- which(pi<lcl); ooc_all <- sort(unique(c(ooc_up,ooc_down)))
  if(!is.null(out_png)){png(out_png,width=width_px,height=height_px,res=res_dpi); on.exit(dev.off(), add=TRUE)}
  op <- par(no.readonly=TRUE); on.exit(par(op), add=TRUE)
  par(mar=c(6,6.2,3.5,2.5),xaxs="i",yaxs="i",cex=1.15,cex.axis=1.05,cex.lab=1.15)
  plot(x, pi, type="l", col=COL_PI, lwd=1.2, ylim=ylim, xaxt="n", xlab="Weeks", ylab=expression(p[i]), main="Additive - General adjusted p-chart")
  axis(1, at=c(1,n), labels=c("1",n)); lines(x,cl,col=COL_CENTER,lwd=1.4,lty=1); lines(x,ucl,col=COL_CONTROL,lwd=1.2,lty=3); lines(x,lcl,col=COL_CONTROL,lwd=1.2,lty=3)
  if(length(ooc_all)>0){points(x[ooc_all], pi[ooc_all], pch=4, col=COL_OOC_CTRL, cex=.65, lwd=1.3)}
  box(); legend("bottom", inset=-.30, xpd=TRUE, horiz=TRUE, bty="n", legend=c("pi","UCL/LCS","CL/LC","LCL/LCI","OOC-Control"), col=c(COL_PI,COL_CONTROL,COL_CENTER,COL_CONTROL,COL_OOC_CTRL), lty=c(1,3,1,3,NA), lwd=c(1.2,1.2,1.4,1.2,NA), pch=c(NA,NA,NA,NA,4), pt.cex=.75)
  invisible(list(total=length(ooc_all), above_ucl=length(ooc_up), below_lcl=length(ooc_down), idx=ooc_all, residuals=n_residuals, n=n))
}

plot_control_and_spec_from_cols <- function(x, pi, ucl, ccl, lcl, usl, csl, lsl, out_png=NULL, title=NULL, width_px=IMG_WIDTH_PX, height_px=IMG_HEIGHT_PX, res_dpi=IMG_DPI){
  ok <- !(is.na(pi)&is.na(ucl)&is.na(ccl)&is.na(lcl)&is.na(usl)&is.na(csl)&is.na(lsl)); n_residuals <- sum(!ok)
  x <- x[ok]; pi <- pi[ok]; ucl <- ucl[ok]; ccl <- ccl[ok]; lcl <- lcl[ok]; usl <- usl[ok]; csl <- csl[ok]; lsl <- lsl[ok]
  n <- length(x); if(n<2) stop("There are not enough data points to plot the window.")
  y_all <- c(pi,ucl,ccl,lcl,usl,csl,lsl); y_min <- min(y_all,na.rm=TRUE); y_max <- max(y_all,na.rm=TRUE); pad <- .02*(y_max-y_min); if(!is.finite(pad)||pad==0) pad <- 1e-6; ylim <- c(y_min-pad,y_max+pad)
  ooc_ctrl_up <- which(pi>ucl); ooc_ctrl_down <- which(pi<lcl); ooc_ctrl_all <- sort(unique(c(ooc_ctrl_up,ooc_ctrl_down)))
  ooc_spec_up <- which(pi>usl); ooc_spec_down <- which(pi<lsl); ooc_spec_all <- sort(unique(c(ooc_spec_up,ooc_spec_down)))
  if(!is.null(out_png)){png(out_png,width=width_px,height=height_px,res=res_dpi); on.exit(dev.off(), add=TRUE)}
  op <- par(no.readonly=TRUE); on.exit(par(op), add=TRUE); layout(matrix(c(1,2), nrow=1), widths=c(4.8,1.8))
  par(mar=c(6,6.2,3.5,1.5), xaxs="i", yaxs="i", cex=1.15, cex.axis=1.05, cex.lab=1.15)
  plot(x,pi,type="l",col=COL_PI,lwd=1.2,ylim=ylim,xaxt="n",xlab="Weeks",ylab=expression(p[i]),main=ifelse(is.null(title),"Window - Control and specification",title)); axis(1,at=c(1,n),labels=c("1",n))
  lines(x,ccl,col=COL_CENTER,lwd=1.4,lty=1); lines(x,ucl,col=COL_CONTROL,lwd=1.2,lty=3); lines(x,lcl,col=COL_CONTROL,lwd=1.2,lty=3)
  lines(x,csl,col=COL_SPEC,lwd=1.3,lty=1); lines(x,usl,col=COL_SPEC,lwd=1.2,lty=2); lines(x,lsl,col=COL_SPEC,lwd=1.2,lty=2)
  if(length(ooc_ctrl_all)>0){points(x[ooc_ctrl_all],pi[ooc_ctrl_all],pch=4,col=COL_OOC_CTRL,cex=.70,lwd=1.4)}
  if(length(ooc_spec_all)>0){points(x[ooc_spec_all],pi[ooc_spec_all],pch=17,col=COL_OOC_SPEC,cex=.45)}
  box(); par(mar=c(2,.5,2,.5)); plot.new(); legend("topleft",bty="n",cex=.98,y.intersp=1.15,legend=c("pi","UCL","CCL","LCL","USL","CEL/CSL","LEL/LSL","OOC-Control","OOC-Specification"), col=c(COL_PI,COL_CONTROL,COL_CENTER,COL_CONTROL,COL_SPEC,COL_SPEC,COL_SPEC,COL_OOC_CTRL,COL_OOC_SPEC), lty=c(1,3,1,3,2,1,2,NA,NA), lwd=c(1.2,1.2,1.4,1.2,1.2,1.3,1.2,NA,NA), pch=c(NA,NA,NA,NA,NA,NA,NA,4,17), pt.cex=c(NA,NA,NA,NA,NA,NA,NA,.90,.55))
  invisible(list(ooc_control_total=length(ooc_ctrl_all), ooc_control_above=length(ooc_ctrl_up), ooc_control_below=length(ooc_ctrl_down), ooc_control_idx=ooc_ctrl_all, ooc_spec_total=length(ooc_spec_all), ooc_spec_above=length(ooc_spec_up), ooc_spec_below=length(ooc_spec_down), ooc_spec_idx=ooc_spec_all, residuals=n_residuals, n=n))
}

if(!"a1" %in% sheets_available){stop("The Excel file does not contain the required sheet 'a1'.")}
df_gw <- readxl::read_excel(XLSX_PATH, sheet="a1"); validate_columns(df_gw, c("Pi","LCS","LC","LCI"), "a1")
GW_PNG <- file.path(OUT_DIR,"Goedhart_Woodall_chart.png")
res_gw <- plot_control_from_cols(x=seq_len(nrow(df_gw)), pi=safe_numeric(df_gw$Pi), ucl=safe_numeric(df_gw$LCS), cl=safe_numeric(df_gw$LC), lcl=safe_numeric(df_gw$LCI), out_png=GW_PNG)
knitr::include_graphics(GW_PNG)

gw_summary <- data.frame(Chart="General additive chart", Observations=res_gw$n, OOC_Control_Total=res_gw$total, Above_UCL_LCS=res_gw$above_ucl, Below_LCL_LCI=res_gw$below_lcl, Residual_Rows_Removed=res_gw$residuals)
knitr::kable(gw_summary, caption="Summary of OOC signals for the general additive chart")

vent_results <- vector("list", N_VENT); names(vent_results) <- paste0("Vent",1:N_VENT)
vent_data_all <- list(); vent_png_files <- character(N_VENT)
for(k in 1:N_VENT){
  sheet_name <- paste0("Vent",k); if(!sheet_name %in% sheets_available){stop(paste0("Missing sheet: ", sheet_name))}
  df_vk <- readxl::read_excel(XLSX_PATH, sheet=sheet_name); validate_columns(df_vk, c("Pi","UCL","CCL","LCL","USL","CEL","LEL"), sheet_name)
  vent_png_files[k] <- file.path(OUT_DIR, sprintf("Vent%02d_Control_Spec_chart.png", k))
  vent_results[[k]] <- plot_control_and_spec_from_cols(x=seq_len(nrow(df_vk)), pi=safe_numeric(df_vk$Pi), ucl=safe_numeric(df_vk$UCL), ccl=safe_numeric(df_vk$CCL), lcl=safe_numeric(df_vk$LCL), usl=safe_numeric(df_vk$USL), csl=safe_numeric(df_vk$CEL), lsl=safe_numeric(df_vk$LEL), out_png=vent_png_files[k], title=paste0("Window ", k, " - Control and specification"))
  tmp <- data.frame(Window=sheet_name, Index=seq_len(nrow(df_vk)), Pi=safe_numeric(df_vk$Pi), UCL=safe_numeric(df_vk$UCL), CCL=safe_numeric(df_vk$CCL), LCL=safe_numeric(df_vk$LCL), USL=safe_numeric(df_vk$USL), CEL=safe_numeric(df_vk$CEL), LEL=safe_numeric(df_vk$LEL))
  tmp <- tmp |> dplyr::mutate(OOC_Control=(Pi>UCL)|(Pi<LCL), OOC_Specification=(Pi>USL)|(Pi<LEL))
  vent_data_all[[k]] <- tmp
}
vent_summary <- dplyr::bind_rows(lapply(seq_len(N_VENT), function(k){r <- vent_results[[k]]; data.frame(Window=paste0("Vent",k), Observations=r$n, OOC_Control_Total=r$ooc_control_total, Above_UCL=r$ooc_control_above, Below_LCL=r$ooc_control_below, OOC_Spec_Total=r$ooc_spec_total, Above_USL=r$ooc_spec_above, Below_LSL=r$ooc_spec_below, Residual_Rows_Removed=r$residuals)}))
vent_data_all_df <- dplyr::bind_rows(vent_data_all)
knitr::kable(vent_summary, caption="Summary of signals by window: control and specification")

vent_long_limits <- vent_data_all_df |> dplyr::select(Window,Index,UCL,CCL,LCL,USL,CEL,LEL) |> tidyr::pivot_longer(cols=c(UCL,CCL,LCL,USL,CEL,LEL), names_to="Limit", values_to="Value")
p_faceted <- ggplot() + geom_line(data=vent_data_all_df, aes(x=Index,y=Pi), color="black", linewidth=.35) + geom_line(data=vent_long_limits, aes(x=Index,y=Value,linetype=Limit,color=Limit), linewidth=.32) + geom_point(data=vent_data_all_df |> dplyr::filter(OOC_Control), aes(x=Index,y=Pi), color=COL_OOC_CTRL, shape=4, size=.9, stroke=.7) + geom_point(data=vent_data_all_df |> dplyr::filter(OOC_Specification), aes(x=Index,y=Pi), color=COL_OOC_SPEC, shape=17, size=.8, alpha=.8) + facet_wrap(~Window, scales="free_x", ncol=3) + scale_color_manual(values=c("UCL"=COL_CONTROL,"CCL"=COL_CENTER,"LCL"=COL_CONTROL,"USL"=COL_SPEC,"CEL"=COL_SPEC,"LEL"=COL_SPEC)) + scale_linetype_manual(values=c("UCL"="dotted","CCL"="solid","LCL"="dotted","USL"="dashed","CEL"="solid","LEL"="dashed")) + labs(title="Sequential windowing: p-charts with control and specification limits", subtitle="Red crosses: OOC-Control; purple triangles: OOC-Specification", x="Local index within the window", y=expression(p[i]), color="Limits", linetype="Limits") + theme_classic(base_size=12) + theme(plot.title=element_text(face="bold",size=16), strip.background=element_rect(fill="#eef3f8",color="#cbd5e1"), strip.text=element_text(face="bold"), legend.position="bottom")
print(p_faceted); FACETED_PNG <- file.path(OUT_DIR,"ALL_WINDOWS_faceted_control_spec.png"); ggsave(FACETED_PNG,p_faceted,width=IMG_FACET_W,height=IMG_FACET_H,dpi=IMG_DPI)

signal_long <- vent_summary |> dplyr::select(Window,OOC_Control_Total,OOC_Spec_Total) |> tidyr::pivot_longer(cols=c(OOC_Control_Total,OOC_Spec_Total), names_to="Signal_Type", values_to="Signals") |> dplyr::mutate(Signal_Type=dplyr::recode(Signal_Type, OOC_Control_Total="OOC-Control", OOC_Spec_Total="OOC-Specification"), Window=factor(Window, levels=paste0("Vent",1:N_VENT)))
p_signal <- ggplot(signal_long, aes(x=Window,y=Signals,fill=Signal_Type)) + geom_col(position="dodge", width=.72) + labs(title="Comparison of signals by window", subtitle="Differentiation between control-limit and specification-limit signals", x="Window", y="Number of signals", fill="Signal type") + theme_classic(base_size=13) + theme(plot.title=element_text(face="bold"), axis.text.x=element_text(angle=45,hjust=1), legend.position="bottom")
print(p_signal); SIGNALS_PNG <- file.path(OUT_DIR,"Signals_by_Window_Control_vs_Spec.png"); ggsave(SIGNALS_PNG,p_signal,width=IMG_BAR_W,height=IMG_BAR_H,dpi=IMG_DPI)