Pregunta 2: Presencias un debate entre colegas. Uno (colega A) planteas que la mortalidad materna en un país puede ser explicada por el nivel de tubercolisis; mientras que el otro (colega B) piensa que es mejor explicada por el nivel de fertilidad. En ambos casos, controlan el nivel de empleo. ¿Si usan los mismos datos, a quién apoyarías con sólo ver ambas tablas de coeficientes? Sustenta tu respuesta.

library(rio)
data1 = import("wbDataMini (2).xlsx")
str(data1)
## 'data.frame':    263 obs. of  6 variables:
##  $ pais              : chr  "MAR" "ABW" "AFG" "AGO" ...
##  $ TasaFertil1mMuje  : num  32.3 24.3 73.1 157.4 20.7 ...
##  $ EmployPerPop      : num  22.4 NA 16.2 69.7 39.5 ...
##  $ Tuberculosis100m  : num  101 14 189 370 16 5.9 NA 0.79 25 47 ...
##  $ MaterMort100m     : num  121 NA 396 477 29 NA 156 6 52 25 ...
##  $ UndernourishPerPop: num  3.5 NA 23 14 4.9 ...

VD: mortalidad materna VI: Tuberculosis VI: Fertilidad VC: enmpleo

Modelo 1: Uno (colega A) planteas que la mortalidad materna en un país puede ser explicada por el nivel de tubercolisis. Mas variable de control: EmployPerPop.

modelo1=formula(MaterMort100m~Tuberculosis100m+EmployPerPop)
reg1=lm(modelo1,data=data1)
summary(reg1)
## 
## Call:
## lm(formula = modelo1, data = data1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -520.80  -99.71  -52.99   35.21 1028.94 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -28.05760   43.65145  -0.643  0.52113    
## Tuberculosis100m   0.75078    0.09557   7.856 2.56e-13 ***
## EmployPerPop       2.33807    0.88920   2.629  0.00923 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 194.2 on 196 degrees of freedom
##   (64 observations deleted due to missingness)
## Multiple R-squared:  0.2864, Adjusted R-squared:  0.2791 
## F-statistic: 39.33 on 2 and 196 DF,  p-value: 4.351e-15
library(modelsummary)
model1=list('apropiacion (I)'=reg1)
modelsummary(model1, title = "Regresion: modelo 1",
             stars = TRUE,
             output = "kableExtra")
Regresion: modelo 1
 apropiacion (I)
(Intercept) −28.058
(43.651)
Tuberculosis100m 0.751***
(0.096)
EmployPerPop 2.338**
(0.889)
Num.Obs. 199
R2 0.286
R2 Adj. 0.279
AIC 2666.8
BIC 2679.9
Log.Lik. −1329.381
RMSE 192.75
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Modelo 2: el otro (colega B) piensa que es mejor explicada por el nivel de fertilidad. Más variable de control: EmployPerPop.

modelo2=formula(MaterMort100m~TasaFertil1mMuje+EmployPerPop)
reg2=lm(modelo2,data=data1)
summary(reg2)
## 
## Call:
## lm(formula = modelo2, data = data1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -320.64  -71.02   -7.28   52.01  876.75 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -89.8159    31.1505  -2.883  0.00432 ** 
## TasaFertil1mMuje   4.3535     0.2488  17.497  < 2e-16 ***
## EmployPerPop       1.0524     0.6383   1.649  0.10057    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 142.5 on 223 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.6037, Adjusted R-squared:  0.6001 
## F-statistic: 169.9 on 2 and 223 DF,  p-value: < 2.2e-16
library(modelsummary)
model2=list('apropiacion (I)'=reg2)
modelsummary(model2, title = "Regresion: modelo 2",
             stars = TRUE,
             output = "kableExtra")
Regresion: modelo 2
 apropiacion (I)
(Intercept) −89.816**
(31.150)
TasaFertil1mMuje 4.353***
(0.249)
EmployPerPop 1.052
(0.638)
Num.Obs. 226
R2 0.604
R2 Adj. 0.600
AIC 2887.9
BIC 2901.6
Log.Lik. −1439.959
RMSE 141.54
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Comparación de modelos:

