Multivariate Analysis of Variance (MANOVA)

Analyzing Iris Flower Morphological Characteristics by Species

1 Introduction

This document demonstrates a complete Multivariate Analysis of Variance (MANOVA) using the classic iris dataset. MANOVA is used when we have multiple continuous dependent variables and one or more categorical independent variables. It tests whether there are differences between groups on a linear combination of dependent variables.

1.1 Research Question

Do the three iris species (setosa, versicolor, virginica) differ significantly in their overall morphological profile across multiple floral characteristics?

2 Required Libraries

The following R packages are required for this analysis:

• car: For MANOVA testing using Manova() function
• biotools: For Box’s M test of covariance matrix homogeneity
• mvnormtest: For multivariate normality testing
• ggplot2: For data visualization
• dplyr: For data manipulation
• tidyr: For data reshaping
• knitr: For formatted table output

library(car)
library(biotools)
library(mvnormtest)
library(ggplot2)
library(dplyr)
library(tidyr)
library(knitr)

library(afcommon)
source("af_common_add_ons.R")

# Load iris data and define basic variables for inline use
data(iris)
n_total <- nrow(iris)
n_species <- length(unique(iris$Species))
n_vars_check <- 4

3 Dataset Overview

The iris dataset contains measurements of 150 iris flowers from 3 species. Each flower has four morphological measurements:

# Basic structure
str(iris)
'data.frame':   150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

# Calculate additional summary statistics for inline use
species_counts <- table(iris$Species)

# Summary statistics by species
summary_stats <- iris %>%
  group_by(Species) %>%
  summarise(
    n = n(),
    mean_sepal_length = round(mean(Sepal.Length), 2),
    mean_sepal_width = round(mean(Sepal.Width), 2),
    mean_petal_length = round(mean(Petal.Length), 2),
    mean_petal_width = round(mean(Petal.Width), 2),
    .groups = 'drop'
  )

summary_stats %>%
  kable(caption = "Descriptive Statistics by Species")

Descriptive Statistics by Species

Species	n	mean_sepal_length	mean_sepal_width	mean_petal_length	mean_petal_width
setosa	50	5.01	3.43	1.46	0.25
versicolor	50	5.94	2.77	4.26	1.33
virginica	50	6.59	2.97	5.55	2.03

3.1 Data Visualization

# Create a pairs plot to visualize relationships
pairs(iris[1:4], 
      col = iris$Species, 
      pch = 19,
      main = "Iris Measurements by Species",
      labels = c("Sepal\nLength", "Sepal\nWidth", "Petal\nLength", "Petal\nWidth"))

# Add legend
legend("top", 
       legend = levels(iris$Species),
       col = 1:3,
       pch = 19,
       ncol = 3,
       inset = c(0, -0.05),
       xpd = TRUE)

4 MANOVA Analysis

4.1 Model Specification

We fit a MANOVA model using the lm() function, treating the four morphological variables as a multivariate response and species as the independent variable.

# Create the dependent variable matrix
Y <- as.matrix(iris[, 1:4])
colnames(Y) <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")

# Fit the MANOVA model using lm()
# The cbind() function combines multiple dependent variables
manova_model <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, 
                   data = iris)

# Display model summary
summary(manova_model)
Response Sepal.Length :

Call:
lm(formula = Sepal.Length ~ Species, data = iris)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6880 -0.3285 -0.0060  0.3120  1.3120 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.0060     0.0728  68.762  < 2e-16 ***
Speciesversicolor   0.9300     0.1030   9.033 8.77e-16 ***
Speciesvirginica    1.5820     0.1030  15.366  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5148 on 147 degrees of freedom
Multiple R-squared:  0.6187,    Adjusted R-squared:  0.6135 
F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16


Response Sepal.Width :

Call:
lm(formula = Sepal.Width ~ Species, data = iris)

Residuals:
   Min     1Q Median     3Q    Max 
