Iris data

Topic of Covariance and the multivariate normal, part of Statistical learning · Also in Estimation and inference for the multivariate normal, Principal component analysis, K-means clustering, Mixtures and the EM algorithm, Statistical classifiers and discriminant analysis, Logistic regression, Networks in practice

Iris data

Type Topic
Taught in Covariance and the multivariate normal · Estimation and inference for the multivariate normal · Principal component analysis · K-means clustering · Mixtures and the EM algorithm · Statistical classifiers and discriminant analysis · Logistic regression · Networks in practice
Part Statistical learning
Questions using it 2

Measurements of sepal length, sepal width, petal length and petal width, in centimetres, on 150 iris flowers, 50 from each of three species: setosa, versicolor and virginica. Edgar Anderson collected them; Fisher’s 1936 paper on discriminant analysis made them famous. The statistical-learning chapters use it throughout, from the multivariate normal and Hotelling’s test to clustering, discriminant analysis, logistic regression and a small network.

In Covariance and the multivariate normal

From Covariance and the multivariate normal

The four iris measurements share centimetre units. Pooling the species changes the question. The calculation below compares the covariance of all 150 flowers with the covariance within setosa, and then checks that the sample conditional-mean slope matches the least-squares slope. A single Gaussian approximation to the pooled data ignores the species structure. In the code, iris_x holds the four measurements and iris_y holds the species labels.

iris_x <- iris[, 1:4]
iris_y <- iris$Species
round(cov(iris_x), 3)
round(cov(subset(iris, Species == "setosa")[, 1:4]), 3)
aggregate(. ~ Species, data = iris, FUN = mean)

setosa <- subset(iris, Species == "setosa")
x <- setosa$Sepal.Length
y <- setosa$Sepal.Width
ols <- coef(lm(y ~ x))
S <- cov(cbind(x, y))
slope <- S[2, 1] / S[1, 1]
gauss <- c(intercept = mean(y) - slope * mean(x), slope = slope)
rbind(ols = unname(ols), gaussian = unname(gauss))
stopifnot(max(abs(unname(ols) - unname(gauss))) < 1e-10)

pairs(iris[, 1:4], col = as.integer(iris_y), pch = 19, cex = 0.55)

The pooled sample covariance, rounded to three decimals, is

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length        0.686      -0.042        1.274       0.516
Sepal.Width        -0.042       0.190       -0.330      -0.122
Petal.Length        1.274      -0.330        3.116       1.296
Petal.Width         0.516      -0.122        1.296       0.581

Within setosa the sample covariance is much smaller, and the sepal and petal blocks are only weakly associated:

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length        0.124       0.099        0.016       0.010
Sepal.Width         0.099       0.144        0.012       0.009
Petal.Length        0.016       0.012        0.030       0.006
Petal.Width         0.010       0.009        0.006       0.011

The species means, in centimetres, are

     Species Sepal.Length Sepal.Width Petal.Length Petal.Width
      setosa        5.006       3.428        1.462       0.246
  versicolor        5.936       2.770        4.260       1.326
   virginica        6.588       2.974        5.552       2.026

On setosa, sepal width against sepal length, the least-squares fit and the slope taken from the sample covariance agree:

             intercept   slope
ols           -0.5694  0.7985
gaussian      -0.5694  0.7985

The scatterplot matrix of the four measurements, coloured by species, shows three clouds. The colours are the observed labels. No model has been fitted to produce them.

The pooled covariance includes both within-species variation and variation between species means. The law of total covariance makes the split explicit. If G denotes the species,

\operatorname{Cov}(\mathbf{X})=\mathbb{E}\bigl\{\operatorname{Cov}(\mathbf{X}\mid G)\bigr\}+\operatorname{Cov}\bigl\{\mathbb{E}(\mathbf{X}\mid G)\bigr\}.

A strong pooled correlation can therefore be a difference between group means. The setosa block shows what remains when that difference is removed.

Key idea

