lunes, 14 de septiembre de 2009

Análisis de Correspondencia en R

Correspondencias3

Correspondencias4

Aplicaciones-Ejemplos
Aplicaciones-Ejemplos


# Correspondence Analysis
#cargar los paquetes necesarios para el análisis
library(ade4)
library(ca)
#cargar datos
datos<-read.table("correspondencias3.txt", header=TRUE)
attach(datos)
#crear una tabla de 2 vías o una tabla de contengencia
i) write.table(T,"T4.txt",sep=";")
mydata<-read.table("T4.txt",sep=";")
mydata<-table(datos,exclude=c(NA,NaN)) #Otra posibilidad ;mydata=datos
ii)mytable<-with(datos,table(categoria,consumo)) #ej:A=categoria, B=consumo; A and B are categorical factors.
prop.table(mytable, 1) # row percentages
prop.table(mytable, 2) # column percentages
> T<-mytable
> dimnames(T)<-NULL
> rownames(T)<-c("Blanca","Negra","Otra")
> colnames(T)<-c("Nor-Este","Sur-Este","Oeste")
#aplicar el análisis de correspondencia
i) acs<-dudi.coa(df=mydata,scann=FALSE,nf=3)
> acs
Duality diagramm
class: coa dudi
$call: dudi.coa(df = datos)

$nf: 2 axis-components saved
$rank: 2
eigen values: 0.04332 0.004576 #valores de la inercia
vector length mode content
1 $cw 3 numeric column weights
2 $lw 4 numeric row weights
3 $eig 2 numeric eigen values

data.frame nrow ncol content
1 $tab 4 3 modified array
2 $li 4 2 row coordinates
3 $l1 4 2 row normed scores
4 $co 3 2 column coordinates
5 $c1 3 2 column normed scores
other elements: N

ii) fit <- ca(mytable,nd=3)
> fit
Principal inertias (eigenvalues):
1 2 3
Value 0.084191 0.015423 0.008126
Percentage 78.14% 14.32% 7.54%

Rows:
1 2 3 4 5
Mass 0.130081 0.126016 0.256098 0.341463 0.146341
ChiDist 0.414477 0.426678 0.264753 0.178265 0.479517
Inertia 0.022347 0.022942 0.017951 0.010851 0.033649
Dim. 1 -1.305799 -1.423420 0.791900 -0.228650 1.534124
Dim. 2 -1.330605 -0.654353 -0.332692 1.326862 -0.767570

Columns:
1 2 3 4
Mass 0.284553 0.203252 0.337398 0.174797
ChiDist 0.291026 0.339951 0.202878 0.514458
Inertia 0.024101 0.023489 0.013887 0.046263
Dim. 1 0.855732 1.032369 -0.465832 -1.694316
Dim. 2 -1.083358 0.602050 1.091590 -1.043475

#

Aplicaciones

Análisis factorial en R

Factorial1b

Factorial 2

Aplicaciones-Ejemplos

 datos<-read.table("factorial.txt",header=TRUE,row.names=1)
attach(datos)
fac<-prcomp(datos, retx=,center=TRUE,scale.=TRUE,tol=NULL)
#En primer lugar, es importante, imprimir el resumen del an´alisis y la gr´afica de los autovalores
#(plot) para determinar el n´umero de factores
summary(fac)
plot(fac)
fac$rotation #cargas factoriales
fac$x

cargas<-matrix(0,7,3) #16=2ºdim(datos) y 4=nºfactores a retener
for (i in 1:3) cargas[,i]<-fac$rotation[,i]*fac$sdev[i]
cargas2<-varimax(cargas,normalize=T)$loadings
print(cargas2,cutoff=0.3)

comunalidad<-matrix(0,7,3)

for (i in 1:7)
{for (j in 1:3)
{comunalidad[i,1]=comunalidad[i,1]+cargas[i,j]^2
comunalidad[i,2]=1-comunalidad[i,1]}}
comunalidad

cargas3<-matrix(0,7,3)
cargas4<-matrix(0,7,3)
for (i in 1:7){
cargas3[i,]<-cargas2[i,]^2
cargas4[i,]<-cargas3[i,]/comunalidad[i,1]}
cargas3
cargas4

ajuste<-cor(datos)-cargas2%*%t(cargas2)
ajuste

biplot(fac)
par(mfrow=c(3,1))
for (i in 1:3){
plot(cargas2[,i],cargas2[,i+1])
text(cargas2[,i],cargas2[,i+1],labels=row.names(datos))}

