CADI Analisis Estadístico Con R

En esta sesión vamos analizar a los admitidos No TEC de Agosto Diciembre 2019 de profesional.

Vamos a utilizar la base de datos (Base_CADI.xlsx) que incluye las siguientes varaiables:

Variable	Descripción
Escuela	Escuela del alumno
Sexo	Hombre, Mujer
PNAi	Promedio de prepa al 2do año
PNA	Promedio de prepa de graduación
PAA	Puntaje del examen de admisión

###Intalación de librerias

##Procesamiento de bases de datos Selección de variables en la base de datos

#Datos1<-filter(Datos,V1=="Categoría")
#Datos1<-filter(Datos,V2>Número)
B0<-select(Datos1,c("Escuela","Sexo","PNA","PNAi","PAA"))

###Identificación y corrección de valores faltantes

colSums(is.na(B0))

## Escuela    Sexo     PNA    PNAi     PAA 
##       0       0       0       0       2

Hay dos valores faltantes en PAA

Podemos quitarlos o reemplazarlos por la media

#Borrar los registros faltantes 
#B0<- na.omit(B0)
#Reemplazarlos por la media de las variables 

B0$PAA <- na.aggregate(B0$PAA,FUN=mean) #se puede usar la media, moda o la mediana
colSums(is.na(B0)) #para validar que funcionó

## Escuela    Sexo     PNA    PNAi     PAA 
##       0       0       0       0       0

#Datos1 <- na.aggregate(B0,FUN=mean) Este funcionaria si no tuviera variables categoricas

###Validar valores atípicos

summary(B0)

##    Escuela              Sexo                PNA              PNAi       
##  Length:2571        Length:2571        Min.   : 72.80   Min.   : 73.30  
##  Class :character   Class :character   1st Qu.: 85.61   1st Qu.: 87.60  
##  Mode  :character   Mode  :character   Median : 91.00   Median : 92.90  
##                                        Mean   : 90.09   Mean   : 91.89  
##                                        3rd Qu.: 95.00   3rd Qu.: 96.90  
##                                        Max.   :110.00   Max.   :100.00  
##       PAA      
##  Min.   :1040  
##  1st Qu.:1260  
##  Median :1340  
##  Mean   :1339  
##  3rd Qu.:1420  
##  Max.   :1590

boxplot(B0$PAA)

boxplot(B0$PNA)

boxplot(B0$PNAi)

Q1=quantile(B0$PNAi,prob=0.25)
Q3=quantile(B0$PNAi,prob=0.75)
#Rango intercuartilico
fs=Q3-Q1 

#observaciones extremas si
#si es mayor a Q3+1.5*fs
#si es menor a Q1-1.5*fs

OLsupPNAi<-which(B0$PNAi>(Q3+1.5*fs))
OLinfPNAi<-which(B0$PNAi<(Q1-1.5*fs))

OLsupPNAi

## integer(0)

OLinfPNAi

## [1] 167

#__________________________________________________________________
Q1=quantile(B0$PAA,prob=0.25)
Q3=quantile(B0$PAA,prob=0.75)
#Rango intercuartilico
fs=Q3-Q1 

#observaciones extremas si
#si es mayor a Q3+1.5*fs
#si es menor a Q1-1.5*fs

OLsupPAA<-which(B0$PAA>(Q3+1.5*fs))
OLinfPAA<-which(B0$PAA<(Q1-1.5*fs))

OLsupPAA

## integer(0)

OLinfPAA

## integer(0)

#__________________________________________________________________
Q1=quantile(B0$PNA,prob=0.25)
Q3=quantile(B0$PNA,prob=0.75)
#Rango intercuartilico
fs=Q3-Q1 

#observaciones extremas si
#si es mayor a Q3+1.5*fs
#si es menor a Q1-1.5*fs

OLsupPNA<-which(B0$PNA>(Q3+1.5*fs))
OLinfPNA<-which(B0$PNA<(Q1-1.5*fs))

