rm(list=ls())
dati <- read.csv("https://www.statforbiology.com/pca/WheatQuality.csv", header=T)
head(dati)

##     Genotype Year   WPH    YB   TKW    PC
## 1 ARCOBALENO 1999 81.67 46.67 44.67 12.71
## 2 ARCOBALENO 2000 82.83 19.67 43.32 11.90
## 3 ARCOBALENO 2001 83.50 38.67 46.78 13.00
## 4 ARCOBALENO 2002 78.60 82.67 43.03 12.40
## 5       BAIO 1999 80.30 41.00 51.83 13.91
## 6       BAIO 2000 81.40 20.00 41.43 12.80

Qual è il principio di fondo?

Come nel caso delle componenti principali, dobbiamo prendere la matrice delle osservazioni originali e ‘ruotarla’, in modo da definire delle nuove variabili, dette variabili canoniche, che abbiano le caratteristiche desiderate. La rotazione è un’operazione di moltiplicazione matriciale, del tipo:

\[ VC = X \times V\]

dove X e VC hanno \(n\) righe e \(p\) colonne (\(n\) soggetti e \(p\) variabili), mentre V è una matrice quadrata \(v \times v\), che definisce la rotazione. Il nostro problema è quindi quello di trovare la matrice \(V\), in modo che le variabili ruotate (VC) permettano di osservare meglio le differenze tra gruppi.

Possiamo iniziare pensando che i gruppi sono facili da discriminare quando gli individui nel gruppo sono molto simili (bassa variabilità entro gruppo) e le medie dei gruppi sono invece molto diverse (alta variabilità tra gruppi). In questa situazione il valore di F nell’ANOVA univariata risulta molto alto e significativo. L’analisi delle variabili canoniche sfrutta questo principio: \(V\) viene scelta in modo che le variabili canoniche assicurino la minima variabilità entro gruppo e la massima variabilità tra gruppi.

Come si calcolano le variabili canoniche?

Per il calcolo, procediamo inizialmente alla standardizzazione dei dati. Anche se ciò non è necessario, ci consente di evitare che una variabile assuma un peso maggiore rispetto ad un’altra; indicheremo con \(Z\) la matrice dei dati standardizzati, della quale mostriamo solo le prime righe:

dataset <- dati[,3:6]
Z <- apply(dataset, 2, scale, center=T, scale=T)
head(Z)

##             WPH           YB         TKW         PC
## [1,]  0.3375814 -0.003873758 -0.37675661 -0.5193726
## [2,]  0.7681267 -1.154020460 -0.59544503 -1.5621410
## [3,]  1.0168038 -0.344657966 -0.03495471 -0.1460358
## [4,] -0.8018791  1.529655178 -0.64242254 -0.9184568
## [5,] -0.1709075 -0.245404565  0.78310197  1.0254693
## [6,]  0.2373683 -1.139963112 -0.90160881 -0.4035095

ATTENZIONE: alcuni programmi statistici (come SPSS) considerano anche un ulteriore tipo di standardizzazione che impiega al posto della deviazione standard totale la deviazione standard entro gruppi. Ne parleremo tra poco.

Indichiamo inoltre con U la matrice ottenuta da Z, sostituendo ad ogni dato in ogni colonna la media del gruppo (varietà) a cui il dato appartiene (matrice dei centroidi).

groups <- as.factor(dati[,1])
maovst <- manova(Z ~ groups)
U <- fitted(maovst)
head(U)

##          WPH           YB         TKW         PC
## 1 0.33015821  0.006775749 -0.41239472 -0.7865016
## 2 0.33015821  0.006775749 -0.41239472 -0.7865016
## 3 0.33015821  0.006775749 -0.41239472 -0.7865016
## 4 0.33015821  0.006775749 -0.41239472 -0.7865016
## 5 0.04529309 -0.337735787 -0.04224433  0.8548930
## 6 0.04529309 -0.337735787 -0.04224433  0.8548930

Con questi ingredienti siamo in grado di iniziare i calcoli. Calcoliamo la matrice delle devianze/codevianze totali (T):

T <- t(Z) %*% Z
T

##            WPH         YB       TKW        PC
## WPH  63.000000 -26.212412 34.638616 -1.293172
## YB  -26.212412  63.000000 -3.271389 -7.052713
## TKW  34.638616  -3.271389 63.000000 30.501366
## PC   -1.293172  -7.052713 30.501366 63.000000