models=list('apropiacion (I)'=reg1,
            'apropiacion (II)'=reg2)
modelsummary(models, title = "Resultados de los modelos II y III",
             stars = TRUE,
             output = "kableExtra")
Resultados de los modelos II y III
 apropiacion (I)  apropiacion (II)
(Intercept) −28.058 −89.816**
(43.651) (31.150)
Tuberculosis100m 0.751***
(0.096)
EmployPerPop 2.338** 1.052
(0.889) (0.638)
TasaFertil1mMuje 4.353***
(0.249)
Num.Obs. 199 226
R2 0.286 0.604
R2 Adj. 0.279 0.600
AIC 2666.8 2887.9
BIC 2679.9 2901.6
Log.Lik. −1329.381 −1439.959
RMSE 192.75 141.54
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Respuesta: Apoyaría al modelo número dos, pues presenta un R2 mayor al del modelo uno. El modelo 1 presenta un R2 de 0.286; mientras que el modelo 2 presenta un R2 de 0.604. Entonces apoyaría al colega B.

Pregunta 3:

Se está estudiando qué factores influyen en que los convictos puestos en libertad vuelvan a la carcel. Los datos puede descargarlos de aqui.

Con esos datos, dos investigadores desean probar algunas hipotesis. El investigador “A” afirma que mientras mayor la edad, y si uno ya está casado, disminuirá el riesgo de volver a la carcel; mientras que el investigador “B” asegura que, mientras menor la edad, y si uno ya está casado disminuye el riesgo de volver a la carcel. Ud, al usar R, corroborá que….

