Análisis de correlación en la Casen 2017
VE-CC-CF
DataIntelligenceViernes 4-06-2021
Abstract
Debemos descubrir correlaciones dentro de las encuesta Casen. La dificultad es que, en general, la información ofrecida es de tipo categórico, por lo que se deben aplicar test no paramétricos. En éste informe y como primer desarrollo los aplicaremos en el estudio de la correlacion ingreso (variable continua) y pertenencia al DAU autonomo, una variable categórica. Queremos comparar si hay diferencias en el valor de una variable dependiente entre tres o más grupos.
1 Introduccion
Comparemos el ingreso que reciben las personas bajo 10 condiciones distintas. ¿Existen diferencias significativas dependiendo de ellas?
2 Analisis exploratorio y tratamiento de los datos
ab <- readRDS(file = "casen_2017_c.rds")# ab <- ab[c(1:1000),]ytotcor: ingreso total corregido dau: Decil autónomo nacional 1 I - 10 X
ingresos <- ab$ytotcor
dau <- ab$dau
nuevo_df <- cbind(ingresos, dau)
head(nuevo_df,10)## ingresos dau
## [1,] 272000 6
## [2,] 159667 6
## [3,] 786932 9
## [4,] 301000 9
## [5,] 220000 4
## [6,] 200000 4
## [7,] NA 4
## [8,] 296634 4
## [9,] 11052 4
## [10,] 500000 4reemplazamos los NA por el promedio
#
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericnuevo_df_sin_NA <- na.aggregate(nuevo_df)
head(nuevo_df_sin_NA,10)## ingresos dau
## [1,] 272000.0 6
## [2,] 159667.0 6
## [3,] 786932.0 9
## [4,] 301000.0 9
## [5,] 220000.0 4
## [6,] 200000.0 4
## [7,] 415297.9 4
## [8,] 296634.0 4
## [9,] 11052.0 4
## [10,] 500000.0 4convertimos a dataframe
nuevo_df_sin_NA <- as.data.frame(nuevo_df_sin_NA)
head(nuevo_df_sin_NA,10)## ingresos dau
## 1 272000.0 6
## 2 159667.0 6
## 3 786932.0 9
## 4 301000.0 9
## 5 220000.0 4
## 6 200000.0 4
## 7 415297.9 4
## 8 296634.0 4
## 9 11052.0 4
## 10 500000.0 4calculamos los promedios por grupo
aggregate(ingresos ~ dau, data = nuevo_df_sin_NA, FUN = mean)## dau ingresos
## 1 1.000000 205108.2
## 2 2.000000 247272.2
## 3 3.000000 277514.0
## 4 4.000000 299047.8
## 5 4.970499 443595.0
## 6 5.000000 327622.3
## 7 6.000000 359488.9
## 8 7.000000 403316.3
## 9 8.000000 480173.7
## 10 9.000000 647311.8
## 11 10.000000 1464206.6calculamos la Desviacion estandar
aggregate(ingresos ~ dau, data = nuevo_df_sin_NA, FUN = sd)## dau ingresos
## 1 1.000000 164055.6
## 2 2.000000 162871.8
## 3 3.000000 159989.2
## 4 4.000000 157984.6
## 5 4.970499 170491.6
## 6 5.000000 167330.3
## 7 6.000000 184688.7
## 8 7.000000 218160.9
## 9 8.000000 278462.3
## 10 9.000000 439739.7
## 11 10.000000 2108267.1Excluimos los outliers
#
Q <- quantile(nuevo_df_sin_NA$ingresos, probs=c(.25, .75), na.rm = T)
iqr <- IQR(nuevo_df_sin_NA$ingresos)
eliminated <- subset(nuevo_df_sin_NA, nuevo_df_sin_NA$ingresos > (Q[1] - 1.5*iqr) & nuevo_df_sin_NA$ingresos < (Q[2]+1.5*iqr))
eliminated <- data.frame(lapply(eliminated, as.character), stringsAsFactors=FALSE)
# despleguemos los primeros 100 registros en pantalla del subset creado:
# eliminated_100 <- eliminated[c(1:100),]#eliminated <- as.double(eliminated)
head(eliminated,10)## ingresos dau
## 1 272000 6
## 2 159667 6
## 3 301000 9
## 4 220000 4
## 5 2e+05 4
## 6 415297.931619351 4
## 7 296634 4
## 8 11052 4
## 9 5e+05 4
## 10 415297.931619351 4convirtamos a matriz:
a = as.double(eliminated$ingresos)
b = as.double(eliminated$dau)
a_b <- cbind(a,b)
head(a_b,10)## a b
## [1,] 272000.0 6
## [2,] 159667.0 6
## [3,] 301000.0 9
## [4,] 220000.0 4
## [5,] 200000.0 4
## [6,] 415297.9 4
## [7,] 296634.0 4
## [8,] 11052.0 4
## [9,] 500000.0 4
## [10,] 415297.9 4boxplot(a ~ b, data=a_b, col="green", cex.axis=0.7,las = 2, ylab="var1", xlab="Treatments", cex.lab=0.75)
ANOVA
Hacemos un modelo ANOVA para nuestra Variable 1 (se haría lo mismo para la Variable 2):
nuestraanova <- lm(ingresos ~ dau, data=eliminated)
nuestromodelo <- anova(nuestraanova)
nuestromodelo## Analysis of Variance Table
##
## Response: ingresos
## Df Sum Sq Mean Sq F value Pr(>F)
## dau 10 8.4733e+14 8.4733e+13 3417.2 < 2.2e-16 ***
## Residuals 198609 4.9248e+15 2.4796e+10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1¿Son nuestros datos normales y homogéneos?
El test ANOVA requiere normalidad de los residuos y homogeneidad en las varianzas.
Representamos visualmente los residuos:
par(mfrow=c(2,2))
plot(nuestraanova)
#shapiro.test(residuals(nuestraanova)) #fligner.test(ingresos ~ deciles, data=eliminated)head(eliminated,10)## ingresos dau
## 1 272000 6
## 2 159667 6
## 3 301000 9
## 4 220000 4
## 5 2e+05 4
## 6 415297.931619351 4
## 7 296634 4
## 8 11052 4
## 9 5e+05 4
## 10 415297.931619351 4eliminated$ingresos <- as.double(eliminated$ingresos)
eliminated$dau <- as.double(eliminated$dau)
head(eliminated,10)## ingresos dau
## 1 272000.0 6
## 2 159667.0 6
## 3 301000.0 9
## 4 220000.0 4
## 5 200000.0 4
## 6 415297.9 4
## 7 296634.0 4
## 8 11052.0 4
## 9 500000.0 4
## 10 415297.9 4Hacemos un test Kruskal-Wallis:
kruskal.test(ingresos~ dau, data = eliminated)##
## Kruskal-Wallis rank sum test
##
## data: ingresos by dau
## Kruskal-Wallis chi-squared = 25013, df = 10, p-value < 2.2e-16pairwise.wilcox.test(eliminated$ingresos,eliminated$dau,
p.adjust.method = c("bonferroni"),paired=F,exact=F)##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: eliminated$ingresos and eliminated$dau
##
## 1 2 3 4 4.97049914212116 5
## 2 < 2e-16 - - - - -
## 3 < 2e-16 < 2e-16 - - - -
## 4 < 2e-16 < 2e-16 < 2e-16 - - -
## 4.97049914212116 < 2e-16 < 2e-16 < 2e-16 < 2e-16 - -
## 5 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.9e-15 -
## 6 < 2e-16 < 2e-16 < 2e-16 < 2e-16 6.1e-10 < 2e-16
## 7 < 2e-16 < 2e-16 < 2e-16 < 2e-16 2.2e-05 < 2e-16
## 8 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1 < 2e-16
## 9 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1 < 2e-16
## 10 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1 < 2e-16
## 6 7 8 9
## 2 - - - -
## 3 - - - -
## 4 - - - -
## 4.97049914212116 - - - -
## 5 - - - -
## 6 - - - -
## 7 < 2e-16 - - -
## 8 < 2e-16 < 2e-16 - -
## 9 < 2e-16 < 2e-16 < 2e-16 -
## 10 < 2e-16 < 2e-16 < 2e-16 < 2e-16
##
## P value adjustment method: bonferronipairwise.wilcox.test(x = eliminated$
ingresos, g = eliminated$dau, p.adjust.method = "holm" )##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: eliminated$ingresos and eliminated$dau
##
## 1 2 3 4 4.97049914212116 5
## 2 < 2e-16 - - - - -
## 3 < 2e-16 < 2e-16 - - - -
## 4 < 2e-16 < 2e-16 < 2e-16 - - -
## 4.97049914212116 < 2e-16 < 2e-16 < 2e-16 < 2e-16 - -
## 5 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 -
## 6 < 2e-16 < 2e-16 < 2e-16 < 2e-16 5.5e-11 < 2e-16
## 7 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.6e-06 < 2e-16
## 8 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.203 < 2e-16
## 9 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.073 < 2e-16
## 10 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.617 < 2e-16
## 6 7 8 9
## 2 - - - -
## 3 - - - -
## 4 - - - -
## 4.97049914212116 - - - -
## 5 - - - -
## 6 - - - -
## 7 < 2e-16 - - -
## 8 < 2e-16 < 2e-16 - -
## 9 < 2e-16 < 2e-16 < 2e-16 -
## 10 < 2e-16 < 2e-16 < 2e-16 < 2e-16
##
## P value adjustment method: holmII Correlaciones entre dos variables categoricas:
ab <- readRDS(file = "casen_2017_c.rds")num_personas <- ab$numper
escolari <- ab$esc
num_escolari <- cbind(num_personas,escolari)
head(num_escolari,10)## num_personas escolari
## [1,] 2 12
## [2,] 2 14
## [3,] 2 12
## [4,] 2 12
## [5,] 3 12
## [6,] 3 12
## [7,] 3 7
## [8,] 5 8
## [9,] 5 9
## [10,] 5 12num_escolari_sin_NA <- na.aggregate(num_escolari)chisq.test(num_escolari_sin_NA)## Warning in chisq.test(num_escolari_sin_NA): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: num_escolari_sin_NA
## X-squared = 219771, df = 216438, p-value = 2.319e-07Referencias:
ANOVA y Kruskal-Wallis
https://rpubs.com/aafernandez1976/anovakw
I Test Kruskal-Wallis
https://rpubs.com/Joaquin_AR/219504
Análisis de variables categóricas con R
https://biocosas.github.io/R/060_analisis_datos_categoricos.html