Skip to content

[Analyse Design] Stepwise analysis #389

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
JohnnyDoorn wants to merge 3 commits into
base: master
Choose a base branch
from
Open

[Analyse Design] Stepwise analysis #389

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
09f1fac
working version of stepwise regression
JohnnyDoorn Sep 22, 2025
07a83d6
improve table aesthetics
JohnnyDoorn Sep 22, 2025
8921661
also for factorial analysis
JohnnyDoorn Sep 24, 2025
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
145 changes: 107 additions & 38 deletions R/doeAnalysis.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,13 +27,15 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
covariates <- options[["covariates"]]
blocks <- options[["blocksFactorial"]]
dependent <- options[["dependentFactorial"]]
stepwiseMethod <- options[["stepwiseMethodFactorial"]]
} else if (options[["designType"]] == "responseSurfaceDesign") {
ready <- length(options[["continuousFactorsResponseSurface"]]) >= 1 && length(options[["dependentResponseSurface"]]) >= 1
discretePredictors <- options[["fixedFactorsResponseSurface"]]
continuousPredictors <- options[["continuousFactorsResponseSurface"]]
covariates <- NULL
blocks <- options[["blocksResponseSurface"]]
dependent <- options[["dependentResponseSurface"]]
stepwiseMethod <- options[["stepwiseMethodResponseSurface"]]
}

if (!ready) {
Expand All @@ -55,7 +57,7 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
jaspResults[[dep]] <- createJaspContainer(title = dep)
}

p <- try(.doeAnalysisMakeState(jaspResults, dataset, options, continuousPredictors, discretePredictors, blocks, covariates, dependent, ready))
p <- try(.doeAnalysisMakeState(jaspResults, dataset, options, continuousPredictors, discretePredictors, blocks, covariates, dependent, stepwiseMethod, ready))

if (isTryError(p)) {
jaspResults[["errorPlot"]] <- createJaspPlot(title = gettext("Error"))
Expand All @@ -70,6 +72,7 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
.doeAnalysisAnovaTable(jaspResults, dependent, options, ready, coded)
.doeAnalysisCoefficientsTable(jaspResults, dependent, options, ready, coded)
.doeAnalysisEquationTable(jaspResults, dependent, options, ready, coded)
.doeAnalysisStepwiseTable(jaspResults, dependent, options, ready, coded, stepwiseMethod)
.doeAnalysisPlotPareto(jaspResults, dependent, options, blocks, covariates, ready)
.doeAnalysisPlotEffectNormalDistribution(jaspResults, dependent, options, blocks, covariates, ready)
.doeAnalysisPlotQQResiduals(jaspResults, dependent, options, ready)
Expand Down Expand Up @@ -120,13 +123,13 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
"highestOrder", "order", "continuousFactorsFactorial", "modelTerms", "blocksFactorial",
"designType", "continuousFactorsResponseSurface", "codeFactors", "rsmPredefinedModel", "fixedFactorsFactorial",
"rsmPredefinedTerms", "dependentFactorial", "covariates", "codeFactorsMethod", "codeFactorsManualTable",
"squaredTerms", "sumOfSquaresType", "squaredTermsCoded")
"squaredTerms", "sumOfSquaresType", "squaredTermsCoded", "stepwiseMethodResponseSurface", "stepwiseMethodFactorial")
return(deps)
}

.scaleDOEvariable <- function(x){2*(x-min(x))/(max(x)-min(x))-1}

.doeAnalysisMakeState <- function(jaspResults, dataset, options, continuousPredictors, discretePredictors, blocks, covariates, dependent, ready) {
.doeAnalysisMakeState <- function(jaspResults, dataset, options, continuousPredictors, discretePredictors, blocks, covariates, dependent, stepwiseMethod, ready) {
if (!ready || jaspResults$getError()) {
return()
}
Expand Down Expand Up @@ -252,8 +255,45 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
covariateString <- paste0(" + ", unlist(covariates), collapse = "")
formulaString <- paste0(formulaString, covariateString)
}

formula <- as.formula(formulaString)

if (stepwiseMethod != "enter") {

stepwiseData <- if (options[["codeFactors"]]) datasetCoded else dataset
fullModel <- lm(formula, data = stepwiseData)

if (stepwiseMethod == "forward") {
initialFormula <- as.formula(paste(dep, "~ 1"))
} else {
initialFormula <- formula
}

stepwiseModel <- step(lm(initialFormula, data = stepwiseData),
direction = stepwiseMethod,
trace = 0,
scope = formula,
keep = function(model, aic)
list(model = model,
aic = aic,
coefficients = coef(model)))

excludedVars <- setdiff(names(fullModel$coefficients), names(stepwiseModel$coefficients))
formula <- as.formula(stepwiseModel)

# Make data frame from stepwise results, with step number, AIC of that model, and formula
stepwiseResult <- data.frame(
stepNr = 1:ncol(stepwiseModel$keep),
AIC = unlist(stepwiseModel$keep[2, ]),
formula = sapply(stepwiseModel$keep[3, ], function(x)
# apply some gsub magic to turn I(A^2) into A^2
paste(gsub("I\\((.*)\\^2\\)", "\\1^2", names(x)[-1]), collapse = " + "))
)

result[["stepwiseResult"]] <- list(excludedVars = gsub("I\\((.*)\\^2\\)", "\\1^2", excludedVars),
stepwiseResult = stepwiseResult)
}

regressionFit <- lm(formula, data = dataset)
regressionFitCoded <- lm(formula, data = datasetCoded)
regressionSummary <- summary(regressionFit)
Expand Down Expand Up @@ -497,40 +537,29 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
resultCoded[["regression"]][["coefficients"]][["p"]] <- rep(NA, length(valid_coefsCoded))
}

