Aporte R

1
Mario A. Morales R.

11 de diciembre de 2006

1 Profesor asistente, departamento de Matemáticas y Estadı́stica, Universidad

de Córdoba
2

Cambios
Este documento contiene, lo que le mandé en la primera oportunidad mas:

Los comandos para los rostros de Chernoff correspondiente al capı́tulo

1.

En el capı́tulo 2 se agrega otra forma de generar datos multinormales

usando la función mvrnorm() de la librerı́a MASS.

En el capı́tulo 2 se agrega otra forma de probar normalidad multiva-

riada usando la función mshapiro.test() de la librerı́a mvnormtest.

Capı́tulo 6 (análisis de factores comunes y únicos)

Capı́tulo 8 (análisis discriminante)

El apéndice: Se escriben todas las ordenes de la sección A.3 en código

de R

De esa forma solo quedan faltando los dos últimos capı́tulos

Acerca del recuadro que rodea el código

Si no le gusta el recuadro que rodea el código, es posible quitarlo totalmente
o usar otro estilo. El recuadro se debe a la opción frame=single del entorno
Verbatim, puede borrarla para que no le salga ningún recuadro o puede usar
una de las siguientes:

frame=bottomline

frame=lines

frame=topline

frame=leftline

frame=none (equivalente a quitarlo).

A continuación unos ejemplos:

frame=none
 3

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

frame=single

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

frame=lines

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

frame=bottomline

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
4

stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

frame=leftline

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

El entorno Verbatim lo proporciona el paquete fancyvrb ası́ que en el preámbu-

lo de su documento debe colocar la lı́nea

\usepackage{fancyvrb}

le adjunto el archivo fancyvrb.sty, de tal forma que si usted no lo tiene

instalado, sólo tiene que colocarlo en el mismo directorio donde tiene los
archivos fuentes .tex
Capı́tulo 1

Conceptos preliminares

Hola profe, para el capı́tulo 1 le propongo los siguientes códigos de R que

realizan la mayorı́a de los gráficos del capı́tulo y las estadı́sticas de resumen,
como el vector de medias, la matriz de covarianzas, la de correlaciones, la
varianza total, la varianza generalizada y por último muestra como se calcula
la tabla 1.3 que muestra la distancia de Mahalanobis y la euclidiana.
(OJO Revisar la tabla 1.3 pag. 32 La segunda columna dice Distancia de
Mahalanobis y según los cálculos que hago con R esa columna contiene es el
cuadrado de las distancias.)

Observación

La lectura de los datos se hace con la función scan(), eso tiene el inconve-
niente que ocupa mucho espacio, otra alternativa es usar read.table() pero
tiene el inconveniente que los datos estarı́an en un archivo plano aparte y
el lector no tendrı́a acceso a ellos a menos que se coloquen en una página
web del libro y tiene el inconveniente que seria mas difı́cil para el lector no
familiarizado con R poner a correr el código, por lo del directorio de trabajo
y esas cosas más avanzadas de lectura de datos externos.
Como las estadı́sticas de resumen que están en el libro las calcula con la tabla
2.1, con estos datos hice lo mismo que con los de la tabla 1.1. Es decir hago
los gráficos y calculo las estadı́sticas.

5
6 CAPÍTULO 1. CONCEPTOS PRELIMINARES

A lo que vinimos
El siguiente código de R lee los datos de la tabla 1.1 y con ellos realiza los
diagramas de tallos y hojas que están inmediatamente debajo, el disperso-
grama de la figura 1.4, el box plot de la figura 1.5, los perfiles de la figura
1.3 y por último un gráfico de estrellas y los rostros de Chernoff.

# Lectura de los datos de la tabla 1.1

datos<-scan()
125 2536 28 86 2505 31 119 2652 32 113 2573 20
101 2382 30 143 2443 30 132 2617 27 106 2556 36
121 2489 34 109 2415 29 88 2434 27 116 2491
24 102 2345 26 75 2350 23 90 2536 24 109 2577
22 104 2464 35 110 2571 24 96 2550 24 101 2437
23 95 2472 36 117 2580 21 115 2436 39 138
2200 41 85 2851 17

datos2<-matrix(datos,ncol=3,byrow=TRUE)
tabla1_1<-data.frame(datos2)
colnames(tabla1_1)<-c("CI","Peso","Edad")

# termina la lectura de datos

Los gráficos que se crean con este código corresponden a los de la sección 1.2

# Diagramas de tallos y hojas

stem(tabla1_1$CI)
stem(tabla1_1$Peso)
stem(tabla1_1$Edad)

# Dispersograma figura 1.4

pairs(tabla1_1)

# En la siguiente linea de código, la función scale()

# estandariza los datos. La función stack()
# convierte la tabla1_1 a un data frame con dos columnas
# la primera contiene los valores y la segunda es un
 7

# factor que identifica a que variable corresponde el valor

datos3<-stack(data.frame(scale(tabla1_1)) )

# Boxplot, figura 1.5 con los datos de la tabla 1.1

plot(datos3$ind,datos3$values)

# El siguiente código crea los perfiles de la matriz de datos