Conditioning a multivariate normal vector on some of its coordinates gives another multivariate normal. The conditional mean is affine in the conditioning values, and the conditional covariance is the Schur complement \Sigma_{AA}-\Sigma_{AB}\Sigma_{BB}^{-1}\Sigma_{AB}^\top, which does not depend on the observed values. Zero covariance is independence inside this family. The iris calculation shows why the population has to be named first: the pooled covariance mixes within-species scatter with the gaps between species means.

In Estimation and inference for the multivariate normal

From Estimation and inference for the multivariate normal

The one-sample test is applied here to the two sepal measurements of Iris setosa, with illustrative null mean (5.0, 3.4) cm. That vector is a chosen comparison point. It is not a biological standard. The probability from the F calibration is exact under independent bivariate normal sampling, with n>p and nonsingular S.

x <- as.matrix(subset(iris, Species == "setosa")[, 1:2])
n <- nrow(x); p <- ncol(x)
d <- colMeans(x) - c(5.0, 3.4)
T2 <- n * drop(t(d) %*% solve(cov(x), d))
Fstat <- (n - p) * T2 / (p * (n - 1))
c(T_squared = T2, F = Fstat,
  df1 = p, df2 = n - p,
  p_value = pf(Fstat, p, n - p, lower.tail = FALSE))
T_squared         F       df1       df2   p_value
    0.4323    0.2117    2.0000   48.0000    0.8099

The computed p-value is 0.81. At the 5\% level this joint null is not rejected. A failure to reject does not establish that the true mean equals (5.0, 3.4). The reference distribution F_{p,n-p} is the exact null law only under the independent normal sampling used to derive the F calibration. The p-value does not check that sampling model.

Key idea

Hotelling’s T^2=n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top S^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0) is n times the squared Mahalanobis distance of the sample mean from the hypothesised mean, in the metric of the sample covariance. It is a strictly decreasing function of the likelihood ratio. Under independent normal sampling with n>p,

\frac{n-p}{p(n-1)}T^2\sim F_{p,\,n-p}

when the null mean is true. When p=1 the statistic is the square of the ordinary t statistic. The two-sample version uses the same geometry with the pooled covariance and the coefficient n_1 n_2/(n_1+n_2).

In Principal component analysis

From Principal component analysis

The fit uses every iris flower, and it is descriptive. Species labels colour the plot after the fit. They do not enter S. Correlation principal components standardise each measurement to unit variance before the eigendecomposition. Covariance principal components keep the original variances, so a measurement with larger sample variance takes a larger share of the first direction. That is the unit dependence in the scaling remark: petal length, whose sample variance in the pooled covariance is about 3.12, dominates a covariance analysis, while standardisation gives each measurement the same marginal scale.

iris_x <- iris[, 1:4]
iris_y <- iris$Species
pc_cov <- prcomp(iris_x, center = TRUE, scale. = FALSE)
pc_cor <- prcomp(iris_x, center = TRUE, scale. = TRUE)
round(rbind(Covariance = pc_cov$sdev^2 / sum(pc_cov$sdev^2),
            Correlation = pc_cor$sdev^2 / sum(pc_cor$sdev^2)), 3)
round(pc_cor$rotation[, 1:2], 3)
stopifnot(max(abs(cov(pc_cor$x) - diag(pc_cor$sdev^2))) < 1e-10)
plot(pc_cor$x[, 1:2], col = as.integer(iris_y), pch = 19,
     xlab = "PC1 score", ylab = "PC2 score")

The proportions of variance are

              PC1   PC2   PC3   PC4
Covariance  0.925 0.053 0.017 0.005
Correlation 0.730 0.229 0.037 0.005

The first two standardised components explain 95.8\% of the total standardised variance: the partial sum of the eigenvalues, divided by the trace, before the row is rounded for display. The leading two correlation loadings, rounded to three decimals, are

                PC1     PC2
Sepal.Length  0.521  -0.377
Sepal.Width  -0.269  -0.923
Petal.Length  0.580  -0.024
Petal.Width   0.565  -0.067