par(mfrow=c(3,2))
for (i in 1:3){
for(j in 1:4){
{plot(fac$x[,i],fac$x[,j])}
text(fac$x[,i],fac$x[,j],labels=row.names(fac$x))}}

facmle<-vector("list",3)
for (i in 1:3) {facmle<-factanal(datos2,i)}
facmle

function (xmat, factors=NULL, cor=TRUE) {
prc <- princomp ( covmat = xmat ,cor = cor )
eig <- prc$sdev^2
if (is.null(factors))
factors <- sum ( eig >= 1 )
loadings <- prc$loadings [ , 1:factors ]
coefficients <- loadings [ , 1:factors ] %*% diag ( prc$sdev[1:factors] )
rotated <- varimax ( coefficients ) $ loadings
fct.ss <- apply( rotated, 2 , function (x) sum (x^2) )
pct.ss <- fct.ss / sum (eig)
cum.ss <- cumsum ( pct.ss )
ss <- t ( cbind ( fct.ss , pct.ss, cum.ss ) )
return ( coefficients , rotated , ss )
}

Análisis Discriminante en R

Act 2 Disc Rim in Ante

Act 2 Disc Rim in Ante 3

Aplicaciones-Ejemplos
Aplicaciones-Ejemplos

#Se puede considerar como un análisis de regresión donde la variable dependiente es categórica y tiene como categorías la etiqueta de cada uno de los grupos,
#y las variables independientes son continuas y determinan a qué grupos pertenecen los objetos.

#OBJETIVOS:
#1)Se pretende encontrar relaciones lineales entre las variables continuas que mejor discriminen en los grupos dados a los objetos.
#2)Construir una regla de decisión que asigne un objeto nuevo, que no sabemos clasificar previamente, a uno de los grupos prefijados con un cierto grado de riesgo.
# 1. Determinar si existen diferencias significativas entre los “perfiles” de un conjunto de variables de dos o m´as grupos definidos a priori.
# 2. Determinar cual de las variables independientes cuantifica mejor las diferencias entre un grupo u otro.
# 3. Establecer un procedimiento para clasificar a un individuo en base a los valores de un conjunto de variables independientes
#Es necesario considerar una serie de restricciones o supuestos
#La discriminaci´on se lleva a cabo estableciendo las ponderaciones del valor te´orico de cada variable, de tal forma que maximicen la varianza entre-grupos frente a la intra-grupos.


##################################################################
# Se carga la librería MASS
library(MASS)
> datos<-read.table("discriminante.txt", header=TRUE)

#Existen 2 clasificadores: el clasificador lineal de Fisher y el clasificador “no param´etrico” k-NN
#el clasificador lineal tiende a quedar consistentemente por encima del k-NN en estas condiciones ideales (Caso normal homoced´astico).

#############I) An´alisis discriminante lineal de Fisher: "lda" se hace un análisis discriminante lineal
> attach(datos)
> discrimi<-lda(discriminante~X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13,prior=c(0.33,0.33,0.34), method="moment", tol=0.001)
#Probabilidades previas: Estos valores se utilizan para la clasificaci´on. Se puede elegir entre:
#1. Todos los grupos iguales: las probabilidades previas ser´an iguales para todos los grupos.
#2. Calcular seg´un tama˜nos de grupos
#En este caso consideramos que las probabilidades de pertenencia a cada grupo son iguales. Si quisiéramos que las probabilidades sean proporcionales al tamaño, no indicaríamos nada.

#Del mismo modo usaremos el m´etodo de los momentos

#Salidas:
#Prior probabilities of groups: tabla con las probabilidades de permanencia a cada grupo.
#Group means: las medias por grupo; con esta opci´on, tenemos que analizar, que para cada variable,
#efectivamente, hay diferencias significativas entre los valores de los tres grupos (se puede realizar
#un an´alisis de la varianza, y utilizar el mismo razonamiento que en la parte de SPSS).
#Coefficients of linear discriminants:
#los coeficientes de la funci´on discriminante lineal de Fisher se extraen mediante
> discrimi$scaling
#Esto significa que la puntuaci´on discriminante de cada caso se obtiene con la fórmula:
y = a1 x1 +a2 x2 +a3 x3 + a4 x4 + a5 x5...
#Las puntuaciones discriminantes por grupos se pueden representar gr´aficamente mediante el comando:
> plot(discrimi) #Las puntuaciones se han trasladado de forma que la media global es igual a cero.
#En muchos casos las puntuaciones van a separar correctamente los datos.
#Hay una zona en la que las dos poblaciones se solapan. Esta es la zona que incrementará el valor de las tasas de error.
#Proportion of trace: proporci´on de varianza explicada por cada eje (el primer eje ser´a mucho m´as discriminante que el segundo).
#Con la orden svd obtendremos el poder discriminate de cada funci´on, es decir, que la primera es mucho m´as discriminante que la segunda.
> discrimi$svd