# (figura 1.3) pero adicionalmente dibuja sobre cada
# histograma la gráfica de la densidad normal. Requiere
# la librerı́a lattice.

library(lattice)
histogram(~values|ind,data=datos3, layout = c(3,1),
type = "density", panel = function(x, ...)
{panel.histogram(x, ...)
panel.mathdensity(dmath = dnorm, col = "black",
args = list(mean=mean(x),sd=sd(x)))
} )

# Gráfico de estrellas
stars(tabla1_1)

# Rostros de Chernoff, requiere la librerı́a aplpack

library(aplpack)
faces(tabla1_1)

Calculo de las estadı́sticas de resumen

# Estadı́sticas de resumen

summary(tabla1_1)

# vector de medias
8 CAPÍTULO 1. CONCEPTOS PRELIMINARES

mean(tabla1_1)

# Matriz de varianzas y covarianzas

cov(tabla1_1)

# Matriz de correlaciones

cor(tabla1_1)

# VT (traza de la matriz de covarianzas)

sum(diag(cov(tabla1_1)))

# VG (determinante de la matriz de covarianzas)

det(cov(tabla1_1))
 9

Lo mismo anterior pero con los datos de la tabla 1.2. Si de mi dependiera

esto serı́a lo que quedarı́a en el libro.

datos<-scan()
1.38 51 4.8 115
1.40 60 5.6 130
1.42 69 5.8 138
1.54 73 6.5 148
1.30 56 5.3 122
1.55 75 7.0 152
1.50 80 8.1 160
1.60 76 7.8 155
1.41 58 5.9 135
1.34 70 6.1 140

datos2<-matrix(datos,ncol=4,byrow=TRUE)
tabla1_2<-data.frame(datos2)
colnames(tabla1_2)<-c("X1","X2","X3","X4")
# termina la lectura de datos

# Diagramas de tallos y hojas

stem(tabla1_2$X1)
stem(tabla1_2$X2)
stem(tabla1_2$X3)
stem(tabla1_2$X3)

# Dispersograma figura 1.4 (con los datos de la tabla 1.2)

pairs(tabla1_2)

# En el siguiente código, la función scale() estandariza

# La función stack() convierte la tabla1_2 a un data frame
# con dos columnas la primera contiene los valores de
# la variable, y la segunda es un factor que identifica
# a qué variable corresponde el valor

datos3<-stack(data.frame(scale(tabla1_2)) )

# Boxplot, figura 1.5 con los datos de la tabla 1.2

10 CAPÍTULO 1. CONCEPTOS PRELIMINARES

plot(datos3$ind,datos3$values)

# El siguiente código crea los perfiles de la matriz de datos

# (figura 1.3) pero adicionalmente dibuja sobre cada
# histograma la gráfica de la densidad normal. Requiere
# la librerı́a lattice.

library(lattice)
histogram(~values|ind,data=datos3, layout = c(3,1),
type = "density", panel = function(x, ...)
{panel.histogram(x, ...)
panel.mathdensity(dmath = dnorm, col = "black",
args = list(mean=mean(x),sd=sd(x)))
} )

# Gráfico de estrellas
stars(tabla1_2)

# Rostros de Chernoff, requiere la librerı́a aplpack

library(aplpack)
faces(tabla1_2)

# Estadı́sticas de resumen

summary(tabla1_2)

# vector de medias

round(mean(tabla1_2),2)

# Matriz de varianzas y covarianzas

round(cov(tabla1_2),3)

# Matriz de correlaciones

round(cor(tabla1_2),3)
 11

# VT (traza de la matriz de covarianzas)

sum(diag(cov(tabla1_2)))

# VG (determinante de la matriz de covarianzas)

det(cov(tabla1_2))

# Cuadrado de la distancia de mahalanobis entre cada

# observación y el vector de medias

S<-cov(tabla1_2)
mahalanobis(tabla1_2,center=mean(tabla1_2),cov=S)

# Distancia de mahalanobis entre cada

# observación y el vector de medias

sqrt( mahalanobis(tabla1_2,center=mean(tabla1_2),cov=S) )

# la distancia euclidiana se obtiene tomando la matriz de

# covarianza igual a la identidad la cual se obtiene con
# el comando diag(1,n) (identidad de orden n)

I4<-diag(1,4)
sqrt(mahalanobis(tabla1_2,center=mean(tabla1_2),cov=I4) )
12 CAPÍTULO 1. CONCEPTOS PRELIMINARES
Capı́tulo 2

Distribuciones multivariantes

Para la generación de muestras a partir de una distribución normal multi-

variante se traduce el código de SAS que usted tiene en el libro a código de
R pero se hace de una manera mucho más interesante; como una función.
Es usuario invoca a la función entregándole n , µ , Σ y opcionalmente una
semilla y la función regresa la matriz de datos generados con esa media y
covarianza.

gmultinorm<-function(n,sigma=diag(1,n),mu=matrix(0,nrow=n),
semilla=NULL){
if(!is.null(semilla)){set.seed(semilla)}
# vector de medias
p<-nrow(sigma) # numero de variables
unos<-matrix(1,nrow=n) #
# repite el vector mu n veces por fila y una vez por columnas
M<-t(mu)%x%unos
G<-chol(sigma) # descomposición de cholesky
Z<-matrix(rnorm(p*n),ncol=p) # genera n vectores $N_p(0,Ip)
# genera n vectores $N_p(mu ,sigma)
Y<-Z%*%G+M
Y
}