The sign of a loading may be flipped without changing the line it spans or the variance it carries. The check cov(pc_cor$x) against diag(pc_cor$sdev^2) confirms that the sample covariance of the scores is diagonal, with the reported component variances on the diagonal. The score plot colours each flower by species after the fit. Separation of species in that plot describes this sample. A prediction experiment has to estimate the centring, the scaling, and the directions inside each training fold, and apply that transformation unchanged to the validation rows.

Key idea

The principal component directions are the eigenvectors of the sample covariance, ordered by eigenvalue. The jth eigenvalue is the sample variance of the jth component, so the proportion of variance explained by the first k components is \sum_{j\le k}\lambda_j\big/\sum_{j\le p}\lambda_j. The same directions are the right singular vectors of the centred data matrix. Scaling changes the matrix being diagonalised, and a principal axis is not a regression line: the axis maximises explanatory variance, while a regression line minimises squared error in a response.

In K-means clustering

From K-means clustering

The iris species labels are hidden while fitting. Standardisation chooses a distance measured in standard deviations. Setting K=3 uses the known number of species. The question is whether three geometric groups resemble those species, rather than whether the algorithm discovers the number three. The cross-tabulation answers that question without pretending that cluster numbers are species names. The default algorithm is used, with 50 random starts, and the objective is checked by recomputing the sum of squared distances to the reported centres.

set.seed(4006)
z <- scale(iris[, 1:4])
km <- kmeans(z, centers = 3, nstart = 50)
table(Cluster = km$cluster, Species = iris$Species)
c(within_SSE = km$tot.withinss,
  adjusted_Rand_index = mclust::adjustedRandIndex(km$cluster, iris$Species))
direct_SSE <- sum((z - km$centers[km$cluster, ])^2)
stopifnot(abs(direct_SSE - km$tot.withinss) < 1e-8)

The cross-tabulation is

Cluster   setosa   versicolor   virginica
      1        0           11          36
      2       50            0           0
      3        0           39          14

and the two summary numbers are

within_SSE   adjusted_Rand_index
  138.8884                0.6201

Setosa is recovered as its own cluster: all 50 setosa flowers sit in cluster 2, and nothing else does. Versicolor and virginica are mixed, 39 against 14 in one cluster and 11 against 36 in the other. The adjusted Rand index compares the clustering with the species partition without requiring cluster 1 to mean species 1. On its chance-corrected scale, identical partitions score 1, and the index can be negative. The value 0.6201 is the agreement on this run. The direct sum of squared distances to the assigned centres reproduces tot.withinss to the tolerance in the check, which is the unweighted criterion \widetilde W.

K-means favours compact groups in the standardised Euclidean metric. A species can have a shape, or an overlap, that this geometry does not represent. The Voronoi picture for three fixed centres in the plane is a partition into three convex cells, with boundaries that are pieces of the perpendicular bisectors. Those cells are what the assignment step builds, in the standardised four-dimensional measurement space as well as in a two-dimensional illustration.

Key idea

K-means turns a geometric preference, that nearby points share a centre, into the unweighted within-cluster sum of squares \widetilde W. Lloyd’s two steps, nearest-centre assignment and mean updates, do not increase that objective, and with the tie convention they stop at a fixed point of the alternation. The fixed point need not be the global minimum. The metric, the number of centres and the starting values each affect the partition. Check all three before interpreting a cluster label as a population.

In Mixtures and the EM algorithm

From Mixtures and the EM algorithm

The fit below is a three-component Gaussian mixture on the two petal measurements, with H held fixed. Fixing H separates the membership probabilities from the choice of how many components to use. Let iris_x be the 150\times 4 matrix of iris measurements in the usual order (sepal length, sepal width, petal length, petal width), and let iris_y be the species labels. Columns 3 and 4 are the petal measurements. The name VVV asks each component for its own unconstrained covariance, which is the normal mixture derived above.

set.seed(5006)
mix <- mclust::Mclust(iris_x[, 3:4], G = 3,
                      modelNames = "VVV", verbose = FALSE)
