rm(list=ls())
library(FactoClass) # carga el paquete

data(admi)
Y<-admi[,8:11] #tomamos las columnas genero, estrato, origen y edad

# estracto de la tabla: múltiplo de 25
sel<-seq(25,425,25)
cat("\nDimensión de Y:",dim(Y)[1],"x",dim(Y)[2])
Y[sel,]

Z<-acm.disjonctif(Y) #función de ade4
cat("\nDimensión de Z:",dim(Z))
names(Z)<-substr(names(Z),6,8)
Z[sel,] #sel<-seq(25,425,25)

cat("\nExtracto de 4 filas de las tablas Y y Z para ver la recodiificación")
Y[sel[1:4],]
Z[sel[1:4],]

Z<-as.matrix(Z)
B<-t(Z)%*%Z
B

cat("\nDp: matriz diagonal con elementos de la tabla de BURT (B)")
Dp<-diag(diag(B))
Dp

n<-nrow(Y) #numero de filas
s<-ncol(Y) #numero de columnas
F<-1/n/s*Z
cat("\nExtracto de 6 columnas de la tabla F")
round(addmargins(head(F[sel,]),2),3)
cat("Sum col",round(colSums(F),3))
cat("\n\nMarginal de cada fila 1/445:",round(1/nrow(Z),5))

head(1/4*Z[sel,],4)

cat("\nDiagonal de la matriz de métrica\n\n")
nj = n*s*1/diag(Dp)
names(nj)=colnames(Z)
t(data.frame(metrica=round(nj,1)))

cat("\nCentro de gravedad en porcentaje\n\n")
t(round(data.frame(centroide=colSums(1/n/s*Z)*100),1))

Y[c(50,100),]

acm<-dudi.acm(Y,scannf=FALSE,nf=3)
inertia(acm)
barplot(acm$eig,main="Barplot eigen-values")
abline(h=sum(acm$eig)/8,col="darkgreen")

plot(acm,Tcol=FALSE,gg=TRUE,roweti=sel,main="Primer plano factorial del ACM")

plot(acm,Trow=FALSE,gg=TRUE,main="Factores 1-2")

ayuacm<-inertia(acm,,T)
names(ayuacm)
round(cbind(acm$co,ayuacm$col.abs,ayuacm$col.rel,ayuacm$col.contrib),2)
round(sum(acm$eig)/12*100,2)

plot(acm,roweti=sel,main="Representación simultanea 1er plano",gg=TRUE)
t(Y[sel,])

plot(acm,2,3,roweti=sel,main="Representación simultanea factores 2-3",gg=TRUE)
t(Y[sel,])

plotfp(acm$cr,main="Factores 1-2, razones de correlación") # función de FactoClass

plotfp(acm$cr,2,3,main="Factores 2-3, razones de correlación")

# Número de admitidos por carreras:
t(data.frame(nCarr=summary(admi$carr)))
# Valor test:
vtQui = sqrt(63*(445-1)/(445-63))*(-0.259)
cat("Valor test de química:",round(vtQui,3))


supCarr<-supqual(acm,admi$carr)

attach(supCarr)
cat("\nCoordenadas y ayudas para la interpretación")
cbind(ncat,round(cbind(dis2,coor,tv,cos2),3))
cat("\nRazones de correlación\n")
round(scr*100,2)
detach(supCarr)

plot(acm,Trow =FALSE,infaxes="no",main="Plano 1-2",
ylim =c( -1.5 ,1.3) , col.col=" black ")
points( supCarr$coor ,col ="darkgreen")
text(supCarr$coor , labels = attributes (admi[,1])$levels ,
col ="darkgreen",pos =1, cex =0.8 , font =3)

plot(acm,2,3,Trow =FALSE,infaxes="no",main="Plano 2-3",
ylim =c( -1.5 ,1.3) , col.col=" black ")
points(supCarr$coor[,c(2,3)] ,col ="darkgreen")
text(supCarr$coor[,c(2,3)] , labels = attributes (admi[,1])$levels ,
col ="darkgreen",pos =1, cex =0.8 , font =3)