-1.128 -0.228  0.026  0.226  0.972 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.42800    0.04804  71.359  < 2e-16 ***
Speciesversicolor -0.65800    0.06794  -9.685  < 2e-16 ***
Speciesvirginica  -0.45400    0.06794  -6.683 4.54e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3397 on 147 degrees of freedom
Multiple R-squared:  0.4008,    Adjusted R-squared:  0.3926 
F-statistic: 49.16 on 2 and 147 DF,  p-value: < 2.2e-16


Response Petal.Length :

Call:
lm(formula = Petal.Length ~ Species, data = iris)

Residuals:
   Min     1Q Median     3Q    Max 
-1.260 -0.258  0.038  0.240  1.348 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.46200    0.06086   24.02   <2e-16 ***
Speciesversicolor  2.79800    0.08607   32.51   <2e-16 ***
Speciesvirginica   4.09000    0.08607   47.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4303 on 147 degrees of freedom
Multiple R-squared:  0.9414,    Adjusted R-squared:  0.9406 
F-statistic:  1180 on 2 and 147 DF,  p-value: < 2.2e-16


Response Petal.Width :

Call:
lm(formula = Petal.Width ~ Species, data = iris)

Residuals:
   Min     1Q Median     3Q    Max 
-0.626 -0.126 -0.026  0.154  0.474 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.24600    0.02894    8.50 1.96e-14 ***
Speciesversicolor  1.08000    0.04093   26.39  < 2e-16 ***
Speciesvirginica   1.78000    0.04093   43.49  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2047 on 147 degrees of freedom
Multiple R-squared:  0.9289,    Adjusted R-squared:  0.9279 
F-statistic:   960 on 2 and 147 DF,  p-value: < 2.2e-16

4.2 MANOVA Test Results

We use the car::Manova() function to obtain multivariate test statistics. By default, it shows only Pillai’s trace, but we can access all four test statistics by specifying the test parameter or examining the summary.

# Perform MANOVA using car package - default shows Pillai's trace 
manova_results <- Manova(manova_model, type = "III")

# Pillai's trace (default)
pillai_result <- Manova(manova_model, type = "III", test = "Pillai")
af_cat(pillai_result)

Type III MANOVA Tests: Pillai test statistic
            Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)  1   0.98304  2086.77      4    144 < 2.2e-16 ***
Species      2   1.19190    53.47      8    290 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Wilks' Lambda
wilks_result <- Manova(manova_model, type = "III", test = "Wilks")
af_cat(wilks_result)

Type III MANOVA Tests: Wilks test statistic
            Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)  1  0.016959  2086.77      4    144 < 2.2e-16 ***
Species      2  0.023439   199.15      8    288 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Hotelling-Lawley trace
hotelling_result <- Manova(manova_model, type = "III", test = "Hotelling-Lawley")
af_cat(hotelling_result)

Type III MANOVA Tests: Hotelling-Lawley test statistic
            Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)  1    57.966  2086.77      4    144 < 2.2e-16 ***
Species      2    32.477   580.53      8    286 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Roy's largest root
roy_result <- Manova(manova_model, type = "III", test = "Roy")
af_cat(roy_result)

Type III MANOVA Tests: Roy test statistic
            Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)  1    57.966   2086.8      4    144 < 2.2e-16 ***
Species      2    32.192   1167.0      4    145 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# COMPREHENSIVE RESULTS EXTRACTION FOR INLINE USE THROUGHOUT DOCUMENT
output <- capture.output(summary(manova_results))

# Find and extract the test statistic lines
species_start <- grep("Multivariate Tests: Species", output)
test_stat_lines <- output[(species_start + 2):(species_start + 5)]

# Extract values using regex
pillai_value <- as.numeric(sub("Pillai\\s+\\d+\\s+([0-9.]+)\\s+.*", "\\1", test_stat_lines[1]))
wilks_value <- as.numeric(sub("Wilks\\s+\\d+\\s+([0-9.]+)\\s+.*", "\\1", test_stat_lines[2]))
hotelling_value <- as.numeric(sub("Hotelling-Lawley\\s+\\d+\\s+([0-9.]+)\\s+.*", "\\1", test_stat_lines[3]))
roy_value <- as.numeric(sub("Roy\\s+\\d+\\s+([0-9.]+)\\s+.*", "\\1", test_stat_lines[4]))
p_value <- as.numeric(sub(".*<\\s+([0-9.e-]+).*", "\\1", test_stat_lines[1]))