table(Component = mix$classification, Species = iris_y)
#>           Species
#> Component setosa versicolor virginica
#>         1      50          0         0
#>         2       0         49         3
#>         3       0          1        47
round(head(mix$z), 3)
#>      [,1] [,2] [,3]
#> [1,]    1     0    0
#> [2,]    1     0    0
#> [3,]    1     0    0
#> [4,]    1     0    0
#> [5,]    1     0    0
#> [6,]    1     0    0
summary(1 - apply(mix$z, 1, max))
#>    Min. 1st Qu. Median   Mean 3rd Qu.   Max.
#>  0.0000  0.0000 0.0003 0.0310  0.0175 0.4596
stopifnot(max(abs(rowSums(mix$z) - 1)) < 1e-10)

Each row of mix$z is the posterior membership at the fitted parameter. The classification table counts 50 setosa in component 1, 49 versicolor and 3 virginica in component 2, and 1 versicolor and 47 virginica in component 3. One minus the largest membership has minimum 0, median 0.0003, mean 0.0310, and maximum 0.4596. That summary is the uncertainty of the assigned component given the fitted parameter and given H=3. It leaves out uncertainty about \Psi and about H. Permuting the component labels leaves the fitted density unchanged, which is why bootstrap components have to be matched before their parameters are averaged.

Key idea

EM replaces a log-of-a-sum by the expected complete-data log-likelihood. The E step computes the soft memberships \tau_{ih}. The M step is the weighted sample proportion, mean, and covariance. An exact step cannot decrease the observed-data log-likelihood. The same unconstrained likelihood is unbounded: a component covariance can collapse onto sample points, so a rising likelihood is checked against the smallest covariance eigenvalues and is not, by itself, a global maximum.

In Statistical classifiers and discriminant analysis

From Statistical classifiers and discriminant analysis

The split holds out 15 flowers of each species and trains on the remaining 35, so both fits below see the same four predictors and the same 105 labelled flowers. The held-out set has 45 flowers. LDA forces one covariance. QDA estimates one per species. That is the only difference between the two rules being compared. The data frames train and test carry Species and the four iris measurements.

lda_fit <- MASS::lda(Species ~ ., data = train)
qda_fit <- MASS::qda(Species ~ ., data = train)
lda_pred <- predict(lda_fit, test)$class
qda_pred <- predict(qda_fit, test)$class
table(Observed = test$Species, Predicted = lda_pred)
#>             Predicted
#> Observed     setosa versicolor virginica
#>   setosa         15          0         0
#>   versicolor      0         13         2
#>   virginica       0          1        14
c(LDA_test_error = mean(lda_pred != test$Species),
  QDA_test_error = mean(qda_pred != test$Species))
#> LDA_test_error QDA_test_error
#>        0.06667        0.06667

LDA misclassifies two versicolor flowers as virginica and one virginica flower as versicolor: 3 mistakes in 45, and 3/45=0.06667. The QDA test error on the same split is also 0.06667. One test flower changes the error by 1/45, about 2.2 percentage points, so a sample of 45 held-out flowers makes a coarse comparison. Equal class counts describe one assessment population. A different prevalence, or a different cost of mistakes, changes which accuracy matters, as in the unequal-cost rule.

Key idea

A statistical classifier estimates the class probabilities \pi_g and the class densities f_g, converts them to posterior probabilities \tau_g(\mathbf{x}), and plugs those probabilities into the Bayes rule. Equal costs choose the largest posterior. Unequal costs choose the class with the smallest expected loss. The Bayes risk is a floor for a fresh observation from the true distribution. Under a shared normal covariance the log posterior odds are linear in \mathbf{x}, which is LDA. Separate covariances leave a quadratic boundary, which is QDA. Fisher’s direction \mathbf{a}\propto S^{-1}(\widehat{\boldsymbol\mu}_1-\widehat{\boldsymbol\mu}_2) is the same vector, obtained by maximising separation of the class means relative to pooled within-class scatter.