OLsupPNA

## [1] 2423

OLinfPNA

## integer(0)

Hay observaciones atípicas en PNA y PNAi

#describe(B0,omit=TRUE)
#t(describe(B0,omit=TRUE)) #t es para trasponer los datos y tener otra vista
round(t(describe(B0,omit=TRUE)),3) #para redondear los valores y que no sea mucho cero

##               PNA     PNAi      PAA
## vars        3.000    4.000    5.000
## n        2571.000 2571.000 2571.000
## mean       90.092   91.891 1339.207
## sd          5.982    5.753   99.540
## median     91.000   92.900 1340.000
## trimmed    90.378   92.311 1339.278
## mad         6.509    6.523  118.608
## min        72.800   73.300 1040.000
## max       110.000  100.000 1590.000
## range      37.200   26.700  550.000
## skew       -0.393   -0.542   -0.022
## kurtosis   -0.574   -0.595   -0.788
## se          0.118    0.113    1.963

Convertir los valores atipicos en NA para su corrección.

B0$PNA[which(B0$PNA>100)] <- NA
B0$PNA <- na.aggregate(B0$PNA,FUN=mean)

summary(B0)

##    Escuela              Sexo                PNA              PNAi       
##  Length:2571        Length:2571        Min.   : 72.80   Min.   : 73.30  
##  Class :character   Class :character   1st Qu.: 85.61   1st Qu.: 87.60  
##  Mode  :character   Mode  :character   Median : 91.00   Median : 92.90  
##                                        Mean   : 90.08   Mean   : 91.89  
##                                        3rd Qu.: 95.00   3rd Qu.: 96.90  
##                                        Max.   :100.00   Max.   :100.00  
##       PAA      
##  Min.   :1040  
##  1st Qu.:1260  
##  Median :1340  
##  Mean   :1339  
##  3rd Qu.:1420  
##  Max.   :1590

Métodos gráficos

####Creación de tablas

#Por escuela
library(dplyr)
Tabla<-B0 %>%
  group_by(Escuela) %>% 
  summarize("n"=n(),
'Porc'=round(100*(n/nrow(B0)),1),
'PNA'=round(mean(PNA,na.rm=TRUE),1),
'PAA'=round(mean(PAA,na.rm=TRUE),0))
Tabla

#pipe(%>%) para la base de datos 0
#los repite por variable categorica
 #cuantos alumnos hay por escuela

#Por Sexo (genero)
library(dplyr)
Tabla1<-B0 %>%
  group_by(Sexo) %>% 
  summarize("n"=n(),
'Porc'=round(100*(n/nrow(B0)),1),
'PNA'=round(mean(PNA,na.rm=TRUE),1),
'PAA'=round(mean(PAA,na.rm=TRUE),0))
Tabla

#Grafico de barras grafico con las n, %de observaciones y donde PAA o del PNA aes(indica qué vas a gráficas ..x=Escuela, y=n) qué tipo de gráfico quieres =geom_col(barras) fill=color geom_text le indico configuraciones para la etiqueta

#Por PNA
ggplot(Tabla,aes(x=Escuela,y=PNA))+
    geom_col(fill = 3)+
  geom_text(aes(label = PNA),vjust=-1)+
    ylim(c(0, 100))+
  scale_x_discrete(limits = c('EIC', 'EN','EAAD','ECSG','EMCS','EHE'))

#Por PAA
ggplot(Tabla,aes(x=Escuela,y=PAA))+
    geom_col(fill = 3)+
  geom_text(aes(label = PAA),vjust=-1)+
    ylim(c(0, 1500))+
  scale_x_discrete(limits = c('EIC', 'EN','EAAD','ECSG','EMCS','EHE'))

#Por PNA
ggplot(Tabla1,aes(x=Sexo,y=PNA))+
    geom_col(fill = 3)+
  geom_text(aes(label = PNA),vjust=-1)+
    ylim(c(0,100))