# el llamado a la función
n<-20 # muestra de tama~ no n=20
# vector de medias

13
14 CAPÍTULO 2. DISTRIBUCIONES MULTIVARIANTES

media<-matrix(c(1,3,0))
# mat de covarianzas
S<-matrix(c(4,2,1,2,3,1,1,1,5),nrow=3,byrow=TRUE)

gmultinorm(20,mu=media,sigma=S,semilla=552154123)

# datos de una normal con media cero y covarianza identidad

gmultinorm(20,mu=matrix(0,nrow=3),sigma=diag(1,3))

Agregado en la ultima versión (Diciembre 11 de 2006 )

Alternativamente, para la generación de datos multinomiales podemos usar
la función mvrnorm() de la librerı́a MASS, de la siguiente forma

library(MASS)
mvrnorm(n, media, S)

# Usando una semilla

set.seed(552154123)
mvrnorm(20, media, S)

Tenga es cuenta que esta función hace la descomposición de la matriz S, via

eigen mientras que la función creada arriba lo hace mediante cholesky.
Para la prueba de multinormalidad de Mardia se traduce el código SAS a
código de R pero como una función, de tal forma que el usuario entrega una
matriz de datos y la función regresa las estadı́sticas

mardia.test<-function(Y){
n<-nrow(Y) # numero de filas de Y
p<-ncol(Y) # numero de columnas de Y
gl_chi<-p*(p+1)*(p+2)/6 # grados de libertad
Q<-diag(n)-(1/n)*matrix(1,n,n) # I_p-(1/n)1_n1’_n
S<-(1/n)*t(Y)%*%Q%*%Y # matriz de covarianzas muestral
# Matriz g_hi de la ecuación 2.12
G_MATRIZ<- Q%*%Y%*%solve(S)%*%t(Y)%*%Q
b_1<-sum(G_MATRIZ^3)/(n^2) # cálculo de la simetı́a
b_2<-sum(diag(G_MATRIZ^2))/n # calculo de la curtosis b_(2,p)
EST_b_1<-n*b_1/6 # calculo de la estadı́stica B1 ec. (2.13a)
 15

# calculo de la estadı́stica B1 ec. (2.13a)

EST_b_2<-(b_2-p*(p+2))/sqrt(8*p*(p+2)/n)
PVAL_ses<-1-pchisq(EST_b_1,gl_chi)
PVAL_cur<-2*(1-pnorm(abs(EST_b_2)))
cat("b_1=",b_1,"b_2=",b_2,"EST_b_1=",EST_b_1,
"EST_b_2=",EST_b_2,"\n")
cat("PVAL_ses=",PVAL_ses,"PVAL_cur=",PVAL_cur,"\n")
}

El llamado a la función considerando los datos de la tabla 3.4

Y1<-scan()
72 66 76 77 60 53 66 63 56 57 64 58 41 29 36 38 32 32 35 36
30 35 34 26 39 39 31 27 42 43 31 25 37 40 31 25 33 29 27 36
32 30 34 28 63 45 74 63 54 46 60 52 47 51 52 43 91 79 100
75 56 68 47 50 79 65 70 61 81 80 68 58 78 55 67 60 46 38
37 38 39 35 34 37 32 30 30 32 60 50 67 54 35 37 48 39 39
36 39 31 50 34 37 40 43 37 39 50 48 54 57 43

Y<-matrix(Y1,ncol=4,byrow=TRUE)

mardia.test(Y)

Agregado en la ultima versión (Diciembre 11 de 2006 )

Test de normalidad multivariada usando la prueba de Shapiro-Wilk

#prueba de Shapiro-Wilk multivariada, requiere la

#librerı́a mvnormtest

library(mvnormtest)
mshapiro.test(t(Y))
16 CAPÍTULO 2. DISTRIBUCIONES MULTIVARIANTES
Capı́tulo 3

Inferencia sobre el vector de

medias

Se traduce el código SAS de la sección 3.6

# lectura de datos ejemplo 3.5

datos<-scan()
1 15 17 24 14 1 17 15 32 26 1 15 14 29 23
1 13 12 10 16 1 20 17 26 28 1 15 21 26 21
1 15 13 26 22 1 13 5 22 22 1 14 7 30 17
1 17 15 30 27 1 17 17 26 20 1 17 20 28 24
1 15 15 29 24 1 18 19 32 28 1 18 18 31 27
1 15 14 26 21 1 18 17 33 26 1 10 14 19 17
1 18 21 30 29 1 18 21 34 26 1 13 17 30 24
1 16 16 16 16 1 11 15 25 23 1 16 13 26 16
1 16 13 23 21 1 18 18 34 24 1 16 15 28 27
1 15 16 29 24 1 18 19 32 23 1 18 16 33 23
1 17 20 21 21 1 19 19 30 28 2 13 14 12 21
2 14 12 14 26 2 12 19 21 21 2 12 13 10 16
2 11 20 16 16 2 12 9 14 18 2 10 13 18 24
2 10 8 13 23 2 12 20 19 23 2 11 10 11 27
2 12 18 25 25 2 14 18 13 26 2 14 10 25 28
2 13 16 8 14 2 14 8 13 25 2 13 16 23 28
2 16 21 26 26 2 14 17 14 14 2 16 16 15 23
2 13 16 23 24 2 2 6 16 21 2 14 16 22 26
2 17 17 22 28 2 16 13 16 14 2 15 14 20 26