In Logistic regression

From Logistic regression

Setosa is easy to separate on these measurements, so a three-class fit would hide the binary difficulty. Restrict training to versicolor and virginica and use sepal length alone. Set Y=1 for virginica, so a positive coefficient means larger modelled odds of virginica. The probability threshold 1/2 is the equal-cost rule. The Brier score below scores the probabilities themselves.

btrain <- droplevels(subset(train, Species != "setosa"))
btest <- droplevels(subset(test, Species != "setosa"))
btrain$y <- as.integer(btrain$Species == "virginica")
logit_fit <- glm(y ~ Sepal.Length, data = btrain, family = binomial())
round(coef(summary(logit_fit)), 4)
#>              Estimate Std. Error z value Pr(>|z|)
#> (Intercept)   -12.446     3.3824  -3.680    2e-04
#> Sepal.Length    1.995     0.5439   3.669    2e-04
exp(coef(logit_fit)["Sepal.Length"])
#> Sepal.Length
#>        7.355
prob <- predict(logit_fit, btest, type = "response")
truth <- as.integer(btest$Species == "virginica")
c(test_error = mean((prob >= 0.5) != truth),
  Brier_score = mean((prob - truth)^2))
#> test_error Brier_score
#>     0.2000      0.1844

The estimated odds multiplier for one additional centimetre of sepal length is 7.35. That number multiplies the odds. It does not add 7.35, or 1.995, to the probability, because the logistic curve is nonlinear. The Brier score averages the squared difference between the fitted probability and the binary label, so two classifiers can make the same hard decisions at the threshold 1/2 and still receive different Brier scores when their probabilities differ. On these held-out flowers the threshold-1/2 error is 0.2000 and the Brier score is 0.1844.

Key idea

Logistic regression models the log odds as \beta^\top x, fits \beta by the Bernoulli score X^\top(y-p)=0, and classifies by the sign of that linear predictor when the threshold is 1/2. A coefficient multiplies the odds. The probability moves by ap/(1-p+ap), which depends on where the curve already sits.

In Networks in practice

From Networks in practice

Use scaled predictors, one hidden layer, three output classes and weight decay. Architecture and decay are fixed for this illustration. A production comparison would tune these choices inside training data, by the validation discipline above. The objects train, test and iris_scaled are the training flowers, the test flowers, and their scaled predictors.

set.seed(16006)
y_onehot <- nnet::class.ind(train$Species)
net <- nnet::nnet(iris_scaled$train, y_onehot,
                  size = 3, softmax = TRUE, decay = 0.1,
                  maxit = 1000, trace = FALSE)
prob <- predict(net, iris_scaled$test, type = "raw")
pred <- factor(colnames(y_onehot)[max.col(prob)],
               levels = levels(train$Species))
c(convergence_code = net$convergence,
  test_error = mean(pred != test$Species),
  max_probability_sum_error = max(abs(rowSums(prob) - 1)))
#>           convergence_code                test_error
#>                    0.000e+00                 8.889e-02
#> max_probability_sum_error
#>                    2.220e-16
stopifnot(net$convergence == 0,
          max(abs(rowSums(prob) - 1)) < 1e-8)

The convergence code reports the optimiser’s stopping status. It is not a proof that the global optimum was found. Weight decay shrinks weights. The species probabilities must sum to one for each flower, which the explicit check verifies. The reported test error is 0.08889, and the largest departure of a probability row from summing to 1 is 2.220\times 10^{-16}.

Key idea

Initialise with small random weights scaled by the fan-in, normalise the inputs, and use cross-entropy with softmax for classification. Choose when to stop, how many hidden units, and how much weight decay by a validation set. Then report performance on a test set that was not used in any of those choices.

Used in

See also

Multivariate normal distribution · Hotelling’s T² test · Binary logistic regression · Joint distribution and density · Mean vector and covariance matrix · Correlation · Squared Mahalanobis distance · Linear images of a normal vector

Glossary and index