## Model formula
## uncoded
coefs <- coef(regressionFit)[!is.na(coef(regressionFit))]
coefs <- round(coefs, .numDecimals)
coefNames <- if (options[["tableAlias"]]) termNamesAliased else termNames
plusOrMin <- sapply(seq_len(length(coefs)), function(x) {
if (coefs[x] > 0) "+" else "-"
})
filledFormula <- sprintf("%s = %s %s %s %s %s", currentDependent, coefs[1], coefNames[1], plusOrMin[2], abs(coefs[2]), coefNames[2])
if (length(coefs) > 2) {
for (i in 3:length(coefs)) {
filledFormula <- sprintf("%s %s %s %s", filledFormula, plusOrMin[i], abs(coefs[i]), coefNames[i])
}
## Model formula for predictors, based on coding/aliasing
if (options[["tableAlias"]]) {
coefNames <- if (options[["codeFactors"]]) termNamesAliasedCoded else termNamesAliased
} else {
coefNames <- if (options[["codeFactors"]]) termNamesCoded else termNames
}

#coded
coefsCoded <- coef(regressionFitCoded)[!is.na(coef(regressionFitCoded))]
coefsCoded <- round(coefsCoded, .numDecimals)
coefNames <- if (options[["tableAlias"]]) termNamesAliasedCoded else termNames
plusOrMin <- sapply(seq_len(length(coefsCoded)), function(x) {
if (coefsCoded[x] > 0) "+" else "-"
})
filledFormulaCoded <- sprintf("%s = %s %s %s %s %s", currentDependent, coefsCoded[1], coefNames[1], plusOrMin[2], abs(coefsCoded[2]), coefNames[2])
if (length(coefsCoded) > 2) {
for (i in 3:length(coefsCoded)) {
filledFormulaCoded <- sprintf("%s %s %s %s", filledFormulaCoded, plusOrMin[i], abs(coefsCoded[i]), coefNames[i])
}
}
coefFormula <- paste(ifelse(sign(coefs[-1,1])==1, " +", " \u2013"),
round(abs(coefs[-1,1]), .numDecimals),
coefNames[-1],
collapse = "", sep = " ")

# Now add dependent name and intercept
filledFormula <- paste0(dep, " = ",
round(coefs[1, 1], .numDecimals),
coefFormula)

result[["regression"]][["filledFormula"]] <- gsubInteractionSymbol(filledFormula)
jaspResults[[currentDependent]][["doeResult"]] <- createJaspState(result)
jaspResults[[currentDependent]][["doeResult"]]$dependOn(options = .doeAnalysisBaseDependencies())


resultCoded[["regression"]][["filledFormula"]] <- gsubInteractionSymbol(filledFormulaCoded)
resultCoded[["regression"]][["filledFormula"]] <- gsubInteractionSymbol(filledFormula)
jaspResults[[currentDependent]][["doeResultCoded"]] <- createJaspState(resultCoded)
jaspResults[[currentDependent]][["doeResultCoded"]]$dependOn(options = .doeAnalysisBaseDependencies())
}
Expand Down Expand Up @@ -1157,9 +1186,9 @@ get_levels <- function(var, num_levels, dataset) {
# calculate model row
modelSS <- sum(anovaFit$`Sum Sq`[-nrow(anovaFit)])
modelDf <- sum(anovaFit$Df[-nrow(anovaFit)])
modelMS <- modelSS / modelDf
modelMS <- if (modelDf != 0) modelSS / modelDf else NA
msError <- anovaFit$`Mean Sq`[nrow(anovaFit)]
modelFValue <- if (msError != 0) modelMS / msError else NA
modelFValue <- if (msError != 0 && modelDf != 0) modelMS / msError else NA
modelPValue <- if (!is.na(modelFValue)) pf(modelFValue, modelDf, anovaFit$Df[nrow(anovaFit)], lower.tail = FALSE) else NA
modelRow <- data.frame(df = modelDf, ss = modelSS, ms = modelMS, f = modelFValue, p = modelPValue)
colnames(modelRow) <- colnames(anovaFit)
Expand Down Expand Up @@ -1326,6 +1355,7 @@ get_levels <- function(var, num_levels, dataset) {
if (!is.null(jaspResults[["tableAnova"]])) {
return()
}

tb <- createJaspTable(gettext("ANOVA"))
tb$addColumnInfo(name = "terms", title = gettext("Source"), type = "string")
tb$addColumnInfo(name = "adjss", title = gettext("Sum of squares"), type = "number")
Expand Down Expand Up @@ -1396,15 +1426,16 @@ get_levels <- function(var, num_levels, dataset) {
tb$dependOn(options = c("codeFactors", "codeFactorsManualTable", "codeFactorsMethod", "tableEquation", "tableAlias", .doeAnalysisBaseDependencies()))
tb2$position <- 4
jaspResults[[dep]][["tableCoefficientsLegend"]] <- tb2

result2 <- jaspResults[[dep]][["doeResult"]]$object[["regression"]][["factorLevels"]]
factorName <- if (options[["tableAlias"]]) result2[["factorNamesAliased"]] else result2[["factorNames"]]
rows2 <- data.frame(
factorName = factorName,
factorLevel = result2[["levels"]]
)
if (length(factorName) > 0) {
rows2 <- data.frame(
factorName = factorName,
factorLevel = result2[["levels"]]
)

tb2$addRows(rows2)
tb2$addRows(rows2)
}
}
}
}
Expand All @@ -1429,6 +1460,44 @@ get_levels <- function(var, num_levels, dataset) {
}
}