17
18 CAPÍTULO 3. INFERENCIA SOBRE EL VECTOR DE MEDIAS

2 12 10 12 9 2 14 17 24 23 2 13 15 18 20
2 11 16 18 28 2 7 7 19 18 2 12 15 7 28
2 6 5 6 13

datos2<-matrix(datos,ncol=5,byrow=TRUE)
ejemp3_5<-data.frame(datos2)
colnames(ejemp3_5)<-c("sexo","X1","X2","X3","X4")

# toma los datos para hombres

hombres<-subset(ejemp3_5,ejemp3_5$sexo==1,select=c(2:5))
# toma los datos para mujeres
mujeres<-subset(ejemp3_5,ejemp3_5$sexo==2,select=c(2:5))
# numero de variables
p<-ncol(hombres)
# numero de hombres
n1<-nrow(hombres)
# numero de mujeres
n2<-nrow(mujeres)
# grados de libertad de Sp
v<-n1+n2-2
# vector de medias de hombres
XMH<-mean(hombres)
# vector de medias de mujeres
XMM<-mean(mujeres)
# Matriz de covarianzas de hombres
SH<-cov(hombres)
# Matriz de covarianzas de mujeres
SM<-cov(mujeres)
# Matris de covarianzas ponderadas
Sp<-1/v*( (n1-1)*SH+(n2-1)*SM )
# estadı́stica T^2
T2<-(n1*n2/(n1+n2))*mahalanobis(XMH, XMM, Sp)
# transformación a la estadı́stica F
F0<-(v-p+1)/(v*p)*T2
# p valor asociado a F0
p_val<-pf(F0,p,n1+n2-p-1,lower.tail=FALSE)
# imprime resultados
 19

cat("\n","T2= ",T2, " P_valor= ",p_val, "\n")

Análisis de varianza multivariado, corresponde a la sección 3.7 con los mismos

datos y genera las salidas del ejemplo 3.9

# lectura de datos elemplo 3.5

datos<-scan()
1 69 75 1 69 70 1 71 73 1 78 82 1 79 81 1 73
75 2 69 70 2 68 74 2 75 80 2 78 85 2 68 68 2
63 68 2 72 74 2 63 66 2 71 76 2 72 78 2 71
73 2 70 73 2 56 59 2 77 83 3 72 79 3 64 65
3 74 74 3 72 75 3 82 84 3 69 68 3 76 76 3
68 65 3 78 79 3 70 71 3 60 61

datos2<-matrix(datos,ncol=3,byrow=TRUE)
ejemp3_9<-data.frame(datos2)
colnames(ejemp3_9)<-c("metodo","matemat","escrit")

# se ubican las columnas y1, y2, en una matriz llamada Mdatos

Mdatos<-as.matrix(ejemp3_9[,-1])
# se define el factor y se llama metodo
ejemp3_9$metodo<- as.factor(ejemp3_9$metodo)

# Análisis de varianza univariado

# para matemat
ajusteM<-lm(matemat~metodo,data=ejemp3_9)
anova(ajusteM)

# Análisis de varianza univariado

# para escritura
ajusteE<-lm(escrit~metodo,data=ejemp3_9)
anova(ajusteE)

# Ajustamos el modelo multivariado de una via

ajuste<-manova(Mdatos~metodo )

# Las diferentes estadisticas

summary(ajuste ,test="Wilks")
20 CAPÍTULO 3. INFERENCIA SOBRE EL VECTOR DE MEDIAS

summary(ajuste,test="Pillai")
summary(ajuste ,test= "Hotelling-Lawley")
summary(ajuste ,test= "Roy")

# Las matrices E y H
M<-summary(ajuste)$SS
H<-M$metodo
E<-M$Residuals
Capı́tulo 4

Inferencia sobre la matriz de

covarianzas

Se implementó el programa de la sección 4.4 en la función sigma.test().

El usuario entrega las matrices Σ0 , S y el valor de n y la función regresa el
valor de la estadı́stica λ∗1 junto con el p−valor

sigma.test<-function(Sigma_0,S,n){
# numero de filas la matriz Sigma
p<-nrow(S)
gl<-(1/2)*p*(p+1) # grados de libertad
# producto entre S y la inversa de Sigma_0
Sisigma_0<-S%*%solve(Sigma_0)
#determinante de sigma_0
D_sigma_0<-det(Sigma_0)
#determinante de S
D_S<-det(S)
# estadı́stica lambda de Ecua 4.6
E_lambda<-(n-1)*(log(D_sigma_0)-log(D_S)+
sum(diag(Sisigma_0))-p)
# estadı́stica lambda de Ecua 4.7
E_lambda1<-(1-(1/(6*(n-1)))*(2*p+1-2/(p+1) ))* E_lambda
p_val<-pchisq(E_lambda1,gl) # pvalor
list(E_lambda1=E_lambda1, P_valor=p_val)
}