#Para clasificar las observaciones se utiliza la sentencia siguiente:
> predclas <- predict(discrimi)$class
#Estudio de los errores
#Para analizarlos se suele utilizar la representaci´on gr´afica de las observaciones,representando las puntuaciones Z discriminantes y buscando los casos mal representados para su estudio.
#La tasa de error aparente se calcula de la siguiente forma:
> tea <- 1 - sum(predclas == datos[,length(datos)])/n #n es el número de casos
#Para calcular la tasa de error por validaci´on cruzada se utiliza la opci´on CV=T al llamar a la funci´on lda:
> predclas.cv <- lda(datos[,-length(datos)],datos[,length(datos)],CV=T)$class #datos[,length(datos)] refiere a la columna de categorías o clases
#CV=TRUE generates jacknifed (i.e., leave one out) predictions.
> tevc <- 1 - sum(predclas.cv == datos[,length(datos)])/n #n es el número de casos

#Validación de los resultados: Tambi´en podemos crear la matriz de clasificaci´on sin m´as que utilizar la siguiente orden:
> table(predict(discrimi)$class, discriminante) #discriminante: columna de categorías o clases de la matriz Datos
#discriminante: permite observar los individuos mal clasificados, ver toda la matriz exc la diagonal

#Realizar una predicción sobre un par de individuos. Por ejemplo:
> Datos2<-rbind(c(5.6,4.2,15.8, 1.7,0.17,1.6,2.3,6.0,4.0,5.27,0.14,216.25,26.66),c(5.45, 3.09,17.82,1.56,0.14,1.62,1.97,5.21,2.68,7.45,0.15,285.71,23.42))
> Datos21<-data.frame(Datos2) #Antes: colnames(Datos2)<-colnames(Datos[,-length(Datos)])
> discrimi2<-predict(discrimi,newdata=Datos21,prior=discrimi$prior,2)
#Obteniendo las siguientes salidas:
#discrimi2$class
#discrimi2$posterior: las probabilidades de que el individuo seleccionado (row) pertenezca a cada grupo (col)
#discrimi2$x
#Otra forma
predict(discrimi,newdata=Datos21)$class ##predice el grupo de pertenencia de los nuevos datos
> grupo <- predict(discrimi,method="plug-in")$class ##Se predicen los datos originales en los grupos segun la function discriminante
> table(grupo,Type) ##Se observa el numero de datos originales bien y mal clasificados


#Representación gráfica:
> plot(discrimi)
> pairs(discrimi) #produce los pares de gr´aficos de los datos frente a las puntuaciones discriminantes
#histogramas de las variables dependientes frente a la variable de grupos.
#En este caso habr´ýa que hacer todas las variables X1,X2, . . . ,X13 frente a la variable discriminante.
#existe un comportamiento muy diferente de la variable dependiendo del grupo al que pertenece.
##Con el fin de identificar las principales caracter´ýsticas de los datos, podemos representar la matriz de diagramas de dispersi´on:
# Panels of histograms and overlayed density plots, for 1st discriminant function
> plot(discrimi, dimen=1, type="both")
## Exploratory Graph for LDA
> library(klaR)
> partimat(formula,data=datos,method="lda") #formula: discriminante~X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13
##You can also produce a scatterplot matrix with color coding by group.Ej. Scatterplot for 3 Group Problem
> pairs(mydata[c("x1","x2","x3")], main="My Title ", pch=22, bg=c("red", "yellow", "blue")[unclass(mydata$G)])

###############II) El m´etodo de k vecinos m´as pr´oximos

###############III) Comparaci´on de los resultados

Análisis cluster en R

Act1 Cluster

Aplicaciones-Ejemplos

#####Como aplicaci´on vamos a realizar un an´alisis cluster al fichero de datos cluster.sav.########3
#Este fichero contiene los datos de veinte variables (nombradas V2,V3,...,V21) de cohesi´on social de 10 pa´ýses europeos.

#cargamos los datos
datos<-read.table("cluster.txt", header=TRUE,row.names=1)
datos