Ora calcoliamo la matrice delle devianze/codevianze tra gruppi (B):

B <- t(U)%*%U
B

##            WPH        YB       TKW         PC
## WPH 20.7760120 -2.707135  7.986179 -0.7945792
## YB  -2.7071352 12.008904  2.639667 -8.5454789
## TKW  7.9861789  2.639667 27.190907 11.6352972
## PC  -0.7945792 -8.545479 11.635297 21.0150194

La matrice delle devianze/codevianze entro gruppo (W) è data da \(T - B = W\):

W <- T - B
W

##             WPH         YB       TKW         PC
## WPH  42.2239880 -23.505277 26.652437 -0.4985925
## YB  -23.5052765  50.991096 -5.911057  1.4927654
## TKW  26.6524371  -5.911057 35.809093 18.8660688
## PC   -0.4985925   1.492765 18.866069 41.9849806

Vediamo bene che esistono forti correlazioni tra le variabili, sia totali, che tra gruppi, che entro gruppo. Insomma, quattro variabili sono ridondanti per la discriminazione tra gruppi, visto che ognuna veicola parte dell’informazione già veicolata da altre. Possiamo calcolare il rapporto tra W e B, che rappresenta il rapporto tra le devianze/codevianze tra gruppi e quelle entro gruppo (una sorta di test F, per semplificare). Utilizzando la moltiplicazione matriciale possiamo scrivere:

WB <- solve(W)%*%B 
WB

##            WPH         YB        TKW         PC
## WPH  1.2427342 -0.4431553 -1.0966307 -0.5766991
## YB   0.4148842  0.1283142 -0.2258146 -0.3764203
## TKW -0.8168720  0.7037502  1.8275515  0.5566547
## PC   0.3481454 -0.5295936 -0.5490796  0.2569372

Gli autovalori della matrice \(WB\) sono:

A <- eigen(WB)$values
A

## [1] 2.1583182 0.7538314 0.4772379 0.0661495

Mentre gli autovettori associati con ciascun autovalore sono:

V1 <- eigen(WB)$vectors
V1

##            [,1]       [,2]       [,3]       [,4]
## [1,] -0.6147150 -0.4677403 0.73129521  0.1977323
## [2,] -0.1543769 -0.5681709 0.06455426  0.7153740
## [3,]  0.7218096 -0.2985051 0.13437782 -0.3705890
## [4,] -0.2780001  0.6076969 0.66556943  0.5583957

La matrice degli autovettori può essere utilizzata per ‘ruotare’ la matrice dei dati originali, analogamente a quanto avviene nell’analisi delle componenti principali. In questo modo otterremo quattro variabili ‘ruotate’ dette variabili canoniche centrate (VCC). Queste variabili hanno:

1. media = 0
2. varianze unitarie
3. correlazioni nulle

VCC <- Z %*% V1
colnames(VCC) <- c("VCC1", "VCC2", "VCC3", "VCC4")
apply(VCC, 2, mean)

##          VCC1          VCC2          VCC3          VCC4 
## -4.770490e-17  2.155394e-16  6.158268e-17  2.595580e-16

cor(VCC)

##              VCC1          VCC2          VCC3          VCC4
## VCC1 1.000000e+00  2.891080e-16  1.626754e-15  3.843447e-16
## VCC2 2.891080e-16  1.000000e+00 -1.021881e-15 -2.740982e-16
## VCC3 1.626754e-15 -1.021881e-15  1.000000e+00 -1.692481e-15
## VCC4 3.843447e-16 -2.740982e-16 -1.692481e-15  1.000000e+00

Se sottoponiamo ad ANOVA le quattro VCC, possiamo calcolare le deviazioni standard dei residui:

zseq <- apply(VCC, 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)
RMSE

## $VCC1
## [1] 0.3116788
## 
## $VCC2
## [1] 0.7872216
## 
## $VCC3
## [1] 0.9997169
## 
## $VCC4
## [1] 0.8114784

Queste deviazioni standard risulteranno utili tra poco.

I coefficienti canonici

Sappiamo che il calcolo degli autovettori di una matrice non ammette un’unica soluzione, ma infinite. In questo caso, gli autovettori indicati in precedenza sono stati scelti in modo che la loro norma fosse unitaria. Infatti:

t(V1) %*% V1

