L’estrazione dei fattori

17. L’estrazione dei fattori#

L’analisi fattoriale mira a descrivere in modo parsimonioso le relazioni tra un grande numero di item. L’obiettivo è identificare un piccolo numero di variabili latenti che, quando controllate, rendono uguali a zero le correlazioni parziali tra gli item. Una volta determinato il numero dei fattori comuni, è possibile stimare le saturazioni fattoriali, che corrispondono alle correlazioni o covarianze tra gli item e i fattori.

In termini matriciali, il modello multifattoriale si scrive

\[ \boldsymbol{\Sigma} =\boldsymbol{\Lambda} \boldsymbol{\Phi} \boldsymbol{\Lambda}^{\mathsf{T}} + \boldsymbol{\Psi} \]

dove \(\boldsymbol{\Phi}\) è la matrice di ordine \(m \times m\) di varianze e covarianze tra i fattori comuni e \(\boldsymbol{\Psi}\) è una matrice diagonale di ordine \(p\) con le unicità delle variabili.

In questo capitolo descriveremo alcuni dei metodi che possono essere usati per stimare \(\boldsymbol{\Lambda}\). Esamineremo il metodo delle componenti principali, il metodo dei fattori principali, il metodo dei fattori principali iterato e il metodo di massima verosimiglianza.

17.1. Metodo delle componenti principali#

L’analisi fattoriale eseguita mediante il metodo delle componenti principali, nonostante il nome, non è un’analisi delle componenti principali. Il metodo delle componenti principali costituisce invece un’applicazione del teorema di scomposizione spettrale di una matrice. Il teorema spettrale afferma che, data la matrice simmetrica \(\textbf{S}_{p \times p}\), è sempre possibile trovare una matrice \(\textbf{C}_{p \times p}\) ortogonale tale che \( \textbf{S} = \textbf{C}\textbf{D}\textbf{C}^{\mathsf{T}} \) con D diagonale. Il teorema specifica inoltre che gli elementi presenti sulla diagonale di D sono gli autovalori di S, mentre le colonne di C rappresentano i rispettivi autovettori normalizzati associati agli autovalori di S.

Facciamo un esempio numerico utilizzando i dati discussi da Rencher(2002). Brown, Williams e Barlow (1984) hanno raccolto le valutazioni di una ragazza dodicenne relativamente a sette persone di sua conoscenza. Ciascuna persona veniva valutata su una scala a nove punti rispetto a cinque variabili: kind, intelligent, happy, likeable e just. La matrice di correlazione per tali variabili è riportata di seguito:

R <- matrix(c(
  1.000, .296, .881, .995, .545,
  .296, 1.000, -.022, .326, .837,
  .881, -.022, 1.000, .867, .130,
  .995, .326, .867, 1.000, .544,
  .545, .837, .130, .544, 1.000
),
ncol = 5, byrow = T, dimnames = list(
  c("K", "I", "H", "L", "J"),
  c("K", "I", "H", "L", "J")
)
)
R

A matrix: 5 × 5 of type dbl
	K	I	H	L	J
K	1.000	0.296	0.881	0.995	0.545
I	0.296	1.000	-0.022	0.326	0.837
H	0.881	-0.022	1.000	0.867	0.130
L	0.995	0.326	0.867	1.000	0.544
J	0.545	0.837	0.130	0.544	1.000

Gli autovalori e gli autovettori si calcolano con la funzione eigen():

e <- eigen(R)
print(e)

eigen() decomposition
$values
[1] 3.2633766259 1.5383820947 0.1679692814 0.0300298228 0.0002421752

$vectors
           [,1]       [,2]        [,3]       [,4]       [,5]
[1,] -0.5367301 -0.1859819 -0.18991841 -0.1247931  0.7910052
[2,] -0.2874672  0.6505666  0.68488713 -0.1198141  0.1034406
[3,] -0.4343545 -0.4736877  0.40694897  0.6136634 -0.2115794
[4,] -0.5373909 -0.1692797 -0.09532905 -0.6293896 -0.5266275
[5,] -0.3896544  0.5377158 -0.56583170  0.4442491 -0.2037363

