# Tabla de contingencia tridimensonal
# ING x OCU X SEX con los datos del archivo Ingresos.rda
load(url("http://pagines.uab.cat/plopez/sites/pagines.uab.cat.plopez/files/Ingresos.rda"))

# Para ejecutar este script en R o RStudio es necesario instalar el paquete Deducer y el paquete vcd
# install.packages("Deducer")
# install.packages("vcd")
library(Deducer)
library(vcd)

# Tabla ING x SEX x OCU con las frecuencias solamente
ING3<-contingency.tables(row.vars=d(ING),col.vars=d(SEX),stratum.var=OCU,data=Ingresos)
print(ING3,prop.r=F,prop.c=F,prop.t=F)

# Tablas ING x OCU y ING x SEX 
ING1<-contingency.tables(row.vars=d(ING),col.vars=d(OCU,SEX),data=Ingresos)
ING1<-add.chi.squared(ING1)
print(ING1,prop.r=T,prop.c=T,prop.t=F,expected.n=T,adj.residuals=T)
# El gráfico se obtiene eliminando los NA con la especificación na.omit(Ingresos)
ggplot() + geom_bar(aes(y = ..count..,x = OCU,fill = ING),data=na.omit(Ingresos),
           position=position_fill()) + scale_y_continuous(labels=scales::percent) + ylab ("%")
ggplot() + geom_bar(aes(y = ..count..,x = SEX,fill = ING),data=na.omit(Ingresos),
           position=position_fill()) + scale_y_continuous(labels=scales::percent) + ylab ("%")
assocstats(table(Ingresos$ING, Ingresos$OCU))
assocstats(table(Ingresos$ING, Ingresos$SEX))

# Tabla OCU x SEX 
ING2<-contingency.tables(row.vars=d(OCU),col.vars=d(SEX),data=Ingresos)
ING2<-add.chi.squared(ING2)
print(ING2,prop.r=T,prop.c=T,prop.t=F,expected.n=T,adj.residuals=T)
ggplot() + geom_bar(aes(y = ..count..,x = SEX,fill = OCU),data=na.omit(Ingresos),
           position=position_fill()) + scale_y_continuous(labels=scales::percent) + ylab ("%")
assocstats(table(Ingresos$OCU, Ingresos$SEX))

# Tabla ING x SEX x OCU
ING3<-add.chi.squared(ING3)
print(ING3,prop.r=T,prop.c=T,prop.t=F,expected.n=T,adj.residuals=T)
ggplot() + geom_bar(aes(y = ..count..,x = SEX,fill = ING),data=na.omit(Ingresos), position=position_fill()) +
           facet_wrap(facets = ~OCU) + ggtitle(label = 'OCU') +
           scale_y_continuous(labels=scales::percent) + ylab ("%")
# Instrucción alternativa para el gráfico en el tratamiento de los valores perdidos
ggplot() + geom_bar(aes(y = ..count..,x = SEX,fill = ING),data=subset(Ingresos,!is.na(ING)&!is.na(OCU)),
           position = position_fill()) + facet_wrap(facets = ~OCU) + ggtitle(label = 'OCU') + 
           scale_y_continuous(labels=scales::percent) + ylab ("%")
summary(subset(Ingresos,subset=OCU=="Bajo"))
OCU1=subset(Ingresos,subset=OCU=="Bajo")
assocstats(table(OCU1$ING, OCU1$SEX))
summary(subset(Ingresos,subset=OCU=="Alto"))
OCU2=subset(Ingresos,subset=OCU=="Alto")
assocstats(table(OCU2$ING, OCU2$SEX))

# Obtención solamente de la V de Cramer de todas las tablas
assocstats(table(Ingresos$ING, Ingresos$OCU))
assocstats(table(Ingresos$ING, Ingresos$SEX))
assocstats(table(Ingresos$OCU, Ingresos$SEX))
OCU1=subset(Ingresos,subset=OCU=="Bajo")
assocstats(table(OCU1$ING, OCU1$SEX))
OCU2=subset(Ingresos,subset=OCU=="Alto")
assocstats(table(OCU2$ING, OCU2$SEX))