# Univariate ANOVA results
univariate_summary <- summary.aov(manova_model)

# Extract F-statistics and degrees of freedom
f_stats <- sapply(1:4, function(i) {
  univariate_summary[[i]][["F value"]][1]
})
df_between <- sapply(1:4, function(i) {
  univariate_summary[[i]][["Df"]][1]
})
df_within <- sapply(1:4, function(i) {
  univariate_summary[[i]][["Df"]][2]
})
names(f_stats) <- names(df_between) <- names(df_within) <- colnames(Y)

# Effect sizes (eta-squared)
eta_squared <- sapply(1:4, function(i) {
  ss_effect <- univariate_summary[[i]][["Sum Sq"]][1]
  ss_total <- sum(univariate_summary[[i]][["Sum Sq"]])
  ss_effect / ss_total
})
names(eta_squared) <- colnames(Y)

# Store key effect size variables
max_eta_var <- names(which.max(eta_squared))
max_eta_value <- max(eta_squared)
min_eta_var <- names(which.min(eta_squared))
min_eta_value <- min(eta_squared)

# Sample size variables
min_group_size <- min(table(iris$Species))
recommended_min <- n_species * n_vars_check + 10
sample_adequate <- min_group_size > n_vars_check && n_total > recommended_min

# Assumption test results
# Multivariate normality
mshapiro_test <- mshapiro.test(t(Y))
mv_normal <- mshapiro_test$p.value > 0.05

# Box's M test
box_m_result <- boxM(Y, iris$Species)
box_m_stat <- box_m_result$statistic
box_m_p <- box_m_result$p.value
covariance_homogeneous <- box_m_p > 0.05

# Linearity assessment
cor_matrix <- cor(Y)
correlations <- cor_matrix[upper.tri(cor_matrix)]
min_correlation <- min(abs(correlations))
max_correlation <- max(abs(correlations))
linearity_adequate <- min_correlation > 0.1 && max_correlation < 0.9

# Multicollinearity
det_cor <- det(cor_matrix)
multicollinearity_ok <- det_cor > 0.00001

# Outliers
center <- colMeans(Y)
cov_matrix <- cov(Y)
mahal_dist <- mahalanobis(Y, center, cov_matrix)
critical_value <- qchisq(0.99, df = n_vars_check)
outliers <- sum(mahal_dist > critical_value)
outlier_proportion <- outliers / nrow(Y)
outliers_acceptable <- outlier_proportion < 0.05

# Species names for profile plot
species_names <- levels(iris$Species)

4.2.1 Interpretation of MANOVA Results

The MANOVA results show four different test statistics, each with specific properties:

1. Pillai’s Trace = 1.192:

Pillai’s Trace (also called Pillai-Bartlett Trace) is a multivariate test statistic that measures the proportion of variance in the combination of dependent variables that is explained by the independent variable(s). It’s calculated as the sum of the squared canonical correlations.

Interpretation:

• Range: 0 to 1 (or can exceed 1 in some cases with multiple groups)
• Value near 0: Little to no group differences (independent variable explains little variance)
• Value near 1: Strong group differences (independent variable explains much variance)

Guidelines:

• 0.00-0.04: Very small effect (negligible group differences)
• 0.04-0.14: Small effect (weak group differences)
• 0.14-0.36: Medium effect (moderate group differences)
• 0.36-1.00: Large effect (strong group differences)

Key Characteristics:

• Most robust to violations of assumptions (especially unequal covariance matrices)
• Conservative - less likely to find false positives
• Recommended when assumptions are questionable or sample sizes are unequal

2. Wilks’ Lambda = 0.02344:

Wilks’ Lambda represents the proportion of variance in the combination of dependent variables that is NOT explained by the independent variable(s). It’s calculated as the ratio of within-group variance to total variance.

Interpretation:

• Range: 0 to 1
• Value near 1: Little to no group differences (most variance is within groups)
• Value near 0: Strong group differences (most variance is between groups)