#Por PAA
ggplot(Tabla1,aes(x=Sexo,y=PAA))+
    geom_col(fill = 3)+
  geom_text(aes(label = PAA),vjust=-1)+
    ylim(c(0,1400))

###Histograma

ggplot(B0,aes(x=PAA))+
  geom_histogram(fill=4,color=1,binwidth = 50)

#bins=30 quiero 30 intervalos

###Grafico de densidad

ggplot(B0, aes(x = PAA)) + 
  geom_density(fill=4, color=1, binwidth = 100,alpha=.3)

## Warning in geom_density(fill = 4, color = 1, binwidth = 100, alpha = 0.3):
## Ignoring unknown parameters: `binwidth`

###Histogramas agrupados

ggplot(B0, aes(x = PAA, fill = Sexo, colour=Sexo)) + 
  geom_histogram(alpha=.4,position="identity", binwidth = 50)

ggplot(B0, aes(x = PAA, fill = Sexo, colour=Sexo)) + 
  geom_density(alpha=.4)

###Caja y bigote

#Sencillo
ggplot(B0, aes(x=PAA)) + 
  geom_boxplot()

#Agrupado
ggplot(B0, aes(x=Sexo, y=PAA)) + 
  geom_boxplot()

#Grafico de caja y bigotes con errores y u observaciones extremas

ggplot(B0, aes(x=Sexo, y=PAA, fill = Sexo)) +
  stat_boxplot(geom = "errorbar",width = 0.25)+
  geom_boxplot()

ggplot(B0, aes(x=Sexo, y=PNA, fill = Sexo)) +
  stat_boxplot(geom = "errorbar",width = 0.25)+
  geom_boxplot()

ggplot(B0, aes(x=Sexo, y=PNAi, fill = Sexo)) +
  stat_boxplot(geom = "errorbar",width = 0.25)+
  geom_boxplot()

##Tablas dinamicas

Tablapivot<-rpivotTable(B0, rows="Escuela", col="Sexo", 
            aggregatorName="Average", vals="PAA",
            rendererName = "Table")
print(Tablapivot)

#ESTADISTICA INFERENCIAL

##PRUEBA DE HIPOTESIS H0: normalidad (asimetria=0y curtosis=3) H1: no normalidad Rechazar H0 a favor de H1 si alpha>=Valor P alpha==0.5

options(scipen = 100, digits = 2)
library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

jarque.bera.test(B0$PAA)

## 
##  Jarque Bera Test
## 
## data:  B0$PAA
## X-squared = 66, df = 2, p-value = 0.000000000000004

Como alfa =0.05 es mayor al valor P se rechaza H0 a favor de H1

No hay normalidad

Como vamos a realizar inferencias comparando hombres y mujeres vs examen de admisión vamos a dividir los datos la base, B0

Mujeres<-filter(B0,Sexo=="M")
Hombres<-filter(B0,Sexo=="H")

Probemos normalidad para hombres y mujeres

jarque.bera.test(Mujeres$PAA)

## 
##  Jarque Bera Test
## 
## data:  Mujeres$PAA
## X-squared = 40, df = 2, p-value = 0.000000002

jarque.bera.test(Hombres$PAA)

## 
##  Jarque Bera Test
## 
## data:  Hombres$PAA
## X-squared = 30, df = 2, p-value = 0.0000004

Como alfa =0.05 es mayor al valor P nos quedamos con la hipótesis de que no hay normalidad para ninguna de las categorias

##Prueba de hipotesis de una muestra para la media

supuesto: normalidad

H0: miu= (cierto valor) media torica del examen de admisión es 1,000 H1: miu=/= 1000 (dos colas)

H1: miu> 1000 (cola superior) para este ejercicio H1: miu< 1000 (cola inferior)

library(stats)
t.test(Mujeres$PAA,alternative="greater",conf.level=0.95,mu=1000)

