Cousas de Multivariante

Esto es más una prueba de fuego que otra cosa. Básicamente voy a poner aquí unas cosas de PCA (análisis de componentes principales) que me gustan bastante. Esta técnica ya de por si me parece muy útil y fácil de visualizar, quizás justamente porque es una técnica de reducción de la dimensionalidad.

Así que ahí a lo bruto van mis scripts anotados de la clase de hoy:

###########################################
#   ANÁLISIS DE COMPONENTES PRINCIPALES   #
###########################################


# Ejemplo. Se ha examinado a 25 alumnos, aspirantes a ingresar en la Facultad de Mate-
# máticas, de 5 materias diferentes: Geometría Diferencial (cuyo resultado se almacena en
# la variable geodif), Análisis Complejo (ancompl), Álgebra (alg), Análisis Real (anreal) y Es-
# tadística (estad). Las puntuaciones obtenidas figuran en la tabla siguiente y se encuentran
# en el fichero aspirantes.txt.


datos<-read.table("aspirantes.txt",header=T)
head(datos)

##   geodif ancompl alg anreal estad
## 1     36      58  43     36    37
## 2     62      54  50     46    52
## 3     31      42  41     40    29
## 4     76      78  69     66    81
## 5     46      56  52     56    40
## 6     12      42  38     38    28

str(datos)

## 'data.frame':    25 obs. of  5 variables:
##  $ geodif : int  36 62 31 76 46 12 39 30 22 9 ...
##  $ ancompl: int  58 54 42 78 56 42 46 51 32 40 ...
##  $ alg    : int  43 50 41 69 52 38 51 54 43 47 ...
##  $ anreal : int  36 46 40 66 56 38 54 52 28 30 ...
##  $ estad  : int  37 52 29 81 40 28 41 32 22 24 ...

# Manualmente: matriz de covarianzas, luego autovectores

n<-nrow(datos);n # Observaciones

## [1] 25

d<-ncol(datos);d # Dimension en el espacio

## [1] 5

Sc<-cov(datos);Sc # matriz de covarianzas corregida

##           geodif   ancompl       alg   anreal    estad
## geodif  348.7733 181.69167 137.54500 176.8850 233.6583
## ancompl 181.6917 145.75000  91.28333 108.4750 142.5000
## alg     137.5450  91.28333  95.39333 106.1383 135.1833
## anreal  176.8850 108.47500 106.13833 166.9567 167.5500
## estad   233.6583 142.50000 135.18333 167.5500 272.6667

S<-Sc*(n-1)/n;S # Matriz de covarianza muestral

##           geodif ancompl      alg   anreal   estad
## geodif  334.8224 174.424 132.0432 169.8096 224.312
## ancompl 174.4240 139.920  87.6320 104.1360 136.800
## alg     132.0432  87.632  91.5776 101.8928 129.776
## anreal  169.8096 104.136 101.8928 160.2784 160.848
## estad   224.3120 136.800 129.7760 160.8480 261.760

aux <- eigen(S) # Autovalores y autovectores de S
lambda<-aux$values;lambda

## [1] 811.66184  81.63974  43.57261  37.62083  13.86339

v<-aux$vectors;v  # Se multiplicaría así ese componente por el dato de la tabla. 

##            [,1]       [,2]       [,3]        [,4]        [,5]
## [1,] -0.5982782  0.6745404  0.1852556 -0.38597894  0.06131111
## [2,] -0.3607532  0.2450733 -0.2490064  0.82871854 -0.24701742
## [3,] -0.3021774 -0.2140882 -0.2114109  0.13484564  0.89441442
## [4,] -0.3890403 -0.3384022 -0.6999921 -0.37537871 -0.32129949
## [5,] -0.5188995 -0.5697232  0.6074477  0.07178665 -0.17892129

# 5 variables en 5 notas dinstintas, si tengo que resumir las notas de un alumno e una sola lo que me dice es que lo mejor
# sería hacer una ponderación delas notas con estos coeficientes. De esta forma el que sea bueno bajará nota y el que la tenga
# mala la mejorará. Las PCS se calculan desde los datos centrados y por eso se tiende a la media, a "centrar" datos.

# Se reduce de 25 datos 5-Dimensionales a 25 datos unidimensionales. ¿Qué varianza tiene?

lambda

## [1] 811.66184  81.63974  43.57261  37.62083  13.86339

# La varianza coincide con el autovalor. La PCA 1º tendrá 811 de varianza, la segunda, de la segunda columna... 0'74 por la primera
#nota, 0'24 por al segunda, -0'21 por la tercera... 


# FORMA AUTOMÁTICA:

pca<-princomp(datos) # Desviaciones típicas de las PCA que sera lo mismo que
sqrt(lambda)

