rm(list = ls())

#Pregunta 1: Logistica

library(rio)
data1 = import('admision.xlsx')
str(data1)
## 'data.frame':    400 obs. of  4 variables:
##  $ admitido : num  0 1 1 1 0 1 1 0 1 0 ...
##  $ letras   : num  380 660 800 640 520 760 560 400 540 700 ...
##  $ ciencias : num  361 367 400 319 293 300 298 308 339 392 ...
##  $ prestigio: num  2 2 4 1 1 3 4 3 2 3 ...
data1$admitido=as.factor(data1$admitido)
data1$admitido=factor(data1$admitido,levels=levels(data1$admitido) ,labels=c("No","Si"))
data1$prestigio=as.factor(data1$prestigio)
data1$prestigio=factor(data1$prestigio,levels=levels(data1$prestigio) ,labels=c("1","2","3","4"))
str(data1)
## 'data.frame':    400 obs. of  4 variables:
##  $ admitido : Factor w/ 2 levels "No","Si": 1 2 2 2 1 2 2 1 2 1 ...
##  $ letras   : num  380 660 800 640 520 760 560 400 540 700 ...
##  $ ciencias : num  361 367 400 319 293 300 298 308 339 392 ...
##  $ prestigio: Factor w/ 4 levels "1","2","3","4": 2 2 4 1 1 3 4 3 2 3 ...
h1 = formula(admitido ~ letras + ciencias + prestigio )
rlog1 =glm(h1, data=data1,family = binomial)
modelsrl=list('Ser admitido (I)'=rlog1)
library(modelsummary)
modelsummary(modelsrl,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "Regresión Logística (Exponenciados)",
             stars = TRUE,
             output = "kableExtra")
Regresión Logística (Exponenciados)
 Ser admitido (I)
(Intercept) 0.004***
[0.000, 0.035]
letras 1.002*
[1.000, 1.004]
ciencias 1.008*
[1.002, 1.015]
prestigio2 1.235
[0.582, 2.743]
prestigio3 2.401*
[1.201, 5.107]
prestigio4 4.718***
[2.125, 11.023]
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

#Pregunta 2: Gauss

data2 = import('limpia1.2.xlsx')
## New names:
## • `` -> `...1`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
data2=data2[complete.cases(data2),]
str(data2)
## 'data.frame':    196 obs. of  9 variables:
##  $ ...1                                                     : chr  "Provincia" "Provincia" "Provincia" "Provincia" ...
##  $ Perú: datos para el planeamiento estratégico, provincial.: chr  "010200" "010300" "010100" "010400" ...
##  $ ...3                                                     : chr  "BAGUA" "BONGARÁ" "CHACHAPOYAS" "CONDORCANQUI" ...
##  $ ...4                                                     : chr  "77637.0" "29552.0" "53637.0" "46520.0" ...
##  $ ...5                                                     : chr  "82193.0" "27085.0" "60419.0" "49800.0" ...
##  $ ...6                                                     : chr  "98110.0" "29335.0" "58773.0" "75888.0" ...
##  $ ...7                                                     : chr  "0.519479274749756" "0.509355187416077" "0.39617508649826" "0.690789043903351" ...
##  $ ...8                                                     : chr  "0.649202856777717" "0.716390890220421" "0.767278998632085" "0.408828965702124" ...
##  $ ...9                                                     : chr  "1050.07668881867" "1288.0695953639" "5186.63721759991" "596.421671278727" ...

#pregunta 3

