Il modello fattoriale

Si consideri il modello fattoriale con \(k=1\) fattore e \(p=3\) variabili: \[\underset{p\times 1}{x} = \underset{p\times k}{\Lambda}\underset{k\times 1}{f} + \underset{p\times 1}{u}\]

Si generino \(n=20\) osservazioni indipendenti e identicamente distribuite da questo modello, salvando ciascuna osservazione \(x\) nella riga \(i\)-sima della matrice \(\underset{n\times p}{X}\). Per generare i dati, impostare il seme generatore dei numeri casuali set.seed(123) e specificare:

  • pesi fattoriali \(\Lambda = (\lambda_1,\lambda_2,\lambda_3)'=(0.5,0.2,-0.1)'\)
  • varianze specifiche \(diag(\Psi) = (\psi_1,\psi_2,\psi_3) = (0.0225, 0.0025, 0.008)\);
  • \(f \sim N(0,1)\)
  • \(u \sim N_3(\underset{3\times 1}{0},\Psi)\)
n <- 20
p <- 3
Lambda <- matrix(c(0.5,0.2,-0.1), ncol=1)
Psi <- diag(c(0.0225, 0.0025, 0.008))
X <- matrix(NA,nrow=n,ncol=p)
set.seed(123)
for (i in 1:n){
f <- rnorm(1)
u <- rnorm(p, mean=0, sd=sqrt(diag(Psi)) )
X[i,] <- Lambda %*% f + u 
}