##              [,1]         [,2]        [,3]        [,4]
## [1,]  1.000000000 -0.009164206 -0.54753700 -0.65471505
## [2,] -0.009164206  1.000000000 -0.01438206 -0.04898404
## [3,] -0.547537002 -0.014382063  1.00000000  0.51263328
## [4,] -0.654715051 -0.048984043  0.51263328  1.00000000

Più comunemente, gli autovettori vengono scalati dividendoli proprio per la deviazione standard entro gruppi della variabile canonica centrata corrispondente. In questo modo si ottengono i coefficienti canonici standardizzati (\(Vst\)):

scaling <- diag(1/unlist(RMSE))
Vst <- V1 %*% scaling
Vst

##            [,1]       [,2]       [,3]       [,4]
## [1,] -1.9722706 -0.5941660 0.73150227  0.2436692
## [2,] -0.4953078 -0.7217420 0.06457254  0.8815688
## [3,]  2.3158763 -0.3791882 0.13441587 -0.4566838
## [4,] -0.8919442  0.7719516 0.66575788  0.6881214

I coefficienti canonici standardizzati sono usualmente utilizzati al posto degli autovettori V1 per ruotare le variabili originali, ottenendo le Variabili Canoniche Standardizzate, dette anche punteggi canonici (canonical scores):

VCst <- Z %*% Vst
colnames(VCst) <- c("VC1", "VC2", "VC3", "VC4")
cor(VCst)

##              VC1           VC2           VC3           VC4
## VC1 1.000000e+00  2.693048e-16  1.612146e-15  3.244759e-16
## VC2 2.693048e-16  1.000000e+00 -1.138091e-15 -2.357008e-16
## VC3 1.612146e-15 -1.138091e-15  1.000000e+00 -1.737112e-15
## VC4 3.244759e-16 -2.357008e-16 -1.737112e-15  1.000000e+00

zseq <- apply(VCst, 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)
RMSE

## $VC1
## [1] 1
## 
## $VC2
## [1] 1
## 
## $VC3
## [1] 1
## 
## $VC4
## [1] 1

MSB <- lapply(zseq, function(mod) anova(mod)[1,2])
MSW <- lapply(zseq, function(mod) anova(mod)[2,2])
unlist(MSB)/unlist(MSW)

##       VC1       VC2       VC3       VC4 
## 2.1583182 0.7538314 0.4772379 0.0661495

Vediamo più sopra che le variabili canoniche standardizzate: 1. non sono correlate tra loro (notare gli zero fuori dalla diagonale principale nella matrice di correlazione) 2. hanno varianze entro gruppi unitarie 3. i rapporti devianza tesi/devianza errore sono pari agli autovalori. La loro somma è pari a 2.158 + 0.754 + 0.477 + 0.066 = 3.455, il che eguaglia esattamente la somma della diagonale principale di W \(^{-1}\) B.

Quest’ultimo aspetto può essere ulteriormente evidenziato se calcoliamo per le variabili canoniche la matrice \(W_y^{-1} B_y\), in analogia alla matrice \(W^{-1} B\) già calcolata per le variabili originali.

maovst2 <- manova(VCst ~ groups)
U2 <- maovst2$fitted
T2 <- t(VCst)%*%VCst
B2 <- t(U2)%*%U2
W2 <- T2 - B2
solve(W2)%*%B2

##               VC1           VC2          VC3           VC4
## VC1  2.158318e+00 -1.259259e-16 5.735773e-16  7.809273e-16
## VC2 -7.139095e-16  7.538314e-01 7.998439e-19  1.625743e-16
## VC3 -2.669248e-15  3.283524e-16 4.772379e-01 -2.328594e-16
## VC4  1.119665e-15  5.061405e-16 5.192494e-16  6.614950e-02

Possiamo osservare che questa matrice conserva la traccia di \(W^{-1} B\), a dimostrazione che il contributo totale delle quattro variabili canoniche alla discriminazione tra gruppi è lo stesso delle quattro variabili originali, ma cambiano i rapporti. Ora la prima variabile canonica da sola garantisce una discriminazione pari al 62.5% del totale delle variabili originali (2.155/3.540 \(\cdot\) 100), mentre due variabili canoniche arrivano quasi all’85% del totale.

** Se quindi usiamo due variabili canoniche invece che le quattro variabili originali, riduciamo la dimensionalità del problema perdendo solo una piccola parte del potere discriminatorio iniziale. E’ un bel guadagno!**

