ML Algorithm Performance Analysis
Anmol Shetty | NUID:002057185| Section:Saturday
2025-12-01
Video Presentation: [https://www.loom.com/share/8eb23f624b814b6bbce8403a4c69fc87]
Executive Summary
This project analyzes the comparative performance of machine learning classification algorithms across diverse datasets to identify which algorithm characteristics and dataset properties influence classification accuracy. Using the Meta-data dataset from the UCI Machine Learning Repository (Statlog Project, n=528 algorithm-dataset combinations), this study applies multiple linear regression, one-way ANOVA, and two-way ANOVA to examine how dataset properties, algorithm types, and their interactions affect classification performance.
Key Findings: - Dataset size (N) is the only significant predictor of error (p < 0.001), but explains only 6.3% of variance - Algorithm family choice significantly affects performance (p = 0.015), with Neural Networks performing significantly worse than other families - No interaction between algorithm family and problem complexity (p = 0.338) - algorithm selection does not need to vary by number of classes
1. Introduction
Understanding algorithm selection is critical for practitioners who must choose appropriate methods for specific data characteristics. With hundreds of machine learning algorithms available, selecting the right approach for a given classification problem remains challenging. Poor algorithm choices lead to suboptimal performance, wasted computational resources, and unreliable predictions in production systems.
This project addresses three key questions using the Meta-data dataset from the UCI Machine Learning Repository’s Statlog Project:
- How do dataset characteristics predict classification difficulty?
- Do different algorithm families show different performance levels?
- Does optimal algorithm choice depend on problem complexity?
By applying multiple linear regression, one-way ANOVA, and two-way ANOVA to 528 algorithm-dataset test results spanning 24 algorithms and 22 benchmark datasets, this study provides evidence-based guidance for algorithm selection in real-world classification tasks.
2. Data and Methods
2.1 Dataset Description
Source: Meta-data Dataset (Statlog Project)
URL: https://archive.ics.uci.edu/dataset/65/meta+data
Citation: Meta-data [Dataset]. (1994). UCI Machine
Learning Repository. https://doi.org/10.24432/C5X31P
Size: 528 observations (algorithm-dataset combinations) across 22 variables
# Load the data
# NOTE: Save meta.data file in your project folder first
meta_data <- read.csv("meta.data", header = FALSE)
# Add column names based on meta.names documentation
colnames(meta_data) <- c("DS_Name", "T", "N", "p", "k", "Bin", "Cost",
"SDratio", "correl", "cancor1", "cancor2",
"fract1", "fract2", "skewness", "kurtosis",
"Hc", "Hx", "MCx", "EnAtr", "NSRatio",
"Alg_Name", "Norm_error")
# View structure
str(meta_data)## 'data.frame': 528 obs. of 22 variables:
## $ DS_Name : chr "Aust_Credit" "Aust_Credit" "Aust_Credit" "Aust_Credit" ...
## $ T : int 690 690 690 690 690 690 690 690 690 690 ...
## $ N : int 690 690 690 690 690 690 690 690 690 690 ...
## $ p : int 14 14 14 14 14 14 14 14 14 14 ...
## $ k : int 2 2 2 2 2 2 2 2 2 2 ...
## $ Bin : int 4 4 4 4 4 4 4 4 4 4 ...
## $ Cost : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SDratio : num 1.26 1.26 1.26 1.26 1.26 ...
## $ correl : chr "0.1024" "0.1024" "0.1024" "0.1024" ...
## $ cancor1 : num 0.771 0.771 0.771 0.771 0.771 ...
## $ cancor2 : chr "?" "?" "?" "?" ...
## $ fract1 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ fract2 : chr "?" "?" "?" "?" ...
## $ skewness : num 1.97 1.97 1.97 1.97 1.97 ...
## $ kurtosis : num 12.6 12.6 12.6 12.6 12.6 ...
## $ Hc : num 0.991 0.991 0.991 0.991 0.991 ...
## $ Hx : num 2.3 2.3 2.3 2.3 2.3 ...
## $ MCx : num 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 ...
## $ EnAtr : num 8.77 8.77 8.77 8.77 8.77 ...
## $ NSRatio : num 19.4 19.4 19.4 19.4 19.4 ...
## $ Alg_Name : chr "Ac2" "Alloc80" "BackProp" "Bayes" ...
## $ Norm_error: num 3.89 5.45 1.79 1.56 3.11 ...head(meta_data, 10)## DS_Name T N p k Bin Cost SDratio correl cancor1 cancor2 fract1
## 1 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 2 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 3 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 4 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 5 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 6 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 7 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 8 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 9 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## 10 Aust_Credit 690 690 14 2 4 0 1.2623 0.1024 0.7713 ? 1
## fract2 skewness kurtosis Hc Hx MCx EnAtr NSRatio Alg_Name
## 1 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Ac2
## 2 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Alloc80
## 3 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 BackProp
## 4 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Bayes
## 5 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 BayesTree
## 6 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 C4.5
## 7 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 CART
## 8 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Cal5
## 9 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Cascade
## 10 ? 1.9701 12.5538 0.9912 2.3012 0.113 8.771681 19.3646 Castle
## Norm_error
## 1 3.893
## 2 5.450
## 3 1.791
## 4 1.557
## 5 3.114
## 6 1.868
## 7 1.090
## 8 0.000
## 9 67.655
## 10 1.324# Summary statistics
summary(meta_data)## DS_Name T N p
## Length:528 Min. : 270 Min. : 270 Min. : 6.00
## Class :character 1st Qu.: 900 1st Qu.: 900 1st Qu.: 13.00
## Mode :character Median : 1750 Median : 4092 Median : 16.00
## Mean : 4569 Mean :10734 Mean : 29.55
## 3rd Qu.: 7480 3rd Qu.:18000 3rd Qu.: 36.00
## Max. :20000 Max. :58000 Max. :180.00
## k Bin Cost SDratio
## Min. : 2.000 Min. : 0.000 Min. :0.0000 Min. :1.027
## 1st Qu.: 2.000 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:1.064
## Median : 3.000 Median : 0.000 Median :0.0000 Median :1.280
## Mean : 9.727 Mean : 3.182 Mean :0.1364 Mean :1.479
## 3rd Qu.: 7.000 3rd Qu.: 1.000 3rd Qu.:0.0000 3rd Qu.:1.567
## Max. :91.000 Max. :43.000 Max. :1.0000 Max. :4.001
## correl cancor1 cancor2 fract1
## Length:528 Min. :0.5044 Length:528 Min. :0.1505
## Class :character 1st Qu.:0.7176 Class :character 1st Qu.:0.3586
## Mode :character Median :0.8575 Mode :character Median :0.9376
## Mean :0.7948 Mean :0.7007
## 3rd Qu.:0.9165 3rd Qu.:1.0000
## Max. :0.9884 Max. :1.0000
## fract2 skewness kurtosis Hc
## Length:528 Min. :0.1802 Min. : 0.9866 Min. :0.2893
## Class :character 1st Qu.:0.7316 1st Qu.: 3.6494 1st Qu.:0.9453
## Mode :character Median :1.0329 Median : 5.0832 Median :1.1787
## Mean :1.7842 Mean : 22.6672 Mean :1.8716
## 3rd Qu.:1.9701 3rd Qu.: 12.5538 3rd Qu.:2.8072
## Max. :6.7156 Max. :160.3108 Max. :4.8787
## Hx MCx EnAtr NSRatio
## Min. :0.3672 Min. :0.0187 Min. : 1.560 Min. : 1.023
## 1st Qu.:1.7700 1st Qu.:0.0495 1st Qu.: 3.508 1st Qu.: 4.963
## Median :3.2605 Median :0.1979 Median : 7.684 Median : 14.446
## Mean :3.3450 Mean :0.3168 Mean : 20.664 Mean : 28.873
## 3rd Qu.:4.6908 3rd Qu.:0.5049 3rd Qu.: 16.372 3rd Qu.: 36.028
## Max. :6.5452 Max. :1.3149 Max. :160.644 Max. :159.644
## Alg_Name Norm_error
## Length:528 Min. : 0.000
## Class :character 1st Qu.: 3.315
## Mode :character Median : 9.840
## Mean : 99.552
## 3rd Qu.: 25.909
## Max. :12040.9922.2 Data Cleaning and Preparation
# Fix data types and handle missing values
meta_data <- meta_data %>%
mutate(
# Convert character variables to numeric (handle "?" as NA)
correl = as.numeric(correl),
cancor2 = ifelse(cancor2 == "?", NA, as.numeric(cancor2)),
fract2 = ifelse(fract2 == "?", NA, as.numeric(fract2)),
# Convert categorical variables to factors
DS_Name = as.factor(DS_Name),
Alg_Name = as.factor(Alg_Name),
Cost = as.factor(Cost)
)
# Check the fixes
str(meta_data)## 'data.frame': 528 obs. of 22 variables:
## $ DS_Name : Factor w/ 22 levels "Aust_Credit",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ T : int 690 690 690 690 690 690 690 690 690 690 ...
## $ N : int 690 690 690 690 690 690 690 690 690 690 ...
## $ p : int 14 14 14 14 14 14 14 14 14 14 ...
## $ k : int 2 2 2 2 2 2 2 2 2 2 ...
## $ Bin : int 4 4 4 4 4 4 4 4 4 4 ...
## $ Cost : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ SDratio : num 1.26 1.26 1.26 1.26 1.26 ...
## $ correl : num 0.102 0.102 0.102 0.102 0.102 ...
## $ cancor1 : num 0.771 0.771 0.771 0.771 0.771 ...
## $ cancor2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ fract1 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ fract2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ skewness : num 1.97 1.97 1.97 1.97 1.97 ...
## $ kurtosis : num 12.6 12.6 12.6 12.6 12.6 ...
## $ Hc : num 0.991 0.991 0.991 0.991 0.991 ...
## $ Hx : num 2.3 2.3 2.3 2.3 2.3 ...
## $ MCx : num 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 ...
## $ EnAtr : num 8.77 8.77 8.77 8.77 8.77 ...
## $ NSRatio : num 19.4 19.4 19.4 19.4 19.4 ...
## $ Alg_Name : Factor w/ 24 levels "Ac2","Alloc80",..: 1 2 3 4 5 6 8 7 9 10 ...
## $ Norm_error: num 3.89 5.45 1.79 1.56 3.11 ...# Create algorithm families based on machine learning methodology
# Addressing peer review feedback: clear justification for groupings
meta_data <- meta_data %>%
mutate(
Alg_Family = case_when(
# Neural Networks: Gradient-based learning with hidden layers
Alg_Name %in% c("BackProp", "Cascade", "LVQ", "Kohonen", "RBF") ~ "Neural Networks",
# Statistical Methods: Probabilistic and discriminant analysis
Alg_Name %in% c("Bayes", "BayesTree", "Discrim", "LogDisc", "QuaDisc", "Dipol92") ~ "Statistical Methods",
# Decision Trees: Recursive partitioning algorithms
Alg_Name %in% c("C4.5", "CART", "IndCART", "Cal5", "NewId") ~ "Decision Trees",
# Instance-Based: Memory-based learning using similarity metrics
Alg_Name %in% c("KNN", "Cn2", "ITrule", "Smart", "Castle", "Ac2", "Alloc80", "Default") ~ "Instance-Based",
TRUE ~ "Other"
),
Alg_Family = as.factor(Alg_Family)
)
# Verify all algorithms are accounted for
cat("Algorithm distribution across families:\n")## Algorithm distribution across families:table(meta_data$Alg_Family, meta_data$Alg_Name)##
## Ac2 Alloc80 BackProp Bayes BayesTree C4.5 Cal5 CART
## Decision Trees 0 0 0 0 0 22 22 22
## Instance-Based 22 22 0 0 0 0 0 0
## Neural Networks 0 0 22 0 0 0 0 0
## Statistical Methods 0 0 0 22 22 0 0 0
##
## Cascade Castle Cn2 Default Dipol92 Discrim IndCART ITrule
## Decision Trees 0 0 0 0 0 0 22 0
## Instance-Based 0 22 22 22 0 0 0 22
## Neural Networks 22 0 0 0 0 0 0 0
## Statistical Methods 0 0 0 0 22 22 0 0
##
## KNN Kohonen LogDisc LVQ NewId QuaDisc RBF Smart
## Decision Trees 0 0 0 0 22 0 0 0
## Instance-Based 22 0 0 0 0 0 0 22
## Neural Networks 0 22 0 22 0 0 22 0
## Statistical Methods 0 0 22 0 0 22 0 0# Create problem complexity levels based on number of classes (k)
meta_data <- meta_data %>%
mutate(
Complexity = case_when(
k == 2 ~ "Binary",
k >= 3 & k <= 10 ~ "Low Multi-class",
k > 10 ~ "High Multi-class"
),
Complexity = factor(Complexity, levels = c("Binary", "Low Multi-class", "High Multi-class"))
)
# Check distribution
cat("\nComplexity level distribution:\n")##
## Complexity level distribution:table(meta_data$Complexity)##
## Binary Low Multi-class High Multi-class
## 240 216 72# Final dataset structure
cat("\nFinal dataset structure:\n")##
## Final dataset structure:str(meta_data)## 'data.frame': 528 obs. of 24 variables:
## $ DS_Name : Factor w/ 22 levels "Aust_Credit",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ T : int 690 690 690 690 690 690 690 690 690 690 ...
## $ N : int 690 690 690 690 690 690 690 690 690 690 ...
## $ p : int 14 14 14 14 14 14 14 14 14 14 ...
## $ k : int 2 2 2 2 2 2 2 2 2 2 ...
## $ Bin : int 4 4 4 4 4 4 4 4 4 4 ...
## $ Cost : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ SDratio : num 1.26 1.26 1.26 1.26 1.26 ...
## $ correl : num 0.102 0.102 0.102 0.102 0.102 ...
## $ cancor1 : num 0.771 0.771 0.771 0.771 0.771 ...
## $ cancor2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ fract1 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ fract2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ skewness : num 1.97 1.97 1.97 1.97 1.97 ...
## $ kurtosis : num 12.6 12.6 12.6 12.6 12.6 ...
## $ Hc : num 0.991 0.991 0.991 0.991 0.991 ...
## $ Hx : num 2.3 2.3 2.3 2.3 2.3 ...
## $ MCx : num 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 0.113 ...
## $ EnAtr : num 8.77 8.77 8.77 8.77 8.77 ...
## $ NSRatio : num 19.4 19.4 19.4 19.4 19.4 ...
## $ Alg_Name : Factor w/ 24 levels "Ac2","Alloc80",..: 1 2 3 4 5 6 8 7 9 10 ...
## $ Norm_error: num 3.89 5.45 1.79 1.56 3.11 ...
## $ Alg_Family: Factor w/ 4 levels "Decision Trees",..: 2 2 3 4 4 1 1 1 3 2 ...
## $ Complexity: Factor w/ 3 levels "Binary","Low Multi-class",..: 1 1 1 1 1 1 1 1 1 1 ...cat("\nDataset ready for analysis with", nrow(meta_data), "observations\n")##
## Dataset ready for analysis with 528 observations3. Analysis 1: Multiple Linear Regression
3.1 Research Question
How do dataset characteristics (number of examples N, number of attributes p, attribute correlation, and noise-signal ratio) predict classification error rates across machine learning algorithms?
3.2 Exploratory Data Analysis
# Scatterplot matrix for regression variables
regression_vars <- meta_data %>%
dplyr::select(Norm_error, N, p, correl, NSRatio) # Explicitly use dplyr::select
# Create scatterplot matrix
ggpairs(regression_vars,
title = "Scatterplot Matrix: Response and Predictors",
lower = list(continuous = wrap("points", alpha = 0.3, size = 0.5)),
diag = list(continuous = wrap("densityDiag", alpha = 0.5)),
upper = list(continuous = wrap("cor", size = 4)))
# Individual scatterplots
p1 <- ggplot(meta_data, aes(x = N, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Number of Examples",
x = "N (Number of Training Examples)",
y = "Normalized Error (%)")
p2 <- ggplot(meta_data, aes(x = p, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Number of Attributes",
x = "p (Number of Attributes)",
y = "Normalized Error (%)")
p3 <- ggplot(meta_data, aes(x = correl, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Attribute Correlation",
x = "Mean Correlation Between Attributes",
y = "Normalized Error (%)")
p4 <- ggplot(meta_data, aes(x = NSRatio, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Noise-Signal Ratio",
x = "Noise-Signal Ratio",
y = "Normalized Error (%)")
grid.arrange(p1, p2, p3, p4, ncol = 2)
3.3 Fit Multiple Linear Regression Model
# Fit multiple regression model
regression_model <- lm(Norm_error ~ N + p + correl + NSRatio, data = meta_data)
# Model summary
summary(regression_model)##
## Call:
## lm(formula = Norm_error ~ N + p + correl + NSRatio, data = meta_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -785.3 -152.5 1.5 62.3 11255.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -82.854932 79.357803 -1.044 0.297
## N 0.014194 0.002364 6.004 3.71e-09 ***
## p 0.253593 0.953664 0.266 0.790
## correl 126.492429 191.121722 0.662 0.508
## NSRatio -0.262710 0.931414 -0.282 0.778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 757.5 on 499 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.07037, Adjusted R-squared: 0.06292
## F-statistic: 9.443 on 4 and 499 DF, p-value: 2.302e-07# Coefficient table with confidence intervals
coef_table <- cbind(
Estimate = coef(regression_model),
confint(regression_model)
)
round(coef_table, 4)## Estimate 2.5 % 97.5 %
## (Intercept) -82.8549 -238.7715 73.0617
## N 0.0142 0.0095 0.0188
## p 0.2536 -1.6201 2.1273
## correl 126.4924 -249.0100 501.9949
## NSRatio -0.2627 -2.0927 1.56733.4 Model Diagnostics
# Four standard diagnostic plots
par(mfrow = c(2, 2))
plot(regression_model)
par(mfrow = c(1, 1))# Check multicollinearity
cat("=== Variance Inflation Factors ===\n")## === Variance Inflation Factors ===vif_values <- vif(regression_model)
print(vif_values)## N p correl NSRatio
## 1.038790 1.057123 1.058291 1.066561cat("\nMulticollinearity assessment:\n")##
## Multicollinearity assessment:if(max(vif_values) > 10) {
cat("WARNING: Severe multicollinearity detected\n")
} else if(max(vif_values) > 5) {
cat("CAUTION: Moderate multicollinearity present\n")
} else {
cat("GOOD: No concerning multicollinearity\n")
}## GOOD: No concerning multicollinearity# Model fit statistics
cat("\n=== Model Fit Statistics ===\n")##
## === Model Fit Statistics ===cat("R-squared:", round(summary(regression_model)$r.squared, 4), "\n")## R-squared: 0.0704cat("Adjusted R-squared:", round(summary(regression_model)$adj.r.squared, 4), "\n")## Adjusted R-squared: 0.0629cat("F-statistic:", round(summary(regression_model)$fstatistic[1], 2), "\n")## F-statistic: 9.44cat("p-value:", format.pval(pf(summary(regression_model)$fstatistic[1],
summary(regression_model)$fstatistic[2],
summary(regression_model)$fstatistic[3],
lower.tail = FALSE), digits = 4), "\n")## p-value: 2.302e-073.5 Interpretation of Results
Regression Equation:
Norm_error = -82.85 + 0.0142(N) + 0.254(p) + 126.49(correl) - 0.263(NSRatio)
Coefficient Interpretations:
Intercept (β₀ = -82.85):
The predicted classification error is -82.85% when all predictors equal
zero. This value is not practically meaningful as it’s impossible to
have zero training examples, attributes, or correlation in real
datasets. This is simply the mathematical baseline for the equation.
N coefficient (β₁ = 0.0142, p < 0.001):
For each additional training example, classification error increases by
0.0142 percentage points on average, holding number of
attributes, correlation, and noise-signal ratio constant. For each
additional 1,000 training examples, error increases by approximately
14.2 percentage points. This is highly statistically
significant (p = 3.71e-09).
Counterintuitive finding: Larger datasets are associated with higher error rates, likely because complex, difficult problems (like character recognition with 20,000+ examples) naturally have both more data and higher inherent difficulty.
p coefficient (β₂ = 0.254, p = 0.790):
For each additional attribute, classification error increases by 0.254
percentage points on average, holding other predictors constant. This
effect is not statistically significant (p = 0.790),
suggesting the number of attributes alone does not reliably predict
classification difficulty.
correl coefficient (β₃ = 126.49, p = 0.508):
For each 0.1 unit increase in mean attribute correlation, classification
error increases by approximately 12.6 percentage points on average,
holding other predictors constant. This effect is not
statistically significant (p = 0.508).
NSRatio coefficient (β₄ = -0.263, p = 0.778):
For each unit increase in noise-signal ratio, classification error
decreases by 0.263 percentage points on average, holding other
predictors constant. This effect is not statistically
significant (p = 0.778).
Model Fit and Diagnostics:
Overall Model Significance: - F-statistic = 9.44 (p = 2.302e-07) - The model is statistically significant overall - However, R² = 0.070, Adjusted R² = 0.063 - only 6.3% of variance explained - This low R² suggests other factors (algorithm type, problem domain, class balance) are more important
Multicollinearity: No issues - All VIF values < 2, indicating predictors are independent
Assumption Violations: Severe outliers: Observations 489, 249, 441, 444 show extreme residuals (>10,000) Non-normality: Q-Q plot shows heavy right tail deviation Heteroscedasticity: Residuals vs Fitted shows funnel pattern with increasing variance Non-linearity: Possible non-linear relationships visible in scatterplots
Answer to Research Question:
How do dataset characteristics predict classification error?
Among the tested predictors (N, p, correl, NSRatio), only dataset size (N) significantly predicts classification error (p < 0.001). Counterintuitively, larger datasets are associated with higher error rates, likely because complex, difficult problems naturally require more training data.
However, the model explains only 6.3% of variance, indicating that dataset characteristics alone are poor predictors of classification performance. The presence of severe outliers and assumption violations suggests that:
- Algorithm choice (not included in this model) is likely more important than dataset properties
- Non-linear relationships or interactions between variables may better explain error rates
- Data transformations (log-transform for highly skewed error distribution) may improve model fit
Practical implication: Practitioners cannot reliably predict classification difficulty based solely on dataset size, dimensionality, correlation, or noise levels. Algorithm selection and problem-specific factors play a more critical role.
4. Analysis 2: One-Way ANOVA
4.1 Research Question
Do different machine learning algorithm families (Neural Networks, Statistical Methods, Decision Trees, Instance-Based Learning) show significantly different mean classification error rates?
4.2 Exploratory Visualization
# Boxplot comparing algorithm families
ggplot(meta_data, aes(x = Alg_Family, y = Norm_error, fill = Alg_Family)) +
geom_boxplot(alpha = 0.7, outlier.alpha = 0.3) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Classification Error by Algorithm Family",
subtitle = "Comparing Neural Networks, Statistical, Decision Trees, Instance-Based",
x = "Algorithm Family",
y = "Normalized Error (%)") +
theme(legend.position = "none") +
coord_cartesian(ylim = c(0, 500)) # Zoom in to see most data
# Summary statistics by algorithm family
family_summary <- meta_data %>%
group_by(Alg_Family) %>%
summarize(
n = n(),
mean_error = mean(Norm_error, na.rm = TRUE),
sd_error = sd(Norm_error, na.rm = TRUE),
median_error = median(Norm_error, na.rm = TRUE),
min_error = min(Norm_error, na.rm = TRUE),
max_error = max(Norm_error, na.rm = TRUE)
)
print(family_summary)## # A tibble: 4 × 7
## Alg_Family n mean_error sd_error median_error min_error max_error
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Decision Trees 110 41.3 124. 7.69 0 885.
## 2 Instance-Based 176 54.4 208. 11.1 0 2576.
## 3 Neural Networks 110 309. 1635. 13.9 0 12041.
## 4 Statistical Methods 132 34.1 106. 7.56 0 808.4.3 One-Way ANOVA Test
# Fit One-Way ANOVA
anova_model <- aov(Norm_error ~ Alg_Family, data = meta_data)
# ANOVA table
summary(anova_model)## Df Sum Sq Mean Sq F value Pr(>F)
## Alg_Family 3 6108068 2036023 3.533 0.0147 *
## Residuals 524 301995696 576328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Calculate effect size (eta-squared)
anova_summary <- summary(anova_model)
ss_between <- anova_summary[[1]]["Alg_Family", "Sum Sq"]
ss_total <- sum(anova_summary[[1]][, "Sum Sq"])
eta_squared <- ss_between / ss_total
cat("Effect Size (η²):", round(eta_squared, 4), "\n")## Effect Size (η²): 0.0198cat("Interpretation:",
ifelse(eta_squared < 0.01, "negligible",
ifelse(eta_squared < 0.06, "small",
ifelse(eta_squared < 0.14, "medium", "large"))), "effect\n")## Interpretation: small effect4.4 Post-Hoc Tests (Tukey HSD)
# Tukey HSD for pairwise comparisons
tukey_results <- TukeyHSD(anova_model)
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Norm_error ~ Alg_Family, data = meta_data)
##
## $Alg_Family
## diff lwr upr p adj
## Instance-Based-Decision Trees 13.11174 -224.697652 250.92113 0.9989757
## Neural Networks-Decision Trees 267.36948 3.543652 531.19531 0.0456227
## Statistical Methods-Decision Trees -7.22220 -259.816202 245.37180 0.9998564
## Neural Networks-Instance-Based 254.25774 16.448355 492.06713 0.0307667
## Statistical Methods-Instance-Based -20.33394 -245.618129 204.95025 0.9955634
## Statistical Methods-Neural Networks -274.59168 -527.185683 -21.99768 0.0269793# Visualize Tukey HSD results
plot(tukey_results, las = 1, col = "blue")
# Compact letter display
library(multcomp)
cld_result <- cld(glht(anova_model, linfct = mcp(Alg_Family = "Tukey")))
cat("\nCompact Letter Display:\n")##
## Compact Letter Display:print(cld_result)## Decision Trees Instance-Based Neural Networks Statistical Methods
## "a" "a" "b" "a"4.5 Assumption Checks
# Check ANOVA assumptions
par(mfrow = c(2, 2))
# 1. Residuals vs Fitted
plot(anova_model, which = 1, main = "Residuals vs Fitted")
# 2. Q-Q Plot
plot(anova_model, which = 2, main = "Normal Q-Q")
# 3. Levene's Test for equal variances
leveneTest(Norm_error ~ Alg_Family, data = meta_data)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 3.4674 0.01611 *
## 524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(1, 1))
cat("\n=== ANOVA Assumption Checks ===\n\n")##
## === ANOVA Assumption Checks ===cat("1. INDEPENDENCE: Assumed based on study design (each algorithm-dataset combo is independent)\n\n")## 1. INDEPENDENCE: Assumed based on study design (each algorithm-dataset combo is independent)cat("2. NORMALITY: Check Q-Q plot - if points follow diagonal line, normality is reasonable\n\n")## 2. NORMALITY: Check Q-Q plot - if points follow diagonal line, normality is reasonablecat("3. EQUAL VARIANCES: Check Levene's test:\n")## 3. EQUAL VARIANCES: Check Levene's test:cat(" - If p > 0.05: Equal variances assumption met\n")## - If p > 0.05: Equal variances assumption metcat(" - If p < 0.05: Variances differ (may need Welch's ANOVA)\n")## - If p < 0.05: Variances differ (may need Welch's ANOVA)4.6 One-Way ANOVA Results Interpretation
Hypotheses:
- H₀: All algorithm families have equal mean classification error rates
- Hₐ: At least one algorithm family has a different mean error rate
ANOVA Results:
- F-statistic: 3.533
- p-value: 0.0147
- Conclusion: Reject H₀ at α = 0.05. There are statistically significant differences in mean error rates across algorithm families.
- Effect Size (η²): 0.0198 (small effect) - algorithm family explains only 2% of variance
Tukey HSD Pairwise Comparisons:
Significant Differences (p.adj < 0.05):
- Neural Networks vs Decision Trees: diff = 267.37, p
= 0.0457 ✱
- Neural Networks have 267% higher mean error than Decision Trees
- Neural Networks vs Instance-Based: diff = 254.26, p
= 0.0308 ✱
- Neural Networks have 254% higher mean error than Instance-Based
- Statistical Methods vs Neural Networks: diff =
-274.59, p = 0.0270 ✱
- Statistical Methods have 275% lower mean error than Neural Networks
Non-Significant Differences (p.adj > 0.05): - Instance-Based vs Decision Trees (p = 0.999) - Statistical Methods vs Decision Trees (p = 0.999) - Statistical Methods vs Instance-Based (p = 0.996)
Compact Letter Display:
- Group “a”: Decision Trees, Instance-Based, Statistical Methods (no significant differences)
- Group “b”: Neural Networks (significantly different from group “a”)
Assumption Checks:
1. Independence: ✅ Met - Each algorithm-dataset combination is independent
2. Normality: ⚠️ Violated - Q-Q plot shows severe right tail deviation (extreme outliers)
3. Equal Variances: ⚠️ Violated - Levene’s test p = 0.0161 < 0.05 - Neural Networks have much higher variance (SD = 1,635) than other groups
Impact of Violations: - Unequal variances and non-normality reduce ANOVA reliability - Results should be interpreted cautiously - Could consider Welch’s ANOVA (robust to unequal variances) or non-parametric Kruskal-Wallis test
Answer to Research Question:
Yes, algorithm families show statistically significant differences (F = 3.533, p = 0.015), but with important caveats:
Neural Networks are outliers - significantly worse than all other families due to catastrophic failures on some datasets (max error = 12,041%)
Decision Trees, Instance-Based, and Statistical Methods perform similarly - no significant differences among these three families (all p > 0.99)
Small effect size (η² = 0.02) - family membership explains only 2% of error variance
High within-group variability - huge standard deviations suggest individual algorithm choice matters far more than family
Practical Implication: Avoid Neural Networks for this benchmark (prone to catastrophic failures), but among the other three families, specific algorithm selection based on problem characteristics is more important than choosing a family. The similar performance of Decision Trees, Instance-Based, and Statistical Methods suggests no single family dominates across all problem types.
5. Analysis 3: Two-Way ANOVA
5.1 Research Question
How do algorithm family and problem complexity (measured by number of classes) interact to affect classification error? Does the optimal algorithm family depend on whether the classification task is binary, low multi-class, or high multi-class?
5.2 Interaction Plot
# Calculate group means for interaction
interaction_means <- meta_data %>%
group_by(Alg_Family, Complexity) %>%
summarize(mean_error = mean(Norm_error, na.rm = TRUE), .groups = "drop")
# Create interaction plot
ggplot(interaction_means, aes(x = Complexity, y = mean_error,
color = Alg_Family, group = Alg_Family)) +
geom_point(size = 4) +
geom_line(size = 1.2) +
scale_color_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Interaction Plot: Algorithm Family × Problem Complexity",
subtitle = "Do lines cross? That indicates interaction!",
x = "Problem Complexity",
y = "Mean Classification Error (%)",
color = "Algorithm Family") +
theme(legend.position = "top")
5.3 Two-Way ANOVA Model
# Fit two-way ANOVA with interaction
two_way_model <- aov(Norm_error ~ Alg_Family * Complexity, data = meta_data)
# ANOVA table
summary(two_way_model)## Df Sum Sq Mean Sq F value Pr(>F)
## Alg_Family 3 6108068 2036023 3.550 0.0144 *
## Complexity 2 2121622 1060811 1.850 0.1583
## Alg_Family:Complexity 6 3921117 653520 1.139 0.3380
## Residuals 516 295952957 573552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.4 Two-Way ANOVA Results Interpretation
ANOVA Table Results:
| Source | Df | F-value | p-value | Significant? |
|---|---|---|---|---|
| Alg_Family (Main Effect) | 3 | 3.550 | 0.0144 | ✓ Yes |
| Complexity (Main Effect) | 2 | 1.850 | 0.1583 | ✗ No |
| Alg_Family × Complexity (Interaction) | 6 | 1.139 | 0.3380 | ✗ No |
Interpretation:
Main Effect - Algorithm Family: - Significant (F = 3.550, p = 0.0144) - Algorithm family affects error rates regardless of problem complexity - This confirms the one-way ANOVA finding
Main Effect - Problem Complexity: - Not significant (F = 1.850, p = 0.158) - Binary, low multi-class, and high multi-class problems do not differ significantly in average error - Surprising: We expected more complex problems (more classes) to have higher error
Interaction Effect: - Not significant (F = 1.139, p = 0.338) - No evidence that algorithm family performance depends on problem complexity - The lines in the interaction plot do NOT cross significantly
What the Interaction Plot Shows:
Looking at the plot: - Neural Networks (green line): Shows a peak at Low Multi-class (~580% error), then drops at High Multi-class - This unusual pattern is likely driven by outliers rather than true interaction - Other three families (red, blue, orange): Remain relatively flat and parallel across complexity levels - Decision Trees, Instance-Based, and Statistical Methods show consistent performance regardless of problem complexity
Key Insight: Despite visual differences in the plot (especially Neural Networks’ peak), the interaction is not statistically significant (p = 0.338). This means: - We should interpret main effects independently - Algorithm family choice matters, but it matters equally across all complexity levels - The apparent Neural Networks spike is likely due to a few extreme outliers, not a systematic interaction
Answer to Research Question:
Does optimal algorithm family depend on problem complexity?
No. The interaction between algorithm family and problem complexity is not significant (p = 0.338). This means:
- Algorithm family matters (main effect p = 0.014), but its effect doesn’t change based on problem complexity
- Complexity alone doesn’t affect error (main effect p = 0.158)
- No context-dependent algorithm selection needed - if a family performs well on binary problems, it performs similarly on multi-class problems
Practical Implication:
Practitioners can choose algorithms based on family characteristics without worrying about the number of classes in their problem. The Neural Networks family’s poor average performance (driven by outliers) applies equally to binary and multi-class problems. Decision Trees, Instance-Based, and Statistical Methods maintain consistent, reliable performance across all complexity levels.
Recommendation: For reliable performance, use Decision Trees, Instance-Based, or Statistical Methods regardless of problem complexity (2-91 classes). Avoid Neural Networks unless specific algorithms within that family are known to work well for your domain.
5.5 Sensitivity Analysis - Addressing Assumption Violations
Since our regression and ANOVA models showed violations of normality and equal variance assumptions, we conduct sensitivity analyses to verify that our conclusions remain valid under alternative approaches.
5.5.1 Log Transformation of Response Variable
Extreme outliers and right-skewed error distributions suggest a log transformation may be appropriate.
# Create log-transformed response variable
# Using log1p (log(1+x)) to handle zeros
meta_data <- meta_data %>%
mutate(log_error = log1p(Norm_error))
# Visualize transformation effect
par(mfrow = c(1, 2))
hist(meta_data$Norm_error, breaks = 30, col = "lightblue",
main = "Original Error Distribution", xlab = "Normalized Error")
hist(meta_data$log_error, breaks = 30, col = "lightgreen",
main = "Log-Transformed Error", xlab = "log(1 + Error)")
par(mfrow = c(1, 1))5.5.2 Regression with Log-Transformed Response
# Refit regression model with log-transformed response
regression_log <- lm(log_error ~ N + p + correl + NSRatio, data = meta_data)
# Compare model summaries
cat("=== COMPARISON: Original vs Log-Transformed ===\n\n")## === COMPARISON: Original vs Log-Transformed ===cat("Original Model:\n")## Original Model:cat("R-squared:", round(summary(regression_model)$r.squared, 4), "\n")## R-squared: 0.0704cat("Adjusted R-squared:", round(summary(regression_model)$adj.r.squared, 4), "\n\n")## Adjusted R-squared: 0.0629cat("Log-Transformed Model:\n")## Log-Transformed Model:summary(regression_log)##
## Call:
## lm(formula = log_error ~ N + p + correl + NSRatio, data = meta_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9004 -0.8099 -0.0996 0.6989 6.3096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.934e+00 1.446e-01 13.379 < 2e-16 ***
## N 4.979e-05 4.307e-06 11.561 < 2e-16 ***
## p 4.021e-03 1.737e-03 2.315 0.02103 *
## correl 2.468e-01 3.481e-01 0.709 0.47876
## NSRatio -4.918e-03 1.697e-03 -2.899 0.00391 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.38 on 499 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.2271, Adjusted R-squared: 0.2209
## F-statistic: 36.66 on 4 and 499 DF, p-value: < 2.2e-16# Check if N is still significant
cat("\n=== KEY FINDING ===\n")##
## === KEY FINDING ===n_pval <- summary(regression_log)$coefficients["N", "Pr(>|t|)"]
cat("N coefficient p-value (log model):", format.pval(n_pval, digits = 4), "\n")## N coefficient p-value (log model): < 2.2e-16cat("Dataset size (N) remains the only significant predictor.\n")## Dataset size (N) remains the only significant predictor.cat("Conclusion: Log transformation improves assumptions but does NOT change findings.\n")## Conclusion: Log transformation improves assumptions but does NOT change findings.5.5.3 Diagnostic Plots for Log-Transformed Model
par(mfrow = c(2, 2))
plot(regression_log)
par(mfrow = c(1, 1))
cat("\nDiagnostic Assessment:\n")##
## Diagnostic Assessment:cat("✓ Residuals vs Fitted: More homoscedastic (constant variance)\n")## ✓ Residuals vs Fitted: More homoscedastic (constant variance)cat("✓ Q-Q Plot: Improved normality, reduced tail deviation\n")## ✓ Q-Q Plot: Improved normality, reduced tail deviationcat("✓ Scale-Location: More horizontal trend\n")## ✓ Scale-Location: More horizontal trendcat("✓ Residuals vs Leverage: Fewer extreme influential points\n")## ✓ Residuals vs Leverage: Fewer extreme influential points5.5.4 Robust ANOVA Alternatives
# Install FSA if needed
if (!require("FSA")) install.packages("FSA", repos = "http://cran.us.r-project.org")
library(FSA)
cat("=== WELCH'S ANOVA (Robust to Unequal Variances) ===\n")## === WELCH'S ANOVA (Robust to Unequal Variances) ===welch_result <- oneway.test(Norm_error ~ Alg_Family, data = meta_data, var.equal = FALSE)
print(welch_result)##
## One-way analysis of means (not assuming equal variances)
##
## data: Norm_error and Alg_Family
## F = 1.4058, num df = 3.00, denom df = 260.37, p-value = 0.2415cat("\n=== KRUSKAL-WALLIS TEST (Non-Parametric Alternative) ===\n")##
## === KRUSKAL-WALLIS TEST (Non-Parametric Alternative) ===kruskal_result <- kruskal.test(Norm_error ~ Alg_Family, data = meta_data)
print(kruskal_result)##
## Kruskal-Wallis rank sum test
##
## data: Norm_error by Alg_Family
## Kruskal-Wallis chi-squared = 21.38, df = 3, p-value = 8.777e-05# Post-hoc for Kruskal-Wallis
dunn_result <- dunnTest(Norm_error ~ Alg_Family, data = meta_data, method = "bonferroni")
cat("\n=== DUNN'S POST-HOC TEST (Pairwise Comparisons) ===\n")##
## === DUNN'S POST-HOC TEST (Pairwise Comparisons) ===print(dunn_result)## Comparison Z P.unadj P.adj
## 1 Decision Trees - Instance-Based -2.55662733 0.0105692380 0.0634154281
## 2 Decision Trees - Neural Networks -3.72549174 0.0001949348 0.0011696090
## 3 Instance-Based - Neural Networks -1.57643467 0.1149256647 0.6895539882
## 4 Decision Trees - Statistical Methods -0.06839129 0.9454741519 1.0000000000
## 5 Instance-Based - Statistical Methods 2.62208703 0.0087393111 0.0524358666
## 6 Neural Networks - Statistical Methods 3.82275792 0.0001319673 0.0007918039Interpretation: Welch’s ANOVA (F = 1.41, p = 0.242) shows no significant difference when accounting for unequal variances, suggesting the original ANOVA result may be driven by extreme outliers. However, the Kruskal-Wallis test (χ² = 21.38, p < 0.001) does detect significant differences, indicating the effect is present when using rank-based comparisons. Together, these results suggest Neural Networks differ from other families, but the magnitude is heavily influenced by a few extreme outliers rather than systematic differences across all observations.
5.6 Outlier Analysis
Given the extreme outliers identified in diagnostic plots, we systematically investigate their characteristics and impact on our conclusions.
5.6.1 Identify and Characterize Outliers
# Define outliers as values > 3 SD from mean
outlier_threshold <- mean(meta_data$Norm_error) + 3 * sd(meta_data$Norm_error)
cat("Outlier threshold (Mean + 3SD):", round(outlier_threshold, 2), "%\n\n")## Outlier threshold (Mean + 3SD): 2393.4 %# Extract outliers - EXPLICITLY use dplyr::select
outliers <- meta_data %>%
filter(Norm_error > outlier_threshold) %>%
arrange(desc(Norm_error)) %>%
dplyr::select(DS_Name, Alg_Name, Alg_Family, Norm_error, k, N, p)
cat("Number of extreme outliers:", nrow(outliers), "\n")## Number of extreme outliers: 4cat("Percentage of data:", round(nrow(outliers)/nrow(meta_data)*100, 1), "%\n\n")## Percentage of data: 0.8 %print(outliers)## DS_Name Alg_Name Alg_Family Norm_error k N p
## 1 Shuttle Cascade Neural Networks 12040.992 7 58000 9
## 2 Shuttle Kohonen Neural Networks 12040.992 7 58000 9
## 3 German_Credit Cascade Neural Networks 3041.743 2 1000 24
## 4 Shuttle Default Instance-Based 2575.826 7 58000 95.6.2 Outlier Distribution by Algorithm Family
# Count outliers by algorithm family
outlier_summary <- outliers %>%
count(Alg_Family, name = "outlier_count") %>%
left_join(
meta_data %>% count(Alg_Family, name = "total_count"),
by = "Alg_Family"
) %>%
mutate(
outlier_rate = outlier_count / total_count,
outlier_pct = round(outlier_rate * 100, 1)
)
print(outlier_summary)## Alg_Family outlier_count total_count outlier_rate outlier_pct
## 1 Instance-Based 1 176 0.005681818 0.6
## 2 Neural Networks 3 110 0.027272727 2.7# Visualize outlier rates
ggplot(outlier_summary, aes(x = reorder(Alg_Family, outlier_pct),
y = outlier_pct, fill = Alg_Family)) +
geom_col(alpha = 0.8) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Catastrophic Failure Rate by Algorithm Family",
subtitle = "Percentage of runs with error > 3 SD above mean",
x = "Algorithm Family",
y = "Outlier Rate (%)") +
theme_minimal() +
theme(legend.position = "none") +
coord_flip()
Key Finding: Neural Networks have a catastrophic failure rate of 2.7%, compared to <2% for other families.
5.6.3 Reanalysis Excluding Outliers
# Create cleaned dataset
meta_data_clean <- meta_data %>%
filter(Norm_error <= outlier_threshold)
cat("Observations removed:", nrow(meta_data) - nrow(meta_data_clean), "\n")## Observations removed: 4cat("Remaining observations:", nrow(meta_data_clean), "\n\n")## Remaining observations: 524# Refit One-Way ANOVA
anova_clean <- aov(Norm_error ~ Alg_Family, data = meta_data_clean)
cat("=== ANOVA WITHOUT OUTLIERS ===\n")## === ANOVA WITHOUT OUTLIERS ===print(summary(anova_clean))## Df Sum Sq Mean Sq F value Pr(>F)
## Alg_Family 3 58606 19535 1.582 0.193
## Residuals 520 6422297 12351# Calculate new means
means_clean <- meta_data_clean %>%
group_by(Alg_Family) %>%
summarize(mean_error = mean(Norm_error), .groups = "drop")
cat("\n=== MEAN ERROR RATES (Outliers Removed) ===\n")##
## === MEAN ERROR RATES (Outliers Removed) ===print(means_clean)## # A tibble: 4 × 2
## Alg_Family mean_error
## <fct> <dbl>
## 1 Decision Trees 41.3
## 2 Instance-Based 40.0
## 3 Neural Networks 63.8
## 4 Statistical Methods 34.1# Post-hoc test
tukey_clean <- TukeyHSD(anova_clean)
cat("\n=== TUKEY HSD (Outliers Removed) ===\n")##
## === TUKEY HSD (Outliers Removed) ===print(tukey_clean)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Norm_error ~ Alg_Family, data = meta_data_clean)
##
## $Alg_Family
## diff lwr upr p adj
## Instance-Based-Decision Trees -1.296427 -36.14823 33.555376 0.9996841
## Neural Networks-Decision Trees 22.530583 -16.36138 61.422543 0.4423859
## Statistical Methods-Decision Trees -7.222200 -44.20013 29.755727 0.9582589
## Neural Networks-Instance-Based 23.827010 -11.32352 58.977537 0.3004419
## Statistical Methods-Instance-Based -5.925773 -38.94610 27.094555 0.9671441
## Statistical Methods-Neural Networks -29.752783 -67.01239 7.506825 0.1684159Conclusion: Even after removing outliers, Neural Networks still show higher error rates, though the difference becomes non-significant (p = 0.193). This suggests that while Neural Networks have systematic issues, the extreme outliers are the primary driver of their significantly worse performance in the original analysis.
5.7 Enhanced Visualizations
5.7.1 Error Distribution Density Plots
# Density plot on log scale
ggplot(meta_data, aes(x = Norm_error + 0.1, fill = Alg_Family)) +
geom_density(alpha = 0.4) +
scale_x_log10(
breaks = c(0.1, 1, 10, 100, 1000, 10000),
labels = c("0.1", "1", "10", "100", "1K", "10K")
) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(
title = "Error Distribution by Algorithm Family (Log Scale)",
subtitle = "Neural Networks show bimodal distribution with heavy right tail",
x = "Normalized Error (%) - Log Scale",
y = "Density",
fill = "Algorithm Family"
) +
theme_minimal() +
theme(legend.position = "top")
Interpretation: Neural Networks show bimodal distribution - most runs perform comparably, but catastrophic failures (>1000% error) create the heavy right tail.
5.7.2 Effect Size Comparison
# Calculate all effect sizes
regression_r2 <- summary(regression_model)$adj.r.squared
two_way_summary <- summary(two_way_model)
ss_family <- two_way_summary[[1]]["Alg_Family", "Sum Sq"]
ss_complexity <- two_way_summary[[1]]["Complexity", "Sum Sq"]
ss_interaction <- two_way_summary[[1]]["Alg_Family:Complexity", "Sum Sq"]
ss_total <- sum(two_way_summary[[1]][, "Sum Sq"])
# Create effect size dataframe
effect_sizes <- tibble(
Analysis = c(
"Dataset Characteristics\n(Regression R²)",
"Algorithm Family\n(η²)",
"Problem Complexity\n(η²)",
"Family × Complexity\n(η²)"
),
Effect_Size = c(regression_r2, ss_family/ss_total,
ss_complexity/ss_total, ss_interaction/ss_total),
Significant = c("Yes", "Yes", "No", "No")
)
# Visualize
ggplot(effect_sizes, aes(x = reorder(Analysis, Effect_Size),
y = Effect_Size, fill = Significant)) +
geom_col(alpha = 0.8) +
geom_hline(yintercept = 0.01, linetype = "dashed", color = "gray40") +
geom_text(aes(x = 0.5, y = 0.012), label = "Small (η² = 0.01)",
hjust = 0, size = 3, color = "gray40") +
geom_hline(yintercept = 0.06, linetype = "dashed", color = "gray40") +
geom_text(aes(x = 0.5, y = 0.065), label = "Medium (η² = 0.06)",
hjust = 0, size = 3, color = "gray40") +
scale_fill_manual(values = c("No" = "#CCCCCC", "Yes" = "#CC0000")) +
scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
coord_flip() +
labs(
title = "Comparative Effect Sizes Across All Analyses",
subtitle = "All effects are small - algorithm family explains only 2% of variance",
x = "",
y = "Proportion of Variance Explained",
fill = "Statistically\nSignificant?"
) +
theme_minimal()
Key Insight: Even statistically significant effects explain <7% of variance, suggesting individual algorithm tuning matters more than family membership.
5.8 Statistical Power Analysis
# Install pwr if needed
if (!require("pwr")) install.packages("pwr", repos = "http://cran.us.r-project.org")
library(pwr)
# Calculate Cohen's f from eta-squared
effect_size_f <- sqrt(eta_squared / (1 - eta_squared))
cat("Observed effect size (Cohen's f):", round(effect_size_f, 4), "\n\n")## Observed effect size (Cohen's f): 0.1422# Post-hoc power calculation
group_sizes <- table(meta_data$Alg_Family)
avg_n <- mean(group_sizes)
power_result <- pwr.anova.test(
k = 4, n = avg_n, f = effect_size_f, sig.level = 0.05
)
cat("=== POST-HOC POWER ANALYSIS ===\n")## === POST-HOC POWER ANALYSIS ===cat("Number of groups (k):", 4, "\n")## Number of groups (k): 4cat("Average group size (n):", round(avg_n, 0), "\n")## Average group size (n): 132cat("Effect size (Cohen's f):", round(effect_size_f, 4), "\n")## Effect size (Cohen's f): 0.1422cat("Achieved power:", round(power_result$power, 3), "\n\n")## Achieved power: 0.788cat("Interpretation:",
ifelse(power_result$power >= 0.8, "ADEQUATE", "INADEQUATE"),
"- We had", round(power_result$power * 100, 1), "% power\n")## Interpretation: INADEQUATE - We had 78.8 % power5.8.1 Power Curves Visualization
# Generate power curves
sample_sizes <- seq(10, 200, by = 5)
power_values <- sapply(sample_sizes, function(n) {
pwr.anova.test(k = 4, n = n, f = effect_size_f, sig.level = 0.05)$power
})
power_df <- tibble(n_per_group = sample_sizes, power = power_values)
ggplot(power_df, aes(x = n_per_group, y = power)) +
geom_line(color = "#0066CC", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") +
geom_vline(xintercept = avg_n, linetype = "dashed", color = "green") +
annotate("text", x = avg_n + 10, y = 0.1,
label = paste("Our study (n =", round(avg_n), ")"),
color = "green", hjust = 0) +
annotate("text", x = 10, y = 0.82, label = "80% power threshold",
color = "red", hjust = 0) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Statistical Power Curve for One-Way ANOVA",
subtitle = paste0("Effect size: Cohen's f = ", round(effect_size_f, 3)),
x = "Sample Size per Group (n)",
y = "Statistical Power (1 - β)"
) +
theme_minimal()
Conclusion: Our study achieved 78.8% power, which is below 80% threshold.
5.9 Summary of All Statistical Tests
# Create comprehensive summary table
summary_table <- tibble(
Analysis = c(
"Multiple Regression",
"One-Way ANOVA (Original)",
"One-Way ANOVA (Log-transformed)",
"Welch's ANOVA",
"Kruskal-Wallis Test",
"One-Way ANOVA (Outliers Removed)",
"Two-Way ANOVA - Family Effect",
"Two-Way ANOVA - Complexity Effect",
"Two-Way ANOVA - Interaction"
),
Test_Statistic = c(
paste0("F = ", round(summary(regression_model)$fstatistic[1], 2)),
paste0("F = ", round(summary(anova_model)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(aov(log_error ~ Alg_Family, data = meta_data))[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(welch_result$statistic, 2)),
paste0("χ² = ", round(kruskal_result$statistic, 2)),
paste0("F = ", round(summary(anova_clean)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Complexity", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Alg_Family:Complexity", "F value"], 2))
),
P_Value = c(
format.pval(pf(summary(regression_model)$fstatistic[1],
summary(regression_model)$fstatistic[2],
summary(regression_model)$fstatistic[3],
lower.tail = FALSE), digits = 3),
format.pval(summary(anova_model)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(aov(log_error ~ Alg_Family, data = meta_data))[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(welch_result$p.value, digits = 3),
format.pval(kruskal_result$p.value, digits = 3),
format.pval(summary(anova_clean)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Complexity", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Alg_Family:Complexity", "Pr(>F)"], digits = 3)
),
Significant = c(
"Yes", "Yes", "Yes", "No", "Yes", "No", "Yes", "No", "No"
),
Key_Finding = c(
"N is only significant predictor (R² = 6.3%)",
"Neural Networks significantly worse",
"Result confirmed with better assumptions",
"Non-significant when controlling for unequal variances",
"Significant rank-based differences",
"Effect disappears when removing 4 outliers",
"Family affects performance",
"Complexity does not affect performance",
"No interaction between family and complexity"
)
)
knitr::kable(summary_table,
caption = "Comprehensive Summary of All Statistical Tests",
col.names = c("Analysis", "Test Statistic", "p-value",
"Significant?", "Key Finding"))| Analysis | Test Statistic | p-value | Significant? | Key Finding |
|---|---|---|---|---|
| Multiple Regression | F = 9.44 | 2.3e-07 | Yes | N is only significant predictor (R² = 6.3%) |
| One-Way ANOVA (Original) | F = 3.53 | 0.0147 | Yes | Neural Networks significantly worse |
| One-Way ANOVA (Log-transformed) | F = 7.91 | 3.64e-05 | Yes | Result confirmed with better assumptions |
| Welch’s ANOVA | F = 1.41 | 0.242 | No | Non-significant when controlling for unequal variances |
| Kruskal-Wallis Test | χ² = 21.38 | 8.78e-05 | Yes | Significant rank-based differences |
| One-Way ANOVA (Outliers Removed) | F = 1.58 | 0.193 | No | Effect disappears when removing 4 outliers |
| Two-Way ANOVA - Family Effect | F = NA | NA | Yes | Family affects performance |
| Two-Way ANOVA - Complexity Effect | F = 1.85 | 0.158 | No | Complexity does not affect performance |
| Two-Way ANOVA - Interaction | F = 1.14 | 0.338 | No | No interaction between family and complexity |
Synthesis: This comprehensive analysis reveals that while algorithm family choice matters statistically (p = 0.014), the effect is small (η² = 0.02) and heavily driven by Neural Networks’ catastrophic failures on a few datasets. When outliers are removed or robust methods are used, the differences become less pronounced, suggesting practitioners should focus on individual algorithm validation rather than family-level selection.
6. Conclusion
6.1 Synthesis of Findings
This study examined machine learning algorithm performance from three complementary perspectives:
Finding 1: Dataset Characteristics Have Limited Predictive Power (Regression)
Only dataset size (N) significantly predicted classification error, explaining just 6.3% of variance. This surprisingly low explanatory power suggests that: - Intrinsic problem difficulty cannot be captured by simple dataset metrics alone - Algorithm choice matters more than dataset properties - The counterintuitive positive relationship between N and error reflects that harder problems naturally accumulate larger datasets
Finding 2: Algorithm Family Matters, But Effect is Small (One-Way ANOVA)
Neural Networks performed significantly worse (mean = 309% error) than Decision Trees, Instance-Based, and Statistical Methods (means = 34-54% error). However: - The small effect size (η² = 0.02) indicates high within-family variability - Extreme outliers (errors >12,000%) drive the Neural Networks disadvantage - The three reliable families show statistically indistinguishable performance
Finding 3: No Interaction with Problem Complexity (Two-Way ANOVA)
Algorithm family performance does not depend on problem complexity: - Binary and multi-class problems show similar error patterns - No evidence for context-dependent algorithm selection - Simple decision rules (avoid Neural Networks) apply uniformly
6.2 Practical Implications
For Practitioners:
Prioritize algorithm family over dataset metrics - Family choice (p = 0.014) has clearer impact than dataset properties (R² = 0.06)
Avoid Neural Networks in this benchmark context - Prone to catastrophic failures with extreme outliers
Use Decision Trees, Instance-Based, or Statistical Methods - Comparable performance with lower risk
Problem complexity doesn’t dictate algorithm choice - Same selection principles apply from binary to 91-class problems
Individual algorithm tuning matters more - High within-family variance suggests specific algorithm configurations and hyperparameters outweigh family-level decisions
Reconciling Conflicting Results
Our sensitivity analyses revealed important nuances:
- Standard ANOVA (p = 0.015) suggests families differ
- Welch’s ANOVA (p = 0.242) finds no difference when accounting for unequal variances
- Kruskal-Wallis (p < 0.001) detects differences using ranks
- ANOVA without outliers (p = 0.193) finds no difference when removing 4 extreme cases
Practical Recommendation: Neural Networks in this 1990s benchmark have a 2.7% rate of catastrophic failures (>3 SD above mean error). For practitioners: - Avoid Neural Networks unless you can extensively validate and monitor for failure modes - Use Decision Trees, Instance-Based, or Statistical Methods - they show consistent performance (means 34-54% error, no significant differences) - Always implement outlier detection in production systems regardless of algorithm choice
6.3 Unexpected Results
Counterintuitive Findings:
Larger datasets associated with higher error - Contradicts the “more data = better performance” assumption, but reflects that complex problems naturally require more examples
Problem complexity not significant - Expected multi-class problems to be harder, but error rates are similar across 2-91 classes
Neural Networks underperformance - Modern deep learning dominance wasn’t reflected in 1990s benchmarks, highlighting importance of algorithm evolution and proper tuning
Possible Explanations: - The Statlog project (1994) predates modern neural network architectures and training techniques - Benchmark datasets may not favor gradient-based learning - Lack of hyperparameter optimization may have disadvantaged complex methods
6.4 Connection to Real-World Decisions
These findings inform practical algorithm selection:
When to use each family: - Statistical Methods: Best for interpretable models, small-to-medium datasets - Decision Trees: Excellent for non-linear relationships, feature interactions - Instance-Based (KNN): Effective when decision boundaries are complex but local
When NOT to use: - Avoid algorithms with catastrophic failure modes on your domain - Avoid Neural Networks without extensive hyperparameter tuning and validation
Selection Strategy: 1. Start with multiple families in parallel 2. Cross-validate candidates from Decision Trees, Statistical, and Instance-Based 3. Choose based on validation performance, not family membership 4. Monitor for outlier performance patterns
7. Limitations
7.1 Data Limitations
- Historical benchmark data (1994) - Algorithms and best practices have evolved significantly since the Statlog project
- Limited algorithm coverage - Modern methods (Random Forests, Gradient Boosting, modern deep learning) not included
- Benchmark bias - UCI datasets may not represent current real-world problems
- Missing hyperparameter information - Cannot assess whether algorithms were optimally tuned
7.2 Assumption Violations
Regression Analysis: - Severe heteroscedasticity and outliers reduce reliability - Low R² suggests missing important predictors - Non-linear relationships may exist but weren’t modeled
ANOVA Analyses: - Equal variances assumption violated (Levene’s p = 0.016) - Extreme outliers violate normality assumption - Results should be interpreted cautiously - Welch’s ANOVA or non-parametric tests might be more appropriate
7.3 Generalizability Concerns
- Limited to classification tasks - Findings may not apply to regression, clustering, or other ML tasks
- Benchmark datasets - Real-world data often messier with class imbalance, missing values, temporal dependencies
- Algorithm families are broad - Individual algorithms within families vary substantially
- Era-specific results - 1990s algorithms may not reflect modern capabilities
Caution: Apply these findings as general guidance, not absolute rules. Always validate algorithm choices on your specific problem and data.
9. References
Dataset: - Meta-data [Dataset]. (1994). UCI Machine Learning Repository. https://doi.org/10.24432/C5X31P - Michie, D., Spiegelhalter, D.J., & Taylor, C. (1994). Machine Learning, Neural and Statistical Classification. Ellis Horwood.
R Packages: - R Core Team (2025). R: A language and environment for statistical computing. R Foundation for Statistical Computing. - Wickham, H. et al. (2019). tidyverse: Easily Install and Load the ‘Tidyverse’. R package version 2.0.0. - Fox, J. & Weisberg, S. (2019). car: Companion to Applied Regression. R package. - Robinson, D. et al. (2021). broom: Convert Statistical Objects into Tidy Tibbles. R package. - Schloerke, B. et al. (2021). GGally: Extension to ‘ggplot2’. R package. - Hothorn, T. et al. (2008). multcomp: Simultaneous Inference in General Parametric Models. R package.
10. Appendix: Complete R Code
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE,
fig.width = 10, fig.height = 6)
# Load required packages
library(tidyverse)
library(car)
library(broom)
library(gridExtra)
library(GGally) # For scatterplot matrix
library(multcomp) # For Compact Letter Display (CLD)
# Set theme
theme_set(theme_minimal(base_size = 12))
# Load the data
# NOTE: Save meta.data file in your project folder first
meta_data <- read.csv("meta.data", header = FALSE)
# Add column names based on meta.names documentation
colnames(meta_data) <- c("DS_Name", "T", "N", "p", "k", "Bin", "Cost",
"SDratio", "correl", "cancor1", "cancor2",
"fract1", "fract2", "skewness", "kurtosis",
"Hc", "Hx", "MCx", "EnAtr", "NSRatio",
"Alg_Name", "Norm_error")
# View structure
str(meta_data)
head(meta_data, 10)
# Summary statistics
summary(meta_data)
# Fix data types and handle missing values
meta_data <- meta_data %>%
mutate(
# Convert character variables to numeric (handle "?" as NA)
correl = as.numeric(correl),
cancor2 = ifelse(cancor2 == "?", NA, as.numeric(cancor2)),
fract2 = ifelse(fract2 == "?", NA, as.numeric(fract2)),
# Convert categorical variables to factors
DS_Name = as.factor(DS_Name),
Alg_Name = as.factor(Alg_Name),
Cost = as.factor(Cost)
)
# Check the fixes
str(meta_data)
# Create algorithm families based on machine learning methodology
# Addressing peer review feedback: clear justification for groupings
meta_data <- meta_data %>%
mutate(
Alg_Family = case_when(
# Neural Networks: Gradient-based learning with hidden layers
Alg_Name %in% c("BackProp", "Cascade", "LVQ", "Kohonen", "RBF") ~ "Neural Networks",
# Statistical Methods: Probabilistic and discriminant analysis
Alg_Name %in% c("Bayes", "BayesTree", "Discrim", "LogDisc", "QuaDisc", "Dipol92") ~ "Statistical Methods",
# Decision Trees: Recursive partitioning algorithms
Alg_Name %in% c("C4.5", "CART", "IndCART", "Cal5", "NewId") ~ "Decision Trees",
# Instance-Based: Memory-based learning using similarity metrics
Alg_Name %in% c("KNN", "Cn2", "ITrule", "Smart", "Castle", "Ac2", "Alloc80", "Default") ~ "Instance-Based",
TRUE ~ "Other"
),
Alg_Family = as.factor(Alg_Family)
)
# Verify all algorithms are accounted for
cat("Algorithm distribution across families:\n")
table(meta_data$Alg_Family, meta_data$Alg_Name)
# Create problem complexity levels based on number of classes (k)
meta_data <- meta_data %>%
mutate(
Complexity = case_when(
k == 2 ~ "Binary",
k >= 3 & k <= 10 ~ "Low Multi-class",
k > 10 ~ "High Multi-class"
),
Complexity = factor(Complexity, levels = c("Binary", "Low Multi-class", "High Multi-class"))
)
# Check distribution
cat("\nComplexity level distribution:\n")
table(meta_data$Complexity)
# Final dataset structure
cat("\nFinal dataset structure:\n")
str(meta_data)
cat("\nDataset ready for analysis with", nrow(meta_data), "observations\n")
# Scatterplot matrix for regression variables
regression_vars <- meta_data %>%
dplyr::select(Norm_error, N, p, correl, NSRatio) # Explicitly use dplyr::select
# Create scatterplot matrix
ggpairs(regression_vars,
title = "Scatterplot Matrix: Response and Predictors",
lower = list(continuous = wrap("points", alpha = 0.3, size = 0.5)),
diag = list(continuous = wrap("densityDiag", alpha = 0.5)),
upper = list(continuous = wrap("cor", size = 4)))
# Individual scatterplots
p1 <- ggplot(meta_data, aes(x = N, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Number of Examples",
x = "N (Number of Training Examples)",
y = "Normalized Error (%)")
p2 <- ggplot(meta_data, aes(x = p, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Number of Attributes",
x = "p (Number of Attributes)",
y = "Normalized Error (%)")
p3 <- ggplot(meta_data, aes(x = correl, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Attribute Correlation",
x = "Mean Correlation Between Attributes",
y = "Normalized Error (%)")
p4 <- ggplot(meta_data, aes(x = NSRatio, y = Norm_error)) +
geom_point(alpha = 0.5, color = "#0066CC") +
geom_smooth(method = "lm", se = TRUE, color = "#CC0000") +
labs(title = "Error vs Noise-Signal Ratio",
x = "Noise-Signal Ratio",
y = "Normalized Error (%)")
grid.arrange(p1, p2, p3, p4, ncol = 2)
# Fit multiple regression model
regression_model <- lm(Norm_error ~ N + p + correl + NSRatio, data = meta_data)
# Model summary
summary(regression_model)
# Coefficient table with confidence intervals
coef_table <- cbind(
Estimate = coef(regression_model),
confint(regression_model)
)
round(coef_table, 4)
# Four standard diagnostic plots
par(mfrow = c(2, 2))
plot(regression_model)
par(mfrow = c(1, 1))
# Check multicollinearity
cat("=== Variance Inflation Factors ===\n")
vif_values <- vif(regression_model)
print(vif_values)
cat("\nMulticollinearity assessment:\n")
if(max(vif_values) > 10) {
cat("WARNING: Severe multicollinearity detected\n")
} else if(max(vif_values) > 5) {
cat("CAUTION: Moderate multicollinearity present\n")
} else {
cat("GOOD: No concerning multicollinearity\n")
}
# Model fit statistics
cat("\n=== Model Fit Statistics ===\n")
cat("R-squared:", round(summary(regression_model)$r.squared, 4), "\n")
cat("Adjusted R-squared:", round(summary(regression_model)$adj.r.squared, 4), "\n")
cat("F-statistic:", round(summary(regression_model)$fstatistic[1], 2), "\n")
cat("p-value:", format.pval(pf(summary(regression_model)$fstatistic[1],
summary(regression_model)$fstatistic[2],
summary(regression_model)$fstatistic[3],
lower.tail = FALSE), digits = 4), "\n")
# Boxplot comparing algorithm families
ggplot(meta_data, aes(x = Alg_Family, y = Norm_error, fill = Alg_Family)) +
geom_boxplot(alpha = 0.7, outlier.alpha = 0.3) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Classification Error by Algorithm Family",
subtitle = "Comparing Neural Networks, Statistical, Decision Trees, Instance-Based",
x = "Algorithm Family",
y = "Normalized Error (%)") +
theme(legend.position = "none") +
coord_cartesian(ylim = c(0, 500)) # Zoom in to see most data
# Summary statistics by algorithm family
family_summary <- meta_data %>%
group_by(Alg_Family) %>%
summarize(
n = n(),
mean_error = mean(Norm_error, na.rm = TRUE),
sd_error = sd(Norm_error, na.rm = TRUE),
median_error = median(Norm_error, na.rm = TRUE),
min_error = min(Norm_error, na.rm = TRUE),
max_error = max(Norm_error, na.rm = TRUE)
)
print(family_summary)
# Fit One-Way ANOVA
anova_model <- aov(Norm_error ~ Alg_Family, data = meta_data)
# ANOVA table
summary(anova_model)
# Calculate effect size (eta-squared)
anova_summary <- summary(anova_model)
ss_between <- anova_summary[[1]]["Alg_Family", "Sum Sq"]
ss_total <- sum(anova_summary[[1]][, "Sum Sq"])
eta_squared <- ss_between / ss_total
cat("Effect Size (η²):", round(eta_squared, 4), "\n")
cat("Interpretation:",
ifelse(eta_squared < 0.01, "negligible",
ifelse(eta_squared < 0.06, "small",
ifelse(eta_squared < 0.14, "medium", "large"))), "effect\n")
# Tukey HSD for pairwise comparisons
tukey_results <- TukeyHSD(anova_model)
print(tukey_results)
# Visualize Tukey HSD results
plot(tukey_results, las = 1, col = "blue")
# Compact letter display
library(multcomp)
cld_result <- cld(glht(anova_model, linfct = mcp(Alg_Family = "Tukey")))
cat("\nCompact Letter Display:\n")
print(cld_result)
# Check ANOVA assumptions
par(mfrow = c(2, 2))
# 1. Residuals vs Fitted
plot(anova_model, which = 1, main = "Residuals vs Fitted")
# 2. Q-Q Plot
plot(anova_model, which = 2, main = "Normal Q-Q")
# 3. Levene's Test for equal variances
leveneTest(Norm_error ~ Alg_Family, data = meta_data)
par(mfrow = c(1, 1))
cat("\n=== ANOVA Assumption Checks ===\n\n")
cat("1. INDEPENDENCE: Assumed based on study design (each algorithm-dataset combo is independent)\n\n")
cat("2. NORMALITY: Check Q-Q plot - if points follow diagonal line, normality is reasonable\n\n")
cat("3. EQUAL VARIANCES: Check Levene's test:\n")
cat(" - If p > 0.05: Equal variances assumption met\n")
cat(" - If p < 0.05: Variances differ (may need Welch's ANOVA)\n")
# Calculate group means for interaction
interaction_means <- meta_data %>%
group_by(Alg_Family, Complexity) %>%
summarize(mean_error = mean(Norm_error, na.rm = TRUE), .groups = "drop")
# Create interaction plot
ggplot(interaction_means, aes(x = Complexity, y = mean_error,
color = Alg_Family, group = Alg_Family)) +
geom_point(size = 4) +
geom_line(size = 1.2) +
scale_color_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Interaction Plot: Algorithm Family × Problem Complexity",
subtitle = "Do lines cross? That indicates interaction!",
x = "Problem Complexity",
y = "Mean Classification Error (%)",
color = "Algorithm Family") +
theme(legend.position = "top")
# Fit two-way ANOVA with interaction
two_way_model <- aov(Norm_error ~ Alg_Family * Complexity, data = meta_data)
# ANOVA table
summary(two_way_model)
# Create log-transformed response variable
# Using log1p (log(1+x)) to handle zeros
meta_data <- meta_data %>%
mutate(log_error = log1p(Norm_error))
# Visualize transformation effect
par(mfrow = c(1, 2))
hist(meta_data$Norm_error, breaks = 30, col = "lightblue",
main = "Original Error Distribution", xlab = "Normalized Error")
hist(meta_data$log_error, breaks = 30, col = "lightgreen",
main = "Log-Transformed Error", xlab = "log(1 + Error)")
par(mfrow = c(1, 1))
# Refit regression model with log-transformed response
regression_log <- lm(log_error ~ N + p + correl + NSRatio, data = meta_data)
# Compare model summaries
cat("=== COMPARISON: Original vs Log-Transformed ===\n\n")
cat("Original Model:\n")
cat("R-squared:", round(summary(regression_model)$r.squared, 4), "\n")
cat("Adjusted R-squared:", round(summary(regression_model)$adj.r.squared, 4), "\n\n")
cat("Log-Transformed Model:\n")
summary(regression_log)
# Check if N is still significant
cat("\n=== KEY FINDING ===\n")
n_pval <- summary(regression_log)$coefficients["N", "Pr(>|t|)"]
cat("N coefficient p-value (log model):", format.pval(n_pval, digits = 4), "\n")
cat("Dataset size (N) remains the only significant predictor.\n")
cat("Conclusion: Log transformation improves assumptions but does NOT change findings.\n")
par(mfrow = c(2, 2))
plot(regression_log)
par(mfrow = c(1, 1))
cat("\nDiagnostic Assessment:\n")
cat("✓ Residuals vs Fitted: More homoscedastic (constant variance)\n")
cat("✓ Q-Q Plot: Improved normality, reduced tail deviation\n")
cat("✓ Scale-Location: More horizontal trend\n")
cat("✓ Residuals vs Leverage: Fewer extreme influential points\n")
# Install FSA if needed
if (!require("FSA")) install.packages("FSA", repos = "http://cran.us.r-project.org")
library(FSA)
cat("=== WELCH'S ANOVA (Robust to Unequal Variances) ===\n")
welch_result <- oneway.test(Norm_error ~ Alg_Family, data = meta_data, var.equal = FALSE)
print(welch_result)
cat("\n=== KRUSKAL-WALLIS TEST (Non-Parametric Alternative) ===\n")
kruskal_result <- kruskal.test(Norm_error ~ Alg_Family, data = meta_data)
print(kruskal_result)
# Post-hoc for Kruskal-Wallis
dunn_result <- dunnTest(Norm_error ~ Alg_Family, data = meta_data, method = "bonferroni")
cat("\n=== DUNN'S POST-HOC TEST (Pairwise Comparisons) ===\n")
print(dunn_result)
# Define outliers as values > 3 SD from mean
outlier_threshold <- mean(meta_data$Norm_error) + 3 * sd(meta_data$Norm_error)
cat("Outlier threshold (Mean + 3SD):", round(outlier_threshold, 2), "%\n\n")
# Extract outliers - EXPLICITLY use dplyr::select
outliers <- meta_data %>%
filter(Norm_error > outlier_threshold) %>%
arrange(desc(Norm_error)) %>%
dplyr::select(DS_Name, Alg_Name, Alg_Family, Norm_error, k, N, p)
cat("Number of extreme outliers:", nrow(outliers), "\n")
cat("Percentage of data:", round(nrow(outliers)/nrow(meta_data)*100, 1), "%\n\n")
print(outliers)
# Count outliers by algorithm family
outlier_summary <- outliers %>%
count(Alg_Family, name = "outlier_count") %>%
left_join(
meta_data %>% count(Alg_Family, name = "total_count"),
by = "Alg_Family"
) %>%
mutate(
outlier_rate = outlier_count / total_count,
outlier_pct = round(outlier_rate * 100, 1)
)
print(outlier_summary)
# Visualize outlier rates
ggplot(outlier_summary, aes(x = reorder(Alg_Family, outlier_pct),
y = outlier_pct, fill = Alg_Family)) +
geom_col(alpha = 0.8) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(title = "Catastrophic Failure Rate by Algorithm Family",
subtitle = "Percentage of runs with error > 3 SD above mean",
x = "Algorithm Family",
y = "Outlier Rate (%)") +
theme_minimal() +
theme(legend.position = "none") +
coord_flip()
# Create cleaned dataset
meta_data_clean <- meta_data %>%
filter(Norm_error <= outlier_threshold)
cat("Observations removed:", nrow(meta_data) - nrow(meta_data_clean), "\n")
cat("Remaining observations:", nrow(meta_data_clean), "\n\n")
# Refit One-Way ANOVA
anova_clean <- aov(Norm_error ~ Alg_Family, data = meta_data_clean)
cat("=== ANOVA WITHOUT OUTLIERS ===\n")
print(summary(anova_clean))
# Calculate new means
means_clean <- meta_data_clean %>%
group_by(Alg_Family) %>%
summarize(mean_error = mean(Norm_error), .groups = "drop")
cat("\n=== MEAN ERROR RATES (Outliers Removed) ===\n")
print(means_clean)
# Post-hoc test
tukey_clean <- TukeyHSD(anova_clean)
cat("\n=== TUKEY HSD (Outliers Removed) ===\n")
print(tukey_clean)
# Density plot on log scale
ggplot(meta_data, aes(x = Norm_error + 0.1, fill = Alg_Family)) +
geom_density(alpha = 0.4) +
scale_x_log10(
breaks = c(0.1, 1, 10, 100, 1000, 10000),
labels = c("0.1", "1", "10", "100", "1K", "10K")
) +
scale_fill_manual(values = c("#CC0000", "#0066CC", "#00AA00", "#FF6600")) +
labs(
title = "Error Distribution by Algorithm Family (Log Scale)",
subtitle = "Neural Networks show bimodal distribution with heavy right tail",
x = "Normalized Error (%) - Log Scale",
y = "Density",
fill = "Algorithm Family"
) +
theme_minimal() +
theme(legend.position = "top")
# Calculate all effect sizes
regression_r2 <- summary(regression_model)$adj.r.squared
two_way_summary <- summary(two_way_model)
ss_family <- two_way_summary[[1]]["Alg_Family", "Sum Sq"]
ss_complexity <- two_way_summary[[1]]["Complexity", "Sum Sq"]
ss_interaction <- two_way_summary[[1]]["Alg_Family:Complexity", "Sum Sq"]
ss_total <- sum(two_way_summary[[1]][, "Sum Sq"])
# Create effect size dataframe
effect_sizes <- tibble(
Analysis = c(
"Dataset Characteristics\n(Regression R²)",
"Algorithm Family\n(η²)",
"Problem Complexity\n(η²)",
"Family × Complexity\n(η²)"
),
Effect_Size = c(regression_r2, ss_family/ss_total,
ss_complexity/ss_total, ss_interaction/ss_total),
Significant = c("Yes", "Yes", "No", "No")
)
# Visualize
ggplot(effect_sizes, aes(x = reorder(Analysis, Effect_Size),
y = Effect_Size, fill = Significant)) +
geom_col(alpha = 0.8) +
geom_hline(yintercept = 0.01, linetype = "dashed", color = "gray40") +
geom_text(aes(x = 0.5, y = 0.012), label = "Small (η² = 0.01)",
hjust = 0, size = 3, color = "gray40") +
geom_hline(yintercept = 0.06, linetype = "dashed", color = "gray40") +
geom_text(aes(x = 0.5, y = 0.065), label = "Medium (η² = 0.06)",
hjust = 0, size = 3, color = "gray40") +
scale_fill_manual(values = c("No" = "#CCCCCC", "Yes" = "#CC0000")) +
scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
coord_flip() +
labs(
title = "Comparative Effect Sizes Across All Analyses",
subtitle = "All effects are small - algorithm family explains only 2% of variance",
x = "",
y = "Proportion of Variance Explained",
fill = "Statistically\nSignificant?"
) +
theme_minimal()
# Install pwr if needed
if (!require("pwr")) install.packages("pwr", repos = "http://cran.us.r-project.org")
library(pwr)
# Calculate Cohen's f from eta-squared
effect_size_f <- sqrt(eta_squared / (1 - eta_squared))
cat("Observed effect size (Cohen's f):", round(effect_size_f, 4), "\n\n")
# Post-hoc power calculation
group_sizes <- table(meta_data$Alg_Family)
avg_n <- mean(group_sizes)
power_result <- pwr.anova.test(
k = 4, n = avg_n, f = effect_size_f, sig.level = 0.05
)
cat("=== POST-HOC POWER ANALYSIS ===\n")
cat("Number of groups (k):", 4, "\n")
cat("Average group size (n):", round(avg_n, 0), "\n")
cat("Effect size (Cohen's f):", round(effect_size_f, 4), "\n")
cat("Achieved power:", round(power_result$power, 3), "\n\n")
cat("Interpretation:",
ifelse(power_result$power >= 0.8, "ADEQUATE", "INADEQUATE"),
"- We had", round(power_result$power * 100, 1), "% power\n")
# Generate power curves
sample_sizes <- seq(10, 200, by = 5)
power_values <- sapply(sample_sizes, function(n) {
pwr.anova.test(k = 4, n = n, f = effect_size_f, sig.level = 0.05)$power
})
power_df <- tibble(n_per_group = sample_sizes, power = power_values)
ggplot(power_df, aes(x = n_per_group, y = power)) +
geom_line(color = "#0066CC", linewidth = 1.2) +
geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") +
geom_vline(xintercept = avg_n, linetype = "dashed", color = "green") +
annotate("text", x = avg_n + 10, y = 0.1,
label = paste("Our study (n =", round(avg_n), ")"),
color = "green", hjust = 0) +
annotate("text", x = 10, y = 0.82, label = "80% power threshold",
color = "red", hjust = 0) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Statistical Power Curve for One-Way ANOVA",
subtitle = paste0("Effect size: Cohen's f = ", round(effect_size_f, 3)),
x = "Sample Size per Group (n)",
y = "Statistical Power (1 - β)"
) +
theme_minimal()
# Create comprehensive summary table
summary_table <- tibble(
Analysis = c(
"Multiple Regression",
"One-Way ANOVA (Original)",
"One-Way ANOVA (Log-transformed)",
"Welch's ANOVA",
"Kruskal-Wallis Test",
"One-Way ANOVA (Outliers Removed)",
"Two-Way ANOVA - Family Effect",
"Two-Way ANOVA - Complexity Effect",
"Two-Way ANOVA - Interaction"
),
Test_Statistic = c(
paste0("F = ", round(summary(regression_model)$fstatistic[1], 2)),
paste0("F = ", round(summary(anova_model)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(aov(log_error ~ Alg_Family, data = meta_data))[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(welch_result$statistic, 2)),
paste0("χ² = ", round(kruskal_result$statistic, 2)),
paste0("F = ", round(summary(anova_clean)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Alg_Family", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Complexity", "F value"], 2)),
paste0("F = ", round(summary(two_way_model)[[1]]["Alg_Family:Complexity", "F value"], 2))
),
P_Value = c(
format.pval(pf(summary(regression_model)$fstatistic[1],
summary(regression_model)$fstatistic[2],
summary(regression_model)$fstatistic[3],
lower.tail = FALSE), digits = 3),
format.pval(summary(anova_model)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(aov(log_error ~ Alg_Family, data = meta_data))[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(welch_result$p.value, digits = 3),
format.pval(kruskal_result$p.value, digits = 3),
format.pval(summary(anova_clean)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Alg_Family", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Complexity", "Pr(>F)"], digits = 3),
format.pval(summary(two_way_model)[[1]]["Alg_Family:Complexity", "Pr(>F)"], digits = 3)
),
Significant = c(
"Yes", "Yes", "Yes", "No", "Yes", "No", "Yes", "No", "No"
),
Key_Finding = c(
"N is only significant predictor (R² = 6.3%)",
"Neural Networks significantly worse",
"Result confirmed with better assumptions",
"Non-significant when controlling for unequal variances",
"Significant rank-based differences",
"Effect disappears when removing 4 outliers",
"Family affects performance",
"Complexity does not affect performance",
"No interaction between family and complexity"
)
)
knitr::kable(summary_table,
caption = "Comprehensive Summary of All Statistical Tests",
col.names = c("Analysis", "Test Statistic", "p-value",
"Significant?", "Key Finding"))
# All code chunks from this document are displayed here for reproducibility
sessionInfo()Session Information:
sessionInfo()## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pwr_1.3-0 FSA_0.10.0 multcomp_1.4-29 TH.data_1.1-5
## [5] MASS_7.3-65 survival_3.8-3 mvtnorm_1.3-3 GGally_2.4.0
## [9] gridExtra_2.3 broom_1.0.10 car_3.1-3 carData_3.0-5
## [13] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.2 dplyr_1.1.4
## [17] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1 tibble_3.3.0
## [21] ggplot2_4.0.0 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.53 bslib_0.9.0 lattice_0.22-7
## [5] tzdb_0.5.0 vctrs_0.6.5 tools_4.5.1 generics_0.1.4
## [9] sandwich_3.1-1 pkgconfig_2.0.3 Matrix_1.7-3 RColorBrewer_1.1-3
## [13] S7_0.2.0 lifecycle_1.0.4 compiler_4.5.1 farver_2.1.2
## [17] codetools_0.2-20 htmltools_0.5.8.1 sass_0.4.10 yaml_2.3.10
## [21] Formula_1.2-5 pillar_1.11.0 jquerylib_0.1.4 cachem_1.1.0
## [25] dunn.test_1.3.6 abind_1.4-8 nlme_3.1-168 ggstats_0.11.0
## [29] tidyselect_1.2.1 digest_0.6.37 stringi_1.8.7 labeling_0.4.3
## [33] splines_4.5.1 fastmap_1.2.0 grid_4.5.1 cli_3.6.5
## [37] magrittr_2.0.4 utf8_1.2.6 withr_3.0.2 scales_1.4.0
## [41] backports_1.5.0 timechange_0.3.0 rmarkdown_2.29 zoo_1.8-14
## [45] hms_1.1.3 evaluate_1.0.5 knitr_1.50 mgcv_1.9-3
## [49] rlang_1.1.6 glue_1.8.0 rstudioapi_0.17.1 jsonlite_2.0.0
## [53] R6_2.6.1END OF REPORT