par(mfrow=c(2,2))
for (varac in 8:11){
    K<-unclass(table(admi$carr,admi[,varac]))
    plotct(K,"row",col=2:5,main=colnames(admi)[varac])
}

rm(list=ls())
# Lectura de los datos, están en formato R
load("../Datos/ninios8a11.rda")
# variables activas
Y <- subset(ninios8a11,select=c(Teat,Libr,Cine,Vide,Radi,Musi))
# variables suplementarias
Ys<-ninios8a11[,c(2,29,30,32,35)]

par(las=1,mfrow=c(2,3),mai=c(0.55,0.5,0.1,0.1))
for(i in 1:6){
  cat<-attributes(Y[,i])$levels;
  per<-tabulate(Y[,i])/nrow(Y)*100;
  pl<-plot(Y[,i],horiz=TRUE,col=gray(seq(1.0,0.9,
           length=length(cat))),ylim=c(0,8),
           xlim=c(0,1400),xlab=colnames(Y)[i]);
  text(800,pl,round(per,1),cex=0.8,pos=4); 
  }  


table(Ys$Edad) #para ver las frecuencias como venían
Ys$Edad=factor(as.numeric(Ys$Edad),labels=paste("a",8:11,sep=""))
table(Ys$Edad) #ver las frecuencias como quedaron              


par(las=1,mfrow=c(2,3),mai=c(0.55,0.5,0.1,0.1))
for(i in 1:5){
  cat<-attributes(Ys[,i])$levels;
  per<-tabulate(Ys[,i])/nrow(Y)*100;
  pl<-plot(Ys[,i],horiz=TRUE,col=gray(seq(1.0,0.9,
           length=length(cat))),ylim=c(0,8),
           xlim=c(0,1400),xlab=colnames(Ys)[i]);
  text(800,pl,round(per,1),cex=0.8,pos=4); 
}  


acm<-dudi.acm(Y,scannf = FALSE,nf=3)
eigtab<-data.frame(valp=acm$eig,porc=acm$eig/sum(acm$eig)*100,
                   pacu=cumsum(acm$eig)/sum(acm$eig)*100)
cbind(eje=1:8,round(eigtab[1:8,],4),eje=9:16,round(eigtab[9:16,],4),
       eje=17:24,round(eigtab[17:24,],4))
#      ,digits=c(0,rep(c(0,3,1,1),3)))
# criterio de Benzecri
s<-6; cat("\n1/s:",round(1/s,4)) 
# --> se calcula tau para los primeros 11 ejes
eig11<-acm$eig[1:11]
tau<-(s/(s-1))^2*(eig11-(1/s))^2
ptau<-tau/sum(tau)*100

barplot(acm$eig,las=2)
abline(h=1/s,col="darkgreen")
cat("\nAltura linea verde = 1/s")

barplot(ptau,las=2)
abline(h=sum(ptau)/11,col="blue")
cat("\nAltura linea azul = promedio de los taus")

plot(acm,Trow=FALSE,gg=TRUE,xlim=c(-1,1),ylim=c(-1,1.3),cframe=1.1,col.col="black",cex.global=0.8)

ayu<-inertia.dudi(acm,,T) 
cat("\ncategoría | peso  |    coordenadas        |     cont.absoluta   |  coseno cuadrado")
round(cbind(peso=acm$cw,acm$co,ayu$col.abs/100,abs(ayu$col.rel)/100),4)

sc<-supqual(acm,Ys)
plotfp(as.data.frame(sc$coor),col="darkgreen",cframe=1,gg=TRUE)


round(data.frame(por=sc$ncat/34.76,dis2=sc$dis2,coor=sc$coor,vt=sc$tv,cos2=sc$cos2),4)

plotfp(acm$cr,xlim=c(0,0.6),cframe=1)