A questi risultati si arriva in modi certamente più ortodossi di quelli qui accennati, ma per questo si rimanda ai testi specializzati (Sadocchi, 1986). A me basta che abbiate compreso i criteri fondamentali di calcolo.

Perchè i coefficienti canonici standardizzati differiscono da un programma all’altro?

In questo tutorial abbiamo iniziato l’analisi standardizzando i dati, cioè sottraendo da ogni dato la media della colonna e dividendo per la sua deviazione standard. Alcuni software statistici (come SPSS o WINSTAT 6.0) eseguono un diverso tipo di standardizzazione, dividendo ogni dato centrato per la deviazione standard entro gruppo. Nel codice sottostante calcoliamo, con quattro ANOVA, la deviazione standard del residuo per le quattro variabili e la utilizziamo per standardizzare i dati, producendo la matrice \(Zw\).

rm(scaling, T, T2, U, U2, W, W2, WB)

#Calculate the RMSE for the original variables
zseq <- apply(dataset, 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)

#Standardise by using the above RMSE
Zw <- t(t(apply(dati[,3:6], 2, function(x) scale(x, center=T, scale=F)) )/unlist(RMSE) )
head(Zw)

##             WPH           YB         TKW         PC
## [1,]  0.3599312 -0.003758425 -0.43619917 -0.5553320
## [2,]  0.8189811 -1.119662034 -0.68939104 -1.6702976
## [3,]  1.0841220 -0.334396532 -0.04046967 -0.1561468
## [4,] -0.8549680  1.484113053 -0.74378040 -0.9820472
## [5,] -0.1822225 -0.238098183  0.90665544  1.0964689
## [6,]  0.2530834 -1.106023212 -1.04385964 -0.4314470

Se iniziamo le analisi da Z\(w\) invece che da \(Z\), i coefficienti canonici standardizzati che otteniamo (\(Vw\)), sono diversi da quelli calcolati in precedenza:

#Calcolo coefficienti canonici a partire da Zw
T <- t(Zw) %*% Zw
maovst <- manova(Zw ~ groups)
U <- fitted(maovst)
B <- t(U)%*%U
W <- T - B
WB <- solve(W)%*%B 
V1w <- eigen(WB)$vectors
VCCw <- Zw %*% V1w
zseq <- apply(VCCw, 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)
scaling <- diag(1/unlist(RMSE))
Vw <- V1w %*% scaling
Vw

##            [,1]       [,2]       [,3]       [,4]
## [1,] -1.8498031 -0.5572715 0.68607988  0.2285387
## [2,] -0.5105070 -0.7438897 0.06655404  0.9086210
## [3,]  2.0002828 -0.3275147 0.11609849 -0.3944497
## [4,] -0.8341882  0.7219655 0.62264811  0.6435635

Possiamo tuttavia notare che il rapporto tra \(Vw\) è \(Vst\) è costante (0.938, 1.031, 0.864 e 0.935) e riflette il rapporto tra la deviazione standard totale e la deviazione standard entro gruppi delle variabili originali.

Vw/Vst

##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.9379053 0.9379053 0.9379053 0.9379053
## [2,] 1.0306864 1.0306864 1.0306864 1.0306864
## [3,] 0.8637261 0.8637261 0.8637261 0.8637261
## [4,] 0.9352471 0.9352471 0.9352471 0.9352471

sdTot <- apply(dataset, 2, sd)
zseq <- apply(dati[3:6], 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)
unlist(RMSE)/sdTot

##       WPH        YB       TKW        PC 
## 0.9379053 1.0306864 0.8637261 0.9352471

Se ruotiamo \(Zw\) utilizzando \(Vw\) arriviamo alle stesse variabili canoniche già determinate in precedenza (anche qui mostriamo solo le prime sei righe).

head(Zw %*% Vw)

##            [,1]       [,2]       [,3]        [,4]
## [1,] -1.0731535 -0.4558524 -0.1497271 -0.10648959
## [2,] -0.9289930 -0.6036012 -0.6326765 -1.63319215
## [3,] -1.7853954 -0.4548743  0.6196159 -0.14060305
## [4,]  0.1553135 -1.0929723 -1.1856243  0.81447716
## [5,]  1.3575325  0.7733359  0.6471100  0.09003153
## [6,] -1.6316285  0.7121127 -0.2898050 -0.81303000

Oltre ai coefficienti canonici standardizzati, abbiamo anche i cosidetti raw coefficients (\(Vr\)) che si ottengono iniziando le analisi dalla matrice dei dati centrati e non standardizzati (\(Zc\)).