Come indicato in precedenza, la matrice R può essere espressa come

\[ \textbf{R} = \textbf{C}\textbf{D}\textbf{C}^{\mathsf{T}} \]

e$vectors %*% diag(e$values) %*% t(e$vectors)

A matrix: 5 × 5 of type dbl
1.000	0.296	0.881	0.995	0.545
0.296	1.000	-0.022	0.326	0.837
0.881	-0.022	1.000	0.867	0.130
0.995	0.326	0.867	1.000	0.544
0.545	0.837	0.130	0.544	1.000

Esaminiamo ora gli autovalori. I primi due autovalori spiegano da soli il 96% della varianza campionaria:

(e$values[1] + e$values[2]) / 5

0.960351744117911

Usando i primi due autovalori e i primi due autovettori è dunque possibile riprodurre in maniera soddisfacente la matrice R operando nel contempo una riduzione di dimensionalità dei dati.

Per fattorizzare \(\textbf{R} = \textbf{C}\textbf{D}\textbf{C}^{\mathsf{T}}\) nella forma \(\hat{\boldsymbol{\Lambda}} \hat{\boldsymbol{\Lambda}}^{\mathsf{T}}\) iniziamo a scrivere

\[ \textbf{D}= \textbf{D}^{1/2} \textbf{D}^{1/2} \]

dove

\[\begin{split} \textbf{D}^{1/2} = \left[ \begin{array}{ c c c c } \sqrt{\theta_1} & 0 & \dots & 0 \\ 0 & \sqrt{\theta_2} & \dots & 0 \\ \dots & \dots & & \dots \\ 0 & 0 & \dots & \sqrt{\theta_p} \end{array} \right] \end{split}\]

Viene qui usata la notazione \(\theta\) per denotare gli autovalori anziché il tradizionale \(\lambda\) per evitare la confusione con la notazione \(\lambda_{jl}\) usata per le saturazioni fattoriali. In questo modo, possiamo scrivere

\[\begin{split} \begin{equation} \begin{aligned} \textbf{R} &= \textbf{C}\textbf{D}\textbf{C}^{\mathsf{T}}\notag\\ &= \textbf{C}\textbf{D}^{1/2}\textbf{D}^{1/2}\textbf{C}^{\mathsf{T}}\notag\\ &= (\textbf{C}\textbf{D}^{1/2}) (\textbf{C}\textbf{D}^{1/2})^{\mathsf{T}}. \end{aligned} \end{equation} \end{split}\]

Non possiamo però limiarci a definire \(\hat{\boldsymbol{\Lambda}}=\textbf{C}\textbf{D}^{1/2}\) in quanto \(\textbf{C}\textbf{D}^{1/2}\) è di ordine \(p \times p\) e non otteniamo quindi una riduzione di dimensioni. Quello che cerchiamo è una matrice \(\hat{\boldsymbol{\Lambda}}\) di ordine \(p \times m\) con \(m < p\). Dunque, definiamo la matrice \(\textbf{D}_1= \text{diag}(\theta_1, \theta_2, \dots, \theta_m)\) come la la matrice contenente gli \(m\) autovalori più grandi di R e \(\textbf{C}_1=( \textbf{c}_1, \textbf{c}_2, \dots, \textbf{c}_m)\) come la matrice contenente i rispettivi autovettori.

Mediante il metodo delle componenti principali, le saturazioni fattoriali \(\hat{\boldsymbol{\Lambda}}\) vengono quindi stimate nel modo seguente:

\[\begin{split} \begin{equation} \begin{aligned} \hat{\boldsymbol{\Lambda}} &= \textbf{C}_1 \textbf{D}_1^{1/2}\notag\\ &= (\sqrt{\theta_1} \textbf{c}_1, \sqrt{\theta_2} \textbf{c}_2, \dots, \sqrt{\theta_m} \textbf{c}_m) \end{aligned} \end{equation} \end{split}\]

Per l’esempio presente, con \(m=2\) e \(p=5\), avremo

