3
$\begingroup$

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.

  1. 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).

  2. 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.

Improve this question
$\endgroup$
9
  • $\begingroup$ I've added the (non-scaled) eigenvectors V to my answer, and they are the same as your two first columns of LDA_V. You can have only 2 but not 4 eigenvectors because there is only 2 discriminant functions. $\endgroup$
    – ttnphns
    Sep 3, 2020 at 23:21
  • 1
    $\begingroup$ The number of discriminants is min(num_of_vars, num_of_classes-1)=2. So you simply disregard 3+th eivenvectors. Then you invert Sw, take its diagonal and multiply by Sw and V, to get the structure matrix. $\endgroup$
    – ttnphns
    Sep 3, 2020 at 23:33
  • $\begingroup$ I've added the Sw and Sb matrices to my old answer so you can compare with your calculations. $\endgroup$
    – ttnphns
    Sep 3, 2020 at 23:40
  • 1
    $\begingroup$ Hey, Pratorum. I see there was a lapse in my formula of the pooled within group correlations here: the square root got somehow lost. I've corrected the formula right now. $\endgroup$
    – ttnphns
    Sep 28, 2020 at 12:27
  • 1
    $\begingroup$ Yes - that has sorted it thanks. I now get the correct value for the structure matrix. I will add a note to Part 2 of my question so that it still makes sense following your correction and I will add an answer with the correct R code for the final step in the calculation of the structure matrix. $\endgroup$
    – Pratorum
    Sep 28, 2020 at 15:03

2 Answers 2

Reset to default
1
$\begingroup$

I learned a ton from this question, thank you very much for posting. Also, I think I may have stumbled upon an answer.

So, the structure matrix is also described as “contain(ing) correlations between predictors and discriminant functions" (Tabachnick & Fidell, 2016, p. 444)…but taken literally, this doesn’t quite workout (as you tested).

It is also described like this: "Mathematically, the matrix of loadings is the pooled within-group correlation matrix multiplied by the matrix of standardized discriminant function coefficients" (p. 444). I eventually found my way to a “pooled within-group correlation matrix”, but my lack of understanding in matrix algebra prevented me from effectively multiplying the two.

However, much latter on, I had a brain blast: What if I just included the discriminant functions in the covariance matrix with the predictors, then converted that to a correlation matrix. lo and behold, the correlations associated with the discriminant functions match those of the example I was using (Tabachnick and Fidell’s first example on the topic).

The code for your example should look something like this (sorry about my R messiness):

> library(MASS)
> 
> iris_lda <- lda(Species ~ ., data = iris) # Create the function
> DAscores <- predict(iris_lda)$x # Get the case sepecific function scores
> 
> iris_DAsocres <- cbind(iris, DAscores) # add the scores to the original dataset
> head(iris_DAsocres)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species      LD1        LD2
1          5.1         3.5          1.4         0.2  setosa 8.061800  0.3004206
2          4.9         3.0          1.4         0.2  setosa 7.128688 -0.7866604
3          4.7         3.2          1.3         0.2  setosa 7.489828 -0.2653845
4          4.6         3.1          1.5         0.2  setosa 6.813201 -0.6706311
5          5.0         3.6          1.4         0.2  setosa 8.132309  0.5144625
6          5.4         3.9          1.7         0.4  setosa 7.701947  1.4617210
> 
> # group specific datasets with just predictors and functions
> table(iris_DAsocres$Species)

    setosa versicolor  virginica 
        50         50         50 
> 
> # setosa
> iris_DAsocres_setosa <- subset(iris_DAsocres, Species=="setosa",
+                              select = c("Sepal.Length", "Sepal.Width", 
+                                         "Petal.Length","Petal.Width",
+                                         "LD1","LD2"))
> cov(iris_DAsocres_setosa) # within group covariance
             Sepal.Length Sepal.Width Petal.Length  Petal.Width         LD1        LD2
