library(rio)
chile=import("MinistrosChile - MinistrosChile (1).csv")
str(chile)
## 'data.frame':    180 obs. of  12 variables:
##  $ id             : int  1 2 4 5 6 7 9 10 11 12 ...
##  $ name           : chr  "Abeliuk Manasevic, René" "Albornoz Pollman, Laura" "Alvarado Constenla, Luis" "Alvear Valenzuela, Soledad" ...
##  $ sex            : chr  "man" "woman" "man" "woman" ...
##  $ age            : int  59 37 49 40 49 43 55 52 46 53 ...
##  $ president      : chr  "Patricio Aylwin" "Michelle Bachelet" "Patricio Aylwin" "Patricio Aylwin" ...
##  $ start_president: chr  "11/03/1990" "11/03/2006" "11/03/1990" "11/03/1990" ...
##  $ end_president  : chr  "11/03/1994" "11/03/2010" "11/03/1994" "11/03/1994" ...
##  $ ministry       : chr  "CORFO" "SERNAM" "Bienes Nacionales" "SERNAM" ...
##  $ start_minister : chr  "11/03/1990" "11/03/2006" "11/03/1990" "03/01/1991" ...
##  $ end_minister   : chr  "11/03/1994" "19/10/2009" "11/03/1994" "11/03/1994" ...
##  $ enministerio   : int  1461 1318 1461 1163 1663 2106 142 1005 2093 125 ...
##  $ evento         : int  0 1 0 0 1 1 0 1 1 1 ...
library(survival)

Creo mi columna survival

chile$survival=with(chile, Surv(time = enministerio, event = evento))
library(magrittr)
chile%>%
  rmarkdown::paged_table()
library(ggplot2)
library(ggfortify)
KM.generico = survfit(survival ~ 1, data = chile)
ejeX='días\n curva cae cuando un ministro es retirado'
ejeY='Probabilidad \n(PERMANECER EN EL MINISTERIO)'
titulo="Curva de Sobrevivencia: permanecer en el ministerio"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)

en el primer día es 100% probable que te retiren del cargo ministerial, pero a los 1000 días de estar ocupando dicho cargo las probabilidades de permacer allí se reducen al 50%. Peor aún al pasar 2000 días ocupando el cargo, las probablilidades de quedarte en el ministerio son menores al 25%.

Comparando según el sexo:

str(chile$sex)
##  chr [1:180] "man" "woman" "man" "woman" "woman" "woman" "woman" "man" ...
chile$sex=as.factor(chile$sex)
str(chile$sex)
##  Factor w/ 2 levels "man","woman": 1 2 1 2 2 2 2 1 1 2 ...
KM_H1=formula(survival ~ sex)

KM.genero = survfit(KM_H1, data = chile)

###
ejeX='días\n curva cae cuando un ministro es retirado'
ejeY="Prob ('seguir en ministerio')"
titulo="Curva de Sobrevivencia: ¿Beneficia el género?"

autoplot(KM.genero,xlab=ejeX,ylab=ejeY, 
         main = titulo,conf.int = F)  + 
        labs(colour = "sexo o genero del ministro") + 
         scale_color_discrete(labels = c("hombre", "mujer"))

LogRank=survdiff(KM_H1, data = chile)
# ver p-valor
LogRank$pvalue
## [1] 0.8969547
autoplot(KM.genero,xlab=ejeX,ylab=ejeY, 
         main = titulo,conf.int = T)+ 
        labs(colour = "sexo o genero del ministro") + 
         scale_color_discrete(labels = c("hombre", "mujer"))

¿crees que la variable tendrá efecto significativo en un potencial modelo? (0.5 puntos)

NO, YA QUE EL P-VALOR ES MAYOR A 0.05, ES DECIR, NO HAY DIFERENCIA SIGNIFICATIVA.

COX_H1= formula(survival~sex+age+president)

#regression
rcox1 <- coxph(COX_H1,data=chile)
modelcox=list('Riesgo - que lo saquen del ministerio'=rcox1,'Riesgo- que lo saquen del ministerio (exponenciado)'=rcox1)

#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
##   backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
## 
## Revert to `kableExtra` for one session:
## 
##   options(modelsummary_factory_default = 'kableExtra')
## 
## Change the default backend persistently:
## 
##   config_modelsummary(factory_default = 'gt')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)
modelsummary(modelcox,
             #fmt=f,
             exponentiate = c(F,T), 
             statistic = 'conf.int',
             title = "Regresión Cox",
             stars = TRUE,
             output = "kableExtra")
Regresión Cox
Riesgo - que lo saquen del ministerio Riesgo- que lo saquen del ministerio (exponenciado)
sexwoman -0.289 0.749
[-0.813, 0.235] [0.444, 1.265]
age -0.001 0.999
[-0.025, 0.023] [0.976, 1.023]
presidentMichelle Bachelet 0.227 1.255
[-0.324, 0.778] [0.723, 2.177]
presidentPatricio Aylwin -1.436** 0.238**
[-2.315, -0.558] [0.099, 0.573]
presidentRicardo Lagos -0.064 0.938
[-0.546, 0.417] [0.579, 1.517]
Num.Obs. 180 180
AIC 825.7 825.7
BIC 841.7 841.7
RMSE 0.72 0.72
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 
library(survminer)
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
ggforest(rcox1, data = chile,main = "¿Quiénes tienen mayor riesgo de ser sacados del ministerio?")

INTERPRETACIONES:

#1) Ninguna variable es significativa, excepto “presidentPatricio Aylwin”. ENTONCES SE PUEDE INFERIR QUE ES MUY POCO PROBABLE QUE EL SEXO Y LA EDAD DE UN MINISTRO TENGAN EFECTO EN LA DECISIÓN DE CESARLOS DEL CARGO

#2) PODEMOS DECIR TAMBIÉN, QUE EN EL GOBIERNO DE PATRICIO AYLWIN DISMINUYE EL RIESGO DE SER SACADO DEL MINISTERIO CON UNA SIGNIFICANCIA DEL 0.01.

(presidentPatricioAylwin=abs(1-exp(coef(rcox1)[4])))
## presidentPatricio Aylwin 
##                0.7621441

#3) Usando el hazard ratio, podemos sostener QUE EL RIESGO DE SER SACADO DE UN MINISTERIO DURANTE EL GOBIERNO DE PATRICIO AYLWIN es en promedio 76% menor a la de los OTROS GOBIERNOS.