## 
##  One Sample t-test
## 
## data:  Mujeres$PAA
## t = 125, df = 1372, p-value <0.0000000000000002
## alternative hypothesis: true mean is greater than 1000
## 95 percent confidence interval:
##  1328  Inf
## sample estimates:
## mean of x 
##      1332

t.test(Hombres$PAA,alternative="greater",conf.level=0.95,mu=1000)

## 
##  One Sample t-test
## 
## data:  Hombres$PAA
## t = 120, df = 1197, p-value <0.0000000000000002
## alternative hypothesis: true mean is greater than 1000
## 95 percent confidence interval:
##  1342  Inf
## sample estimates:
## mean of x 
##      1347

Se rechaza H0 a favor de H1, la media es mayor a 1300 puntos

##Pureba de hipotesis NO parametricos de una muestra para la media

H0: miu= (cierto valor) media torica del examen de admisión es 1,000 H1: miu=/= 1000 (dos colas)

H1: miu> 1000 (cola superior) para este ejercicio H1: miu< 1000 (cola inferior)

library("ggpubr")
wilcox.test(Mujeres$PAA, mu = 1000, alternative = "greater")

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Mujeres$PAA
## V = 943251, p-value <0.0000000000000002
## alternative hypothesis: true location is greater than 1000

wilcox.test(Hombres$PAA, mu = 1000, alternative = "greater")

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Hombres$PAA
## V = 718201, p-value <0.0000000000000002
## alternative hypothesis: true location is greater than 1000

Se rechaza h0 a favor de h1, la media es mayor a 1300 puntos.

##Pureba de hipotesis para la difernecia de dos medias con muestras independientes.

Supuestos:

Normalidad en las poblaciones

Muestras sean independientes (que la muestra A es indiependiente de B, en este caso mujer y hombre son independientes) Variazas poblacionales desconocidas pero asumidas iguales

Normalidad no se cumple Muestras independientes entre sí

##Prueba de hipotesis para varianzas

H0: varianza 1 / varianza 2 ==1, (creencia a priori) H1: varianza 1 / varianza 2 =/=1, (creencia a priori)

var.test(x = Hombres$PAA, y = Mujeres$PAA, null.value=1, alternative="two.sided",
                conf.level=0.95)

## 
##  F test to compare two variances
## 
## data:  Hombres$PAA and Mujeres$PAA
## F = 1, num df = 1197, denom df = 1372, p-value = 0.6
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.92 1.15
## sample estimates:
## ratio of variances 
##                  1

Como alfa es menor al valor p, nos queamos con H0, las varianzas son iguales

Normalidad no se cumple(por teorema de limite central si, pero probaremos pruebas no parametricas) Muestras independientes sí se cumple Varianza constantes sí

Ho: MediaH-MediaM=0, (Δ=0 creencia a priori) Ha: MediaH-MediaM>0 (prueba cola superir)

t.test(x = Hombres$PAA, y = Mujeres$PAA,
alternative="greater",
           mu=0, paired=FALSE, var.equal=TRUE, conf.level = 0.95)

## 
##  Two Sample t-test
## 
## data:  Hombres$PAA and Mujeres$PAA
## t = 4, df = 2569, p-value = 0.00009
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  8.2 Inf
## sample estimates:
## mean of x mean of y 
##      1347      1332

Como alfa > valor P, se rechaza H0 a favor de H1 al PAA promedio de los hombres es mayor al de las mujeres

###Pureba de hipotesisi no parametrica para comparar dos medias Supuesto no requiere, solo que las muestras sean independientes.

wilcox.test(x=Hombres$PAA,y=Mujeres$PAA,
            alternative = "greater",paired = FALSE,conf.level = 0.95)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Hombres$PAA and Mujeres$PAA
## W = 885586, p-value = 0.0004
## alternative hypothesis: true location shift is greater than 0

Misma conclusión parametrica o no parametricas

Como alfa > valor P, se rechaza H0 a favor de H1 al PAA promedio de los hombres es mayor al de las mujeres

Análisis de Varianza (ANOVA)