Sepal.Length   0.12424898 0.099216327  0.016355102  0.010330612  0.19025929 0.23183940
Sepal.Width    0.09921633 0.143689796  0.011697959  0.009297959  0.25089470 0.32890802
Petal.Length   0.01635510 0.011697959  0.030159184  0.006069388 -0.05192976 0.01484082
Petal.Width    0.01033061 0.009297959  0.006069388  0.011106122 -0.02173788 0.04625080
LD1            0.19025929 0.250894696 -0.051929761 -0.021737877  0.71818979 0.53432908
LD2            0.23183940 0.328908018  0.014840823  0.046250797  0.53432908 0.83500044
> 
> # versicolor
> iris_DAsocres_versicolor <- subset(iris_DAsocres, Species=="versicolor",
+                                select = c("Sepal.Length", "Sepal.Width", 
+                                           "Petal.Length","Petal.Width",
+                                           "LD1","LD2"))
> cov(iris_DAsocres_versicolor) # within group covariance
             Sepal.Length Sepal.Width Petal.Length Petal.Width         LD1        LD2
Sepal.Length   0.26643265  0.08518367   0.18289796  0.05577959 -0.20767811  0.1787257
Sepal.Width    0.08518367  0.09846939   0.08265306  0.04120408 -0.07599126  0.2551522
Petal.Length   0.18289796  0.08265306   0.22081633  0.07310204 -0.41299348  0.1850795
Petal.Width    0.05577959  0.04120408   0.07310204  0.03910612 -0.16133037  0.1334358
LD1           -0.20767811 -0.07599126  -0.41299348 -0.16133037  1.07364854 -0.2426600
LD2            0.17872572  0.25515218   0.18507955  0.13343580 -0.24266002  0.7629597
> 
> # virginica
> iris_DAsocres_virginica <- subset(iris_DAsocres, Species=="virginica",
+                                    select = c("Sepal.Length", "Sepal.Width", 
+                                               "Petal.Length","Petal.Width",
+                                               "LD1","LD2"))
> cov(iris_DAsocres_virginica) # within group covariance
             Sepal.Length Sepal.Width Petal.Length Petal.Width         LD1         LD2
Sepal.Length   0.40434286  0.09376327   0.30328980  0.04909388 -0.32635130  0.06944266
Sepal.Width    0.09376327  0.10400408   0.07137959  0.04762857 -0.05362318  0.29608525
Petal.Length   0.30328980  0.07137959   0.30458776  0.04882449 -0.44660957  0.01658269
Petal.Width    0.04909388  0.04762857   0.04882449  0.07543265 -0.20567139  0.27294322
LD1           -0.32635130 -0.05362318  -0.44660957 -0.20567139  1.20816167 -0.29166906
LD2            0.06944266  0.29608525   0.01658269  0.27294322 -0.29166906  1.40203983
> 
> # create pooled within-group covariance matrix (sorry I suck with functions)
> # different n per group (just in case)
> pooled_cov_matrix_LDs <-
+   (((length(iris_DAsocres_setosa$Sepal.Length)-1)*(cov(iris_DAsocres_setosa))) + 
+      ((length(iris_DAsocres_versicolor$Sepal.Length)-1)*(cov(iris_DAsocres_versicolor)))+
+      ((length(iris_DAsocres_virginica$Sepal.Length)-1)*(cov(iris_DAsocres_virginica))))/((length(iris_DAsocres_setosa$Sepal.Length)-1)+(length(iris_DAsocres_versicolor$Sepal.Length)-1)+(length(iris_DAsocres_virginica$Sepal.Length)-1))
> pooled_cov_matrix_LDs
             Sepal.Length Sepal.Width Petal.Length Petal.Width           LD1           LD2
Sepal.Length   0.26500816  0.09272109   0.16751429  0.03840136 -1.145900e-01  1.600026e-01
Sepal.Width    0.09272109  0.11538776   0.05524354  0.03271020  4.042675e-02  2.933818e-01
Petal.Length   0.16751429  0.05524354   0.18518776  0.04266531 -3.038443e-01  7.216769e-02
Petal.Width    0.03840136  0.03271020   0.04266531  0.04188163 -1.295799e-01  1.508766e-01
LD1           -0.11459004  0.04042675  -0.30384427 -0.12957988  1.000000e+00 -5.921189e-16
LD2            0.16000259  0.29338181   0.07216769  0.15087661 -5.921189e-16  1.000000e+00
> 
> # convert pooled within-group covariance matrix to pooled within-group correlation matrix
> pooled_cor_matrix_LDs <- cov2cor(pooled_cov_matrix_LDs)
> pooled_cor_matrix_LDs
             Sepal.Length Sepal.Width Petal.Length Petal.Width           LD1           LD2