Guidelines:

• 0.90-1.00: Very small effect (negligible group differences)
• 0.70-0.89: Small effect (weak group differences)
• 0.50-0.69: Medium effect (moderate group differences)
• 0.00-0.49: Large effect (strong group differences)

Key Characteristics:

• Most commonly reported in literature (historical preference)
• Assumes equal covariance matrices across groups
• Less robust than Pillai’s Trace when assumptions are violated
• More powerful than Pillai’s Trace when assumptions are perfectly met

3. Hotelling-Lawley Trace = 32.477:

Hotelling-Lawley Trace represents the ratio of between-group variance to within-group variance across all discriminant functions. It’s calculated as the sum of the eigenvalues from the hypothesis matrix relative to the error matrix.

Interpretation:

Range: 0 to ∞ (theoretically unbounded) Value near 0: Little to no group differences Larger values: Stronger group differences

Guidelines:

• 0.00-0.10: Very small effect
• 0.10-0.25: Small effect
• 0.25-0.50: Medium effect
• 0.50+: Large effect

Key Characteristics:

• More powerful than both Pillai’s Trace and Wilks’ Lambda when assumptions are perfectly met
• Sensitive to violations of assumptions (especially outliers and unequal sample sizes)
• Can become unstable with small sample sizes or many dependent variables
• Not bounded between 0 and 1 (can exceed 1)

When to Use:

• When assumptions are perfectly met
• With large, balanced sample sizes
• When seeking maximum statistical power
• When you want to emphasize group separation magnitude

4. Roy’s Largest Root = 32.192:

Roy’s Largest Root represents the proportion of variance explained by the first and strongest discriminant function (canonical variate). It focuses only on the largest eigenvalue, capturing the maximum possible separation between groups along one dimension.

Interpretation:

• Range: 0 to 1 (or can exceed 1 in some cases)
• Value near 0: Little to no group differences
• Value near 1: Very strong group differences along one primary dimension

Guidelines:

• 0.00-0.04: Very small effect
• 0.04-0.14: Small effect
• 0.14-0.36: Medium effect
• 0.36-1.00: Large effect

Key Characteristics:

• Most powerful test when group differences exist primarily along one dimension
• Least robust to violations of assumptions
• Most liberal test - most likely to find significant results
• Focuses only on the strongest discriminant function (ignores others)
• Can be misleading if group differences are multidimensional

When to Use:

• When you have strong theoretical reasons to expect one primary dimension of group differences
• When you want the most powerful test for detecting any group differences
• In exploratory analyses where you want to maximize sensitivity
• When other tests show borderline significance

Important Caveat:

Roy’s Largest Root can be overly optimistic if group differences are actually spread across multiple dimensions. It may detect significance even when other more conservative tests do not, but this doesn’t necessarily mean the effect is stronger - just that it’s concentrated.

Roy’s test is like using a spotlight that only illuminates the brightest part of the pattern, potentially missing subtler but still important differences in other dimensions.

All test statistics are highly significant (p < <0.001), providing strong convergent evidence that iris species differ significantly in their overall morphological profile across the four measurements.

Effect Size Interpretation:

• Pillai’s trace of 1.192 is very large (maximum possible ≈ min(groups-1, variables) = 2)
• Wilks’ Lambda close to 0 indicates strong group separation (ranges from 0 to 1, with smaller values indicating larger effects)
• The consistency across all four statistics provides robust evidence for significant multivariate group differences.

4.3 Univariate Follow-up Tests

Since the MANOVA is significant, we examine univariate ANOVAs for each dependent variable to understand which specific variables contribute to the group differences.

# Display univariate ANOVA results (values already extracted above)
print(univariate_summary)
 Response Sepal.Length :
             Df Sum Sq Mean Sq F value    Pr(>F)    
Species       2 63.212  31.606  119.26 < 2.2e-16 ***
Residuals   147 38.956   0.265                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response Sepal.Width :
             Df Sum Sq Mean Sq F value    Pr(>F)    