Estudia el efecto de una variable categórica (con más de una categória) sobre una variable numérica.

y=PAA(Puntaje examen de admisión) x=Escuela (escuela de inscripción)

H0: m1=m2=m3=m4=m5=m6 = m H1: al menos una es diferente ~ = “en función a”

Estimacion<-aov(PAA ~ Escuela,B0)
summary(Estimacion)

##               Df   Sum Sq Mean Sq F value              Pr(>F)    
## Escuela        5  1044939  208988    21.9 <0.0000000000000002 ***
## Residuals   2565 24419012    9520                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: como alfa es mayor al valor p, se rechaza H0, al menos una media es diferente.

Validación de problemas en la estimación y supuestos

Observaciones extremas e influyentes

Residual = Observado - Estimado …Obtenemos la t para cada error

Observaciones extremas

B0$RS<-rstudent(Estimacion)

Obs_extremas<-which(abs(B0$RS) >= 2)
Obs_extremas

##   30   88  180  222  251  277  327  332  341  347  420  444  445  455  468  510 
##   30   88  180  222  251  277  327  332  341  347  420  444  445  455  468  510 
##  529  558  625  638  640  641  720  776  803  849  965  975 1001 1026 1036 1162 
##  529  558  625  638  640  641  720  776  803  849  965  975 1001 1026 1036 1162 
## 1176 1229 1240 1249 1257 1334 1393 1474 1507 1594 1655 1662 1688 1717 1725 1735 
## 1176 1229 1240 1249 1257 1334 1393 1474 1507 1594 1655 1662 1688 1717 1725 1735 
## 1799 1812 1829 1847 1853 1861 1874 1878 2084 2143 2282 2286 2293 2313 2387 2412 
## 1799 1812 1829 1847 1853 1861 1874 1878 2084 2143 2282 2286 2293 2313 2387 2412 
## 2433 2470 2479 2489 2522 2536 2554 2561 
## 2433 2470 2479 2489 2522 2536 2554 2561

Observaciones influyentes

B0$DC<-cooks.distance(Estimacion)

#Valor F de tablas
Fc<-qf(0.5, (6+1), (nrow(B0)-(6+1)))

#k = 6 escuelas

Obs_influyentes<-which(B0$DC > Fc)
Obs_influyentes

## named integer(0)

No hay influyentes

Quitar observaciones influyentes

B0<-B0[-Obs_extremas,]

Supuesto de normalidad

H0: Normalidad

H1: No normalida d

jarque.bera.test(resid(Estimacion))

## 
##  Jarque Bera Test
## 
## data:  resid(Estimacion)
## X-squared = 49, df = 2, p-value = 0.00000000002

No hay normalidad.

Alternativa

1.-Transformación Box-cox (busca que la variable tenga una distribución normal, con una potencia que asegure normalidad)

2.-Prueba no parametrica de ANOVA

3.-Teorema del limite central (tengo una muestra grande)

Transformación Box-Cox lambda= potencia

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

BC<-boxcox(Estimacion,lambda=seq(-100,100,.01))

lambda<-BC$x[which.max(BC$y)]
lambda

## [1] 1.5

El valor de lambda de 1.5 maximiza la probabilidad de que tengamos normalidad

with= transformar variables

B0$PAAT<-with(B0,PAA^lambda)

EstimacionT<-aov(PAAT ~ Escuela,B0)
summary(EstimacionT)

##               Df      Sum Sq    Mean Sq F value              Pr(>F)    
## Escuela        5  5107285276 1021457055    34.5 <0.0000000000000002 ***
## Residuals   2493 73879406429   29634740                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

jarque.bera.test(resid(EstimacionT))

## 
##  Jarque Bera Test
## 
## data:  resid(EstimacionT)
## X-squared = 72, df = 2, p-value = 0.0000000000000002

No se lora corregir la normalidad, tenemos dos opciones más

ANOVA no parametrica
Teorema del limite central