\[\begin{split} \left[ \begin{array}{ c c } \hat{\lambda}_{11} & \hat{\lambda}_{12} \\ \hat{\lambda}_{21} & \hat{\lambda}_{22} \\ \hat{\lambda}_{31} & \hat{\lambda}_{32} \\ \hat{\lambda}_{41} & \hat{\lambda}_{42} \\ \hat{\lambda}_{51} & \hat{\lambda}_{52} \end{array} \right] = \left[ \begin{array}{ c c } c_{11} & c_{12} \\ c_{21} & c_{22} \\ c_{31} & c_{32} \\ c_{41} & c_{42} \\ c_{51} & c_{52} \end{array} \right] \left[ \begin{array}{ c c } \sqrt{\theta_1} & 0\\ 0 &\sqrt{\theta_2} \end{array} \right] \end{split}\]

Le saturazioni fattoriali stimate sono dunque uguali a

\[\begin{split} \left[ \begin{array}{ c c } \sqrt{\theta_1}c_{11} & \sqrt{\theta_2}c_{12} \\ \sqrt{\theta_1}c_{21} & \sqrt{\theta_2}c_{22} \\ \sqrt{\theta_1}c_{31} & \sqrt{\theta_2}c_{32} \\ \sqrt{\theta_1}c_{41} & \sqrt{\theta_2}c_{42} \\ \sqrt{\theta_1}c_{51} & \sqrt{\theta_2}c_{52} \end{array} \right] \end{split}\]

Svolgendo i calcoli con \(\textsf{R}\) otteniamo:

L <- cbind(
  e$vectors[, 1] * sqrt(e$values[1]),
  e$vectors[, 2] * sqrt(e$values[2])
)

round(L, 3)

A matrix: 5 × 2 of type dbl
-0.970	-0.231
-0.519	0.807
-0.785	-0.588
-0.971	-0.210
-0.704	0.667

La matrice di correlazione riprodotta (con le comunalità sulla diagonale principale) diventa

R_hat <- round(L %*% t(L), 3)
R_hat

A matrix: 5 × 5 of type dbl
0.993	0.317	0.896	0.990	0.529
0.317	0.921	-0.067	0.335	0.904
0.896	-0.067	0.961	0.885	0.160
0.990	0.335	0.885	0.987	0.543
0.529	0.904	0.160	0.543	0.940

La matrice di correlazioni residue (con le specificità sulla diagonale principale) è la seguente.

R - R_hat

A matrix: 5 × 5 of type dbl
	K	I	H	L	J
K	0.007	-0.021	-0.015	0.005	0.016
I	-0.021	0.079	0.045	-0.009	-0.067
H	-0.015	0.045	0.039	-0.018	-0.030
L	0.005	-0.009	-0.018	0.013	0.001
J	0.016	-0.067	-0.030	0.001	0.060

Possiamo ora capire il motivo del nome metodo delle componenti principali: le saturazioni fattoriali sono proporzionali agli autovettori di \(\textbf{R}\). Tuttavia, l’interpretazione è diversa da quella che viene assegnata ai risultati dell’analisi delle componenti principali.

È possibile condurre l’analisi fattoriale con il metodo delle componenti principali sia utilizzando la matrice \(\textbf{S}\) di varianze-covarianze sia la matrice \(\textbf{R}\) delle correlazioni. Tuttavia, le soluzioni ottenute usando \(\textbf{S}\) o \(\textbf{R}\) non sono legate da una relazione algebrica: il metodo delle componenti principali non è invariante rispetto ai cambiamenti di scala delle osservazioni. Un altro svantaggio del metodo delle componenti principali è che non fornisce un test di bontà di adattamento. Tale test può essere invece svolto quando la soluzione viene trovata con il metodo della massima verosimiglianza.

17.2. Metodo dei fattori principali#

Il metodo dei fattori principali (principal factor method, anche detto principal axis method) è uno dei metodi maggiormente usati per la stima delle saturazioni fattoriali e delle comunalità. Il metodo delle componenti principali che abbiamo usato in precedenza trascura la specificità \(\boldsymbol{\Psi}\) e si limita a fattorializzare le covarianze di S o le correlazioni di R. Il metodo dei fattori principali affronta questo problema ricorrendo ad una procedura simile al metodo delle componenti principali, nella quale però viene utilizzata una matrice ridotta di varianze-covarianze \(\textbf{S} - \hat{\boldsymbol{\Psi}}\) in cui una stima delle comunalità viene sostituita alle varianze presenti sulla diagonale principale della matrice S o R.