Species       2 11.345  5.6725   49.16 < 2.2e-16 ***
Residuals   147 16.962  0.1154                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response Petal.Length :
             Df Sum Sq Mean Sq F value    Pr(>F)    
Species       2 437.10 218.551  1180.2 < 2.2e-16 ***
Residuals   147  27.22   0.185                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response Petal.Width :
             Df Sum Sq Mean Sq F value    Pr(>F)    
Species       2 80.413  40.207  960.01 < 2.2e-16 ***
Residuals   147  6.157   0.042                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.1 Interpretation of Univariate Results

All four morphological variables show statistically significant species differences:

• Sepal.Length: F(2, 147) = 119.3, p < 0.001
• Sepal.Width: F(2, 147) = 49.2, p < 0.001
  
• Petal.Length: F(2, 147) = 1180.2, p < 0.001
• Petal.Width: F(2, 147) = 960, p < 0.001

Petal measurements show the largest F-statistics, indicating that petal characteristics are the primary drivers of species differentiation, with petal length showing the strongest discrimination (F = 1180.2).

5 Assumption Testing

MANOVA has several important assumptions that must be verified for valid interpretation.

5.1 1. Sample Size Adequacy

# Display sample size information (values already computed above)
cat("Total sample size:", n_total, "\n")
Total sample size: 150 

cat("Number of groups:", n_species, "\n")
Number of groups: 3 

cat("Number of dependent variables:", n_vars_check, "\n")
Number of dependent variables: 4 

cat("Group sizes:\n")
Group sizes:

print(table(iris$Species))

    setosa versicolor  virginica 
        50         50         50 

cat("\nMinimum group size:", min_group_size, "\n")

Minimum group size: 50 

cat("Recommended minimum (groups × variables + 10):", recommended_min, "\n")
Recommended minimum (groups × variables + 10): 22 

cat("Sample size adequate:", sample_adequate, "\n")
Sample size adequate: TRUE 

Result: Sample size is adequate. Each group (n=50) exceeds the number of variables (4), and total sample size (150) exceeds the recommended minimum (22).

5.2 2. Multivariate Normality

# Display multivariate normality test (values already computed above)
cat("Multivariate Shapiro-Wilk Test:\n")
Multivariate Shapiro-Wilk Test:

cat("Statistic (W):", round(mshapiro_test$statistic, 4), "\n")
Statistic (W): 0.9793 

cat("p-value:", format(mshapiro_test$p.value, scientific = TRUE), "\n")
p-value: 2.342118e-02 

cat("Assumption met (p > 0.05):", mv_normal, "\n")
Assumption met (p > 0.05): FALSE 

Result: Multivariate normality is violated (W = 0.9793, p = 2.342118e-02). However, MANOVA is generally robust to moderate violations of normality, especially with balanced designs and adequate sample sizes like we have here.

5.3 3. Homogeneity of Covariance Matrices (Box’s M Test)

# Display Box's M test (values already computed above)
print(box_m_result)

    Box's M-test for Homogeneity of Covariance Matrices

data:  Y
Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16

cat("\nInterpretation:\n")

Interpretation:

cat("Chi-square statistic:", round(box_m_stat, 2), "\n")
Chi-square statistic: 140.94 

cat("p-value:", format(box_m_p, scientific = TRUE), "\n")
p-value: 3.352034e-20 

cat("Assumption met (p > 0.05):", covariance_homogeneous, "\n")
Assumption met (p > 0.05): FALSE 

Result: The assumption of equal covariance matrices is violated (χ² = 140.94, p = 3.352034e-20). This suggests that the three species have different patterns of relationships among the morphological variables. However, this violation can be addressed by using Pillai’s trace, which is robust to covariance matrix heterogeneity.

5.4 4. Linearity Assessment

# Display correlation matrix (values already computed above)
print(round(cor_matrix, 3))
             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length        1.000      -0.118        0.872       0.818
Sepal.Width        -0.118       1.000       -0.428      -0.366
Petal.Length        0.872      -0.428        1.000       0.963
Petal.Width         0.818      -0.366        0.963       1.000