La geometria del modello fattoriale può essere rappresentata graficamente come segue: i punti neri sono le osservazioni \(x\) generate come \(\Lambda f\) (punti rossi) e aggiungendo l’errore \(u\) (linee grigie). La retta rossa rappresenta il sottospazio dimensionale lungo il quale devono cadere tutti i valori di \(f\) mostrato in rosso.

  1. Calcolare la matrice di varianze/covarianze \(\Sigma = \Lambda \Lambda' + \Psi\) e confrontarla con la sua stima \(S\) ottenuta dalla matrice dei dati \(X\)
Sigma = Lambda %*% t(Lambda) + Psi
round(Sigma,2)
      [,1]  [,2]  [,3]
[1,]  0.27  0.10 -0.05
[2,]  0.10  0.04 -0.02
[3,] -0.05 -0.02  0.02
S = var(X) * ((n-1)/n)
round(S,2)
      [,1]  [,2]  [,3]
[1,]  0.18  0.06 -0.04
[2,]  0.06  0.03 -0.02
[3,] -0.04 -0.02  0.02
  1. Calcolare la matrice di correlazione \(\mathbb{C}\mathrm{orr}(x)=D^{-1/2}\Sigma D^{-1/2}\), con \(D = diag(\Sigma)=diag(\sigma_{11},\ldots,\sigma_{pp})\), e confrontarla con la sua stima \(R\) ottenuta dalla matrice dei dati \(X\);
D = diag(diag(Sigma)^(-1/2))
Corr_x = D %*% Sigma %*% D
round(Corr_x,2)
      [,1]  [,2]  [,3]
[1,]  1.00  0.93 -0.71
[2,]  0.93  1.00 -0.72
[3,] -0.71 -0.72  1.00
R = cor(X) 
round(R,2)
      [,1]  [,2]  [,3]
[1,]  1.00  0.90 -0.72
[2,]  0.90  1.00 -0.73
[3,] -0.72 -0.73  1.00
  1. Calcolare i pesi fattoriali \(\Lambda_z = D^{-1/2} \Lambda\) e le varianze specifiche \(\Psi_z = D^{-1/2} \Psi (D^{-1/2})'\) del modello fattoriale relativo alla standardizzazione \[\underset{p\times 1}{z} = D^{-1/2} \underset{p\times 1}{y} = \Lambda_zf + u_z\]
# pesi fattoriali z
Lambda_z = D %*% Lambda
round(Lambda_z,2)
      [,1]
[1,]  0.96
[2,]  0.97
[3,] -0.75
# varianze specifiche z
Psi_z = D %*% Psi %*% t(D)
diag(Psi_z)
[1] 0.08256881 0.05882353 0.44444444
# modello fattoriale per dati standardizzati
Lambda_z %*% t(Lambda_z) + Psi_z
           [,1]       [,2]       [,3]
[1,]  1.0000000  0.9292280 -0.7139216
[2,]  0.9292280  1.0000000 -0.7231015
[3,] -0.7139216 -0.7231015  1.0000000
Corr_x
           [,1]       [,2]       [,3]
[1,]  1.0000000  0.9292280 -0.7139216
[2,]  0.9292280  1.0000000 -0.7231015
[3,] -0.7139216 -0.7231015  1.0000000
  1. Stimare con il metodo della massima verosimiglianza il modello fattoriale con \(k=1\) fattore (funzione factanal(), argomenti factors=1 e rotation="none"), ed ottenere
  • la stima della matrice dei pesi fattoriali \(\hat{\Lambda}_z\);
  • la stima delle comunalità;
  • la stima delle varianze specifiche, i.e. gli elementi diagonali di \(\hat{\Psi}_z\);
af <- factanal(x=X, factors=1, rotation="none", scores = "regression") 
af

Call:
factanal(x = X, factors = 1, scores = "regression", rotation = "none")

Uniquenesses:
[1] 0.122 0.082 0.415

Loadings:
     Factor1
[1,]  0.937 
[2,]  0.958 
[3,] -0.765 

               Factor1
SS loadings      2.381
Proportion Var   0.794

The degrees of freedom for the model is 0 and the fit was 0 
# stima dei pesi fattoriali
hatLambda_z <- af$loadings[,]
round(hatLambda_z,3)
[1]  0.937  0.958 -0.765
# veri valori
round(Lambda_z,3)
       [,1]
[1,]  0.958
[2,]  0.970
[3,] -0.745
# stima delle comunalità
hatH2 <- apply(af$loadings^2,1,sum)
round(hatH2,3)
[1] 0.878 0.918 0.585
# veri valori
H2 <- apply(Lambda_z^2,1,sum)
round(H2,3)
[1] 0.917 0.941 0.556
# stima delle varianze specifiche
hatPsi_z <- diag(af$uniquenesses)
round(diag(hatPsi_z),3)
[1] 0.122 0.082 0.415
# veri valori
round(diag(Psi_z),3)
[1] 0.083 0.059 0.444
  1. Calcolare la differenza tra \(\mathbb{C}\mathrm{orr}(x)\) e \(\hat{\Lambda}_z \hat{\Lambda}_z' + \hat{\Psi}_z\)
# differenza vera correlazione e stima
round(Corr_x - ( hatLambda_z%*%t(hatLambda_z) + hatPsi_z ),5)
        [,1]    [,2]    [,3]
[1,] 0.00000 0.03147 0.00268
[2,] 0.03147 0.00000 0.00938
[3,] 0.00268 0.00938 0.00000
  1. Calcolare le stime \(\hat\Lambda = D^{1/2}\hat \Lambda_z\) e \(\hat\Psi = D^{1/2}\hat\Psi_z D^{1/2}\)
A = diag(diag(Sigma))
hatLambda = A^(1/2) %*% hatLambda_z
round(hatLambda, 4)
        [,1]
[1,]  0.4892
[2,]  0.1975
[3,] -0.1026
# vero valore
Lambda
     [,1]
[1,]  0.5
[2,]  0.2
[3,] -0.1
hatPsi <- A^(1/2) %*% hatPsi_z %*% A^(1/2)
round(hatPsi, 4)
       [,1]   [,2]   [,3]
[1,] 0.0332 0.0000 0.0000
[2,] 0.0000 0.0035 0.0000
[3,] 0.0000 0.0000 0.0075
# vero valore
Psi
       [,1]   [,2]  [,3]
[1,] 0.0225 0.0000 0.000
[2,] 0.0000 0.0025 0.000
[3,] 0.0000 0.0000 0.008
  1. Calcolare i punteggi fattoriali con il metodo di Thompson \[\hat f_i = \hat \Lambda' S^{-1} x_i\] e confrontarli con le vere realizzazioni di \(f\) tramite un diagramma di dispersione.
Xtilde = scale(X, center=TRUE, scale=FALSE)
hatf_thompson = apply(Xtilde,1, function(x) t(hatLambda) %*% solve(S) %*% x) 
plot(hatf_thompson, f)
abline(a=0,b=1)

  1. Calcolare i punteggi fattoriali con il metodo di Bartlett \[\hat f_i = (\hat \Lambda' \hat{\Psi}^{-1} \hat \Lambda) \hat \Lambda' \hat{\Psi}^{-1} x_i\] e confrontarli con le vere realizzazioni di \(f\) tramite un diagramma di dispersione.
hatf_bartlett = apply(Xtilde,1, function(x) solve(t(hatLambda) %*% solve(hatPsi) %*% hatLambda) %*% t(hatLambda) %*% solve(hatPsi)%*% x) 
# uguale a 
hatf = af$scores
sum(hatf_bartlett- hatf)
[1] -6.106227e-16
plot(hatf_bartlett, f)
abline(a=0,b=1, "dashed")

Dati Esami

Si consideri la seguente matrice di correlazione

rm(list=ls())
R = matrix(c(
 1.000, 0.439, 0.410, 0.288, 0.329, 0.248,
 0.439, 1.000, 0.351, 0.354, 0.320, 0.329,
 0.410, 0.351, 1.000, 0.164, 0.190, 0.181,
 0.288, 0.354, 0.164, 1.000, 0.595, 0.470,
 0.329, 0.320, 0.190, 0.595, 1.000, 0.464,
 0.248, 0.329, 0.181, 0.470, 0.464, 1.000
), byrow=T, ncol=6)
colnames(R) = c("Gaelic", "English", "History", "Arithmetic", "Algebra", "Geometry")
R
     Gaelic English History Arithmetic Algebra Geometry
[1,]  1.000   0.439   0.410      0.288   0.329    0.248
[2,]  0.439   1.000   0.351      0.354   0.320    0.329
[3,]  0.410   0.351   1.000      0.164   0.190    0.181
[4,]  0.288   0.354   0.164      1.000   0.595    0.470
[5,]  0.329   0.320   0.190      0.595   1.000    0.464
[6,]  0.248   0.329   0.181      0.470   0.464    1.000

relativa ai voti di \(n=220\) studenti maschi nelle materie Gaelic, English, History, Arithmetic, Algebra e Geometry (\(p=6\)). Interpretare la matrice di correlazione.

library(corrplot)
corrplot(R, diag=FALSE, type = 'upper')

I voti tra le materie sono correlati positivamente, ovvero gli studenti che vanno bene in una specifica materia vanno bene anche nelle altre. Le correlazioni più forti si trovano nel gruppo delle materie matematiche, e in misura leggermente inferiore, nel gruppo delle materie umanistiche, mentre c’è meno correlazione tra i due gruppi di materie.

  1. Stimare con il metodo della massima verosimiglianza il modello fattoriale con \(k=2\) fattori (funzione factanal(), argomenti factors=2 e rotation="none"), ottenendo la stima della matrice dei pesi fattoriali, interpretando i due fattori ottenuti;
n = 220
# modello fattoriale con 2 fattori senza rotazione
af <- factanal(covmat=R, factors=2, rotation="none", n.obs=n) 
af

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

Uniquenesses:
    Gaelic    English    History Arithmetic    Algebra   Geometry 
     0.510      0.594      0.644      0.377      0.431      0.628 

Loadings:
     Factor1 Factor2
[1,]  0.553   0.429 
[2,]  0.568   0.288 
[3,]  0.392   0.450 
[4,]  0.740  -0.273 
[5,]  0.724  -0.211 
[6,]  0.595  -0.132 

               Factor1 Factor2
SS loadings      2.209   0.606
Proportion Var   0.368   0.101
Cumulative Var   0.368   0.469

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 2.33 on 4 degrees of freedom.
The p-value is 0.674 
# stima dei pesi fattoriali
Lambda <- af$loadings[,]
round(Lambda,3)
     Factor1 Factor2
[1,]   0.553   0.429
[2,]   0.568   0.288
[3,]   0.392   0.450
[4,]   0.740  -0.273
[5,]   0.724  -0.211
[6,]   0.595  -0.132
# primo fattore "andare bene a scuola"
# secondo fattore "math vs non-math"
  1. Stimare il modello fattoriale con \(k=2\) fattori eseguendo la rotazione varimax della matrice dei pesi fattoriali (argomento rotation="varimax"). Confrontare le stime delle varianze specifiche e della matrice dei pesi fattoriali con la soluzione senza rotazione.
# rotazione varimax
af2 <- factanal(covmat = R, factors=2, rotation="varimax", n.obs=n)
af2

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

Uniquenesses:
    Gaelic    English    History Arithmetic    Algebra   Geometry 
     0.510      0.594      0.644      0.377      0.431      0.628 

Loadings:
     Factor1 Factor2
[1,] 0.235   0.659  
[2,] 0.323   0.549  
[3,]         0.590  
[4,] 0.771   0.170  
[5,] 0.724   0.213  
[6,] 0.572   0.210  

               Factor1 Factor2
SS loadings      1.612   1.203
Proportion Var   0.269   0.201
Cumulative Var   0.269   0.469

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 2.33 on 4 degrees of freedom.
The p-value is 0.674 
# matrice di rotazione
af2$rotmat
           [,1]      [,2]
[1,]  0.8420397 0.5394156
[2,] -0.5394156 0.8420397
# la matrice dei pesi fattoriali cambia
af2$loadings

Loadings:
     Factor1 Factor2
[1,] 0.235   0.659  
[2,] 0.323   0.549  
[3,]         0.590  
[4,] 0.771   0.170  
[5,] 0.724   0.213  
[6,] 0.572   0.210  

               Factor1 Factor2
SS loadings      1.612   1.203
Proportion Var   0.269   0.201
Cumulative Var   0.269   0.469
# primo fattore "math"
# secondo fattore "non-math"

# rappresentazione grafica della rotazione
plot(af$loadings,xlim=c(-1,1), ylim=c(-1,1), pch="", asp=1)
text(af$loadings, colnames(R))
abline(h=0)
abline(v=0)
af2 <- varimax(af$loadings,normalize=FALSE)
abline(0, af2$rotmat[2,1]/af2$rotmat[1,1], lty=2)
abline(0, af2$rotmat[2,2]/af2$rotmat[1,2], lty=2)

  1. Si riporti il \(p\)-value del test rapporto di verosimiglianza per verificare l’ipotesi nulla \(H_0:\) `\(k=1\) fattore è sufficiente’, e si concluda ad un livello di significatività del \(5\%\) se è opportuno utilizzare il modello fattoriale ad 1 fattore.
factanal(covmat=R, factors=1, n.obs = n)$PVAL
   objective 
4.528823e-08
  1. Si calcoli la statistica test rapporto di verosimiglianza \[T= n \log (|\hat{\Lambda} \hat{\Lambda}' + \hat{\Psi}|/|R|)\] e il corrispondente \(p\)-value.
k=1
p=6
af = factanal(covmat=R, factors=k, n.obs = n)
hatLambda = af$loadings[,]
hatPsi = diag(af$uniqueness)
fit = hatLambda %*% t(hatLambda) + hatPsi
t = n*log(det(fit)/det(R))
t
[1] 53.08106
gdl = ((p-k)^2 - p - k)/2
pchisq(t, lower.tail=FALSE, df=gdl)
[1] 2.821254e-08

Dati Open/Closed Book Examination Data

Caricare il dati scor presenti nella libreria bootstrap e calcolare la matrice di correlazione \(\underset{p\times p}{R}\).

library(bootstrap)

# carico i dati:
data(scor)
X <- scor
n <- nrow(X)
p <- ncol(X)

# matrice di correlazione
R <- cor(X)
round(R, 2)
     mec  vec  alg  ana  sta
mec 1.00 0.55 0.55 0.41 0.39
vec 0.55 1.00 0.61 0.49 0.44
alg 0.55 0.61 1.00 0.71 0.66
ana 0.41 0.49 0.71 1.00 0.61
sta 0.39 0.44 0.66 0.61 1.00
  1. Utilizzando la funzione factanal(), stimare con il metodo della massima verosimiglianza il modello fattoriale con \(k=2\) fattori (argomento factors=2) eseguendo la rotazione varimax della matrice dei pesi fattoriali (argomento rotation="varimax").
    Stimare i punteggi fattoriali con il metodo di Thompson e rappresentarli graficamente.
af <- factanal(X, factors=2, rotation = "varimax", scores="regression")
af$loadings

Loadings:
    Factor1 Factor2
mec 0.265   0.681  
vec 0.356   0.674  
alg 0.740   0.514  
ana 0.738   0.322  
sta 0.696   0.290  

               Factor1 Factor2
SS loadings      1.774   1.370
Proportion Var   0.355   0.274
Cumulative Var   0.355   0.629
# primo fattore "open book"
# secondo fattore "closed book"
puntfat <- af$scores
# punteggi fattoriali con il metodo di thompson
plot(puntfat,pch="")
text(puntfat, labels=c(1:88))
abline(h=0)
abline(v=0)

# studente 66
X[66,]
   mec vec alg ana sta
66  59  53  37  22  19

Stima del modello con il metodo dei fattori principali

Caricare il dati scor presenti nella libreria bootstrap e calcolare la matrice di correlazione \(\underset{p\times p}{R}\).

library(bootstrap)

# carico i dati:
data(scor)
X <- scor
n <- nrow(X)
p <- ncol(X)

# matrice di correlazione
R <- cor(X)
round(R, 2)
     mec  vec  alg  ana  sta
mec 1.00 0.55 0.55 0.41 0.39
vec 0.55 1.00 0.61 0.49 0.44
alg 0.55 0.61 1.00 0.71 0.66
ana 0.41 0.49 0.71 1.00 0.61
sta 0.39 0.44 0.66 0.61 1.00
  1. Calcolare la stima iniziale delle comunalità \(\hat{h}^2_i\) come \[\hat{h}^2_i = \max_{j\neq i}|r_{ij}|\] per \(i=1,\ldots,p\), dove \(r_{ij}\) è l’elemento di posizione \((i,j)\) della matrice \(R\). Ottenere la matrice di correlazione ridotta \[R^*= R - \hat{\Psi}\] dove \(\hat{\Psi} = diag(\hat{\psi}_1,\ldots,\hat{\psi}_p)\) con \(\hat{\psi}_i = 1 - \hat{h}^2_i\).
# sostituisco 1 sulla diagonale di R con 0
R0 <- R - diag(rep(1,p))
# calcolo la stima iniziale della comunalità
h2 <- apply(abs(R0), 2, max)
h2 
      mec       vec       alg       ana       sta 
0.5534052 0.6096447 0.7108059 0.7108059 0.6647357 
# calcolo la stima di Psi
Psi = diag(1-h2)
Psi
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.4465948 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.3903553 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.2891941 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.2891941 0.0000000
[5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.3352643
# calcolo la matrice di correlazione ridotta R* 
Rstar <- R - Psi 
round(Rstar,2)
     mec  vec  alg  ana  sta
mec 0.55 0.55 0.55 0.41 0.39
vec 0.55 0.61 0.61 0.49 0.44
alg 0.55 0.61 0.71 0.71 0.66
ana 0.41 0.49 0.71 0.71 0.61
sta 0.39 0.44 0.66 0.61 0.66
  1. Considerando il modello fattoriale con \(k=2\) fattori

  • calcolare la stima della matrice dei pesi fattoriali \[\hat{\Lambda} = V_k L^{1/2}_k\] dopo aver ottenuto la decomposizione spettrale di \(R^* = V L V'\) e dove \(\underset{p \times k}{V_k}\) contiene le prime \(k\) colonne di \(V\) e \(\underset{k \times k}{L_k} = diag(l_1,\ldots,l_k)\).

  • aggiornare la stima delle comunalità \[\hat{h}^2_i = \sum_{j=1}^{k} \hat{\lambda}^2_{ij}\] dove \(\hat{\lambda}_{ij}\) è l’elemento di posizione \((i,j)\) della matrice \(\hat{\Lambda}\) ottenuta al passo precedente

  • aggiornare la stima della matrice di correlazione ridotta \[R^*= R - \hat{\Psi}\] dove \(\hat{\Psi} = diag(\hat{\psi}_1,\ldots,\hat{\psi}_p)\) con \(\hat{\psi}_i = 1 - \hat{h}^2_i\) dove \(\hat{h}^2_i\) è la stima ottenuta al passo precedente

# decomposizione spettrale di R*
eigen <- eigen(Rstar)

#  k=2
k = 2

# stima di Lambda
Lambda <- eigen$vectors[,1:k] %*% diag(eigen$values[1:k]^{1/2})
Lambda
           [,1]        [,2]
[1,] -0.6459890  0.35354561
[2,] -0.7125528  0.30300833
[3,] -0.8639382 -0.05134392
[4,] -0.7864129 -0.24856797
[5,] -0.7419269 -0.27558100
# nuova stima comunalità 
h2.new = apply(Lambda^2, 1, sum)
h2.new
[1] 0.5422963 0.5995455 0.7490254 0.6802313 0.6264004
# nuova stima di Psi
Psi.new <- diag(1-h2.new)
Psi.new
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.4577037 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.4004545 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.2509746 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.3197687 0.0000000
[5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.3735996
# nuova stima di R*
Rstar.new = R - Psi.new
Rstar.new
          mec       vec       alg       ana       sta
mec 0.5422963 0.5534052 0.5467511 0.4093920 0.3890993
vec 0.5534052 0.5995455 0.6096447 0.4850813 0.4364487
alg 0.5467511 0.6096447 0.7490254 0.7108059 0.6647357
ana 0.4093920 0.4850813 0.7108059 0.6802313 0.6071743
sta 0.3890993 0.4364487 0.6647357 0.6071743 0.6264004
  1. Iterare per 100 volte la procedura descritta al punto 3, ottenendo le stime finali per \(\hat{\Lambda}\) e \(\hat{\Psi}\). Calcolare la differenza tra \(R\) e \(\hat{\Lambda} \hat{\Lambda}' + \hat{\Psi}\) e commentare.
for (i in 1:100){
  h2 <- apply(Lambda^2, 1, sum)
  Rstar <- R0 + diag(h2)
  eigen <- eigen(Rstar)
  Lambda <- eigen$vectors[,1:k] %*% diag(eigen$values[1:k]^{1/2})
}

# stima finale per Lambda
Lambda
           [,1]        [,2]
[1,] -0.6427483  0.34238878
[2,] -0.7081925  0.28687387
[3,] -0.8966056 -0.08718706
[4,] -0.7710248 -0.23504626
[5,] -0.7179803 -0.22818566
# stima finale per le comunalità
h2 <- apply(Lambda^2, 1, sum)
h2
[1] 0.5303554 0.5838332 0.8115032 0.6497260 0.5675644
# stima finale per Psi
Psi = diag(1-h2)

# differenza
fit = Lambda%*%t(Lambda) + Psi
round( R - fit, 4)
        mec     vec     alg     ana     sta
mec  0.0000  0.0000  0.0003 -0.0057  0.0057
vec  0.0000  0.0000 -0.0003  0.0065 -0.0066
alg  0.0003 -0.0003  0.0000 -0.0010  0.0011
ana -0.0057  0.0065 -0.0010  0.0000  0.0000
sta  0.0057 -0.0066  0.0011  0.0000  0.0000