##Haga Click Aqui para ver Certificado Machine Learning MIT https://www.credential.net/4dd365ea-ea5a-46a2-a72e-539e70545c6e

##Haga Click Aqui para ver Certificado Columbia Python for Managers https://certificates.emeritus.org/0a2e1de7-add2-4710-ad49-417d1dadfb61#gs.4a92hv

Machine Learning Aplicado/ Columbia University: https://certificates.emeritus.org/4f62c75a-2e43-48f8-bfdd-586ccfc5acec#gs.dp16fk

Algunos Dashboards elaborados son:

Para el Estado del Clima https://rchang.shinyapps.io/rchang-app_clima_ho/

Para Machine Learning https://rchang.shinyapps.io/rchang-app/

Para Empresariales e Industriales https://rchang.shinyapps.io/rchang-app_final_emp/

Prueba de hipótesis

En este notebook se muestran las funciones que hay disponibles en R para realizar prueba de hipótesis para:

  1. la media μ,

  2. la proporción p,

  3. la varianza σ2,

  4. la diferencia de medias,μ1−μ2 para muestras independientes y dependientes (o pareadas),

  5. la diferencia de proporciones p1− p2 , y

  6. la razón de varianzas σ21 / σ22

  7. Prueba de hipótesis para μ de una población normal

Para realizar este tipo de prueba se puede usar la función t.test que tiene la siguiente estructura.

t.test(x, y = NULL, alternative = c(“two.sided”, “less”, “greater”), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …) Los argumentos a definir dentro de t.test para hacer la prueba son:

x: vector numérico con los datos. alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la hipótesis alterna es
≠ , “less” para el caso
< y “greater” para
> . mu: valor de referencia de la prueba. conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).

Ejemplo

Para verificar si el proceso de llenado de bolsas de café con 500 gramos está operando correctamente se toman aleatoriamente muestras de tamaño diez cada cuatro horas. Una muestra de bolsas está compuesta por las siguientes observaciones: 502, 501, 497, 491, 496, 501, 502, 500, 489, 490.

¿Está el proceso llenando bolsas conforme lo dice la envoltura? Use un nivel de significancia del 5%.

Solución

Lo primero es explorar si la muestra proviene de una distribución normal, para eso ingresamos los datos y aplicamos la prueba Anderson-Darling por medio de la función ad.test disponible en el paquete nortest (Gross and Ligges 2015) como se muestra a continuación.

cafe <- c(510, 492, 494, 498, 492, 496, 502, 491, 507, 496)

require(nortest) # Se debe haber instalado antes nortest ad.test(cafe) ## ## Anderson-Darling normality test ## ## data: contenido ## A = 0.49161, p-value = 0.1665

#install.packages("nortest")
cafe <- c(510, 492, 494, 498, 492,
               496, 502, 491, 507, 496)
library(nortest)
ad.test(cafe)
## 
##  Anderson-Darling normality test
## 
## data:  cafe
## A = 0.49161, p-value = 0.1665
shapiro.test(cafe)
## 
##  Shapiro-Wilk normality test
## 
## data:  cafe
## W = 0.88468, p-value = 0.1476

Como el p-value es > 0.05, se puede asumir que la muestra proviene de una población normal.

Luego de haber explorado la normalidad retornamos al problema de interés que se puede resumir así:

H0: μ = 500gr

H1: μ ≠ 500gr

La prueba de hipótesis se puede realizar usando la función t.test por medio del siguiente código.

t.test(cafe, alternative='two.sided',
       conf.level=0.95, mu=500)
## 
##  One Sample t-test
## 
## data:  cafe
## t = -1.0629, df = 9, p-value = 0.3155
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.1176 502.4824
## sample estimates:
## mean of x 
##     497.8

Ejemplo 2 Proporciones

Un fabricante de un quitamanchas afirma que su producto quita 90% de todas las manchas. Para poner a prueba esta afirmación se toman 200 camisetas manchadas de las cuales a solo 174 les desapareció la mancha. Pruebe la afirmación del fabricante a un nivel
α = 0.05

Solución

En este problema interesa probar lo siguiente:

H0: p = 0.90

H1: p < 0.90

Del anterior conjunto de hipótesis se observa que el valor de referencia de la prueba es
p0 = 0.90 De la información inicial se tiene que de las n = 200 pruebas se observó que en
x = 174 la mancha desapareció, con esta información se puede calcular el estadístico

z así:

z <- (174/200 - 0.90) / sqrt(0.90 * (1 - 0.90) / 200)
z  # Para obtener el valor del estadístico
## [1] -1.414214

Para obtener el valor-P de la prueba debemos tener en cuenta el sentido en la hipótesis alternativa H1: p < 0.90, por esa razón el valor-P será P (Z < z) y para obtenerlo usamos el siguiente código