Nel caso della matrice ridotta di correlazioni \(\textbf{R} - \hat{\boldsymbol{\Psi}}\), per la comunalità \(i\)-esima \(\sum_{j}\lambda_{ij}^2\) si sceglie il quadrato del coefficiente di correlazione multipla tra \(Y_i\) e tutte le altre \(p-1\) variabili. Tale valore si può trovare nel modo seguente

\[ \hat{h}^2_i=R^2_i=1-\frac{1}{r^{ii}}, \]

dove \(r^{ii}\) è l’elemento diagonale \(i\)-esimo di \(\textbf{R}^{-1}\). Nel caso di una matrice ridotta di varianze-covarianze \(\textbf{S} - \hat{\boldsymbol{\Psi}}\), le comunalità possono essere stimate calcolando

\[ \hat{h}_i^2=s_{ii}-\frac{1}{r^{ii}}, \]

dove \(s_{ii}\) è l’elemento diagonale \(i\)-esimo di \(\textbf{S}\).

Affinché le stime comunalità possano essere calcolate come descritto sopra, la matrice \(\textbf{R}\) non deve essere singolare. Se \(\textbf{R}\) è singolare, per la stima della comunalità \(i\)-esima, \(\hat{h}^2_i\), si utilizza il valore assoluto del più elevato coefficiente di correlazione lineare tra \(Y_i\) e le altre variabili.

Scelta la stima della comunalità, la matrice ridotta di varianze-covarianze si ottiene sostituendo le stime delle comunalità alle varianze sulla diagonale principale della matrice \(\textbf{S}\):

\[\begin{split} \textbf{S} - \hat{\boldsymbol{\Psi}} = \left[ \begin{array}{ c c c c } \hat{h}^2_1 & s_{12} & \dots & s_{1p} \\ s_{21} & \hat{h}^2_2 & \dots & s_{2p} \\ \dots & \dots & & \dots\\ s_{p1} & s_{p2} & \dots & \hat{h}^2_p \end{array} \right] \end{split}\]

In maniera equivalente, la matrice ridotta di correlazioni si ottiene nel modo seguente:

\[\begin{split} \textbf{R} - \hat{\boldsymbol{\Psi}} = \left[ \begin{array}{ c c c c } \hat{h}^2_1 & r_{12} & \dots & r_{1p} \\ r_{21} & \hat{h}^2_2 & \dots & r_{2p} \\ \dots & \dots & & \dots\\ r_{p1} & r_{p2} & \dots & \hat{h}^2_p \end{array} \right] \end{split}\]

Facciamo ora un esempio concreto usando la matrice di correlazione dell’esempio precedente. Quale stima della comunalità \(i\)-esima, useremo il valore assoluto più elevato nella riga \(i\)-esima della matrice R. Per i dati dell’esempio, le stime delle comunalità sono dunque pari a \(0.995\), \(0.837\), \(0.881\), \(0.995\) e \(0.837\).

Inserendo tali valori nella diagonale principale, otteniamo la matrice ridotta delle correlazioni \(\textbf{R} - \hat{\boldsymbol{\Psi}}\):

R1 <- R
h.hat <- c(.995, .837, .881, .995, .837)
R1[cbind(1:5,1:5)] <- h.hat
R1

A matrix: 5 × 5 of type dbl
	K	I	H	L	J
K	0.995	0.296	0.881	0.995	0.545
I	0.296	0.837	-0.022	0.326	0.837
H	0.881	-0.022	0.881	0.867	0.130
L	0.995	0.326	0.867	0.995	0.544
J	0.545	0.837	0.130	0.544	0.837

Gli autovalori della matrice ridotta di correlazioni \(\textbf{R} - \hat{\boldsymbol{\Psi}}\) sono:

ee <- eigen(R1)
print(round(ee$values, 3))

[1]  3.202  1.394  0.029  0.000 -0.080

