Coisas de Análise Multivariada

Isto é mais uma prova de fogo do que outra coisa. Basicamente, vou colocar aqui algumas coisas sobre PCA (análise de componentes principais) de que gosto bastante. Esta técnica, por si só, já me parece muito útil e fácil de visualizar, talvez precisamente por ser uma técnica de redução da dimensionalidade.

Por isso, aqui vão, sem grandes cerimónias, os meus scripts anotados da aula de hoje:

###########################################
#   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 ﬁguran en la tabla siguiente y se encuentran
# en el ﬁchero 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

É giro.

< Voltar ao início