pnorm(q=z, lower.tail=TRUE)  # Para obtener el valor-P
## [1] 0.0786496

Como el valor-P obtenido fue mayor que el nivel de significancia α = 0.05 se concluye que no hay evidencias suficientes para rechazar la hipótesis nula.

Prueba Chi 2 de Pearson Para realizar la prueba Chi2 de Pearson se usa la función prop.test que tiene la siguiente estructura.

prop.test(x, n, p = NULL, alternative = c(“two.sided”, “less”, “greater”), conf.level = 0.95, correct = TRUE) Los argumentos a definir dentro de prop.test para hacer la prueba son:

x: número de éxitos en la muestra. n: número de observaciones en la muestra. alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la alterna es
≠ , “less” para el caso < y “greater” para > . p: valor de referencia de la prueba. correct: valor lógico para indicar si se usa la corrección. conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).

Ejemplo

Un fabricante de un quitamanchas afirma que su producto quita 90% de todas las manchas. Para poner a prueba esta afirmación se toman 200 camisetas manchadas de las cuales a solo 174 les desapareció la mancha. Pruebe la afirmación del fabricante a un nivel
α = 0.05

Solución

En este problema interesa probar lo siguiente:

H0 : p = 0.90

H1 : p < 0.90

La forma de usar la función prop.test para realizar la prueba se muestra a continuación.