#Calcolo coefficienti canonici a partire da Zc
rm(scaling, T, U, W, WB)
Zc <- apply(dati[,3:6], 2, function(x) scale(x, center=T, scale=F))
T <- t(Zc) %*% Zc
maovst <- manova(Zc ~ groups)
U <- fitted(maovst)
B <- t(U)%*%U
W <- T - B
WB <- solve(W)%*%B 
V1r <- eigen(WB)$vectors
VCCr <- Zc %*% V1r
zseq <- apply(VCCr, 2, function(x) lm(x ~ Genotype, data = dati))
RMSE <- lapply(zseq, function(mod) summary(mod)$sigma)
scaling <- diag(1/unlist(RMSE))
Vr <- V1r %*% scaling
Vr

##             [,1]        [,2]        [,3]        [,4]
## [1,] -0.73202750 -0.22053053 0.271504213  0.09044022
## [2,] -0.02109913 -0.03074479 0.002750663  0.03755309
## [3,]  0.37515210 -0.06142523 0.021774217 -0.07397885
## [4,] -1.14826063  0.99378600 0.857075589  0.88586568

Anche in questo caso, se ruotiamo \(Zc\) utilizzando \(Vr\), arriviamo agli stessi punteggi canonici.

head(Zc %*% Vr)

##            [,1]       [,2]       [,3]        [,4]
## [1,] -1.0731535 -0.4558524 -0.1497271 -0.10648959
## [2,] -0.9289930 -0.6036012 -0.6326765 -1.63319215
## [3,] -1.7853954 -0.4548743  0.6196159 -0.14060305
## [4,]  0.1553135 -1.0929723 -1.1856243  0.81447716
## [5,]  1.3575325  0.7733359  0.6471100  0.09003153
## [6,] -1.6316285  0.7121127 -0.2898050 -0.81303000

Insomma i coefficienti canonici possono cambiare perchè la matrice di partenza può essere diversa. Tuttavia le variabili canoniche non cambiano: esse possono essere ottenute (indifferentemente):

1. moltiplicando \(Z\) per \(Vst\) (come abbiamo fatto più sopra);
2. moltiplicando \(Zw\) per \(Vw\);
3. moltiplicando \(Zc\) per \(Vr\) (“raw coefficients”).

Punteggi canonici dei centroidi

Riepilogando, i coefficienti canonici sono un metodo per riorganizzare le nostre osservazioni sperimentali in uno spazio diverso, con meno dimensioni, ma senza perdere la capacità di differenziare i gruppi (ricordiamo che nel caso della PCA si cercava invece di non perdere la distanza tra individui!).

In realtà, più che le variabili canoniche, ci interessa calcolare i punteggi canonici dei centroidi di ogni gruppo (le medie delle variabili originali), cioè la loro posizione nel nuovo spazio di riferimento dopo la rotazione. Le nuove coordinate dei centroidi (punteggi canonici dei centroidi) possono essere ottenute allo stesso modo delle variabili canoniche, moltiplicando la matrice delle medie delle variabili standardizzate (con 16 elementi, cioè le sedici medie per ogni combinazione varietà x anno) per la matrice dei coefficienti canonici standardizzati (Vst). I punteggi, ottenuti con R, sono i seguenti:

medie <- aov(Z ~ groups-1)$coefficients
scores <- medie %*% Vst
scores

##                         [,1]       [,2]       [,3]        [,4]
## groupsARCOBALENO  -0.9080571 -0.6518251 -0.3371030 -0.26645191
## groupsBAIO        -0.7823965  0.8928011  0.5747970  0.32086167
## groupsCLAUDIO     -0.5496314 -1.3845288  0.3590569  0.16423169
## groupsCOLORADO    -1.1481765  1.4231696 -0.5535939 -0.13201115
## groupsCOLOSSEO     1.0654126 -1.3108147 -0.1201515 -0.05709275
## groupsCRESO        0.3070820 -0.3947049  0.6469757 -0.09565658
## groupsDUILIO       0.2985911 -0.3083817  0.3302409  0.39103654
## groupsGIOTTO       0.5551808 -0.1276729 -0.6353720 -0.14515354
## groupsGRAZIA      -2.7043586  0.1540981  0.2141260 -0.09434961
## groupsIRIDE       -0.9975478 -0.6185305 -0.7408252 -0.13364342
## groupsMERIDIANO    0.7232796  0.1328308 -0.1189112  0.31997541
## groupsSANCARLO    -0.6322898  0.7065535  0.9005722 -0.01388980
## groupsSIMETO       3.1209198  0.4928771  0.5853148 -0.18719341
## groupsSOLEX       -0.3851960 -0.1854706  0.1043439  0.02715539
## groupsTORREBIANCA  0.9167304  0.5916795  0.2310080 -0.37993393
## groupsVERDI        1.1204575  0.5879195 -1.4404785  0.28211539