link="https://docs.google.com/spreadsheets/d/e/2PACX-1vQSlGaMI8Q8qlXI0Bp3m7BQcEh8ZLzaP7RymVtRYkg3ah1sZVlCi6-HmeKCic1RjfuH3gL_wrbMms88/pub?gid=1573532387&single=true&output=csv"
carcel=read.csv(link, stringsAsFactors = T)
str(carcel)
## 'data.frame':    432 obs. of  10 variables:
##  $ semanasLibre       : int  20 17 25 52 52 52 23 52 52 52 ...
##  $ fueArrestado       : int  1 1 1 0 0 0 1 0 0 0 ...
##  $ tuvoApoyoDinero    : int  0 0 0 1 0 0 0 1 0 0 ...
##  $ edad               : int  27 18 19 23 19 24 25 21 22 20 ...
##  $ esNegro            : int  1 1 0 1 0 1 1 1 1 1 ...
##  $ expLaboralPrevia   : int  0 0 1 1 1 1 1 1 0 1 ...
##  $ casado             : int  0 0 0 1 0 0 1 0 0 0 ...
##  $ libertadCondicional: int  1 1 1 1 1 0 1 1 0 0 ...
##  $ vecesEnCarcel      : int  3 8 13 1 3 2 0 4 6 0 ...
##  $ nivelEduca         : int  2 3 2 4 2 3 3 2 2 4 ...
carcel[,c(2,3,5,6,7,8)]=lapply(carcel[,c(2,3,5,6,7,8)], as.factor)
carcel$nivelEduca=as.ordered(carcel$nivelEduca)
str(carcel)
## 'data.frame':    432 obs. of  10 variables:
##  $ semanasLibre       : int  20 17 25 52 52 52 23 52 52 52 ...
##  $ fueArrestado       : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 2 1 1 1 ...
##  $ tuvoApoyoDinero    : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 2 1 1 ...
##  $ edad               : int  27 18 19 23 19 24 25 21 22 20 ...
##  $ esNegro            : Factor w/ 2 levels "0","1": 2 2 1 2 1 2 2 2 2 2 ...
##  $ expLaboralPrevia   : Factor w/ 2 levels "0","1": 1 1 2 2 2 2 2 2 1 2 ...
##  $ casado             : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 2 1 1 1 ...
##  $ libertadCondicional: Factor w/ 2 levels "0","1": 2 2 2 2 2 1 2 2 1 1 ...
##  $ vecesEnCarcel      : int  3 8 13 1 3 2 0 4 6 0 ...
##  $ nivelEduca         : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 2 4 2 3 3 2 2 4 ...
summary(carcel$semanasLibre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   50.00   52.00   45.85   52.00   52.00
table(carcel$fueArrestado)
## 
##   0   1 
## 318 114

H1: El investigador “A” afirma que mientras mayor la edad, y si uno ya está casado, disminuirá el riesgo de volver a la carceL

H2: El investigador “B” asegura que, mientras menor la edad, y si uno ya está casado disminuye el riesgo de volver a la carcel

library(survival)
# note que necesito el factor como numérico
carcel$survival=with(carcel,Surv(time = semanasLibre,event =  as.numeric(fueArrestado)))
# que es:

library(magrittr) # needed for pipe %>% 
carcel%>%
    rmarkdown::paged_table()
COX_H1= formula(survival~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)

#regression
rcox1 <- coxph(COX_H1,data=carcel)
modelcox=list('Riesgo - Re arrestado'=rcox1,'Riesgo- Re arrestado (exponenciado)'=rcox1)

#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
modelsummary(modelcox,
             #fmt=f,
             exponentiate = c(F,T), 
             statistic = 'conf.int',
             title = "Regresión Cox",
             stars = TRUE,
             output = "kableExtra")
Regresión Cox
Riesgo - Re arrestado Riesgo- Re arrestado (exponenciado)
tuvoApoyoDinero1 −0.415* 0.660*
[−0.789, −0.042] [0.454, 0.959]
nivelEduca.L −0.558 0.573
[−1.976, 0.861] [0.139, 2.366]
nivelEduca.Q −0.694 0.500
[−1.907, 0.519] [0.149, 1.680]
nivelEduca.C 0.396 1.486
[−0.506, 1.298] [0.603, 3.662]
nivelEduca^4 −0.008 0.992
[−0.582, 0.566] [0.559, 1.762]
vecesEnCarcel 0.087** 1.091**
[0.032, 0.142] [1.032, 1.152]
Num.Obs. 432 432
AIC 1338.1 1338.1
BIC 1362.6 1362.6
RMSE 0.51 0.51
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Respuesta: El investigador A está en lo cierto.

Pregunta 3: Se está estudiando los factores que influencian la admisión a los colegios elite del Perú. Para ello se tienen datos de los postulantes procedentes de diferentes escuelas primarias. Use la data de este LINK (note que el documento incluye la metadata).

Si usas todas las variables disponibles como predictoras de la admisión, explica cada una de las variables en la tabla de resultados. Haz los cálculos que necesites para interpretar los coeficientes obtenidos, de ser necesario.

library(rio)
data3 = import("https://docs.google.com/spreadsheets/d/1qrRFEAB-jZBUUjn9XwEm8mi-l0l4mJnCWUPX6aXtpz4/edit#gid=0")
str(data3)
## 'data.frame':    400 obs. of  4 variables:
##  $ admitido : int  0 1 1 1 0 1 1 0 1 0 ...
##  $ letras   : int  380 660 800 640 520 760 560 400 540 700 ...
##  $ ciencias : int  361 367 400 319 293 300 298 308 339 392 ...
##  $ prestigio: int  2 2 4 1 1 3 4 3 2 3 ...
### first hypothesis
h1=formula(admitido~letras)

#regression
rlog1=glm(h1, data=data3,family = binomial)
modelrl=list('Ser Admitido (I)'=rlog1)

library(modelsummary)
modelsummary(modelrl,
             title = "Regresión Logística",
             stars = TRUE,
             output = "kableExtra")
Regresión Logística
 Ser Admitido (I)
(Intercept) −2.901***
(0.606)
letras 0.004***
(0.001)
Num.Obs. 400
AIC 490.1
BIC 498.0
Log.Lik. −243.028
F 13.199
RMSE 0.46
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Respuesta: Ser admitido en letras es significativo para el modelo.