21
22CAPÍTULO 4. INFERENCIA SOBRE LA MATRIZ DE COVARIANZAS

# llamado de la función

Sigma_0<-matrix(c(4,3,2,3,6,5,2,5,10),nrow=3 )
S<-matrix(c(3.42,2.60,1.89,2.60,8,6.51,1.89,6.51,9.62),nrow=3)
sigma.test(Sigma_0,S,n=20)

La siguiente función implementa la prueba de igualdad de matrices de va-

rianza covarianza para dos o mas poblaciones, está basada en el estadı́stico
λ1 que usted muestra en la ecuación 4.11, pero se usan unas estadı́sticas que
muestra Rencher en la sección 7.3.2 que permiten una mejor aproximación
de la distribución de las estadı́sticas. Con mas tiempito se puede adaptar
completamente para que use solo las ecuaciones del libro. Tiene la ventaja
que el usuario entrega la matriz de datos y un factor que clasifica los datos
y la función se encarga de calcular todo lo demás y regresa las estadı́sticas
en mención.

var.igual<-function(datos,grupos){
# matriz de ceros para guardar la suma de v_i*S_i
St<-matrix(0,nrow=ncol(datos),ncol=ncol(datos))
# matrices de covarianzas.
Covs<-by(datos,grupos,cov)
# numero de obs en cada grupo
ni<-as.vector( by(datos,grupos,nrow ) )
vi<-ni-1
lev<-levels(grupos) # niveles del factor
# un vector para guardar los logaritmos de los determinantes.
LnSi<-numeric(1)
# este ciclo for calcula los logaritmos de los determinantes
# y la matriz para calcular Sp
for(i in 1:length(lev)){
Si<- Covs[[i]]
LnSi[i]<-log(det(Si))
St<-St+vi[i]*Si
}
# matriz Sp
Sp<- St/sum(vi)
# se implementa la ec 7.23 Rencher
LnM<-0.5*sum(vi* LnSi) - 0.5*sum(vi)*log(det(Sp))
 23

p<-ncol(Sp)
k<-length(ni) # numero de grupos
# ec. 7.22 Rencher
c1<-(sum(1/vi)-1/sum(vi))*( (2*p^2+3*p-1)/(6*(p+1)*(k-1)) )
# ec. 7.24 Rencher
c2<-(sum(1/vi^2)-1/sum(vi)^2)*( (p-1)*(p+2)/( 6*(k-1) ) )
a1<-0.5*(k-1)*p*(p+1)
a2<-(a1+2)/(abs(c2-c1^2))
b1<-(1-c1-a1/a2)/a1
b2<- (1-c1-2/a2)/a2
# ec. 7.21
u<- -2*(1-c1)*LnM
if(c2>=c1^2){
F<- -2*b1*LnM # ec. 7.25 Rencher
pvalor<-pf(F,a1,a2,lower.tail=FALSE)
}else{
F<- -( a2*b2*LnM ) /(a1*(1+2*b2*LnM)) # ec. 7.26 Rencher
pvalor<-pf(F,a1,a2,lower.tail=FALSE)
}
Resp<-matrix(c(LnM,u,c1,c2,a1,a2,b1,b2,F,pvalor,-2*LnM),
nrow=1)
colnames(Resp)<-c("LnM","u","c1","c2","a1","a2","b1","b2","F",
"pvalor","-2*LnM")
print(Resp)
}

Para probar la hipótesis que las matrices de varianzas y covarianzas de los

tratamientos del ejemplo 3.9 son iguales, usamos el comando: (después de
haber definido la función)

var.igual(Mdatos,ejemp3_9[,1] )
24CAPÍTULO 4. INFERENCIA SOBRE LA MATRIZ DE COVARIANZAS
Capı́tulo 5

Análisis de componentes
principales

Con el siguiente código se realizan los cálculos del ejemplo 5.1

# lectura de datos datos<-scan()

156.0 245.0 31.6 18.5 20.5
154.0 240.0 30.4 17.9 19.6

insertar aquı́ el resto de datos

162.0 245.0 32.5 18.5 21.1

164.0 248.0 32.3 18.8 20.9

datos2<-matrix(datos,ncol=5,byrow=TRUE)
ejemp5_1<-data.frame(datos2)

# matriz de varianza covarianza

cov(ejemp5_1)
# vector de medias
mean(ejemp5_1)
# desviaciones estándar
sqrt(diag(cov(ejemp5_1) ))
# matriz de coreelaciones
cor(ejemp5_1)

25
26 CAPÍTULO 5. ANÁLISIS DE COMPONENTES PRINCIPALES

# análisis de componentes principales desde la matriz

# de correlaciones
acp<-princomp(ejemp5_1,cor=TRUE)
summary(acp)
# gráfico scree
plot(acp)
# la desviación estándar de cada componente principal
# es decir la raiz de los valores propios de la matriz
acp$sdev
# Matriz con los vectores propios
acp$loadings
# la media de las variables originales con la que se corrigen
# las obs
acp$center
# numero de observaciones
acp$n.obs
# las coordenadas factoriales
acp$scores