La struttura canonica

Immaginiamo di calcolare, oltre alle variabili canoniche, anche i punteggi canonici relativi alla matrice delle medie delle variabili standardizzate (U) e quella dei residui (differenza tra i dati standardizzati e le loro medie, cioè Z - U). Otteniamo rispettivamente \(Vu\) e \(Vzu\)

VU <- U%*%Vst
head(VU)

##        [,1]       [,2]        [,3]      [,4]
## 1 -7.183974 -0.1496154 -0.08796465  1.099194
## 2 -7.183974 -0.1496154 -0.08796465  1.099194
## 3 -7.183974 -0.1496154 -0.08796465  1.099194
## 4 -7.183974 -0.1496154 -0.08796465  1.099194
## 5  2.490094  6.2612890 -0.01564151 -6.383678
## 6  2.490094  6.2612890 -0.01564151 -6.383678

VZU  <- (Z-U)%*%Vst
head(VZU)

##           [,1]       [,2]        [,3]       [,4]
## [1,]  6.110820 -0.3062371 -0.06176242 -1.2056833
## [2,]  6.254981 -0.4539858 -0.54471185 -2.7323858
## [3,]  5.398578 -0.3052589  0.70758051 -1.2397967
## [4,]  7.339287 -0.9433569 -1.09765964 -0.2847165
## [5,] -1.132562 -5.4879531  0.66275151  6.4737092
## [6,] -4.121723 -5.5491763 -0.27416353  5.5706477

Le nuove variabili così ottenute, non sono non sono correlate tra loro, ma hanno un certo grado di correlazione con le variabili originali. La matrice di correlazione delle variabili canoniche (\(VCst\)) con le variabili originali (Z) prende il nome di Struttura Canonica Totale.

cor(Z, VCst)

##            VC1        VC2        VC3         VC4
## WPH -0.3059273 -0.4482850  0.7209655 -0.43088021
## YB   0.1965463 -0.4682288 -0.3028492  0.80648226
## TKW  0.5320720 -0.2548998  0.8064667 -0.03920395
## PC   0.2096499  0.5894317  0.6679134  0.40312384

Allo stesso modo, la correlazione tra U e \(Vu\) prende il nome di Struttura Canonica tra gruppi (‘between’):

cor(U, VU)

##           [,1]        [,2]      [,3]       [,4]
## WPH  0.1113207 -0.02298862 0.7867065 -0.2081955
## YB  -0.3336955 -0.96964998 0.2556828  0.9840392
## TKW  0.8358955 -0.34964079 0.7717052 -0.0294356
## PC   0.7895205  0.45777788 0.1090482 -0.6304306

Ancora, la correlazione tra i residui (Z - U) e i loro punteggi canonici (\(Vzu\)) prende il nome di Struttura Canonica entro gruppi (‘within’):

cor((Z-U), VZU)

##           [,1]        [,2]      [,3]        [,4]
## WPH  0.1308540 -0.01011711 0.6889033 -0.18997474
## YB  -0.4177260 -0.96838132 0.4180063  0.98299197
## TKW  0.7929795 -0.34729954 0.6239354 -0.02966873
## PC  -0.1111922 -0.00194042 0.4295819  0.14763277

La struttura canonica è estremamente importante perchè rappresenta il collegamento tra le variabili canoniche e le variabili originali e quindi può aiutarci ad interpretare i risultati dell’analisi.

---

Esame dei punteggi canonici

L’esame dei punteggi canonici mostra che SIMETO e GRAZIA sono più distanti sulla prima variabile canonica, mentre CLAUDIO e COLORADO sono molto distanti sulla seconda. Inoltre, alcune varietà sono caratterizzate da punteggi positivi o negativi su entrambe le variabili canoniche, mentre altre sono caratterizzate da situazioni opposte (positive su una variabile e negative sull’altra). In ogni caso, emergono chiaramente le differenze tra varietà, sotto forma di distanze su una o entrambe le variabili canoniche (anche se le distanze euclidea tra centroidi non rappresentano al meglio le distanze euclidee tra individui, come nel caso delle componenti principali).