La somma degli autovalori è uguale a

print(sum(ee$values))

[1] 4.545

I primi due autovalori di \(\textbf{R} - \hat{\boldsymbol{\Psi}}\) sono:

print(round(ee$vectors[, 1:2], 3))

       [,1]   [,2]
[1,] -0.548 -0.177
[2,] -0.272  0.656
[3,] -0.431 -0.461
[4,] -0.549 -0.159
[5,] -0.373  0.549

Moltiplicando tali valori per la radice quadrata dei rispettivi autovalori si ottengono le stime delle saturazioni fattoriali:

round(ee$vectors[,1:2] %*% sqrt(diag(ee$values[1:2])), 3)

A matrix: 5 × 2 of type dbl
-0.981	-0.209
-0.487	0.774
-0.772	-0.544
-0.982	-0.187
-0.667	0.648

Tale risultato replica quello riportato da Rencher (2002).

17.3. Metodo dei fattori principali iterato#

Solitamente, per migliorare la stima della comunalità, la diagonale della matrice \(\textbf{S} - \hat{\boldsymbol{\Psi}}\) o \(\textbf{R} - \hat{\boldsymbol{\Psi}}\) viene ottenuta per iterazione. Dopo avere trovato \(\hat{\boldsymbol{\Lambda}}\) a partire da \(\textbf{S} - \hat{\boldsymbol{\Psi}}\) o \(\textbf{R} - \hat{\boldsymbol{\Psi}}\) come indicato in precedenza, utilizzando le stime delle saturazioni fattoriali così ottenute possiamo stimare le comunalità nel modo seguente:

\[\hat{h}^2_i = \sum_{i=1}^m \hat{\lambda}_{ij}^2.\]

I valori di \(\hat{h}^2_i\) così trovati vengono quindi sostituiti nella diagonale della matrice ridotta \(\textbf{S} - \hat{\boldsymbol{\Psi}}\) o \(\textbf{R} - \hat{\boldsymbol{\Psi}}\). A partire da questa nuova matrice, usando il metodo descritto in precedenza, possiamo ottenere una nuova stima delle saturazioni fattoriali \(\hat{\boldsymbol{\Lambda}}\). Mediante questa nuova stima di \(\hat{\boldsymbol{\Lambda}}\), possiamo procedere ad una nuova stima delle comunalità. Tale processo continua iterativamente sino alla convergenza. Gli autovalori e gli autovettori della versione finale di \(\textbf{S} - \hat{\boldsymbol{\Psi}}\) o \(\textbf{R} - \hat{\boldsymbol{\Psi}}\) vengono infine usati per stimare i pesi fattoriali. Il metodo dei fattori principali iterato e il metodo delle componenti principali producono risultati molto simili quando \(m\) assume un piccolo valore (questo si verifica quando le correlazioni sono alte) e quando \(p\) (il numero delle variabili) è grande.

17.3.1. Casi di Heywood#

Tra gli inconvenienti del metodo dei fattori principali iterato vi è il fatto che può talvolta portare a soluzioni inammissibili (quando viene fattorizzata la matrice R) caratterizzate da valori di comunalità maggiori di uno (caso di Heywood). Se \(\hat{h}^2_i > 1\) allora \(\hat{\psi}_i < 0\) il che è chiaramente assurdo in quanto una varianza non può assumere un valore negativo. Solitamente, quando la stima di una comunalità è maggiore di uno, il processo iterativo viene interrotto e il programma riporta che non può essere trovata una soluzione.

Per fare un esempio, usiamo la funzione fa() contenuta nel pacchetto psych. Tale funzione per trovare la soluzione dei fattori principali mediante il metodo iterativo.

pa <- fa(R, nfactors = 2, rotate = "none", fm = "pa")
pa

Warning message in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
“The estimated weights for the factor scores are probably incorrect.  Try a different factor score estimation method.”

Warning message in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, :
“An ultra-Heywood case was detected.  Examine the results carefully”

Factor Analysis using method =  pa
Call: fa(r = R, nfactors = 2, rotate = "none", fm = "pa")
Standardized loadings (pattern matrix) based upon correlation matrix
   PA1   PA2   h2     u2 com