cat("\nCorrelation range:", round(min_correlation, 3), "to", round(max_correlation, 3), "\n")

Correlation range: 0.118 to 0.963 

cat("Adequate linearity (0.1 < |r| < 0.9):", linearity_adequate, "\n")
Adequate linearity (0.1 < |r| < 0.9): FALSE 

Result: Linear relationships appear problematic. Correlations range from 0.118 to 0.963, suggesting possible multicollinearity issues.

5.5 5. Multicollinearity Assessment

# Display multicollinearity assessment (values already computed above)
cat("Correlation matrix determinant:", format(det_cor, scientific = TRUE), "\n")
Correlation matrix determinant: 8.109611e-03 

cat("No extreme multicollinearity (det > 0.00001):", multicollinearity_ok, "\n")
No extreme multicollinearity (det > 0.00001): TRUE 

Result: No extreme multicollinearity detected. The correlation matrix determinant (8.109611e-03) is well above the critical threshold, indicating that variables provide relatively independent information.

5.6 6. Multivariate Outliers

# Display outlier analysis (values already computed above)
cat("Critical value (99th percentile):", round(critical_value, 2), "\n")
Critical value (99th percentile): 13.28 

cat("Number of outliers:", outliers, "\n")
Number of outliers: 0 

cat("Proportion of outliers:", round(outlier_proportion * 100, 1), "%\n")
Proportion of outliers: 0 %

cat("Acceptable level (< 5%):", outliers_acceptable, "\n")
Acceptable level (< 5%): TRUE 

# Identify outlier cases if any
if(outliers > 0) {
  outlier_cases <- which(mahal_dist > critical_value)
  cat("Outlier cases:", outlier_cases, "\n")
}

Result: Very few outliers detected (0% of cases), which is within acceptable limits. These cases represent legitimate extreme observations rather than problematic data points.

6 Effect Size and Practical Significance

# Display effect sizes (values already computed above)
cat("Eta-squared (proportion of variance explained by species):\n")
Eta-squared (proportion of variance explained by species):

for(i in 1:4) {
  cat(names(eta_squared)[i], ":", round(eta_squared[i], 3), 
      "(", round(eta_squared[i] * 100, 1), "%)\n")
}
Sepal.Length : 0.619 ( 61.9 %)
Sepal.Width : 0.401 ( 40.1 %)
Petal.Length : 0.941 ( 94.1 %)
Petal.Width : 0.929 ( 92.9 %)

cat("\nEffect size interpretation:\n")

Effect size interpretation:

cat("Small effect: η² ≈ 0.01, Medium effect: η² ≈ 0.06, Large effect: η² ≈ 0.14\n")
Small effect: η² ≈ 0.01, Medium effect: η² ≈ 0.06, Large effect: η² ≈ 0.14

Effect Sizes: All variables show large effect sizes, with species explaining: - Petal.Length: 94.1% of variance (extremely large effect) - Petal.Width: 92.9% of variance (extremely large effect)
- Sepal.Length: 61.9% of variance (large effect) - Sepal.Width: 40.1% of variance (large effect)

These results indicate not only statistical significance but also strong practical significance, with Petal.Length showing the strongest effect (η² = 0.941).

7 Profile Plot Analysis

# Create means for profile plot
species_means <- iris %>%
  group_by(Species) %>%
  summarise(
    Sepal.Length = mean(Sepal.Length),
    Sepal.Width = mean(Sepal.Width),
    Petal.Length = mean(Petal.Length),
    Petal.Width = mean(Petal.Width),
    .groups = 'drop'
  ) %>%
  tidyr::pivot_longer(cols = -Species, names_to = "Variable", values_to = "Mean")

# Store species names for inline reference
species_names <- levels(iris$Species)