tres<-read.table("Libro1.txt", header=FALSE)

dist_tres<-dist(t(tres),"eucl")
dist_tres
V1 V2 V3 V4 V5 V6
V2 17.857491
V3 32.336048 44.004659
V4 15.479018 23.482973 23.298498
V5 9.099451 8.825531 37.714188 17.837601
V6 12.323149 10.307764 39.729083 18.878559 7.005712
V7 28.680830 33.971459 20.152667 17.021457 30.292738 32.227783
Z<-dist_tres^2

#En este caso nos ocurre igual que en la secci´on anterior, los datos son muy dispares y es
#conveniente tipificarlos:
datos2<-scale(datos)

par(mar=c(8,2.5,.7,.7),lwd=3,bty="l",cex=1.7)
boxplot(data.frame(scale(datos)),pch=20,cex=.5,las=2,cex.axis=.7)

#tenemos que seleccionar una de las medidas, en este caso aplicaremos la distancia euclídea:
distancia<-dist(datos2,"eucl")
distancia<-distancia^2

#vamos a usar cuatro métodos: el linkaje simple, completo, el de la media y el de Ward.
cluster1<-hclust(distancia,method="complete")
cluster2<-hclust(distancia,method="single")
cluster3<-hclust(distancia,method="ave")
cluster4<-hclust(distancia,method="ward")

distancia->nombre
as.matrix(nombre)->Z
rownames(Z,do.NULL=T,prefix="p")->rownames(Z)

opar <- par(mfrow = c(2, 2))
plot(cluster,main="Average",ylab="Rescaled Distance Cluster Combine", xlab="Label",hang=1,cex.axis=1,cex=.8,cex.lab=.8,labels=x)
rect.hclust(cluster,k=5)
plot(cluster1,main="Complete",ylab="Rescaled Distance Cluster Combine", xlab="Label",hang=1,cex.axis=1,cex=.8,cex.lab=.8,labels=x)
rect.hclust(cluster1,k=5)
plot(cluster2,main="Single",ylab="Rescaled Distance Cluster Combine", xlab="Label",hang=1,cex.axis=1,cex=.8,cex.lab=.8,labels=x)
rect.hclust(cluster2,k=5)
plot(cluster4,main="Ward",ylab="Rescaled Distance Cluster Combine", xlab="Label",hang=1,cex.axis=1,cex=.8,cex.lab=.8,labels=x)
rect.hclust(cluster4,k=5)

#Para terminar solo tendr´iamos que graficar los dendogramas de cada m´etodo, y analizar
#los diferentes resultados para intentar ver la consistencia de los resultados (en este caso solo
#haremos el primero de ellos):
plot(cluster1)

#Obteniendo los mismos resultados que con SPSS. Finalmente calcularemos los respectivos
#coeficientes cofen´eticos para decidir cual de los m´etodos es mejor. Para ello, primero calcularemos
#las matrices cofen´eticas de los cuatro m´etodos y posteriormente la correlaci´on de estas
#con la matriz de distancias, eucl´idea en este caso.
co11<-cophenetic(cluster1)
co12<-cophenetic(cluster2)
co13<-cophenetic(cluster3)
co14<-cophenetic(cluster4)

cor(distancia,co11)
cor(distancia,co12)
cor(distancia,co13)
cor(distancia,co14)

# Por lo que concluiríamos que el mejor método es el “cluster3”, con un coeficiente cofenético
#de 0.795, por lo que le mejor m´etodo es el de la media aunque existen pocas diferencias entre ambos m´etodos.

Procesos markovianos

Aplicaciones

Análisis de componentes Principales (ACP) en R

ACP1

ACP2


 datos<-read.table("comprincipales.txt",header=T,row.names=1)
attach(datos)
datos
#Evidentemente, antes de realizar cualquier an´alisis de componentes principales, podemos realizar un resumen estad´ýstico o gr´aficos descriptivos bidimensionales,
plot(datos) #boxplt(datos,cex.axis=0.75) para gráfico de cajas
summary(datos) #print(summary(datos)) para estadísticos básicos y print(apply(datos,2,sd)) para desviaciones estándar
#Boxplot: la diferencia esos dos grupos de variables y justifica la normalizaci¶on de los datos antes de realizar el ACP (ACP normado).
#Del mismo modo, es interesante, estudiar la matriz de correlaciones, y ver, que estas sean en general altas, ya que esta es una de las hip´otesis para el an´alisis de componentes principales.
#El espacio de las variables del ACP normado es una imagen de la matriz de correlaciones, que se obtiene en R con la funci¶on cor, cuando se aplica sobre un objeto de tipo data.frame o matrix de datos num¶ericos. La funci¶on plot sobre el mismo data.frame, cafe produce la ¯gura 2.
cor(datos)
#matriz de correlaciones y matriz de gr¶aficas de dispersi¶on: print(cor(datos),2) y plot(datos)