K 0.98 -0.21 1.01 -0.008 1.1
I 0.48  0.74 0.77  0.230 1.7
H 0.78 -0.56 0.92  0.085 1.8
L 0.98 -0.19 0.99  0.010 1.1
J 0.69  0.69 0.95  0.049 2.0

                       PA1  PA2
SS loadings           3.22 1.41
Proportion Var        0.64 0.28
Cumulative Var        0.64 0.93
Proportion Explained  0.70 0.30
Cumulative Proportion 0.70 1.00

Mean item complexity =  1.5
Test of the hypothesis that 2 factors are sufficient.

df null model =  10  with the objective function =  12
df of  the model are 1  and the objective function was  5.6 

The root mean square of the residuals (RMSR) is  0.01 
The df corrected root mean square of the residuals is  0.04 

Fit based upon off diagonal values = 1

Si noti che, in questo caso, le unicità assumono valori negativi, il che suggerisce che la soluzione è impropria.

17.4. Metodo di massima verosimiglianza#

Il metodo di massima verosimiglianza è indicato quando si può assumere che le variabili manifeste seguano una distribuzione normale multivariata. In queste condizioni, il metodo produce stime dei pesi fattoriali che più verosimilmente hanno prodotto le correlazioni osservate. Gli stimatori di massima verosimiglianza sono preferibili a quelli ottenuti con altri metodi, a condizione che le premesse siano pienamente realizzate.

La funzione \(F\) che viene minimizzata dal metodo di massima verosimiglianza rappresenta una misura di “distanza” tra la matrice di covarianza osservata e quella predetta dal modello. Uguagliando a zero le derivate di \(F\) rispetto ai parametri del modello \(\boldsymbol{\Lambda}\) e \(\boldsymbol{\Psi}\) si ottengono le equazioni per le stime di massima verosimiglianza di \(\hat{\boldsymbol{\Lambda}}\) e \(\hat{\boldsymbol{\Psi}}\). Risolvendo tali equazioni rispetto alle incognite \(\hat{\boldsymbol{\Lambda}}\) e \(\hat{\boldsymbol{\Psi}}\) si ricavano le stime di massima verosimiglianza.

Poiché non esiste una soluzione analitica per le equazioni che stimano i parametri \(\boldsymbol{\Lambda}\) e \(\boldsymbol{\Psi}\), si utilizzano metodi numerici iterativi per minimizzare la differenza tra la matrice di covarianze (o correlazioni) osservata e quella predetta dal modello. Tuttavia, le stime di massima verosimiglianza possono presentare problemi di convergenza e casi di Heywood, in cui le stime di comunalità sono superiori a 1. Nonostante ciò, la soluzione è indipendente dall’unità di misura delle variabili manifeste e si ottiene la stessa soluzione sia analizzando la matrice delle varianze e covarianze sia quella delle correlazioni.

La stima di massima verosimiglianza si ottiene usando factanal. È inoltre il metodo di stima di default di lavaan.

Consideriamo nuovamente i dati dell’esempio precedente. La stima di massima verosimiglianza dei parametri \(\boldsymbol{\Lambda}\) e \(\boldsymbol{\Psi}\) si ottiene nel modo seguente. i

factanal(covmat=R, factors=2, rotation="none", n.obs=225)

Call:
factanal(factors = 2, covmat = R, n.obs = 225, rotation = "none")

Uniquenesses:
    K     I     H     L     J 
0.005 0.268 0.055 0.008 0.005 

Loadings:
  Factor1 Factor2
K  0.955  -0.289 
I  0.528   0.673 
H  0.720  -0.653 
L  0.954  -0.287 
J  0.764   0.642 

               Factor1 Factor2
SS loadings      3.203   1.457
Proportion Var   0.641   0.291
Cumulative Var   0.641   0.932

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 648.09 on 1 degree of freedom.
The p-value is 5.81e-143 

Si noti che il risultato è molto simile a quello ottenuto in precedenza. Si noti inoltre che le stime di massima verosimiglianza consentono un test di bontà di adattamento del modello ai dati (test chi quadrato).