Skip to content

Process Capability: Adjust parameter estimation precision settings #462

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JTPetter wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Process Capability: Adjust parameter estimation precision settings #462

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ce3e67f
adjust precision settings weibull and fix tests
JTPetter Mar 31, 2026
53adf41
fix typo
JTPetter Mar 31, 2026
a1dd302
Adjust prob. table to use common dist. function
JTPetter Apr 13, 2026
2f640e1
Remove other precision settings where applicable
JTPetter Apr 17, 2026
cf5e770
New unit test
JTPetter Apr 17, 2026
fec976b
update test plots
JTPetter Jun 18, 2026
2848874
Copilot feedback
JTPetter Jun 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 0 additions & 5 deletions R/commonQualityControl.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1456,12 +1456,8 @@ KnownControlStats.RS <- function(N, sigma = 3) {
theta <- fix.arg[["scale"]] # scale
} else {
fitWeibull <- try(fitdistrplus::fitdist(data, "weibull", method = "mle",
control = list(maxit = 10000),
fix.arg = fix.arg))




if (jaspBase::isTryError(fitWeibull))
stop(estimationErrorMessage, call. = FALSE)

Expand Down Expand Up @@ -1510,7 +1506,6 @@ KnownControlStats.RS <- function(N, sigma = 3) {

# Estimate parameters using fitdistrplus, because it can keep values fixed
lnorm3Fit <- try(fitdistrplus::fitdist(data, "lnorm3Temp", method = "mle",
control = list(maxit = 10000),
start = lnorm3startList,
fix.arg = fix.arg))
if (jaspBase::isTryError(lnorm3Fit))
Expand Down
101 changes: 52 additions & 49 deletions R/processCapabilityStudies.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1131,20 +1131,20 @@ processCapabilityStudies <- function(jaspResults, dataset, options) {
if (options[["target"]] && options[["processCapabilityPlotSpecificationLimits"]])
xLimits <- range(xLimits, options[["targetValue"]])

# Addition to consider that if distributions are set historically, they may fall outside the usual limits
if (distribution == "normal" && options[["historicalMean"]])
xLimits <- range(xLimits, options[["historicalMeanValue"]] - 1.5 * sdo, options[["historicalMeanValue"]] + 1.5 * sdo)
if (distribution == "weibull" || distribution == "3ParameterWeibull" && options[["historicalScale"]])
xLimits <- range(xLimits, options[["historicalScaleValue"]] - 1.5 * sdo, options[["historicalScaleValue"]] + 1.5 * sdo)
if (distribution == "lognormal" || distribution == "3ParameterLognormal" && options[["historicalLogMean"]])
xLimits <- range(xLimits, exp(options[["historicalLogMeanValue"]]) - 1.5 * sdo,
exp(options[["historicalLogMeanValue"]]) + 1.5 * sdo)
if (distribution == "exponential")
xLimits <- range(xLimits, 0)
if (distribution == "gamma" && options[["historicalShape"]] && options[["historicalScale"]])
xLimits <- range(xLimits, options[["historicalScaleValue"]] * (options[["historicalShapeValue"]] - 1))
if ((distribution == "logistic" || distribution == "loglogistic") && options[["historicalLocation"]])
xLimits <- range(xLimits, options[["historicalLocationValue"]])
# Addition to consider that if distributions are set historically, they may fall outside the usual limits
if (distribution == "normal" && options[["historicalMean"]])
xLimits <- range(xLimits, options[["historicalMeanValue"]] - 1.5 * sdo, options[["historicalMeanValue"]] + 1.5 * sdo)
if ((distribution == "weibull" || distribution == "3ParameterWeibull") && options[["historicalScale"]])
xLimits <- range(xLimits, options[["historicalScaleValue"]] - 1.5 * sdo, options[["historicalScaleValue"]] + 1.5 * sdo)
if ((distribution == "lognormal" || distribution == "3ParameterLognormal") && options[["historicalLogMean"]])
xLimits <- range(xLimits, exp(options[["historicalLogMeanValue"]]) - 1.5 * sdo,
exp(options[["historicalLogMeanValue"]]) + 1.5 * sdo)
if (distribution == "exponential")
xLimits <- range(xLimits, 0)
if (distribution == "gamma" && options[["historicalShape"]] && options[["historicalScale"]])
xLimits <- range(xLimits, options[["historicalScaleValue"]] * (options[["historicalShapeValue"]] - 1))
if ((distribution == "logistic" || distribution == "loglogistic") && options[["historicalLocation"]])
xLimits <- range(xLimits, options[["historicalLocationValue"]])

xBreaks <- jaspGraphs::getPrettyAxisBreaks(c(allData))
xStep <- diff(xBreaks)[1]
Expand Down Expand Up @@ -3093,16 +3093,16 @@ processCapabilityStudies <- function(jaspResults, dataset, options) {
}
table$addColumnInfo(name = "n", title = gettext("N"), type = "integer")
if (options[["nullDistribution"]] == "normal") {
table$addColumnInfo(name = "mean", title = gettextf("Mean (%1$s)", "\u03BC"), type = "number")
table$addColumnInfo(name = "sd", title = gettextf("Std. dev. (%1$s)", "\u03C3"), type = "number")
table$addColumnInfo(name = "parameterEstimate1", title = gettextf("Mean (%1$s)", "\u03BC"), type = "number")
table$addColumnInfo(name = "parameterEstimate2", title = gettextf("Std. dev. (%1$s)", "\u03C3"), type = "number")
} else if (options[["nullDistribution"]] == "lognormal") {
table$addColumnInfo(name = "mean", title = gettextf("Log mean (%1$s)", "\u03BC"), type = "number")
table$addColumnInfo(name = "sd", title = gettextf("Log std.dev (%1$s)", "\u03C3"), type = "number")
table$addColumnInfo(name = "parameterEstimate1", title = gettextf("Log mean (%1$s)", "\u03BC"), type = "number")
table$addColumnInfo(name = "parameterEstimate2", title = gettextf("Log std.dev (%1$s)", "\u03C3"), type = "number")
} else if (options[["nullDistribution"]] == "weibull") {
table$addColumnInfo(name = "mean", title = gettextf("Shape (%1$s)", "\u03B2"), type = "number")
table$addColumnInfo(name = "sd", title = gettextf("Scale (%1$s)", "\u03B8"), type = "number")
table$addColumnInfo(name = "parameterEstimate1", title = gettextf("Shape (%1$s)", "\u03B2"), type = "number")
table$addColumnInfo(name = "parameterEstimate2", title = gettextf("Scale (%1$s)", "\u03B8"), type = "number")
}
table$addColumnInfo(name = "ad", title = gettext("AD"), type = "integer")
table$addColumnInfo(name = "ad", title = gettext("AD"), type = "number")
table$addColumnInfo(name = "p", title = gettext("<i>p</i>-value"), type = "pvalue")
table$addFootnote(gettextf("The Anderson-Darling statistic A<i>D</i> is calculated against the %1$s distribution.", distributionTitle))
table$addFootnote(gettextf("Red dotted lines in the probability plot below represent a 95%% confidence interval."))
Expand All @@ -3115,47 +3115,50 @@ processCapabilityStudies <- function(jaspResults, dataset, options) {
container[["probabilityTable"]] <- table
return()
}
tableColNames <- c("mean", "sd", "n", "ad", "p")
tableColNames <- c("parameterEstimate1", "parameterEstimate2", "n", "ad", "p")
if (nStages > 1)
tableColNames <- c("stage", tableColNames)
tableDf <- data.frame(matrix(ncol = length(tableColNames), nrow = 0))
colnames(tableDf) <- tableColNames
for (i in seq_len(nStages)) {
stage <- unique(dataset[[stages]])[i]
dataCurrentStage <- dataset[which(dataset[[stages]] == stage), ][!names(dataset) %in% stages]
values <- as.vector(na.omit(unlist(dataCurrentStage[measurements]))) # distribution fitting function complains if this is not explicitly a vector
stageValues <- as.vector(na.omit(unlist(dataCurrentStage[measurements]))) # distribution fitting function complains if this is not explicitly a vector

if (options[["nullDistribution"]] == "normal") {
meanx <- mean(values)
sdx <- sd(values)
test <- goftest::ad.test(x = values, "norm", mean = meanx, sd = sdx)
parameterEstimate1 <- mean(stageValues)
parameterEstimate2 <- stats::sd(stageValues)
andersonDarlingTest <- goftest::ad.test(x = stageValues, "norm", mean = parameterEstimate1, sd = parameterEstimate2)
} else if (options[["nullDistribution"]] == "lognormal") {
fit <- EnvStats::elnorm(values)
meanx <- fit$parameters[1]
sdx <- fit$parameters[2]
test <- goftest::ad.test(x = values, "plnorm", meanlog = meanx, sdlog = sdx)
fitPars <- .distributionParameters(stageValues, distribution = "lognormal")
parameterEstimate1 <- fitPars$beta
parameterEstimate2 <- fitPars$theta
andersonDarlingTest <- goftest::ad.test(x = stageValues, "plnorm", meanlog = parameterEstimate1, sdlog = parameterEstimate2)
} else if (options[["nullDistribution"]] == "weibull") {
fit <- fitdistrplus::fitdist(values, 'weibull',
control = list(maxit = 10000))
meanx <- fit$estimate[1]
sdx <- fit$estimate[2]
test <- goftest::ad.test(x = values, "pweibull", shape = meanx, scale = sdx)
}
n <- length(values)
ad <- test$statistic
adStar <- ad*(1 + (0.75/n) + (2.25/(n^2)))
if(ad >= 0.6) {
p <- exp(1.2937 - (5.709 * adStar) + 0.0186 * (adStar^2))
} else if(adStar < 0.6 && adStar > 0.34) {
p <- exp(0.9177 - (4.279 * adStar) - 1.38 * (adStar^2))
} else if(adStar < 0.34 && adStar > 0.2) {
p <- 1 - exp(-8.318 + (42.796 * adStar) - 59.938 * (adStar^2))
} else if(adStar <= 0.2) {
p <- 1 - exp(-13.436 + (101.14 * adStar) - 223.73 * (adStar^2)) #Jaentschi & Bolboaca (2018)
fitPars <- .distributionParameters(stageValues, distribution = "weibull")
parameterEstimate1 <- fitPars$beta
parameterEstimate2 <- fitPars$theta
andersonDarlingTest <- goftest::ad.test(x = stageValues, "pweibull", shape = parameterEstimate1, scale = parameterEstimate2)
}
sampleSize <- length(stageValues)
adStatistic <- andersonDarlingTest$statistic
adStatisticAdjusted <- adStatistic * (1 + (0.75 / sampleSize) + (2.25 / (sampleSize^2)))
if (adStatisticAdjusted >= 0.6) {
pValue <- exp(1.2937 - (5.709 * adStatisticAdjusted) + 0.0186 * (adStatisticAdjusted^2))
} else if (adStatisticAdjusted < 0.6 && adStatisticAdjusted > 0.34) {
pValue <- exp(0.9177 - (4.279 * adStatisticAdjusted) - 1.38 * (adStatisticAdjusted^2))
} else if (adStatisticAdjusted < 0.34 && adStatisticAdjusted > 0.2) {
Comment on lines +3144 to +3150

Copilot AI Apr 17, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The p-value approximation branches mix adjusted and unadjusted AD statistics: the first threshold compares adStatistic to 0.6, while subsequent branches compare adStatisticAdjusted. Since the approximation formulas use the adjusted statistic, the thresholding should be based on the adjusted value consistently (otherwise you can select the wrong branch for small/medium sample sizes).

Copilot uses AI. Check for mistakes.
pValue <- 1 - exp(-8.318 + (42.796 * adStatisticAdjusted) - 59.938 * (adStatisticAdjusted^2))
} else if (adStatisticAdjusted <= 0.2) {
pValue <- 1 - exp(-13.436 + (101.14 * adStatisticAdjusted) - 223.73 * (adStatisticAdjusted^2)) #Jaentschi & Bolboaca (2018)
} else {
p <- test$p.value
pValue <- andersonDarlingTest$p.value
}
tableDfCurrentStage <- data.frame(mean = meanx, sd = sdx, n = n, ad = round(ad, .numDecimals), p = round(p, .numDecimals))
tableDfCurrentStage <- data.frame(parameterEstimate1 = parameterEstimate1,
parameterEstimate2 = parameterEstimate2,
n = sampleSize,
ad = round(adStatistic, .numDecimals),
p = round(pValue, .numDecimals))
if (i == 1 && nStages > 1)
baseLineDf <- tableDfCurrentStage
if (i > 1) {
Expand Down
Loading
Loading