EXAMEN PARCIAL
2023-10-13
PREGUNTA 1 –> REGRESION COX
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….
library(readxl)
dataCarcel <- read_excel("dataCarcel.xlsx")
View(dataCarcel)
str(dataCarcel)## tibble [432 × 10] (S3: tbl_df/tbl/data.frame)
## $ semanasLibre : num [1:432] 20 17 25 52 52 52 23 52 52 52 ...
## $ fueArrestado : num [1:432] 1 1 1 0 0 0 1 0 0 0 ...
## $ tuvoApoyoDinero : num [1:432] 0 0 0 1 0 0 0 1 0 0 ...
## $ edad : num [1:432] 27 18 19 23 19 24 25 21 22 20 ...
## $ esNegro : num [1:432] 1 1 0 1 0 1 1 1 1 1 ...
## $ expLaboralPrevia : num [1:432] 0 0 1 1 1 1 1 1 0 1 ...
## $ casado : num [1:432] 0 0 0 1 0 0 1 0 0 0 ...
## $ libertadCondicional: num [1:432] 1 1 1 1 1 0 1 1 0 0 ...
## $ vecesEnCarcel : num [1:432] 3 8 13 1 3 2 0 4 6 0 ...
## $ nivelEduca : num [1:432] 2 3 2 4 2 3 3 2 2 4 ...dataCarcel[,c(2,3,5,6,7,8)]=lapply(dataCarcel[,c(2,3,5,6,7,8)], as.factor)
dataCarcel$nivelEduca=as.ordered(dataCarcel$nivelEduca)
str(dataCarcel)## tibble [432 × 10] (S3: tbl_df/tbl/data.frame)
## $ semanasLibre : num [1:432] 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 : num [1:432] 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 : num [1:432] 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 ...Creacion de la variable survival
library(survival)
dataCarcel$survival=with(dataCarcel,Surv(time = semanasLibre,event = as.numeric(fueArrestado)))
# que es:
library(magrittr) # needed for pipe %>%
dataCarcel%>%
rmarkdown::paged_table()Analisis Kaplan Meier
library(ggplot2)
library(ggfortify)
#aqui el generico
KM.generico = survfit(survival ~ 1, data = dataCarcel)
###graficando
ejeX='SEMANAS\n curva cae cuando alguien es arrestado'
ejeY='Probabilidad \n(PERMANECER LIBRE)'
titulo="Curva de Sobrevivencia: permanecer libre"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)
HIPOTESIS 1 –> LO QUE DICE EL INVESTIGADOR A
El investigador “A” afirma que mientras mayor la edad, y si uno ya está casado, disminuirá el riesgo de volver a la carcel
COX_H1= formula(survival~edad+casado)
#regression
rcox1 <- coxph(COX_H1,data=dataCarcel)
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")| Riesgo - Re arrestado | Riesgo- Re arrestado (exponenciado) | |
|---|---|---|
| edad | -0.067** | 0.935** |
| [-0.108, -0.026] | [0.898, 0.974] | |
| casado1 | -0.493 | 0.611 |
| [-1.224, 0.237] | [0.294, 1.268] | |
| Num.Obs. | 432 | 432 |
| AIC | 1337.5 | 1337.5 |
| BIC | 1345.6 | 1345.6 |
| RMSE | 0.51 | 0.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
HIPOTESIS 2 –> LO QUE DICE EL INVESTIGADOR B
“B” asegura que, mientras menor la edad, y si uno ya está casado disminuye el riesgo de volver a la carcel
#Teniendo en cuenta que ambos investigadores poseen hipotesis con las mismas variables, podriamos decir, con el primer cuadro, que si bien la edad es una variable significativa al 0.05 y que influye en el riesgo de ser arrestado, la variable casado no posee significancia en el modelo, por lo que ninguno de los investigadores habia realizado una hipotesis acertada.
PREGUNTA 2 –> REGRESION LOGISTICA
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(readxl)
admision <- read_excel("admision.xlsx")
View(admision)
str(admision)## tibble [400 × 4] (S3: tbl_df/tbl/data.frame)
## $ admitido : num [1:400] 0 1 1 1 0 1 1 0 1 0 ...
## $ letras : num [1:400] 380 660 800 640 520 760 560 400 540 700 ...
## $ ciencias : num [1:400] 361 367 400 319 293 300 298 308 339 392 ...
## $ prestigio: num [1:400] 2 2 4 1 1 3 4 3 2 3 ...admision$admitido=as.factor(admision$admitido)
admision$admitido=factor(admision$admitido,levels=levels(admision$admitido) ,labels=c("No","Si"))
admision$prestigio=as.factor(admision$prestigio)
admision$prestigio=factor(admision$prestigio,levels=levels(admision$prestigio),labels=c("1","2", "3", "4"), ordered = T)
str(admision)## tibble [400 × 4] (S3: tbl_df/tbl/data.frame)
## $ admitido : Factor w/ 2 levels "No","Si": 1 2 2 2 1 2 2 1 2 1 ...
## $ letras : num [1:400] 380 660 800 640 520 760 560 400 540 700 ...
## $ ciencias : num [1:400] 361 367 400 319 293 300 298 308 339 392 ...
## $ prestigio: Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 2 2 4 1 1 3 4 3 2 3 ...set.seed(1)
modeloA = glm(admitido ~ letras, family = binomial, data = admision)
modeloB = glm(admitido ~ letras + ciencias, family = binomial, data = admision)
modelo1= glm(admitido ~ letras + ciencias + prestigio, family = binomial, data = admision)
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(modeloA, modeloB,modelo1, type = "text")##
## ===============================================
## Dependent variable:
## -----------------------------
## admitido
## (1) (2) (3)
## -----------------------------------------------
## letras 0.004*** 0.003** 0.002**
## (0.001) (0.001) (0.001)
##
## ciencias 0.008** 0.008**
## (0.003) (0.003)
##
## prestigio.L 1.189***
## (0.287)
##
## prestigio.Q 0.232
## (0.252)
##
## prestigio.C -0.099
## (0.212)
##
## Constant -2.901*** -4.949*** -4.882***
## (0.606) (1.075) (1.113)
##
## -----------------------------------------------
## Observations 400 400 400
## Log Likelihood -243.028 -240.172 -229.259
## Akaike Inf. Crit. 490.056 486.344 470.517
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01library(modelsummary)
modelsummary(modelo1,
exponentiate = T,
statistic = 'conf.int',
title = "Regresión Logística (Coeficientes Exponenciados)",
stars = TRUE,
output = "kableExtra")| (1) | |
|---|---|
| (Intercept) | 0.008*** |
| [0.001, 0.064] | |
| letras | 1.002* |
| [1.000, 1.004] | |
| ciencias | 1.008* |
| [1.002, 1.015] | |
| prestigio.L | 3.285*** |
| [1.897, 5.879] | |
| prestigio.Q | 1.261 |
| [0.764, 2.057] | |
| prestigio.C | 0.906 |
| [0.595, 1.369] | |
| Num.Obs. | 400 |
| AIC | 470.5 |
| BIC | 494.5 |
| Log.Lik. | -229.259 |
| F | 7.228 |
| RMSE | 0.44 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
#Se puede decir que todas las variables independientes poseen significancia. En el caso de las variables letras y ciencias, se obtiene una significancia al 0.1 y en el caso de la variable prestigio, posee una significancia al 0.01. Hay 3 coeficientes mayores a 1, lo que nos indica que podria existir un efecto positivo y que, ademas, todas son significativas para el modelo.
PREGUNTA 3 –> REGRESION GAUSSIANA
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(readxl)
madres <- read_excel("wbDataMini.xlsx")
View(madres)
str(madres)## tibble [263 × 6] (S3: tbl_df/tbl/data.frame)
## $ pais : chr [1:263] "MAR" "ABW" "AFG" "AGO" ...
## $ TasaFertil1mMuje : num [1:263] 32.3 24.3 73.1 157.4 20.7 ...
## $ EmployPerPop : num [1:263] 22.4 NA 16.2 69.7 39.5 ...
## $ Tuberculosis100m : num [1:263] 101 14 189 370 16 5.9 NA 0.79 25 47 ...
## $ MaterMort100m : num [1:263] 121 NA 396 477 29 NA 156 6 52 25 ...
## $ UndernourishPerPop: num [1:263] 3.5 NA 23 14 4.9 ...Nos quedamos solo con los casos completos
madres= madres[complete.cases(madres),]HIPOTESIS 1
h1=formula(MaterMort100m~ Tuberculosis100m + EmployPerPop)
regre1 = lm(h1, data = madres)HIPOTESIS 2
h2=formula(MaterMort100m~ TasaFertil1mMuje + EmployPerPop)
regre2 = lm(h2, data = madres)library(stargazer)
stargazer(regre1, regre2, type = "text")##
## ===========================================================
## Dependent variable:
## ----------------------------
## MaterMort100m
## (1) (2)
## -----------------------------------------------------------
## Tuberculosis100m 0.735***
## (0.094)
##
## TasaFertil1mMuje 4.314***
## (0.261)
##
## EmployPerPop 2.187** 0.412
## (0.927) (0.689)
##
## Constant -29.995 -70.858**
## (45.622) (33.366)
##
## -----------------------------------------------------------
## Observations 182 182
## R2 0.293 0.626
## Adjusted R2 0.286 0.621
## Residual Std. Error (df = 179) 185.875 135.301
## F Statistic (df = 2; 179) 37.176*** 149.575***
## ===========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Comparacion de modelos
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(magrittr)
library(knitr)
tanova=anova(regre1,regre2)
tanova## Analysis of Variance Table
##
## Model 1: MaterMort100m ~ Tuberculosis100m + EmployPerPop
## Model 2: MaterMort100m ~ TasaFertil1mMuje + EmployPerPop
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 179 6184330
## 2 179 3276819 0 2907512kable(tanova,
caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 179 | 6184330 | NA | NA | NA | NA |
| 179 | 3276819 | 0 | 2907512 | NA | NA |
```