#biplots
par(mfrow=c(2,2))
# primer plano factorial
biplot(acp)
#segundo plano factorial
biplot(acp,choices = c(1,3))
# tercer plano factorial
biplot(acp,choices = c(2,3))

# Acp desde la matriz de covarianzas

# ( opción por defecto )
acpCov<-princomp(ejemp5_1)
Capı́tulo 6

Análisis de factores comunes y

únicos

Se realiza el ejemplo 6.1 usando el método de la componente principal y el

método de máxima verosimilitud.

# Matriz R del ejemplo 6.1

R<-matrix(c(1,.02,.96,.42,.01,
.02,1,.13,.71,.85,
.96,.13,1,.5,.11,
.42,.71,.50,1,.79,
.01,.85,.11,.79,1 ),nrow=5)

# Valores y vectores propios

eig<-eigen(R)
# Matriz P (vectores propios)
P<-eig$vectors
P
# valores propios
eig$values
# ponderaciones factoriales
f<-P[,1:2]%*%D
f
# Comunalidad
comun<-matrix(rowSums(f^2))
comun

27
28 CAPÍTULO 6. ANÁLISIS DE FACTORES COMUNES Y ÚNICOS

# varianza especı́fica
1-comun

# Rotación varimax
varimax(f)

# Para realizar la rotación cuartimax se requiere instalar la

# librerı́a GPArotation

library(GPArotation)
TV <- GPFoblq(f, method="quartimax", normalize=TRUE)
print(TV)
summary(TV)

#La libreria GPArotation realiza varios tipos de rotación

#para otros métodos pida ayuda mediante el comando
#help(’GPFoblq’)

#Análisis de factores usando lo función factanal()

#tenga en cuenta que esta función realiza el análisis de
#factores usando el método de maxima verosimilitud.

fac<-factanal(factors=2, covmat=R,rotation="none")
# los pesos se obtienen ası́:
fac$loadings
# Varianza especı́fica
fac$uniquenesses
# Comunalidades
1-fac$uniquenesses
# Matriz de correlaciones usada
fac$correlation
# Análisis de factores con rotación varimax
fac<-factanal(factors=2, covmat=R,rotation="varimax")
# Como "varimax" es la opción por defecto para rotation,
# el siguiente comando produce el mismo resultado.
fac<-factanal(factors=2, covmat=R)
Capı́tulo 7

Análisis de conglomerados

Si ilustra la conformación de conglomerados con R usando la matriz de dis-

tancias que aparece en la página 272. Se realiza el análisis usando varios
métodos.

x<-matrix(c(0,3,7,11,10,
3,0,6,10,9,
7,6,0,5,6,
11,10,5,0,4,
10,9,6,4,0),ncol=5 )

dimnames(x)<-list(paste("O",1:5,sep=""),
paste("O",1:5,sep="") )
y<-as.dist(x)

# Enlace completo (la vecina más lejana): opción por defecto

cl<-hclust(y)
plot(cl,hang = -1)
abline(h=4.5,lty=2)

# Enlace simple (la vecina más cercana)

cl<-hclust(y, method = "single")
plot(cl,hang = -1)
abline(h=4.5,lty=2)

# Union mediante el promedio

29
30 CAPÍTULO 7. ANÁLISIS DE CONGLOMERADOS

cl<-hclust(y, method="average")
plot(cl,hang = -1)
abline(h=4.5,lty=2)

# método de ward
cl<-hclust(y, method="ward")
plot(cl,hang = -1)
abline(h=4.5,lty=2)

Si se tienen los datos, en lugar de la matriz de distancias, se puede calcular

esta mediate la función dist(), a continuación se ilustra como hacerlo.

x<-c(1,2,4,5,3,3)
y<-c(2,1,1,4,5,3)
datos<-data.frame(x=x,y=y,row.names=LETTERS[1:length(x)])
# distancia máxima entre dos componentes de X y Y
dd<-dist(datos,method="maximum")
# distancia de manhattan
dd<-dist(datos,method="manhattan")
# distancia de canberra
dd<-dist(datos,method="canberra")
# distancia de minkowski
dd<-dist(datos,method="minkowski")
# distancia euclidiana
dd<-dist(datos)
# distancia euclidiana
dd<-dist(datos,method="euclidean")

# conformación de los conglomerados

cl<-hclust(dd)
plot(cl,hang = -1)
Capı́tulo 8

Análisis discriminate

Se presenta el ejemplo 8.2, como son muchos datos (172 filas) la lectura se
hace desde un archivo externo usando la función (read.table())

# lectura de los datos

ejemp8_2<-read.table("ejemplo8_2.txt",header=TRUE)
# transformación mediante raı́z cuadrada
ejemp8_2$DOWN_A<-sqrt(ejemp8_2$DOWN_A)
ejemp8_2$RIGHT_A<-sqrt(ejemp8_2$RIGHT_A)
# definición del factor
ejemp8_2$GRP<-factor(ejemp8_2$GRP)
head(ejemp8_2)
# requiere la librerı́a MASS
library(MASS)
# análisis discriminante lineal
z<-lda(GRP ~.,ejemp8_2,prior =c(36/172,36/172,50/172,50/172))