# Create profile plot
ggplot(species_means, aes(x = Variable, y = Mean, color = Species, group = Species)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  theme_minimal() +
  labs(
    title = "Profile Plot: Mean Morphological Characteristics by Species",
    subtitle = "Iris Dataset MANOVA Results",
    x = "Morphological Variables",
    y = "Mean Measurement (cm)",
    color = "Species"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "bottom"
  ) +
  scale_color_viridis_d(option = "plasma", end = 0.8)

Profile Plot Interpretation: The plot reveals distinct morphological profiles for each species:

• I. setosa: Characterized by wide sepals but very small petals
• I. versicolor: Intermediate measurements across all variables
• I. virginica: Large flowers overall, especially long petals and sepals

The non-parallel lines confirm significant interaction effects and support the MANOVA findings of species differences.

8 Post-hoc Pairwise Comparisons

# Display pairwise comparisons
cat("Pairwise t-tests with Bonferroni correction:\n\n")
Pairwise t-tests with Bonferroni correction:

variables <- colnames(Y)

for(var in variables) {
  cat("Variable:", var, "\n")
  pairwise_result <- pairwise.t.test(iris[[var]], iris$Species, 
                                     p.adjust.method = "bonferroni")
  print(pairwise_result)
  cat("\n")
}
Variable: Sepal.Length 

    Pairwise comparisons using t tests with pooled SD 

data:  iris[[var]] and iris$Species 

           setosa  versicolor
versicolor 2.6e-15 -         
virginica  < 2e-16 8.3e-09   

P value adjustment method: bonferroni 

Variable: Sepal.Width 

    Pairwise comparisons using t tests with pooled SD 

data:  iris[[var]] and iris$Species 

           setosa  versicolor
versicolor < 2e-16 -         
virginica  1.4e-09 0.0094    

P value adjustment method: bonferroni 

Variable: Petal.Length 

    Pairwise comparisons using t tests with pooled SD 

data:  iris[[var]] and iris$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: bonferroni 

Variable: Petal.Width 

    Pairwise comparisons using t tests with pooled SD 

data:  iris[[var]] and iris$Species 

           setosa versicolor
versicolor <2e-16 -         
virginica  <2e-16 <2e-16    

P value adjustment method: bonferroni 

# extract and store pairwise p-values in the extract_all_results chunk
pairwise_results <- list()
for(var in colnames(Y)) {
  pairwise_results[[var]] <- pairwise.t.test(iris[[var]], iris$Species, 
                                             p.adjust.method = "bonferroni")$p.value
}

# Then make dynamic in interpretation
all_significant <- all(sapply(pairwise_results, function(x) all(x < 0.05, na.rm = TRUE)))

Post-hoc Results: Pairwise t-tests with Bonferroni correction show that all comparisons are highly significant (p < 0.001), indicating that each species differs significantly from every other species on the morphological variables.

9 Conclusions

9.1 Statistical Conclusions

1. Primary Finding: MANOVA revealed highly significant differences between iris species in their overall morphological profile (Pillai’s trace = 1.192, p = 2.22e-16).

2. Univariate Effects: All four morphological variables contribute significantly to species differentiation, with petal length showing the strongest discrimination power (F = 1180.2).

3. Effect Magnitude: Species differences are not only statistically significant but also practically meaningful, with species explaining 40-94% of variance in morphological traits (large to very large effects).

4. Assumption Evaluation: While multivariate normality was violated and covariance homogeneity was violated, the analysis remains valid due to MANOVA’s robustness properties, balanced design (n = 50 per group), and adequate sample sizes. The use of Pillai’s trace provides additional protection against assumption violations.

9.2 Biological Interpretation

Petal characteristics emerge as the primary drivers of species classification (explaining 94.1% of variance in Petal.Length)

The morphological profiles suggest: - Evolutionary divergence in floral architecture among species - Adaptive specialization potentially related to different pollination strategies - Clear taxonomic boundaries supporting the species-level classification

9.3 Methodological Notes

This analysis demonstrates proper MANOVA methodology: - Comprehensive assumption testing with appropriate interpretation - Use of robust test statistics (Pillai’s trace) when assumptions are violated - Integration of multivariate and univariate perspectives - Assessment of both statistical and practical significance - Clear documentation of analytical decisions and limitations

The significant results, combined with large effect sizes (η² range: 0.401 - 0.941) and biological plausibility, provide compelling evidence for genuine morphological differentiation among iris species rather than statistical artifacts.