Varianza Constante

H0: Varianza Constante

H1: Varianza no constante

Estimacion<- aov(PAA ~ Escuela, B0)
summary(Estimacion)

##               Df   Sum Sq Mean Sq F value              Pr(>F)    
## Escuela        5  1434540  286908      34 <0.0000000000000002 ***
## Residuals   2493 21056071    8446                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(lmtest)
bptest(Estimacion)

## 
##  studentized Breusch-Pagan test
## 
## data:  Estimacion
## BP = 5, df = 5, p-value = 0.4

como alfa < valor p, nos quedamos con Ho la varianza es constante

Comparación de medias

Debido a que nos quedamos con que las medias son parametricas ahora hay que probar cuales son firerentes a iguales.

Tukey<-TukeyHSD(Estimacion)
Tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = PAA ~ Escuela, data = B0)
## 
## $Escuela
##             diff    lwr   upr p adj
## ECSG-EAAD  -0.93  -20.0  18.1  1.00
## EHE-EAAD    4.81  -16.1  25.8  0.99
## EIC-EAAD   26.88   10.4  43.3  0.00
## EMCS-EAAD  72.50   51.3  93.7  0.00
## EN-EAAD    -9.52  -26.3   7.3  0.59
## EHE-ECSG    5.74  -16.0  27.4  0.97
## EIC-ECSG   27.81   10.4  45.2  0.00
## EMCS-ECSG  73.43   51.5  95.4  0.00
## EN-ECSG    -8.59  -26.3   9.1  0.74
## EIC-EHE    22.08    2.6  41.5  0.02
## EMCS-EHE   67.70   44.1  91.3  0.00
## EN-EHE    -14.33  -34.1   5.4  0.30
## EMCS-EIC   45.62   25.9  65.4  0.00
## EN-EIC    -36.40  -51.3 -21.5  0.00
## EN-EMCS   -82.02 -102.0 -62.0  0.00

plot(Tukey,las=2)

Donde alfa > al valor p ajustados, la diferencia en las medias será estadisticamente diferente de cero.

Medicina Seguida de Ingenieria

Cualquier intervalo que no cruce cero, tiene media mayor

Prueba de hipotesis para no parametricas para comparar más de dos medias

supuestos, solo requiere que las muestras sean indepenties

KT<-kruskal.test(PAA~Escuela,B0)
KT

## 
##  Kruskal-Wallis rank sum test
## 
## data:  PAA by Escuela
## Kruskal-Wallis chi-squared = 156, df = 5, p-value <0.0000000000000002

Como el ANOVA, las medias son diferentes alfa > valor p

Prueba para decir cuales medias son diferentes alternativa de tukey para no parametricas

B0$Escuela<-as.factor(B0$Escuela)

library(PMCMRplus)
kwAllPairsConoverTest(PAA~Escuela,B0, p.adjust.method ="holm")

## Warning in kwAllPairsConoverTest.default(c(1240, 1370, 1390, 1370, 1180, : Ties
## are present. Quantiles were corrected for ties.

## 
##  Pairwise comparisons using Conover's all-pairs test

## data: PAA by Escuela

##      EAAD                 ECSG                 EHE              
## ECSG 1.000                -                    -                
## EHE  1.000                1.000                -                
## EIC  0.000073089730191    0.000073089730191    0.012            
## EMCS < 0.0000000000000002 < 0.0000000000000002 0.000000000000011
## EN   0.462                0.688                0.231            
##      EIC               EMCS                
## ECSG -                 -                   
## EHE  -                 -                   
## EIC  -                 -                   
## EMCS 0.000000000509954 -                   
## EN   0.000000000114425 < 0.0000000000000002

## 
## P value adjustment method: holm

en la tabla viene el valor p de la comparación de medias , si alfa > valor p se concluye que las medias son diferentes

Terminan las pruebas de la estadística inferencial _______________________________________________________________________________________

Miércoles_Sesió1

Miriam Alonso

2023-06-28