## [1] 28.489680  9.035471  6.600955  6.133582  3.723358

summary(pca)

## Importance of components:
##                            Comp.1     Comp.2     Comp.3     Comp.4     Comp.5
## Standard deviation     28.4896795 9.03547104 6.60095491 6.13358179 3.72335754
## Proportion of Variance  0.8212222 0.08260135 0.04408584 0.03806395 0.01402668
## Cumulative Proportion   0.8212222 0.90382353 0.94790936 0.98597332 1.00000000

lambda[1]/sum(lambda)

## [1] 0.8212222

lambda[2]/sum(lambda)

## [1] 0.08260135

pca$sdev^2 # Standard Desviation

##    Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
## 811.66184  81.63974  43.57261  37.62083  13.86339

screeplot(pca) # PCA's representados gráficamente

windows()
screeplot(pca,type="lines")

names(pca)

## [1] "sdev"     "loadings" "center"   "scale"    "n.obs"    "scores"   "call"

plot(pca$scores[,1],pca$scores[,2])
text(pca$scores[,1],pca$scores[,2],1:n)

# Coeficientes de los componentes o autovectores:

loadings(pca)

## 
## Loadings:
##         Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## geodif   0.598  0.675  0.185  0.386       
## ancompl  0.361  0.245 -0.249 -0.829 -0.247
## alg      0.302 -0.214 -0.211 -0.135  0.894
## anreal   0.389 -0.338 -0.700  0.375 -0.321
## estad    0.519 -0.570  0.607        -0.179
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings       1.0    1.0    1.0    1.0    1.0
## Proportion Var    0.2    0.2    0.2    0.2    0.2
## Cumulative Var    0.2    0.4    0.6    0.8    1.0

barplot(loadings(pca), beside = TRUE) # Representación coeficientes

pca$scores # Puntuaciones de los individuos en las componentes

##            Comp.1      Comp.2      Comp.3      Comp.4     Comp.5
##  [1,]  -7.5403215  10.2167650   2.5374713  -8.6708997 -4.3011164
##  [2,]  20.3610372  13.3460340   8.9820585   6.4124949 -1.3548711
##  [3,] -19.5031539   6.5552439  -1.6414327   5.0042015 -2.2980498
##  [4,]  65.9652730   1.3136646   5.1988159  -5.2093505 -1.0457668
##  [5,]   9.7780565   6.0680143  -9.1921582   2.9249303 -2.1069944
##  [6,] -33.0739529  -4.3722312  -3.7345190  -2.6038323 -5.3246839
##  [7,]   1.4212177  -0.7833317  -5.7800406   7.8225646 -0.4967347
##  [8,]  -6.7011638  -0.4667798 -13.3936497  -0.3040531  2.6525120
##  [9,] -36.1916160   5.6543609   2.9062814   5.5458471  6.5171961
## [10,] -38.0586178  -3.8266815  -2.5248244  -6.0338259  6.3212284
## [11,]  -1.6837296  -5.9262762  10.1237093  -4.9410716  4.5102357
## [12,]   6.4961788   3.9922267   5.0805842 -10.8805481 -1.3208797
## [13,]  48.6656614   2.5248899  -7.4227003  -6.1976282  3.0745497
## [14,]  12.8648106 -20.9375616   1.3328311  13.8459667  1.2286291
## [15,]  -5.1279619 -14.2990082  -2.9194088   2.7778854 -9.4299722
## [16,]  14.7627376   1.9542134   2.1019732   6.8802981  0.7262472
## [17,]   0.5206666   0.3328352  12.2272438   5.5083908  5.1001830
## [18,] -55.8465116   0.7024956   1.3349377  -4.9981377 -2.2873123
## [19,]   1.2809636  24.1913016   1.8618184   5.1875656 -3.1248129
## [20,]  44.0731110  -1.0133086  -8.2094067  -5.5695807  3.6879677
## [21,]   2.7282968 -15.4819784  13.1612877  -5.6870830 -3.1343155
## [22,] -22.3031608  -2.8413762  -5.9443642  -4.0929718  2.9545412
## [23,]  10.4526967  -7.6936163  -3.2010012  -3.0864389 -0.6469043
## [24,]  28.7147101  -2.0746379  -2.7048913   5.2002364 -0.7510021
## [25,] -42.0552279   2.8647426  -0.1806153   1.1650401  0.8501260

biplot(pca)

# Estos resultados se pueden obtener de forma manual:

n <- nrow(datos)
s <- cov(datos) * (n - 1)/n
auto <- eigen(s)
lambda = auto$values
plot(lambda, type = "l")

v <- auto$vectors

Ta guapo.