I am trying to use R to replicate the more detailed output from a Linear Discriminant Analysis that is produced by SPSS.
The R output lacks several of the statistics which are given with SPSS; however, it should be possible to calculate these from the available information. I have been using the Iris data set (https://en.wikipedia.org/wiki/Iris_flower_data_set). Having read through previous answers on this issue, I can see that @ttnphns here gave a detailed comparison of the SPSS and R output, as well as instructions on how to calculate the various statistics here. This is also complemented by the question and answer by @Keaton Wilson here.
However, I am still having difficulty replicating the structure matrix produced by SPSS in R.
My question has two parts which I will summarise here before explaining the details:
Firstly, I can produce a structure matrix using R; however, it does not match the one given by SPSS. I am interested in what the matrix that R produces is and whether it is a useful measure to describe the results of the Linear Discriminant Analysis.
Secondly, I have tried calculating the structure matrix more directly, but end up with a matrix that doesn't match either the R output or the SPSS output, so I suspect I have made a mistake somewhere.
Here are the iris data:
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
7 4.6 3.4 1.4 0.3 setosa
8 5.0 3.4 1.5 0.2 setosa
9 4.4 2.9 1.4 0.2 setosa
10 4.9 3.1 1.5 0.1 setosa
11 5.4 3.7 1.5 0.2 setosa
12 4.8 3.4 1.6 0.2 setosa
13 4.8 3.0 1.4 0.1 setosa
14 4.3 3.0 1.1 0.1 setosa
15 5.8 4.0 1.2 0.2 setosa
16 5.7 4.4 1.5 0.4 setosa
17 5.4 3.9 1.3 0.4 setosa
18 5.1 3.5 1.4 0.3 setosa
19 5.7 3.8 1.7 0.3 setosa
20 5.1 3.8 1.5 0.3 setosa
21 5.4 3.4 1.7 0.2 setosa
22 5.1 3.7 1.5 0.4 setosa
23 4.6 3.6 1.0 0.2 setosa
24 5.1 3.3 1.7 0.5 setosa
25 4.8 3.4 1.9 0.2 setosa
26 5.0 3.0 1.6 0.2 setosa
27 5.0 3.4 1.6 0.4 setosa
28 5.2 3.5 1.5 0.2 setosa
29 5.2 3.4 1.4 0.2 setosa
30 4.7 3.2 1.6 0.2 setosa
31 4.8 3.1 1.6 0.2 setosa
32 5.4 3.4 1.5 0.4 setosa
33 5.2 4.1 1.5 0.1 setosa
34 5.5 4.2 1.4 0.2 setosa
35 4.9 3.1 1.5 0.2 setosa
36 5.0 3.2 1.2 0.2 setosa
37 5.5 3.5 1.3 0.2 setosa
38 4.9 3.6 1.4 0.1 setosa
39 4.4 3.0 1.3 0.2 setosa
40 5.1 3.4 1.5 0.2 setosa
41 5.0 3.5 1.3 0.3 setosa
42 4.5 2.3 1.3 0.3 setosa
43 4.4 3.2 1.3 0.2 setosa
44 5.0 3.5 1.6 0.6 setosa
45 5.1 3.8 1.9 0.4 setosa
46 4.8 3.0 1.4 0.3 setosa
47 5.1 3.8 1.6 0.2 setosa
48 4.6 3.2 1.4 0.2 setosa
49 5.3 3.7 1.5 0.2 setosa
50 5.0 3.3 1.4 0.2 setosa
51 7.0 3.2 4.7 1.4 versicolor
52 6.4 3.2 4.5 1.5 versicolor
53 6.9 3.1 4.9 1.5 versicolor
54 5.5 2.3 4.0 1.3 versicolor
55 6.5 2.8 4.6 1.5 versicolor
56 5.7 2.8 4.5 1.3 versicolor
57 6.3 3.3 4.7 1.6 versicolor
58 4.9 2.4 3.3 1.0 versicolor
59 6.6 2.9 4.6 1.3 versicolor
60 5.2 2.7 3.9 1.4 versicolor
61 5.0 2.0 3.5 1.0 versicolor
62 5.9 3.0 4.2 1.5 versicolor
63 6.0 2.2 4.0 1.0 versicolor
64 6.1 2.9 4.7 1.4 versicolor
65 5.6 2.9 3.6 1.3 versicolor
66 6.7 3.1 4.4 1.4 versicolor
67 5.6 3.0 4.5 1.5 versicolor
68 5.8 2.7 4.1 1.0 versicolor
69 6.2 2.2 4.5 1.5 versicolor
70 5.6 2.5 3.9 1.1 versicolor
71 5.9 3.2 4.8 1.8 versicolor
72 6.1 2.8 4.0 1.3 versicolor
73 6.3 2.5 4.9 1.5 versicolor
74 6.1 2.8 4.7 1.2 versicolor
75 6.4 2.9 4.3 1.3 versicolor
76 6.6 3.0 4.4 1.4 versicolor
77 6.8 2.8 4.8 1.4 versicolor
78 6.7 3.0 5.0 1.7 versicolor
79 6.0 2.9 4.5 1.5 versicolor
80 5.7 2.6 3.5 1.0 versicolor
81 5.5 2.4 3.8 1.1 versicolor
82 5.5 2.4 3.7 1.0 versicolor
83 5.8 2.7 3.9 1.2 versicolor
84 6.0 2.7 5.1 1.6 versicolor
85 5.4 3.0 4.5 1.5 versicolor
86 6.0 3.4 4.5 1.6 versicolor
87 6.7 3.1 4.7 1.5 versicolor
88 6.3 2.3 4.4 1.3 versicolor
89 5.6 3.0 4.1 1.3 versicolor
90 5.5 2.5 4.0 1.3 versicolor
91 5.5 2.6 4.4 1.2 versicolor
92 6.1 3.0 4.6 1.4 versicolor
93 5.8 2.6 4.0 1.2 versicolor
94 5.0 2.3 3.3 1.0 versicolor
95 5.6 2.7 4.2 1.3 versicolor
96 5.7 3.0 4.2 1.2 versicolor
97 5.7 2.9 4.2 1.3 versicolor
98 6.2 2.9 4.3 1.3 versicolor
99 5.1 2.5 3.0 1.1 versicolor
100 5.7 2.8 4.1 1.3 versicolor
101 6.3 3.3 6.0 2.5 virginica
102 5.8 2.7 5.1 1.9 virginica
103 7.1 3.0 5.9 2.1 virginica
104 6.3 2.9 5.6 1.8 virginica
105 6.5 3.0 5.8 2.2 virginica
106 7.6 3.0 6.6 2.1 virginica
107 4.9 2.5 4.5 1.7 virginica
108 7.3 2.9 6.3 1.8 virginica
109 6.7 2.5 5.8 1.8 virginica
110 7.2 3.6 6.1 2.5 virginica
111 6.5 3.2 5.1 2.0 virginica
112 6.4 2.7 5.3 1.9 virginica
113 6.8 3.0 5.5 2.1 virginica
114 5.7 2.5 5.0 2.0 virginica
115 5.8 2.8 5.1 2.4 virginica
116 6.4 3.2 5.3 2.3 virginica
117 6.5 3.0 5.5 1.8 virginica
118 7.7 3.8 6.7 2.2 virginica
119 7.7 2.6 6.9 2.3 virginica
120 6.0 2.2 5.0 1.5 virginica
121 6.9 3.2 5.7 2.3 virginica
122 5.6 2.8 4.9 2.0 virginica
123 7.7 2.8 6.7 2.0 virginica
124 6.3 2.7 4.9 1.8 virginica
125 6.7 3.3 5.7 2.1 virginica
126 7.2 3.2 6.0 1.8 virginica
127 6.2 2.8 4.8 1.8 virginica
128 6.1 3.0 4.9 1.8 virginica
129 6.4 2.8 5.6 2.1 virginica
130 7.2 3.0 5.8 1.6 virginica
131 7.4 2.8 6.1 1.9 virginica
132 7.9 3.8 6.4 2.0 virginica
133 6.4 2.8 5.6 2.2 virginica
134 6.3 2.8 5.1 1.5 virginica
135 6.1 2.6 5.6 1.4 virginica
136 7.7 3.0 6.1 2.3 virginica
137 6.3 3.4 5.6 2.4 virginica
138 6.4 3.1 5.5 1.8 virginica
139 6.0 3.0 4.8 1.8 virginica
140 6.9 3.1 5.4 2.1 virginica
141 6.7 3.1 5.6 2.4 virginica
142 6.9 3.1 5.1 2.3 virginica
143 5.8 2.7 5.1 1.9 virginica
144 6.8 3.2 5.9 2.3 virginica
145 6.7 3.3 5.7 2.5 virginica
146 6.7 3.0 5.2 2.3 virginica
147 6.3 2.5 5.0 1.9 virginica
148 6.5 3.0 5.2 2.0 virginica
149 6.2 3.4 5.4 2.3 virginica
150 5.9 3.0 5.1 1.8 virginica
In R, the lda can be performed using:
library(MASS)
iris_lda <- lda(Species ~ ., data = iris)
The unstandardised discriminant coefficents and discriminant scores match those in the SPSS output and can be obtained using:
#Unstandardised discriminant coefficients
iris_lda$scaling
LD1 LD2
Sepal.Length 0.8293776 0.02410215
Sepal.Width 1.5344731 2.16452123
Petal.Length -2.2012117 -0.93192121
Petal.Width -2.8104603 2.83918785
#Discriminant scores
predict(iris_lda)$x
LD1 LD2
1 8.0617998 0.300420621
2 7.1286877 -0.786660426
3 7.4898280 -0.265384488
4 6.8132006 -0.670631068
5 8.1323093 0.514462530
6 7.7019467 1.461720967
7 7.2126176 0.355836209
8 7.6052935 -0.011633838
9 6.5605516 -1.015163624
10 7.3430599 -0.947319209
...etc
Additional outputs can be obtained using the package candisc mentioned in this helpful post by @Keaton Wilson here.
library(candisc)
#Run the lda
man1 <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data = iris)
can_lda <- candisc(man1)
#E.g. Standardised discriminant coefficients:
can_lda$coeffs.std
Can1 Can2
Sepal.Length -0.4269548 0.01240753
Sepal.Width -0.5212417 0.73526131
Petal.Length 0.9472572 -0.40103782
Petal.Width 0.5751608 0.58103986
Part 1
The structure matrix from candisc (which I believe is the same as the pooled within-groups correlations, i.e. as mentioned here) doesn't match the SPSS output:
In R:
can_lda$structure
Can1 Can2
Sepal.Length 0.7918878 0.21759312
Sepal.Width -0.5307590 0.75798931
Petal.Length 0.9849513 0.04603709
Petal.Width 0.9728120 0.22290236
And the SPSS output (copied from @ttnphns answer). A friend was also able to replicate this same output for me in SPSS.
Pooled within-groups correlations between variables and discriminants
Dis1 Dis2
SLength .2225959415 .3108117231
SWidth -.1190115149 .8636809224
PLength .7060653811 .1677013843
PWidth .6331779262 .7372420588
It should be possible to calculate the structure matrix between variables and discriminants by calculating the covariance between standardised discriminant scores and the original variables, so I tried this:
#Store the scores
dfs <- predict(iris_lda)$x
#Then we standardise these
z_dfs <- apply(dfs, 2, FUN = function (x) {(x - mean(x)) / sd(x)})
#Then we calculate the covariance between these and the original variables,
#divided by the standard deviation of the original variables
apply(iris[,-5], 2, FUN = function (x) {cov(x, z_dfs) / sd (x)})
Sepal.Length Sepal.Width Petal.Length Petal.Width
[1,] -0.7918878 0.5307590 -0.98495127 -0.9728120
[2,] 0.2175931 0.7579893 0.04603709 0.2229024
However, this gives an identical structure matrix to that obtained from the candisc package (I note that some of the signs have reversed, but this doesn't seem to be an issue). So although I have apparently calculated something useful here, it still doesn't match the SPSS output.
Does this matrix produced by R have a use in interpreting the discriminant loadings, and how is it related to the SPSS output?
Part 2.
I am also interested in whether I can calculate the structure matrix from the original data.
To do this I am following the detailed guidelines provided by @ttnphns here, which have been very helpful in replicating the analysis. This says we require two bits of information.
The matrix $\mathbf {S_w}$, described as "the pooled within-group scatter matrix (i.e. the sum of $\mathbf k$
p x p
SSCP matrices of the variables, centered about the respective groups' centroid)". k are the number of groups (here species).The discriminant eigenvectors $ \mathbf V$, obtained using $\mathbf {S_w}$, the total scatter matrix $\mathbf {S_t}$ and the between group scatter matrix $\mathbf {S_b} = \mathbf{S_t} - \mathbf{S_w}$. I think neither the lda function in MASS nor candisc seem to output the eigenvectors directly.
Calculating $\mathbf{S_w}$:
#Group centering the dataset by columns
gc_iris_set <- apply(iris[which(iris$Species == "setosa"), 1:4], 2, function(x) {x - mean (x)})
gc_iris_ver <- apply(iris[which(iris$Species == "versicolor"), 1:4],2, function(x) {x - mean (x)})
gc_iris_vir <- apply(iris[which(iris$Species == "virginica"), 1:4], 2, function(x) {x - mean (x)})
#Calculating an SSCP matrix (see: https://stats.stackexchange.com/a/22520) for each group
SSCP_set_gc <- crossprod(gc_iris_set)
SSCP_ver_gc <- crossprod(gc_iris_ver)
SSCP_vir_gc <- crossprod(gc_iris_vir)
#Taking the sum of these to give Sw
Sw <- SSCP_set_gc + SSCP_ver_gc + SSCP_vir_gc
Sw
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 38.9562 13.6300 24.6246 5.6450
Sepal.Width 13.6300 16.9620 8.1208 4.8084
Petal.Length 24.6246 8.1208 27.2226 6.2718
Petal.Width 5.6450 4.8084 6.2718 6.1566
Calculating the discriminant eigenvectors $\mathbf V$:
#Centering the iris data to calculate the total scatter matrix
c_iris <- apply(iris[,1:4], 2, FUN = function(x) {(x - mean(x))})
#Calculating the total scatter matrix
St <- crossprod(c_iris)
#And the between group scatter matrix
Sb <- St - Sw
#The cholesky root of Sw
U <- chol(Sw)
#Calculation of the eigenvectors of the LDA
LDA_V <- solve(U) %*% eigen(t(solve(U)) %*% Sb %*% solve(U))$vectors
#The eigenvectors
LDA_V
[,1] [,2] [,3] [,4]
Sepal.Length -0.06840592 -0.001987912 0.1824441 0.18919900
Sepal.Width -0.12656121 -0.178526702 -0.2192389 -0.02956174
Petal.Length 0.18155288 0.076863566 -0.2478258 -0.01788111
Petal.Width 0.23180286 -0.234172267 0.3513745 -0.13460680
The structure matrix should be calculated using $\mathbf R = diag(\mathbf {S_w})^{-1}\mathbf {S_w}\mathbf V$, so:
solve(diag(diag(Sw))) %*% Sw %*% LDA_V[,c(1,2)]
[,1] [,2]
[1,] 0.03566391 -0.04979768
[2,] -0.02889685 -0.20970790
[3,] 0.13532565 -0.03214192
[4,] 0.25518509 -0.29712530
#I was initially unsure whether to take the inverse of Sw before creating the
#diagonal matrix or do this the other way round; however this was confirmed in a
#comment by @ttnphns below.
#Neither approach gives results which match either the R or SPSS output
This matches neither of the outputs produced above. I would welcome any assistance in identifying what is wrong with my calculation here.
I assume that it is the final calculation $\mathbf R = diag(\mathbf {S_w})^{-1}\mathbf {S_w}\mathbf V$ that is the issue, because I am fairly sure I have correct values for $\mathbf V$ and $\mathbf{S_w}$. I can use these values to correctly produce other statistics from the LDA - for instance the standardized discriminant coefficients:
sqrt(diag(Sw)) * LDA_V[,1:2]
[,1] [,2]
[1,] -0.4269548 -0.01240753
[2,] -0.5212417 -0.73526131
[3,] 0.9472572 0.40103782
[4,] 0.5751608 -0.58103986
#which match
can_lda$coeffs.std
Can1 Can2
Sepal.Length -0.4269548 0.01240753
Sepal.Width -0.5212417 0.73526131
Petal.Length 0.9472572 -0.40103782
Petal.Width 0.5751608 0.58103986
NB. Following the comments below by @ttnphns, there was a missing square root from the final equation. This has now been corrected in @ttnphns answer here and I have added an answer below detailing this final step in R.