# clasificación de las observaciones por medio de la regla

# obtenida.
clasif<-predict(z,ejemp8_2[,-9])$class
clasif

# tabla de clasificación
addmargins(table(ejemp8_2$GRP,clasif))

Suponga que se tiene la observación futura

31
32 CAPÍTULO 8. ANÁLISIS DISCRIMINATE

DOWN A DOWN P DOWN L DOWN B RIGHT A RIGHT P RIGHT L RIGHT B

54.36 220.07 90.08 44.84 53.86 220.17 90.27 43.53

y queremos aplicar la regla de discriminación sobre ella, en R procedemos ası́:

nuevo<-data.frame(DOWN_A=54.36,DOWN_P=220.07,DOWN_L=90.08,
DOWN_B=44.84,RIGHT_A=53.86,RIGHT_P=220.17,
RIGHT_L=90.27,RIGHT_B=43.53)

predict(z,nuevo)$class

Clasificación usando la función de discriminación cuadrática

zq<-qda(GRP ~.,ejemp8_2,prior =c(36/172,36/172,50/172,50/172))

clasifq<-predict(zq,ejemp8_2[,-9])$class
clasifq
tabla<-table(ejemp8_2$GRP,clasif)
addmargins(tabla)

Profe no se (hasta el momento) como obtener la tabla de clasificación usando

validación cruzada, en realidad no he encontrado una función que lo haga
directamente, ası́ que escribı́ el código que sigue, el cual reproduce la tabla
8.4 del libro.

# Estimación de la probabilidad de clasificación

# errónea por validación cruzada
clasifq<-numeric(nrow(ejemp8_2))
for(i in 1:nrow( ejemp8_2)){
zq<-qda(GRP~.,ejemp8_2[-i,],prior=c(36,36,50,50)/172)
clasifq[i]<-as.numeric(predict(zq,ejemp8_2[i,-9])$class)
}
tabla8_3<-table(ejemp8_2$GRP,clasifq)
addmargins(tabla8_3)

A partir de ese código se puede escribir una función que lo haga para cualquier
situación
Capı́tulo 9

Análisis de correlación canónica

Se realiza el ejemplo 9.2 , hay que revisar porque no están coincidiendo

las correlaciones, que aparecen en el libro.

# lectura de datos
datosY<-scan()
500 43 98 17 800 20 92 32 570 28 98 26 550 28 98 26
550 28 98 26 380 15 99 28 930 21 99 28 650 10 101 27
600 10 101 27 1500 19 99 23 1750 22 101 27 2000
58 100 18 2500 34 102 16 2000 21 105 20 7850 42 84 5
10500 50 81 -12

datosY<-matrix(datos,ncol=4,byrow=TRUE)
colnames(datosY)<-paste("Y",1:4,sep="")

datosX<-scan()
0 3 22 57 17 1 0 16 20 38 13 13 0 6 28 46
17 3 0 4 19 47 27 3 0 1 8 50 35 6 0 2 19
44 32 3 0 0 15 50 27 8 10 21 40 25 4 0 14
26 32 28 0 0 0 1 6 80 12 1 1 4 34 33 22
6 0 7 14 66 13 0 0 9 15 47 21 8 3 7 17
32 27 14 0 5 7 84 4 0 0 3 1 94 4 0

datosX<-matrix(datosX,ncol=6,byrow=TRUE)
colnames(datosX)<-paste("X",1:6,sep="")

33
34 CAPÍTULO 9. ANÁLISIS DE CORRELACIÓN CANÓNICA

ejemp9_2<-data.frame(cbind(datosX,datosY))

# matriz de correlaciones de los datos

round(cor(ejemp9_2),3)

# análisis de correlación canónica

corcan<-cancor(datosX,datosY)
# correlaciones
corcan$cor
# coeficientes estimados para las variables X
corcan$xcoef
# coeficientes estimados para las variables Y
corcan$ycoef
# valores usados para ajustar X
corcan$xcenter
# valores usados para ajustar Y
corcan$ycenter
Capı́tulo 10

Escalamiento multidimensional

No tengo nada de eso

35
36 CAPÍTULO 10. ESCALAMIENTO MULTIDIMENSIONAL
Capı́tulo 11

Análisis de correspondencia

aún no tengo nada

37
38 CAPÍTULO 11. ANÁLISIS DE CORRESPONDENCIA
Apéndice

Rutinas para vectores y matrices

Conformación de matrices
Las instrucciones
u<-matrix(c(2,3,-1,1),nrow=1)
v<-matrix(c(0,2,1),ncol=1)
A<-matrix(c(3,2,1,5),nrow=2)
B<-matrix(c(2,0,4,1),nrow=2)
C<-matrix(c(1,3,4,3,2,1,4,1,3),nrow=3)
Clase<-matrix(c("Pedro","Olga","Pilar","Carlos"),ncol=1 )

Producen los siguientes vectores y matrices

 
0    