prop.test(x=174, n=200, p=0.9, alternative='less',
          conf.level=0.95, correct=FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  174 out of 200, null probability 0.9
## X-squared = 2, df = 1, p-value = 0.07865
## alternative hypothesis: true p is less than 0.9
## 95 percent confidence interval:
##  0.0000000 0.9042273
## sample estimates:
##    p 
## 0.87

Prueba de hipótesis para la diferencia de medias μ1−μ2 con varianzas iguales

Para realizar este tipo de prueba se puede usar la función t.test que tiene la siguiente estructura.

t.test(x, y = NULL, alternative = c(“two.sided”, “less”, “greater”), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …)

Ejemplo:

T1 presentan la productividad de una linea de producción de una empresa antes de la inversión en maquinaria y equipo, y T2 presenta la reformulación del flujo de trabajo con la inversión en maquinaria y equipo. Un directivo menciona que no hay diferencias entre T1 y T2

T1 <- c(76, 85, 74,78, 82, 75, 82) 
T2 <- c(57, 67, 55, 64, 61, 63, 63)
bartlett.test(list(T1, T2))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(T1, T2)
## Bartlett's K-squared = 0.00016629, df = 1, p-value = 0.9897
datos <- data.frame(tiempo=c(T1, T2), trat=rep(1:2, each=7))
boxplot(tiempo ~ trat, data=datos, las=1,
        xlab='Tratamiento', ylab='Tiempo (min)')

este problema interesa estudiar el siguiente conjunto de hipótesis.

H0 : μ1 − μ2 = 0 H1 : μ1 − μ2 ≠ 0

t.test(x=T1, y=T2, alternative="two.sided", mu=0, 
       paired=FALSE, var.equal=TRUE, conf.level=0.97)
## 
##  Two Sample t-test
## 
## data:  T1 and T2
## t = 7.8209, df = 12, p-value = 4.737e-06
## alternative hypothesis: true difference in means is not equal to 0
## 97 percent confidence interval:
##  11.94503 22.91212
## sample estimates:
## mean of x mean of y 
##  78.85714  61.42857

Prueba de hipótesis para la diferencia de medias pareadas

Ejemplo

Diez empleados fueron seleccionados para participar en un programa de capacitación empresarial orientado a incrementar las ventas corporativas.

Las ventas de los empleados medidas antes y después de haber participado en el programa y los datos en lempiras aparecen abajo. ¿Hay evidencia que soporten la afirmación de el programa aumentó las ventas? Usar nivel de significancia del 5%.

Empleado 001 002 003 004 005 006 007 008 009 010 Después 195 213 247 201 187 210 215 246 294 310 Antes 187 195 221 190 175 197 199 221 278 285

despu <- c(195, 213, 247, 201, 187, 210, 215, 246, 294, 310)
antes <- c(187, 195, 221, 190, 175, 197, 199, 221, 278, 285)
dif <- antes - despu
qqnorm(dif, pch=19, main='')
qqline(dif)

Se puede aplicar la prueba de normalidad Kolmogorov-Smirnov para estudiar si las diferencias dif provienen de una población normal, esto se puede realizar por medio del siguiente código.

require(nortest) lillie.test(dif)

require(nortest)
lillie.test(dif)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dif
## D = 0.19393, p-value = 0.3552
t.test(x=despu, y=antes, alternative="greater", mu=0, 
       paired=TRUE, conf.level=0.95)
## 
##  Paired t-test
## 
## data:  despu and antes
## t = 8.3843, df = 9, p-value = 7.593e-06
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  13.2832     Inf
## sample estimates:
## mean of the differences 
##                      17

Prueba de hipótesis para la diferencia de proporciones
p1 − p2

Para realizar pruebas de hipótesis para la proporción se usa la función prop.test.

mplo Se quiere determinar si un cambio en el método de fabricación de una piezas ha sido efectivo o no. Para esta comparación se tomaron 2 muestras, una antes y otra después del cambio en el proceso y los resultados obtenidos son los siguientes.

Num piezas Antes Después Defectuosas 75 80 Analizadas 1500 2000 Realizar una prueba de hipótesis con un nivel de significancia del 10%.

En este problema interesa estudiar el siguiente conjunto de hipótesis.

H0 : pantes -p despues = 0 H1 : pantes -p despues > 0

Para realizar la prueba se usa la función prop.test como se muestra a continuación.

prop.test(x=c(75, 80), n=c(1500, 2000),
          alternative='greater', conf.level=0.90)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(75, 80) out of c(1500, 2000)
## X-squared = 1.7958, df = 1, p-value = 0.09011
## alternative hypothesis: greater
## 90 percent confidence interval:
##  0.0002765293 1.0000000000
## sample estimates:
## prop 1 prop 2 
##   0.05   0.04
library(datarium)
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata

y nota del examen

x nota de la tarea

y<-c(95,80,0,0,79,77,72,66,98,90,0,95,35,50,72,55,75,66)
x<-c(96,77,0,0,78,64,89,47,90,93,18,86,0,30,59,77,74,67)
n<-18

Queremos explicar la nota del examen en terminos de la nota de la tarea

cor(y,x) # coeficiente de correlacion de pearson
## [1] 0.91041

modelos lineales con la función lm

modelo<-lm(y ~ x)

summary(modelo)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.4345  -8.8437   0.3528   9.6466  24.2731 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.72691    6.61733   1.621    0.125    
## x            0.87265    0.09914   8.802 1.57e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.85 on 16 degrees of freedom
## Multiple R-squared:  0.8288, Adjusted R-squared:  0.8181 
## F-statistic: 77.48 on 1 and 16 DF,  p-value: 1.571e-07
modelo2<-lm(y ~ 0 + x)

#Y = Beta X #Y = 1.01242 X

summary(modelo2)
## 
## Call:
## lm(formula = y ~ 0 + x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.956  -2.102   0.056  11.137  35.000 
## 
## Coefficients:
##   Estimate Std. Error t value Pr(>|t|)    
## x  1.01242    0.05121   19.77 3.62e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.5 on 17 degrees of freedom
## Multiple R-squared:  0.9583, Adjusted R-squared:  0.9559 
## F-statistic: 390.8 on 1 and 17 DF,  p-value: 3.617e-13
modelo<-lm(y ~ x)
summary(modelo)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.4345  -8.8437   0.3528   9.6466  24.2731 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.72691    6.61733   1.621    0.125    
## x            0.87265    0.09914   8.802 1.57e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.85 on 16 degrees of freedom
## Multiple R-squared:  0.8288, Adjusted R-squared:  0.8181 
## F-statistic: 77.48 on 1 and 16 DF,  p-value: 1.571e-07

Ho: Beta_i=0: H1: Beta_i=0:

plot(x,y)
abline(a=modelo$coefficients[1],b=modelo$coefficients[2])

confint(modelo)
##                  2.5 %    97.5 %
## (Intercept) -3.3012022 24.755021
## x            0.6624861  1.082807

Evaluación del modelo

res <- residuals( modelo ) # residuos 
pre <- predict(modelo) #predicciones

BIC(modelo)
## [1] 152.2632
AIC(modelo)
## [1] 149.5921

como saber cual modelo se ajusta mejor

R cuadrado y Rcuadrado ajustado - mas grande mejor (>70)

AIC y BIC - mas pequeño mejor

diagnóstico del modelo

una vez tenemos el modelo ajustado procedemos con su

#diagnostico, que se realiza a traves de analisis de residuos

homogeneidad de varianzas: homocedasticidad

##library(“lmtest”)

#install.packages("lmtest")

#H0:Existe homogeneidad en las varianzas #H1:No existe homogeneidad en las varianzas

Test de Breusch-Pagan

library("lmtest")
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 3.9853, df = 1, p-value = 0.0459

no existe homogeneidad en las varianzas

test de normalidad en los residuos

residuos<-residuals(modelo)
hist(residuos) # histograma de los residuos estandarizados 

boxplot(residuos) # diagrama de cajas de los residuos estandarizados 

qqnorm(residuos) # gráfico de cuantiles de los residuos estandarizados 
qqline(residuos)  

shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96484, p-value = 0.697

Autocorrelacion o independecia en los residuos

#supuesto de no correlacion # test de Durbin-Watson

#H0: No hay autocorrelacion #H1: Si hay autocorrelacion

dwtest(modelo, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.9597, p-value = 0.9112
## alternative hypothesis: true autocorrelation is not 0

regresion lineal multiple/ modelo saturado

marketing<-datarium::marketing
modelo <- lm(sales ~ youtube + facebook + newspaper, data = marketing)
summary(modelo)
## 
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5932  -1.0690   0.2902   1.4272   3.3951 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.526667   0.374290   9.422   <2e-16 ***
## youtube      0.045765   0.001395  32.809   <2e-16 ***
## facebook     0.188530   0.008611  21.893   <2e-16 ***
## newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.023 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16
cor(marketing)
##              youtube   facebook  newspaper     sales
## youtube   1.00000000 0.05480866 0.05664787 0.7822244
## facebook  0.05480866 1.00000000 0.35410375 0.5762226
## newspaper 0.05664787 0.35410375 1.00000000 0.2282990
## sales     0.78222442 0.57622257 0.22829903 1.0000000
plot(marketing)

modelo2  <- lm(sales ~ youtube + facebook, data = marketing)

summary(modelo2)
## 
## Call:
## lm(formula = sales ~ youtube + facebook, data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5572  -1.0502   0.2906   1.4049   3.3994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.50532    0.35339   9.919   <2e-16 ***
## youtube      0.04575    0.00139  32.909   <2e-16 ***
## facebook     0.18799    0.00804  23.382   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8962 
## F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16
confint(modelo2)
##                  2.5 %     97.5 %
## (Intercept) 2.80841159 4.20222820
## youtube     0.04301292 0.04849671
## facebook    0.17213877 0.20384969

Regresión Logística

Concepto de ODDS o razón de probabilidad, ratio de ODDS y logaritmo de ODDS

En la regresión lineal simple, se modela el valor de la variable dependiente Y en función del valor de la variable independiente X. Sin embargo, en la regresión logística, tal como se ha descrito en la sección anterior, se modela la probabilidad de que la variable respuesta Y pertenezca al nivel de referencia 1 en función del valor que adquieran los predictores, mediante el uso de LOG of ODDs.

Supóngase que la probabilidad de que un evento sea verdadero es de 0.8, por lo que la probabilidad de evento falso es de 1 - 0.8 = 0.2. Los ODDs o razón de probabilidad de verdadero se definen como el ratio entre la probabilidad de evento verdadero y la probabilidad de evento falso pq.

En este caso los ODDs de verdadero son 0.8 / 0.2 = 4, lo que equivale a decir que se esperan 4 eventos verdaderos por cada evento falso.

La trasformación de probabilidades a ODDs es monotónica, si la probabilidad aumenta también lo hacen los ODDs, y viceversa. El rango de valores que pueden tomar los ODDs es de \[0,∞\]. Dado que el valor de una probabilidad está acotado entre \[0,1\] se recurre a una trasformación logit (existen otras) que consiste en el logaritmo natural de los ODDs. Esto permite convertir el rango de probabilidad previamente limitado a \[0,1\] a \[−∞,+∞\].

Condiciones

Independencia: las observaciones tienen que ser independientes unas de otras.

Relación lineal entre el logaritmo natural de odds y la variable continua: patrones en forma de U son una clara violación de esta condición.

La regresión logística no precisa de una distribución normal de la variable continua independiente.

Número de observaciones: no existe una norma establecida al respecto, pero se recomienda entre 50 a 100 observaciones.

Regresion logistica simple

#El modelo logit es un modelo de regresión típico, Y=f(X+ε) #en el que la variable respuesta (variable aleatoria Y) es dicotómica o #binaria (toma dos valores: 0 y 1), habitualmente sobre si el individuo tiene #una característica (1) o no (0), y las variables predictivas (vector aleatorio X) #son continuas. El modelo logit es un caso particular de los llamados modelos lineales #generalizados (GLMs, Generalized Linear Model).

#Un estudio quiere establecer un modelo que permita calcular la probabilidad de obtener una mencion #de honor al final del bachillerato en función de la nota que se ha obtenido en matemáticas. #La variable matrícula está codificada como 0 si no se tiene mencion y 1 si se tiene.

notas <- read.csv("../entradas/notas.csv")
notas
##       X matricula matematicas
## 1     1         0          41
## 2     2         0          53
## 3     3         0          54
## 4     4         0          47
## 5     5         0          57
## 6     6         0          51
## 7     7         0          42
## 8     8         0          45
## 9     9         0          54
## 10   10         0          52
## 11   11         0          51
## 12   12         1          51
## 13   13         0          71
## 14   14         1          57
## 15   15         0          50
## 16   16         0          43
## 17   17         0          51
## 18   18         0          60
## 19   19         1          62
## 20   20         0          57
## 21   21         0          35
## 22   22         1          75
## 23   23         0          45
## 24   24         0          57
## 25   25         0          45
## 26   26         0          46
## 27   27         1          66
## 28   28         0          57
## 29   29         0          49
## 30   30         0          49
## 31   31         0          57
## 32   32         0          64
## 33   33         1          63
## 34   34         0          57
## 35   35         0          50
## 36   36         1          58
## 37   37         0          75
## 38   38         1          68
## 39   39         0          44
## 40   40         0          40
## 41   41         0          41
## 42   42         0          62
## 43   43         0          57
## 44   44         0          43
## 45   45         1          48
## 46   46         0          63
## 47   47         0          39
## 48   48         0          70
## 49   49         0          63
## 50   50         0          59
## 51   51         1          61
## 52   52         0          38
## 53   53         0          61
## 54   54         0          49
## 55   55         1          73
## 56   56         0          44
## 57   57         0          42
## 58   58         0          39
## 59   59         0          55
## 60   60         0          52
## 61   61         0          45
## 62   62         1          61
## 63   63         0          39
## 64   64         0          41
## 65   65         0          50
## 66   66         0          40
## 67   67         0          60
## 68   68         0          47
## 69   69         0          59
## 70   70         0          49
## 71   71         0          46
## 72   72         0          58
## 73   73         1          71
## 74   74         0          58
## 75   75         0          46
## 76   76         0          43
## 77   77         1          54
## 78   78         0          56
## 79   79         0          46
## 80   80         0          54
## 81   81         0          57
## 82   82         0          54
## 83   83         0          71
## 84   84         1          48
## 85   85         0          40
## 86   86         1          64
## 87   87         0          51
## 88   88         0          39
## 89   89         0          40
## 90   90         0          61
## 91   91         1          66
## 92   92         0          49
## 93   93         1          65
## 94   94         0          52
## 95   95         0          46
## 96   96         1          61
## 97   97         1          72
## 98   98         1          71
## 99   99         0          40
## 100 100         1          69
## 101 101         0          64
## 102 102         0          56
## 103 103         0          49
## 104 104         0          54
## 105 105         0          53
## 106 106         0          66
## 107 107         1          67
## 108 108         0          40
## 109 109         0          46
## 110 110         1          69
## 111 111         0          40
## 112 112         0          41
## 113 113         0          57
## 114 114         1          58
## 115 115         1          57
## 116 116         0          37
## 117 117         0          55
## 118 118         1          62
## 119 119         0          64
## 120 120         0          40
## 121 121         0          50
## 122 122         0          46
## 123 123         0          53
## 124 124         0          52
## 125 125         1          45
## 126 126         0          56
## 127 127         0          45
## 128 128         0          54
## 129 129         0          56
## 130 130         0          41
## 131 131         0          54
## 132 132         1          72
## 133 133         1          56
## 134 134         0          47
## 135 135         0          49
## 136 136         1          60
## 137 137         0          54
## 138 138         0          55
## 139 139         0          33
## 140 140         0          49
## 141 141         0          43
## 142 142         0          50
## 143 143         0          52
## 144 144         0          48
## 145 145         0          58
## 146 146         0          43
## 147 147         1          41
## 148 148         0          43
## 149 149         0          46
## 150 150         0          44
## 151 151         0          43
## 152 152         0          61
## 153 153         0          40
## 154 154         0          49
## 155 155         1          56
## 156 156         0          61
## 157 157         0          50
## 158 158         0          51
## 159 159         0          42
## 160 160         1          67
## 161 161         1          53
## 162 162         0          50
## 163 163         1          51
## 164 164         1          72
## 165 165         0          48
## 166 166         0          40
## 167 167         0          53
## 168 168         0          39
## 169 169         1          63
## 170 170         0          51
## 171 171         0          45
## 172 172         0          39
## 173 173         0          42
## 174 174         0          62
## 175 175         0          44
## 176 176         0          65
## 177 177         1          63
## 178 178         0          54
## 179 179         0          45
## 180 180         1          60
## 181 181         1          49
## 182 182         0          48
## 183 183         1          57
## 184 184         1          55
## 185 185         1          66
## 186 186         1          64
## 187 187         0          55
## 188 188         0          42
## 189 189         1          56
## 190 190         0          53
## 191 191         0          41
## 192 192         0          42
## 193 193         0          53
## 194 194         0          42
## 195 195         1          60
## 196 196         0          52
## 197 197         0          38
## 198 198         0          57
## 199 199         1          58
## 200 200         1          65
table(notas$matricula)
## 
##   0   1 
## 151  49
boxplot(matematicas~matricula,notas)

modelo <- glm(matricula ~ matematicas, data = notas, family = "binomial")
summary(modelo)
## 
## Call:
## glm(formula = matricula ~ matematicas, family = "binomial", data = notas)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0332  -0.6785  -0.3506  -0.1565   2.6143  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.79394    1.48174  -6.610 3.85e-11 ***
## matematicas  0.15634    0.02561   6.105 1.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 222.71  on 199  degrees of freedom
## Residual deviance: 167.07  on 198  degrees of freedom
## AIC: 171.07
## 
## Number of Fisher Scoring iterations: 5
confint(object = modelo, level = 0.95 )
## Waiting for profiling to be done...
##                   2.5 %     97.5 %
## (Intercept) -12.9375208 -7.0938806
## matematicas   0.1093783  0.2103937

Tenga en cuenta que la salida del modelo logístico está en escala de enlace (logit); por lo tanto, la salida numérica del modelo corresponde a las probabilidades logarítmicas. Por conveniencia, anotemos el modelo de regresión logística con la intersección y el coeficiente calculados.

P(X) = probabilidad de que sea matemáticas con mención honorífica

exp(coefficients(modelo)[2])
## matematicas 
##    1.169224

P(X) = probabilidad de que sea matemáticas y que no sea mención honorífica

exp(coefficients(modelo)[1])
##  (Intercept) 
## 5.578854e-05
confint.default(modelo)[2,]
##     2.5 %    97.5 % 
## 0.1061467 0.2065340
exp(-9.79394)
## [1] 5.578866e-05

Acorde al modelo, el logaritmo de los odds de que un estudiante tenga matrícula está positivamente relacionado con la puntuación obtenida en matemáticas (coeficiente de regresión = 0.1563404).

Esto significa que, por cada unidad que se incrementa la variable matemáticas, se espera que el logaritmo de odds de la variable matrícula se incremente en promedio 0.1563404 unidades.

Aplicando la inversa del logaritmo natural e (0.1563404)=1.169) se obtiene que, por cada unidad que se incrementa la variable matemáticas, los odds de obtener matrícula se incremente en promedio 1.169 unidades. No hay que confundir esto último con que la probabilidad de matrícula se incremente un 1.169 %.

exp(0.15634)
## [1] 1.169224

A diferencia de la regresión lineal en la que β1 se corresponde con el cambio promedio en la variable dependiente Y debido al incremento en una unidad del predictor X, en regresión logística, β1 indica el cambio en el logaritmo de odds debido al incremento de una unidad de X, o lo que es lo mismo, multiplica los odds por eβ1. Dado que la relación entre p(Y) y X no es lineal, β1 no se corresponde con el cambio en la probabilidad de Y asociado con el incremento de una unidad de X. Cuánto se incremente la probabilidad de Y por unidad de X depende del valor de X, es decir, de la posición en la curva logística en la que se encuentre.

Además del valor de las estimaciones de los coeficientes parciales de correlación del modelo, es conveniente calcular sus correspondientes intervalos de confianza.

confint(object = modelo, level = 0.95 )
## Waiting for profiling to be done...
##                   2.5 %     97.5 %
## (Intercept) -12.9375208 -7.0938806
## matematicas   0.1093783  0.2103937
plot(matricula ~ matematicas, notas, col = "darkblue",
     main = "Modelo regresión logística",
     ylab = "P(matrícula=1|matemáticas)",
     xlab = "matemáticas", pch = "I")
curve(predict(modelo, data.frame(matematicas = x), type = "response"),
      col = "blue", lwd = 2.5, add = TRUE)

#si alguien tiene una nota en matematica de 40, cual es la #probabilidad de tener una meción honorífica

P(X) = e (Bo + B1 * X) / 1 +P(X) = e (Bo + B1 * X)

notaMat<-70

exp(-9.79394+0.15634*(notaMat))/(1+exp(-9.79394+0.15634*(notaMat)))*100
## [1] 75.94853

#regresion logistica multiple

#Un estudio considera que existe relación entre el hecho de que un estudiante asista a clases de repaso de lectura (sí = 1, no = 0), #la nota que obtiene en un examen de lectura estándar (realizado antes de iniciar las clases de repaso) #y el sexo (hombre = 1, mujer = 0). Se quiere generar un modelo en el que a partir de las variables puntuación del examen y sexo, #prediga la probabilidad de que el estudiante tenga que asistir a clases de repaso.

clases<-read.csv("../entradas/clases.csv")
clases
##       X   sexo examen_lectura clases_repaso
## 1     1 hombre           91.0             0
## 2     2 hombre           77.5             0
## 3     3  mujer           52.5             0
## 4     4  mujer           54.0             0
## 5     5  mujer           53.5             0
## 6     6 hombre           62.0             0
## 7     7  mujer           59.0             0
## 8     8 hombre           51.5             0
## 9     9  mujer           61.5             0
## 10   10  mujer           56.5             0
## 11   11 hombre           47.5             0
## 12   12 hombre           75.0             0
## 13   13 hombre           47.5             0
## 14   14 hombre           53.5             0
## 15   15  mujer           50.0             0
## 16   16  mujer           50.0             0
## 17   17 hombre           49.0             0
## 18   18  mujer           59.0             0
## 19   19 hombre           60.0             0
## 20   20  mujer           60.0             0
## 21   21 hombre           60.5             0
## 22   22  mujer           50.0             0
## 23   23  mujer          101.0             0
## 24   24 hombre           60.0             0
## 25   25 hombre           60.0             0
## 26   26  mujer           83.5             0
## 27   27  mujer           61.0             0
## 28   28  mujer           75.0             0
## 29   29 hombre           84.0             0
## 30   30 hombre           56.5             0
## 31   31 hombre           56.5             0
## 32   32  mujer           45.0             0
## 33   33 hombre           60.5             0
## 34   34  mujer           77.5             0
## 35   35 hombre           62.5             0
## 36   36  mujer           70.0             0
## 37   37  mujer           69.0             0
## 38   38  mujer           62.0             0
## 39   39  mujer          107.5             0
## 40   40  mujer           54.5             0
## 41   41 hombre           92.5             0
## 42   42  mujer           94.5             0
## 43   43 hombre           65.0             0
## 44   44  mujer           80.0             0
## 45   45  mujer           45.0             0
## 46   46  mujer           45.0             0
## 47   47  mujer           66.0             0
## 48   48 hombre           66.0             0
## 49   49  mujer           57.5             0
## 50   50 hombre           42.5             0
## 51   51  mujer           60.0             0
## 52   52 hombre           64.0             0
## 53   53  mujer           65.0             0
## 54   54  mujer           47.5             0
## 55   55 hombre           57.5             0
## 56   56 hombre           55.0             0
## 57   57 hombre           55.0             0
## 58   58 hombre           76.5             0
## 59   59 hombre           51.5             0
## 60   60 hombre           59.5             0
## 61   61 hombre           59.5             0
## 62   62 hombre           59.5             0
## 63   63 hombre           55.0             0
## 64   64  mujer           70.0             0
## 65   65 hombre           66.5             0
## 66   66 hombre           84.5             0
## 67   67 hombre           57.5             0
## 68   68 hombre          125.0             0
## 69   69  mujer           70.5             0
## 70   70 hombre           79.0             0
## 71   71  mujer           56.0             0
## 72   72 hombre           75.0             0
## 73   73 hombre           57.5             0
## 74   74 hombre           56.0             0
## 75   75  mujer           67.5             0
## 76   76 hombre          114.5             0
## 77   77  mujer           70.0             0
## 78   78  mujer           67.0             0
## 79   79 hombre           60.5             0
## 80   80  mujer           95.0             0
## 81   81  mujer           65.5             0
## 82   82  mujer           85.0             0
## 83   83 hombre           55.0             0
## 84   84 hombre           63.5             0
## 85   85 hombre           61.5             0
## 86   86 hombre           60.0             0
## 87   87 hombre           52.5             0
## 88   88  mujer           65.0             0
## 89   89  mujer           87.5             0
## 90   90  mujer           62.5             0
## 91   91  mujer           66.5             0
## 92   92 hombre           67.0             0
## 93   93  mujer          117.5             0
## 94   94  mujer           47.5             0
## 95   95  mujer           67.5             0
## 96   96  mujer           67.5             0
## 97   97  mujer           77.0             0
## 98   98  mujer           73.5             0
## 99   99  mujer           73.5             0
## 100 100  mujer           68.5             0
## 101 101  mujer           55.0             0
## 102 102  mujer           92.0             0
## 103 103 hombre           55.0             0
## 104 104  mujer           55.0             0
## 105 105 hombre           60.0             0
## 106 106 hombre          120.5             0
## 107 107  mujer           56.0             0
## 108 108  mujer           84.5             0
## 109 109  mujer           60.0             0
## 110 110 hombre           85.0             0
## 111 111  mujer           93.0             0
## 112 112 hombre           60.0             0
## 113 113  mujer           65.0             0
## 114 114  mujer           58.5             0
## 115 115  mujer           85.0             0
## 116 116 hombre           67.0             0
## 117 117  mujer           67.5             0
## 118 118 hombre           65.0             0
## 119 119  mujer           60.0             0
## 120 120 hombre           47.5             0
## 121 121  mujer           79.0             0
## 122 122 hombre           80.0             0
## 123 123  mujer           57.5             0
## 124 124  mujer           64.5             0
## 125 125  mujer           65.0             0
## 126 126  mujer           60.0             0
## 127 127  mujer           85.0             0
## 128 128  mujer           60.0             0
## 129 129  mujer           58.0             0
## 130 130  mujer           61.5             0
## 131 131 hombre           60.0             1
## 132 132  mujer           65.0             1
## 133 133 hombre           93.5             1
## 134 134 hombre           52.5             1
## 135 135 hombre           42.5             1
## 136 136 hombre           75.0             1
## 137 137 hombre           48.5             1
## 138 138 hombre           64.0             1
## 139 139 hombre           66.0             1
## 140 140  mujer           82.5             1
## 141 141  mujer           52.5             1
## 142 142  mujer           45.5             1
## 143 143 hombre           57.5             1
## 144 144 hombre           65.0             1
## 145 145  mujer           46.0             1
## 146 146  mujer           75.0             1
## 147 147 hombre          100.0             1
## 148 148  mujer           77.5             1
## 149 149 hombre           51.5             1
## 150 150 hombre           62.5             1
## 151 151 hombre           44.5             1
## 152 152  mujer           51.0             1
## 153 153  mujer           56.0             1
## 154 154  mujer           58.5             1
## 155 155  mujer           69.0             1
## 156 156 hombre           65.0             1
## 157 157 hombre           60.0             1
## 158 158  mujer           65.0             1
## 159 159 hombre           65.0             1
## 160 160 hombre           40.0             1
## 161 161  mujer           55.0             1
## 162 162 hombre           52.5             1
## 163 163 hombre           54.5             1
## 164 164 hombre           74.0             1
## 165 165 hombre           55.0             1
## 166 166  mujer           60.5             1
## 167 167 hombre           50.0             1
## 168 168 hombre           48.0             1
## 169 169  mujer           51.0             1
## 170 170  mujer           55.0             1
## 171 171 hombre           93.5             1
## 172 172 hombre           61.0             1
## 173 173 hombre           52.5             1
## 174 174 hombre           57.5             1
## 175 175 hombre           60.0             1
## 176 176  mujer           71.0             1
## 177 177  mujer           65.0             1
## 178 178  mujer           60.0             1
## 179 179  mujer           55.0             1
## 180 180 hombre           60.0             1
## 181 181 hombre           77.0             1
## 182 182 hombre           52.5             1
## 183 183  mujer           95.0             1
## 184 184 hombre           50.0             1
## 185 185  mujer           47.5             1
## 186 186 hombre           50.0             1
## 187 187 hombre           47.0             1
## 188 188 hombre           71.0             1
## 189 189  mujer           65.0             1

#Las tablas de frecuencia así, como representaciones gráficas de las observaciones son útiles para intuir #si las variables independientes escogidas están relacionadas con la variable respuesta y por lo tanto #ser buenos predictores.

tabla <- table(clases$clases_repaso, clases$sexo,
               dnn = c("clases de repaso","sexo"))
tabla
##                 sexo
## clases de repaso hombre mujer
##                0     57    73
##                1     36    23
addmargins(tabla)
##                 sexo
## clases de repaso hombre mujer Sum
##              0       57    73 130
##              1       36    23  59
##              Sum     93    96 189
addmargins(prop.table(tabla))
##                 sexo
## clases de repaso    hombre     mujer       Sum
##              0   0.3015873 0.3862434 0.6878307
##              1   0.1904762 0.1216931 0.3121693
##              Sum 0.4920635 0.5079365 1.0000000
mosaicplot(tabla)

library("ggplot2")
## Warning: package 'ggplot2' was built under R version 4.1.3
ggplot(clases, aes(x=factor(clases_repaso),y= examen_lectura, fill=factor(sexo))) +
  geom_boxplot()

modelo <- glm(clases_repaso ~ examen_lectura + sexo, data = clases,
              family = "binomial")
summary(modelo)
## 
## Call:
## glm(formula = clases_repaso ~ examen_lectura + sexo, family = "binomial", 
##     data = clases)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2079  -0.8954  -0.7243   1.2592   2.0412  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     1.18365    0.78559   1.507   0.1319  
## examen_lectura -0.02617    0.01223  -2.139   0.0324 *
## sexomujer      -0.64749    0.32484  -1.993   0.0462 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 224.64  on 186  degrees of freedom
## AIC: 230.64
## 
## Number of Fisher Scoring iterations: 4

Comparación de las predicciones con las observaciones

#Para este estudio se va a emplear un umbral de 0.5. Si la probabilidad predicha de asistir a #clases de repaso es superior a 0.5 se asigna al nivel 1 (sí asiste), si es menor se asigna al #nivel 0 (no clases de repaso).

predicciones <- ifelse(modelo$fitted.values > 0.5, 1, 0)

matriz_confusion <- table(modelo$model$clases_repaso, predicciones,
                          dnn = c("observaciones", "predicciones"))

addmargins(matriz_confusion)
##              predicciones
## observaciones   0   1 Sum
##           0   129   1 130
##           1    56   3  59
##           Sum 185   4 189
mosaicplot(matriz_confusion)

Ccorrecta<-(129+3)/189*100
Ccorrecta
## [1] 69.84127

##Cómo explica el modelo en que porcentaje

Ccorrecta<-(3551+12)/3921*100
Ccorrecta
## [1] 90.86968