Esame dei coefficienti canonici e della struttura canonica

Le considerazioni appena esposte non hanno molto senso se non si riesce a comprendere come ogni variabile canonica si rapporta con le variabili originali. Per questo è necessario considerare, oltre ai punteggi canonici delle medie, anche la struttura canonica o i coefficienti canonici.

In particolare, questi ultimi costituiscono il contributo unico (peso) che ogni variabile originale ha dato nella creazione dei punteggi canonici delle medie (si veda in precedenza l’esempio d’uso delle funzioni canoniche). Per esempio, considerando il primo vettore dei coefficienti canonici (quello che crea la prima variabile canonica) possiamo notare che l’alto coefficiente relativo al peso di 1000 cariossidi implica che le varietà con valori elevati per questa caratteristica avranno probabilmente un elevato valore positivo sulla prima variabile canonica. Allo stesso modo, il valore negativo per il coefficiente del peso ettolitrico implica che le varietà con un elevato peso ettolitrico avranno probabilmente valori negativi sulla prima variabile canonica.

Considerazioni analoghe potrebbero essere fatte osservando la struttura canonica.

Visualizzazione geometrica dei risultati dell’analisi

Come già detto più volte, il calcolo delle variabili canoniche ci permette di riorganizzare nello spazio le nostre osservazioni, in modo più agevole da comprendere e spiegare. Ovviamente, per poter rappresentare in un grafico il risultato di questa riorganizzazione dobbiamo limitarci alle prime due variabili canoniche, anche se abbiamo già visto come con esse, in questo esempio, si riesca a conservare un potere discriminante pari all’85% di quello originale.

Possiamo quindi utilizzare un biplot, nel quale vengono riportati i punteggi canonici delle medie, insieme ai coefficienti canonici standardizzati, in forma di vettori geometrici (uno per ogni variabile originale). Questo biplot (CVA biplot) non ha la stessa interpretazione geometrica di un PCA biplot, perchè è frutto di elaborazioni algebriche/geometriche diverse. Per la sua interpretazione dovremo considerare che:

1 - La distanza euclidea tra i soggetti (le varietà) nel biplot approssima la loro distanza di Mahalanobis nella scala orginale ed è quindi la più adatta per valutare le differenze tra gruppi. 2 - I vettori geometrici ‘tirano’ le medie lungo la loro direzione in modo tanto più intenso quanto più è alto il loro valore originale della caratteristica rilevata. In questo modo si può facilmente comprendere quali caratteristiche contribuiscono più delle altre alla discriminazione.

Si può osservare ad esempio che Colorado, Baio e Sancarlo sono distinguibili per il loro alto contenuto in proteine (il loro punteggio canonico giace lungo la direzione del vettore corrispondente) associato ad un basso grado di bianconatura (il vettore relativo a questa caratteristica punta in direzione opposta). D’altra parte, Claudio e Colosseo sono discriminate dal loro basso contenuto in proteine e da un’alta percentuale di bianconatura.

Grazia ha mostrato un contenuto di proteine abbastanza alto, associato ad un buon peso ettolitrico (la sua posizione in figura risulta dalla combinazione dei due vettori corrispondenti allecaratteristiche sopracitate) ed un basso peso dei 1000 semi.

Arcobaleno è Iride hanno invece mostrato un alto peso ettolitrico, un alto grado di bianconatura, un basso contenuto di proteine e un basso peso dei 1000 semi.

Si deve comunque notare che alcune situazioni sono poco “chiare”; in particolare alcune varietà si trovano vicine all’origine degli assi, il che significa che non hanno mostrato caratteristiche di spicco (in senso negativo o positivo) e si sono mantenute vicino alla media. Inoltre, Verdi mostra un punteggio abbastanza poco chiaro, dato che la sua posizione sul biplot non evidenzia l’alto grado di bianconatura che era invece chiaro dai dati originali (non dimentichiamoci che parte dell’informazione iniziale è comunque persa nel processo di combinazione lineare delle variabili originali).

Insomma, anche se l’analisi canonica riveste un’indubbia utilità per riassumere i dati, essa non potrà mai sostituire un’accurata analisi puntiforme dei risultati dell’esperimento, altrimenti il rischio di perdere o confondere le informazioni può essere abbastanza serio.