.doeAnalysisStepwiseTable <- function(jaspResults, dependent, options, ready, coded, stepwiseMethod) {
for (dep in dependent) {
if (stepwiseMethod == "enter" || !is.null(jaspResults[[dep]][["doeResult"]])) {
return()
}

tb <- createJaspTable(gettext("Stepwise results"))

tb$addColumnInfo(name = "stepNr", title = gettext("Step"), type = "string", combine = TRUE)
tb$addColumnInfo(name = "AIC", title = "AIC", type = "number")
tb$addColumnInfo(name = "formula", title = gettext("Formula"), type = "string")

tb$dependOn(options = .doeAnalysisBaseDependencies())
tb$position <- 4

jaspResults[[dep]][["stepwiseTable"]] <- tb

if (!ready || is.null(jaspResults[[dep]][["doeResult"]]) || jaspResults[[dep]]$getError()) {
return()
}

result <- jaspResults[[dep]][["doeResult"]]$object[["stepwiseResult"]]

tb$setData(result[["stepwiseResult"]])

excludedVars <- result[["excludedVars"]]
if (length(excludedVars) > 0) {
message <- sprintf(ngettext(length(excludedVars),
"The following covariate was considered but not included: %s.",
"The following covariates were considered but not included: %s."),
paste(excludedVars, collapse=", "))
tb$addFootnote(message)
}

}
}


.doeAnalysisPlotPareto <- function(jaspResults, dependent, options, blocks, covariates, ready) {
for (dep in dependent) {
if (!is.null(jaspResults[[dep]][["plotPareto"]]) || !options[["plotPareto"]]) {
Expand Down
24 changes: 24 additions & 0 deletions inst/qml/doeAnalysis.qml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,18 @@ Form
label: qsTr("Responses")
height: 50 * preferencesModel.uiScale
}

DropDown
{
name: "stepwiseMethodFactorial"
label: qsTr("Method")
info: qsTr("Specify the order in which the predictors are entered into the model. A block of one or more predictors represents one step in the hierarchy. Note that the present release does not allow for more than one block.")
values: [
{ label: qsTr("Enter"), info: qsTr("All predictors are entered into the model simultaneously.") ,value: "enter"},
{ label: qsTr("Backward"), info: qsTr("Starting with the full model, predictors are removed sequentially based on AIC."), value: "backward"},
{ label: qsTr("Forward"), info: qsTr("Starting with the intercept-only model, predictors are entered sequentially based on AIC.") , value: "forward"}
]
}

AssignedVariablesList
{
Expand Down Expand Up @@ -89,6 +101,18 @@ Form
label: qsTr("Responses")
height: 50 * preferencesModel.uiScale
}

DropDown
{
name: "stepwiseMethodResponseSurface"
label: qsTr("Method")
info: qsTr("Specify the order in which the predictors are entered into the model. A block of one or more predictors represents one step in the hierarchy. Note that the present release does not allow for more than one block.")
values: [
{ label: qsTr("Enter"), info: qsTr("All predictors are entered into the model simultaneously.") ,value: "enter"},
{ label: qsTr("Backward"), info: qsTr("Starting with the full model, predictors are removed sequentially based on AIC."), value: "backward"},
{ label: qsTr("Forward"), info: qsTr("Starting with the intercept-only model, predictors are entered sequentially based on AIC.") , value: "forward"}
]
}

AssignedVariablesList
{
Expand Down
Loading