Sepal.Length    1.0000000   0.5302358    0.7561642   0.3645064 -2.225959e-01  3.108117e-01
Sepal.Width     0.5302358   1.0000000    0.3779162   0.4705346  1.190115e-01  8.636809e-01
Petal.Length    0.7561642   0.3779162    1.0000000   0.4844589 -7.060654e-01  1.677014e-01
Petal.Width     0.3645064   0.4705346    0.4844589   1.0000000 -6.331779e-01  7.372421e-01
LD1            -0.2225959   0.1190115   -0.7060654  -0.6331779  1.000000e+00 -5.921189e-16
LD2             0.3108117   0.8636809    0.1677014   0.7372421 -5.921189e-16  1.000000e+00
> round(pooled_cor_matrix_LDs, digits = 8)
             Sepal.Length Sepal.Width Petal.Length Petal.Width        LD1       LD2
Sepal.Length    1.0000000   0.5302358    0.7561642   0.3645064 -0.2225959 0.3108117
Sepal.Width     0.5302358   1.0000000    0.3779162   0.4705346  0.1190115 0.8636809
Petal.Length    0.7561642   0.3779162    1.0000000   0.4844589 -0.7060654 0.1677014
Petal.Width     0.3645064   0.4705346    0.4844589   1.0000000 -0.6331779 0.7372421
LD1            -0.2225959   0.1190115   -0.7060654  -0.6331779  1.0000000 0.0000000
LD2             0.3108117   0.8636809    0.1677014   0.7372421  0.0000000 1.0000000
> 
> # Check out the last 2 columns and the first 4 rows: These are the exact values noted in the SPSS output.
> # However, the signs in the first LD are reversed for some reason.
> # The same thing happened in another example I was using, which was a bit worrying.
```
Improve this answer
$\endgroup$
7
  • $\begingroup$ This answer can appear to be correct and valuable. But we cannot say until you post the results themselves - the figures that your code computes. $\endgroup$
    – ttnphns
    Sep 3, 2020 at 23:46
  • $\begingroup$ Thank you a bunch for the updates above and on that other thread. I'm trying going through the steps to see if they match. Anyway, I just wanted to update this "answer" with the computed output. $\endgroup$
    – PeteSugar
    Sep 8, 2020 at 21:15
  • $\begingroup$ This is super - thank you very much for your answer and comments. I don't think I would have thought of the approach you used. I feel like there may be a way of doing it without outputting the unwanted information in the matrix, but I don't know what that would be! $\endgroup$
    – Pratorum
    Sep 28, 2020 at 11:00
  • 1
    $\begingroup$ To summarise my understanding from your answer: The 'structure matrix' that the candisc package in R gives is the overall correlation matrix between the discriminant scores and the original variables. The 'structure matrix' that SPSS gives is the pooled within group correlation matrix between the discriminant scores and the original variables. $\endgroup$
    – Pratorum
    Sep 28, 2020 at 11:04
  • 1
    $\begingroup$ It doesn't matter that the signs are reversed for the first LD using your method. As noted here, the decision on whether these are positive or negative seems to be arbitrary. $\endgroup$
    – Pratorum
    Sep 28, 2020 at 11:35
2
$\begingroup$

In answer to Part 2 of my question, there was a missing square root from the final equation. This has now been corrected in the guide I was following here. I have included the R script giving the correct answer for the structure matrix here in case this is useful for anyone wanting a complete example of how to do this. I'm not sure if an answer is the best place for this; however, I didn't want to edit the question as it then wouldn't make sense.

The correct equation is:

$\bf R= {\it \sqrt{diag \bf (S_w)}} ^{-1} \bf S_w V$

In R this gives:

solve(sqrt(diag(diag(Sw)))) %*% Sw %*% LDA_V[,c(1,2)]

           [,1]       [,2]
[1,]  0.2225959 -0.3108117
[2,] -0.1190115 -0.8636809
[3,]  0.7060654 -0.1677014
[4,]  0.6331779 -0.7372421

Matching the structure matrix produced by SPSS. There is some sign reversal, though as noted above, this is not of concern.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.