3 1 2 4
u = 2 3 −1 1 , v =  1  , A = ,B=
2 5 0 1
2
 
  P edro
1 3 4  Olga 
C =  3 2 1  y Clase =   P ilar 

4 1 3
Carlos
Con las instrucciones
M_UNOS<-matrix(rep(1,6),nrow=2)
M_CEROS <-matrix(rep(0,6),nrow=2)
V_UNOS<-matrix(rep(1,4),nrow=1)
I_2<-diag(1,2)

39
40 CAPÍTULO 11. ANÁLISIS DE CORRESPONDENCIA

se genera una matriz de tamaño 2 × 3 de unos, una matriz nula, un vector

de unos y una matriz identidad de tamaño 2 × 2
   
1 1 1 0 0 0
M _U N OS = , M _CEROS =
1 1 1 0 0 0
 

1 0
V _U N OS = 1 1 1 1 e I_2 =
0 1

Traspaso de un archivo de datos a una matriz

Podemos crear una matriz usando las columnas de un marco de datos (data
frame). Para el ejemplo usaremos el marco de datos women, de los datos de
ejemplo de R1

data(women)
W <- as.matrix(women)

Operaciones y transformaciones sobre matrices

Entre las matrices A y B, señaladas arriba, se desarrollan las operaciones y
transformaciones respectivas mediante la sintaxis de R, la cual se describe a
continuación.  
5 5
C<-A+B; produce la suma entre A y B, .
2 6
 
1 −3
D<-A-B; Resta B de A, ..
2 4
 
1 −3
E<-A%*%B; Producto entre A y B, .
2 4
 
6 4
F_1<- A*B; Producto entre elementos correspondientes de A y B, .
0 5
 
9 1
F_2<- A^B; Cada elemento de A elevado al correspondiente de B, .
1 5
1
Para tener un listado y una corta descripción de los marcos de datos disponibles
en el paquete (librerı́a) datasets de R, use el comando data(), si quiere un listado de
los marcos de datos de ejemplo de todas las librerı́as instaladas use data(package =
.packages(all.available = TRUE))
 41
 
9 1
G_1<-A^2; Las entradas G_1 son el cuadrado de las entradas de A, .
4 25
 
3 49 51
G_2<-A%*%A%*%A; Matriz A multiplicada por si misma tres veces (A ), .
102 151
 
3 1
 2 5 
H<-rbind(A,B); Dispone la matriz B debajo de la matriz A , H =   2 4 .


 0 1 
3 1 2 4
H<-cbind(A,B); Dispone la matriz B al lado de la matriz A , H = .
2 5 0 1
 
6 12 2 4
 0 3 0 1 
K<-A%x%B; Producto directo (Kronecker) entre A y B , K =   4 8 10 20 .


0 2 0 5 
0,67 4,00
L<-B/A; Divide cada elemento de B por el respectivo de A , L = .
0,00 0,20

m<-6:10; Genera el vector m con valores enteros entre 6 y 10, m = 6 7 8
9 10 .
n<-8:5; Genera el vector n con valores enteros entre 8 y 5, m = 8 7 6 5 .
library(Matrix);
O<-as.matrix(bdiag(A,B));
 Matriz en bloques; A y B en la diagonal O =
3 1 0 0
 2 5 0 0 
 0 0 2 4 .
 

0 0 0 1
Det_A<-det(A); Determinante de A , Det_A = 13.
diag(A); Vector

de tipo numérico cuyas componentes son la diagonal de A,
d= 3 5 .
DG_A<-diag(diag(A));
  Transforma la matriz A en una matriz diagonal,
3 0
DG_A = .
0 5
 
0,3846 −0,0769
Inv_A<-solve(A); Inversa de la matriz A, Inv_A = .
 −0,1538
 0,2308
3 2
T_A<-t(A); Produce la transpuesta de A, T _A = .
1 5
Vap_C<-eigen(C)$values Obtiene
los valores propios de C,
V ap_C = 7,470 1,399 −2,870
Vep_C<-eigen(C)$vectors Obtiene los vectores propios de C,
42 CAPÍTULO 11. ANÁLISIS DE CORRESPONDENCIA
 
−0,61 −0,04 0,79
V ep_C =  −0,45 −0,80 −0,39 
−0,65 0,60 −0,47
Filas_C<-nrow(C); Cuenta el numero de filas de C , F ilas_C = 3.
Colum_C<-ncol(C); Cuenta el numero de columnas de C , Colum_C = 3.
Tra_C<-sum(diag(C)); Calcula la traza de la matriz C , T ra_C = 6.

b<-matrix(c(3,2),ncol=1)
X<-solve(A,b)
    
3x1 + x2 = 3 3 1 3
resuelve el sistema con A = y b = .
  2x 1 + 5x 2 = 2 2 5 2
1
Solución X =
0

svd(A)$u
svd(A)$v
diag(svd(A)$d)

Encuentra
 la descompocición
 singular de la matriz
 A (A2.29)
 esta es: U 
=
−0,416 −0,909 −0,526 −0,851 5,834 0,000
,V = y∆=
−0,909 0,416 −0,851 0,526 0,000 2,228