data3 = import('wbDataMini (1).xlsx')
str(data3)
## '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 ...
modeloa=formula(MaterMort100m ~Tuberculosis100m )
reg3=lm(modeloa,data=data3)
summary(reg3)
## 
## Call:
## lm(formula = modeloa, data = data3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -551.14  -88.25  -68.68   21.30 1048.92 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      78.95383   17.80788   4.434 1.53e-05 ***
## Tuberculosis100m  0.75612    0.09361   8.077 6.15e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 197.9 on 200 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.246,  Adjusted R-squared:  0.2422 
## F-statistic: 65.24 on 1 and 200 DF,  p-value: 6.145e-14
modelob=formula(MaterMort100m ~ TasaFertil1mMuje )
reg4=lm(modelob,data=data3)
summary(reg4)
## 
## Call:
## lm(formula = modelob, data = data3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -324.14  -67.34   -0.36   49.59  879.17 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -44.3893    15.1407  -2.932  0.00371 ** 
## TasaFertil1mMuje   4.4385     0.2411  18.411  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 142.3 on 227 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.5989, Adjusted R-squared:  0.5972 
## F-statistic:   339 on 1 and 227 DF,  p-value: < 2.2e-16

estandarizamos

library(modelsummary)

modelo3<-lm(MaterMort100m~Tuberculosis100m + TasaFertil1mMuje, data3)

modelo3=formula(scale(MaterMort100m)~scale(Tuberculosis100m)+scale(TasaFertil1mMuje))
reg3a=lm(modelo3,data3)
model3=list(reg3a)
modelsummary(model3, title = "Regresion: modelo con \ncoeficientes estandarizados",
             stars = TRUE,
             output = "kableExtra")
Regresion: modelo con coeficientes estandarizados
 (1)
(Intercept) -0.045
(0.045)
scale(Tuberculosis100m) 0.210***
(0.048)
scale(TasaFertil1mMuje) 0.655***
(0.047)
Num.Obs. 202
R2 0.615
R2 Adj. 0.611
AIC 393.4
BIC 406.6
Log.Lik. -192.690
RMSE 0.63
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

#Pregunta 4 Cox

carcel = import('dataCarcel.xlsx')
str(carcel)
## 'data.frame':    432 obs. of  10 variables:
##  $ semanasLibre       : num  20 17 25 52 52 52 23 52 52 52 ...
##  $ fueArrestado       : num  1 1 1 0 0 0 1 0 0 0 ...
##  $ tuvoApoyoDinero    : num  0 0 0 1 0 0 0 1 0 0 ...
##  $ edad               : num  27 18 19 23 19 24 25 21 22 20 ...
##  $ esNegro            : num  1 1 0 1 0 1 1 1 1 1 ...
##  $ expLaboralPrevia   : num  0 0 1 1 1 1 1 1 0 1 ...
##  $ casado             : num  0 0 0 1 0 0 1 0 0 0 ...
##  $ libertadCondicional: num  1 1 1 1 1 0 1 1 0 0 ...
##  $ vecesEnCarcel      : num  3 8 13 1 3 2 0 4 6 0 ...
##  $ nivelEduca         : num  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)
str(carcel)
## 'data.frame':    432 obs. of  10 variables:
##  $ semanasLibre       : num  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  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  3 8 13 1 3 2 0 4 6 0 ...
##  $ nivelEduca         : num  2 3 2 4 2 3 3 2 2 4 ...
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()
library(ggplot2)
library(ggfortify)

#aqui el generico
KM.generico = survfit(survival ~ 1, data = carcel)

###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)

KM_H1=formula(survival ~ casado)

KM.fondos = survfit(KM_H1, data = carcel)

###
ejeX='SEMANAS\n curva cae cuando alguien es arrestado'
ejeY="Prob ('seguir libre')"
titulo="Curva de Sobrevivencia: ¿Beneficia el estar casado?"

autoplot(KM.fondos,xlab=ejeX,ylab=ejeY, 
         main = titulo,conf.int = F)  + 
        labs(colour = "casado") + 
         scale_color_discrete(labels = c("No", "Sí"))

LogRank=survdiff(KM_H1, data = carcel)
# ver p-valor
LogRank$pvalue
## [1] 0.04722271
COX_H1= formula(survival~edad+casado)
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)
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