#Se realiza el ACP normado del ejemplo caf¶e con todas las variables como activas. Las funciones utilizados de ade4 son dudi.pca, que realiza el ACP; inertia.dudi, que calcula las contribuciones absolutas y relativas a la inercia.
library(ade4)
acp<-dudi.pca(df=datos,scannf=T,nf=2) #no es scannf=F??? ver "analmultir.pdf"
#Select the number of axes: 2
#De esta manera generaremos el an´alisis de componentes principales y a su vez obtenemos la representaci´on de las gr´afica de los autovalores.
#Para ver la importancia (contribuciones) absolutas y relativas, vamos a usar la funci´on
acpi<-inertia.dudi(acp, row.inertia=T, col.inertia=T) #contiene las ayudas a la interpretación del ACP

#A continuaci´on vamos a ir analizando las salidas que nos proporciona R.
# En primer lugar vemos los resultados para las filas. En este caso, obtendremos la representaci´on de cada fila en el espacio bidimensional (normalizado y sin normalizar).
acp$l1
acp$li
#Estos puntos los podemos representar con la orden s.label de la forma:
s.label(acp$li)
#Del mismo modo la representaci´on de las columnas ser´a:
acp$co
acp$c1
s.label(acp$co)
#Tambi´en podemos analizar las contribuciones a la inercia de la filas (o columnas) de modo que:
acpi
#Tambi´en podemos obtener la representaci´on conjunta de filas y columnas sin m´as que:
biplot(acp$co,acp$li)
#Podemos obtener, finalmente, una representaci´on de las correlaciones de las variables, con la orden:
s.corcircle(acp$li); s.corcircle(acp$co)

#NOTAS: para imprimir todas las salidas
cat("\n Estadísticas básicas \n")
print(summary(datos))
cat("\n Desviaciones estándar \n")
print(apply(datos,2,sd))
par(ask=TRUE) # para pausa antes de la gr¶afica
boxplot(datos,cex.axis=0.75)
acp <- dudi.pca(cafe,scannf=F,nf=2) # realiza el ACP
acpI <- inertia.dudi(acp,row.inertia=T,col.inertia=T) # acpI contiene las ayudas a la interpretaci¶on del ACP
impresión de los objetos de acp y acpi con títulos:
cat("\n Valores propios \n")
print(acpI$TOT,2)
plot(acp$eig) #produce la grafica de los valores propios, que sirve como com-plemento a la tabla de correlaciones.
cat("\n Vectores propios \n")
print(acp$c1)
cat("\n Coordenadas de las columnas \n")
print(acp$co)
cat("\n Contribuciones de las columnas a los ejes \n")
print(acpI$col.abs/100)
cat("\n Calidad de representaci¶on de las columnas \n")
print(acpI$col.rel/100)
cat("\n Calidad de representaci¶on de las columnas en el plano \n")
print(acpI$col.cum/100)
cat("\n Coordenadas de las filas \n")
print(acp$li)
cat("\n Contribuciones de las filas a los ejes \n")
print(acpI$row.abs/100)
cat("\n Calidad de representaci¶on de las filas en los ejes \n")
print(acpI$row.rel/100)
cat("\n Calidad de representaci¶on de las filas en el plano \n")
print(acpI$row.cum/100)
# gr¶aficas del acp
par(mfrow=c(1,2)) # para 2 gr¶aficas simult¶aneas
s.corcircle(acp$co,sub="Cafe - C¶³rculo de correlaciones",possub= "bottomright") #presenta el círculo de correlaciones; muestra el efecto "tamaño" del primer factor; es la clave para la lectura del primer plano factorial de los cafés
iden <- c("EC","C4M","C4C","C2M","C2C","EO","O4M","O4C","O2M","O2C") # se define iden para tener etiquetas m¶as cortas en las gr¶aficas
s.label(acp$li,label=iden,sub="Preparaciones de caf¶e",possub= "bottomright") # s.label es de ade4 y coloca puntos con etiquetas en planos factoriales #primer plano factorial de lso cafés