Alumna: Romina Ariana Loayza Gaitán

Código: 20211323

Regresión Gaussiana

library(rio)
ceplancito=import("Ceplan.xlsx")
str(ceplancito)

## 'data.frame':    228 obs. of  8 variables:
##  $ TIPO     : chr  "Nacional" "Nacional_esp" "Región" "Zona" ...
##  $ UBIGEO   : chr  "NACIONAL" "NACIONAL SIN LA PROV. DE LIMA" "150000" "ZONA NORTE" ...
##  $ PROVINCIA: chr  NA NA "LIMA" NA ...
##  $ IVIA     : chr  "0.49689161300659201" "-" "0.26664412021637002" "-" ...
##  $ IDE      : chr  "0.75900264069939904" "-" "0.82746140000000001" "-" ...
##  $ ABASTOS  : chr  "2612" "1490" "1232" "-" ...
##  $ IDH      : chr  "0.58576377555328896" "-" "0.70733610728959495" "-" ...
##  $ POBREZA  : num  30.1 31.7 26.2 31.8 33 ...

Cambiar a numérica:

ceplancito$IVIA = as.numeric(ceplancito$IVIA)

## Warning: NAs introducidos por coerción

ceplancito$IDE = as.numeric(ceplancito$IDE)

## Warning: NAs introducidos por coerción

ceplancito$ABASTOS = as.numeric(ceplancito$ABASTOS)

## Warning: NAs introducidos por coerción

ceplancito$IDH = as.numeric(ceplancito$IDH)

## Warning: NAs introducidos por coerción

ceplancito$POBREZA = as.numeric(ceplancito$POBREZA)

str(ceplancito)

## 'data.frame':    228 obs. of  8 variables:
##  $ TIPO     : chr  "Nacional" "Nacional_esp" "Región" "Zona" ...
##  $ UBIGEO   : chr  "NACIONAL" "NACIONAL SIN LA PROV. DE LIMA" "150000" "ZONA NORTE" ...
##  $ PROVINCIA: chr  NA NA "LIMA" NA ...
##  $ IVIA     : num  0.497 NA 0.267 NA 0.868 ...
##  $ IDE      : num  0.759 NA 0.827 NA 0.663 ...
##  $ ABASTOS  : num  2612 1490 1232 NA 26 ...
##  $ IDH      : num  0.586 NA 0.707 NA 0.418 ...
##  $ POBREZA  : num  30.1 31.7 26.2 31.8 33 ...

Completamos los casos:

ceplancito=ceplancito[complete.cases(ceplancito),]

str(ceplancito)

## 'data.frame':    209 obs. of  8 variables:
##  $ TIPO     : chr  "Región" "Región" "Provincia" "Provincia" ...
##  $ UBIGEO   : chr  "150000" "010000" "010200" "010300" ...
##  $ PROVINCIA: chr  "LIMA" "AMAZONAS" "BAGUA" "BONGARÁ" ...
##  $ IVIA     : num  0.267 0.868 0.519 0.509 0.396 ...
##  $ IDE      : num  0.827 0.663 0.649 0.716 0.767 ...
##  $ ABASTOS  : num  1232 26 6 2 4 ...
##  $ IDH      : num  0.707 0.418 0.461 0.413 0.543 ...
##  $ POBREZA  : num  26.2 33 34.8 33.2 22 ...

ceplancito<- subset(ceplancito, TIPO == "Provincia")

Hipótesis: El Índice de Vulnerabilidad a la Inseguridad Alimentaria (IVIA) de las provincias de Perú está determinado por el número de mercados de abastos (ABASTOS), el Índice de Desarrollo Humano (IDH), el Índice de Densidad del Estado (IDE) y el porcentaje de la población en pobreza total (POBREZA).

modelaso1=formula(IVIA~ABASTOS+IDH+IDE+POBREZA)

reg1=lm(modelaso1,data=ceplancito)
summary(reg1)

## 
## Call:
## lm(formula = modelaso1, data = ceplancito)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.120806 -0.032349  0.003728  0.032799  0.145702 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.0379787  0.0584848  17.748  < 2e-16 ***
## ABASTOS     -0.0007544  0.0002477  -3.045 0.002678 ** 
## IDH         -0.9959502  0.0768319 -12.963  < 2e-16 ***
## IDE         -0.2360260  0.0621235  -3.799 0.000199 ***
## POBREZA      0.0021360  0.0005498   3.885 0.000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05284 on 178 degrees of freedom
## Multiple R-squared:  0.9024, Adjusted R-squared:  0.9003 
## F-statistic: 411.7 on 4 and 178 DF,  p-value: < 2.2e-16

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)

modelito1=list('IVIA'=reg1)
modelsummary(modelito1, title = "Regresion: modelo 1",
             stars = TRUE,
             output = "kableExtra")

Regresion: modelo 1

	IVIA
(Intercept)	1.038***
	(0.058)
ABASTOS	-0.001**
	(0.000)
IDH	-0.996***
	(0.077)
IDE	-0.236***
	(0.062)
POBREZA	0.002***
	(0.001)
Num.Obs.	183
R2	0.902
R2 Adj.	0.900
AIC	-549.9
BIC	-530.7
Log.Lik.	280.975
F	411.659
RMSE	0.05
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Interpretación de los coeficientes: El modelo en general es válido, pues tiene un pvalue menor a 0.05 * Todas las variables tiene un efecto estadísitcamente significativo * Las variables ‘ABASTOS’, ‘IDH’ e ‘IDE’ poseen un signo negativo, por lo cual poseen una relación inversa con la variable dependiente (IVIA). Mientras que, ‘POBREZA’ tiene una relación positiva con la VD (IVIA). - La magnitud de ‘ABASTOS’ es de -0.001, lo que significa que cuando ‘ABASTOS’ aumenta en 1 punto, ‘IVIA’ disminuye en 0.001. - La magnitud de ‘IDH’ es de -0.996, lo que significa que cuando ‘IDH’ aumenta en 1 punto, ‘IVIA’ disminuye en 0.996. - La magnitud de ‘IDE’ es de -0.236, lo que significa que cuando ‘IDE’ aumenta en 1 punto, ‘IVIA’ disminuye en 0.236 - La magnitud de ‘POBREZA’ es de 0.002, lo que significa que cuando ‘POBREZA’ aumenta en 1 punto, ‘IVIA’ aumenta en 0.002. * El R2 ajustado nos indica que el modelo explica un 90% de la variabilidad de ‘IVIA’

Diagnóstico ## Linealidad:

plot(reg1,1)

mean(reg1$residuals)

## [1] 3.396889e-18

En base al gráfico y a la diferencia entre el valor esperado y el valor obtenido, se puede evidenciar que sí hay linealidad en el modelo. Esto se debe a que el promedio de los residuos que se obtiene es de 3.396889e-18, el cual es un valor menor a 0.05.

Homocedasticidad:

plot(reg1, 3)

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(reg1)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg1
## BP = 18.107, df = 4, p-value = 0.001176

El valor obtenido, 0.001176, es menor a 0.05. Lo que significa que hay Heterocedasticidad en el modelo, por lo que a medida en la que se encuentren más variables, habrá un mayor margen de error.

Normalidad de residuos:

plot(reg1, 2)

shapiro.test(reg1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  reg1$residuals
## W = 0.99142, p-value = 0.3492

Al no obtener un p-value menor a 0.05, se concluye que los residuos no se distribuyen de una manera normal. Por tanto, esto limita la posibilidad de hacer inferencias con el modelo.

No multicolinealidad:

library(DescTools)

## 
## Attaching package: 'DescTools'

## The following objects are masked from 'package:modelsummary':
## 
##     Format, Mean, Median, N, SD, Var

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, union

VIF(reg1)

##  ABASTOS      IDH      IDE  POBREZA 
## 1.252023 4.971380 1.562895 4.183900

Los resultados obtenidos demuestran que no hay correlación entre los predictores, pues todos los valores son menores a 5. Sin embargo, tanto IDH como POBREZA se acercan a dicho valor.

Valores influyentes:

plot(reg1, 5)

valorsitos = as.data.frame(influence.measures(reg1)$is.inf)
valorsitos[valorsitos$cook.d & valorsitos$hat,]

## [1] dfb.1_   dfb.ABAS dfb.IDH  dfb.IDE  dfb.POBR dffit    cov.r    cook.d  
## [9] hat     
## <0 rows> (or 0-length row.names)

Los resultados obtenidos demuestran que no hay presencia de valores influyentes. Por ende, no hay valores que afecten el resultado de la regresión.

Regresión Poisson

ceplansote=import("Ceplan_2.xlsx")
str(ceplansote)

## 'data.frame':    228 obs. of  9 variables:
##  $ TIPO        : chr  "Nacional" "Nacional_esp" "Región" "Zona" ...
##  $ Ubigeo      : chr  "NACIONAL" "NACIONAL SIN LA PROV. DE LIMA" "150000" "ZONA NORTE" ...
##  $ PROVINCIA   : chr  NA NA "LIMA" NA ...
##  $ POBLACIÓN   : num  31237385 22075063 10135009 8987709 417365 ...
##  $ DESNUTRICION: num  12.07 14.86 5.09 15.74 17.47 ...
##  $ IDE         : chr  "0.75900264069939904" "-" "0.82746140000000001" "-" ...
##  $ SENATI      : num  70 61 12 22 2 0 0 1 0 0 ...
##  $ ABASTOS     : chr  "2612" "1490" "1232" "-" ...
##  $ POBREZA     : num  10352427 7673663 2956030 3171351 155851 ...

Convertir a numéricas:

ceplansote$IDE = as.numeric(ceplansote$IDE)

## Warning: NAs introducidos por coerción

ceplansote$ABASTOS = as.numeric(ceplansote$ABASTOS)

## Warning: NAs introducidos por coerción

ceplansote=ceplansote[complete.cases(ceplansote),]
str(ceplansote)

## 'data.frame':    209 obs. of  9 variables:
##  $ TIPO        : chr  "Región" "Región" "Provincia" "Provincia" ...
##  $ Ubigeo      : chr  "150000" "010000" "010200" "010300" ...
##  $ PROVINCIA   : chr  "LIMA" "AMAZONAS" "BAGUA" "BONGARÁ" ...
##  $ POBLACIÓN   : num  10135009 417365 82193 27085 60419 ...
##  $ DESNUTRICION: num  5.09 17.47 24.9 20.09 13.69 ...
##  $ IDE         : num  0.827 0.663 0.649 0.716 0.767 ...
##  $ SENATI      : num  12 2 0 0 1 0 0 0 1 4 ...
##  $ ABASTOS     : num  1232 26 6 2 4 ...
##  $ POBREZA     : num  2956030 155851 34184 9751 12943 ...

ceplansote<- subset(ceplansote, TIPO == "Provincia")

Hipótesis: El número de habitantes en situación de pobreza (POBREZA) de las provincias de Perú está determinado por el número de mercados de abastos (ABASTOS), el Índice de Densidad del Estado (IDE), el número de sedes de SENATI (SENATI) y el porcentaje de desnutrición crónica en niños menores de 5 años (DESNUTRICION).

hipopotamo=formula(POBREZA~ABASTOS+IDE+SENATI+DESNUTRICION)

names(ceplansote)=gsub(pattern = "POBLACIÓN",
                           replacement = "POBLA",
                           x = names(ceplansote))

ceplansote$POBREZA <- as.integer(ceplansote$POBREZA) #Se lo cambia a entero, ya que la poisson es de conteos y r lo estaba leyendo como número con decimales

str(ceplansote)

## 'data.frame':    183 obs. of  9 variables:
##  $ TIPO        : chr  "Provincia" "Provincia" "Provincia" "Provincia" ...
##  $ Ubigeo      : chr  "010200" "010300" "010100" "010400" ...
##  $ PROVINCIA   : chr  "BAGUA" "BONGARÁ" "CHACHAPOYAS" "CONDORCANQUI" ...
##  $ POBLA       : num  82193 27085 60419 49800 48140 ...
##  $ DESNUTRICION: num  24.9 20.1 13.7 41.5 19.4 ...
##  $ IDE         : num  0.649 0.716 0.767 0.409 0.705 ...
##  $ SENATI      : num  0 0 1 0 0 0 1 0 2 0 ...
##  $ ABASTOS     : num  6 2 4 1 3 2 8 1 14 9 ...
##  $ POBREZA     : int  34184 9751 12942 43128 24472 11073 37887 47331 136834 56750 ...

rp1=glm(hipopotamo, data = ceplansote, 
        offset=log(POBLA), 
        family = poisson(link = "log"))

summary(rp1)

## 
## Call:
## glm(formula = hipopotamo, family = poisson(link = "log"), data = ceplansote, 
##     offset = log(POBLA))
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.052e-01  5.233e-03  -58.34   <2e-16 ***
## ABASTOS      -5.139e-04  1.528e-05  -33.64   <2e-16 ***
## IDE          -2.112e+00  6.693e-03 -315.55   <2e-16 ***
## SENATI       -1.605e-01  6.002e-04 -267.38   <2e-16 ***
## DESNUTRICION  3.582e-02  6.479e-05  552.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2151531  on 182  degrees of freedom
## Residual deviance:  736313  on 178  degrees of freedom
## AIC: 738467
## 
## Number of Fisher Scoring iterations: 4

modelslmpoi=list('POBREZA'=rp1)

modelsummary(modelslmpoi, title = "Regresion Poisson",
             stars = TRUE,
             output = "kableExtra")

Regresion Poisson

	POBREZA
(Intercept)	-0.305***
	(0.005)
ABASTOS	-0.001***
	(0.000)
IDE	-2.112***
	(0.007)
SENATI	-0.160***
	(0.001)
DESNUTRICION	0.036***
	(0.000)
Num.Obs.	183
AIC	738466.7
BIC	738482.8
Log.Lik.	-369228.368
F	380557.753
RMSE	14554.43
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Coeficientes exponenciados:

cbind(exp(coef(rp1)),exp(confint(rp1)))

## Waiting for profiling to be done...

##                            2.5 %    97.5 %
## (Intercept)  0.7369444 0.7294243 0.7445401
## ABASTOS      0.9994862 0.9994562 0.9995161
## IDE          0.1210089 0.1194322 0.1226069
## SENATI       0.8517430 0.8507415 0.8527453
## DESNUTRICION 1.0364642 1.0363326 1.0365958

Justificación del OFFset Considero necesario usar el offset, debido a que el número de habitantes en situación de pobreza guarda relación con la población, asimismo es una manera de controlar la VD. Por ejemplo, no es lo mismo tener 1000 pobres en una población de 5000 a tenerlos en una población de 2000.

Desripción de resultados Todos los predictores son significativos. ## Impacto de ABASTOS

exp_abas = 1-0.9994862
exp_abas = exp_abas*-1
exp_abas*100

## [1] -0.05138

Por cada unidad que aumente ABASTOS, la cantidad esperada de pobreza disminuye en 0.05%.

Impacto de IDE

exp_ide = 1-0.1210089
exp_ide = exp_ide*-1
exp_ide*100

## [1] -87.89911

Por cada unidad que aumente IDE, la cantidad esperada de personas en situación de pobreza disminuye en 87.9%

Impacto de SENATI

exp_senati = 1-0.8517430
exp_senati = exp_senati*-1
exp_senati*100

## [1] -14.8257

Por cada unidad que aumente SENATI, la cantidad esperada de personas en situación de pobreza disminuye en 14.8%

Impacto de DESNUTRICIÓN

exp_desnutrición = 1-1.0364642
exp_desnutrición = exp_desnutrición*-1
exp_desnutrición*100

## [1] 3.64642

Por cada unidad que aumente desnutrición, la cantidad esperada de personas en situación de pobreza aumenta en 3.64%.

Supuestos

Equidispersión

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(magrittr)
overdispersion=AER::dispersiontest(rp1,alternative='greater')$ p.value<0.05
underdispersion=AER::dispersiontest(rp1,alternative='less')$ p.value<0.05
# tabla
testResult=as.data.frame(rbind(overdispersion,underdispersion))
names(testResult)='Es probable?'
testResult%>%kable(caption = "Test de Equidispersión")%>%kableExtra::kable_styling()

Test de Equidispersión

	Es probable?
overdispersion	TRUE
underdispersion	FALSE

En base a los resultados de la tabla, el modelo presenta sobredispersión.

Modelo Alternativo

Quasi Poisson

laquasi = glm(hipopotamo, data = ceplansote,
          offset=log(POBLA),
          family = quasipoisson(link = "log"))

modelsQPexp=list('QuasiPoisson asegurados (II) exponenciado'=laquasi)


modelsummary(modelsQPexp,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "EXP() de la Regresión Quasi Poisson (II) para Interpretación",
             stars = TRUE,
             output = "kableExtra")

EXP() de la Regresión Quasi Poisson (II) para Interpretación

	QuasiPoisson asegurados (II) exponenciado
(Intercept)	0.737
	[0.389, 1.382]
ABASTOS	0.999
	[0.998, 1.001]
IDE	0.121***
	[0.054, 0.274]
SENATI	0.852***
	[0.791, 0.915]
DESNUTRICION	1.036***
	[1.028, 1.045]
Num.Obs.	183
F	99.766
RMSE	14554.43
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Binomial negativa

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

hipoff=formula(POBREZA~ABASTOS+IDE+SENATI+DESNUTRICION+offset(log(POBLA)))
negativa=glm.nb(hipoff,data=ceplansote)

modelsQP_BN=list('Poisson asegurados (II)'=rp1,
                 'QuasiPoisson asegurados (II)'=laquasi,
                 'Binomial Negativa asegurados (II)'=negativa)


modelsummary(modelsQP_BN,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "EXP() de la Regresiones Poisson, Quasi Poisson  y Binomial Negativa",
             stars = TRUE,
             output = "kableExtra")

EXP() de la Regresiones Poisson, Quasi Poisson y Binomial Negativa

	Poisson asegurados (II)	QuasiPoisson asegurados (II)	Binomial Negativa asegurados (II)
(Intercept)	0.737***	0.737	0.605
	[0.729, 0.745]	[0.389, 1.382]	[0.300, 1.241]
ABASTOS	0.999***	0.999	0.999
	[0.999, 1.000]	[0.998, 1.001]	[0.996, 1.003]
IDE	0.121***	0.121***	0.142***
	[0.119, 0.123]	[0.054, 0.274]	[0.055, 0.360]
SENATI	0.852***	0.852***	0.864**
	[0.851, 0.853]	[0.791, 0.915]	[0.784, 0.956]
DESNUTRICION	1.036***	1.036***	1.043***
	[1.036, 1.037]	[1.028, 1.045]	[1.033, 1.052]
Num.Obs.	183	183	183
AIC	738466.7		3824.2
BIC	738482.8		3843.5
Log.Lik.	-369228.368		-1906.115
F	380557.753	99.766	58.083
RMSE	14554.43	14554.43	15019.97
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Se verifica la presencia de la sobredispersión:

performance::check_overdispersion(rp1)

## # Overdispersion test
## 
##        dispersion ratio =   3814.455
##   Pearson's Chi-Squared = 678973.049
##                 p-value =    < 0.001

## Overdispersion detected.

performance::check_overdispersion(laquasi)

## # Overdispersion test
## 
##        dispersion ratio =   3814.455
##   Pearson's Chi-Squared = 678973.049
##                 p-value =    < 0.001

## Overdispersion detected.

performance::check_overdispersion(negativa)

## # Overdispersion test
## 
##  dispersion ratio = 0.623
##           p-value =  0.08

## No overdispersion detected.

En base a los resultados de la tabla y los resultados de sobredispersión, se identifica que el mejor modelo es la regresión Binomial Negativa. Asimismo, en el caso de los AIC, se evidencia un mejor resultado para la BN.

Regresión Logísitca

ceplancito$Vulnerable <- ifelse(ceplancito$IVIA < 0.5, 0, 1)

str(ceplancito)

## 'data.frame':    183 obs. of  9 variables:
##  $ TIPO      : chr  "Provincia" "Provincia" "Provincia" "Provincia" ...
##  $ UBIGEO    : chr  "010200" "010300" "010100" "010400" ...
##  $ PROVINCIA : chr  "BAGUA" "BONGARÁ" "CHACHAPOYAS" "CONDORCANQUI" ...
##  $ IVIA      : num  0.519 0.509 0.396 0.691 0.582 ...
##  $ IDE       : num  0.649 0.716 0.767 0.409 0.705 ...
##  $ ABASTOS   : num  6 2 4 1 3 2 8 1 14 9 ...
##  $ IDH       : num  0.461 0.413 0.543 0.253 0.341 ...
##  $ POBREZA   : num  34.8 33.2 22 56.8 48.1 ...
##  $ Vulnerable: num  1 1 0 1 1 1 1 1 0 1 ...

hiphip=formula(Vulnerable~ABASTOS+IDH+IDE+POBREZA)
snooplog=glm(hiphip, data=ceplancito,family = binomial)
modelrl=list('Ser Vulnerable (I)'=snooplog)
library(modelsummary)
modelsummary(modelrl,
             title = "Regresión Logística",
             stars = TRUE,
             output = "kableExtra")

Regresión Logística

	Ser Vulnerable (I)
(Intercept)	23.958**
	(7.708)
ABASTOS	0.030
	(0.074)
IDH	-43.278***
	(11.131)
IDE	-10.077
	(7.008)
POBREZA	0.073
	(0.052)
Num.Obs.	183
AIC	67.0
BIC	83.0
Log.Lik.	-28.483
F	6.794
RMSE	0.22
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Evaluación del modelo

ceplancito$predictedProbs <- predict(snooplog, ceplancito,type = "response")
head(ceplancito$predictedProbs)

## [1] 0.546590568 0.795312840 0.003943837 0.999997851 0.996658123 0.977952660

ceplancito$predicted=ifelse(ceplancito$predictedProbs>0.5,"1","0")

head(ceplancito[,c("predicted","Vulnerable")],10)

##    predicted Vulnerable
## 6          1          1
## 7          1          1
## 8          0          0
## 9          1          1
## 10         1          1
## 11         1          1
## 12         0          1
## 14         1          1
## 15         0          0
## 16         1          1

library(cvms)
ConfMatrix=confusion_matrix(predictions =  ceplancito$predicted,
                      targets= ceplancito$Vulnerable)
library(cvms)
plot_confusion_matrix(ConfMatrix,
                      class_order=c("0","1"))

## Warning in plot_confusion_matrix(ConfMatrix, class_order = c("0", "1")):
## 'ggimage' is missing. Will not plot arrows and zero-shading.

## Warning in plot_confusion_matrix(ConfMatrix, class_order = c("0", "1")): 'rsvg'
## is missing. Will not plot arrows and zero-shading.

Sensitividad

ConfMatrix$Sensitivity

## [1] 0.950495

Especificidad

ConfMatrix$Specificity

## [1] 0.9268293

Kapa

ConfMatrix$Kappa

## [1] 0.8783318

El resultado del Kappa es 0.87, el cual es un valor cercano a 1. Por lo que, se evidencia que hay un alto grado de predicción de variables.

library(margins)
library(dplyr)
library(kableExtra)
marginalsData=summary(margins(snooplog))
marginalsData%>% kable(caption = "Efectos Marginales Promedio (AME)- Modelo III") %>%kableExtra::kable_styling(full_width = T)

Efectos Marginales Promedio (AME)- Modelo III

factor	AME	SE	z	p	lower	upper
ABASTOS	0.0014139	0.0035122	0.4025778	0.6872588	-0.0054699	0.0082977
IDE	-0.4791767	0.3241918	-1.4780653	0.1393903	-1.1145811	0.1562276
IDH	-2.0578376	0.3592293	-5.7284796	0.0000000	-2.7619140	-1.3537611
POBREZA	0.0034730	0.0024163	1.4372931	0.1506347	-0.0012629	0.0082089

Interpretación de los coeficientes De las cuatro variables usadas como predictores, se evidencia que solo el IDH tiene un efecto estadisticamente significativo respecto a la ser vulnerable de inseguridad alimentaria. Asimismo, esto se confirma al evaluar los margins, pues su intervalo de confianza no incluye al 0. Entonces, por cada punto que IDH aumente, la probabilidad de ser VULNERABLE ante la IVIA disminuye en 2.05.

library(ggplot2)
base= ggplot(marginalsData,aes(x=factor, y=AME)) + geom_point()
base +  geom_errorbar(aes(ymin=lower, ymax=upper))

# Regresión Cox

chilenitos=import("MinistrosChile.csv")
str(chilenitos)

## 'data.frame':    180 obs. of  12 variables:
##  $ id             : int  1 2 4 5 6 7 9 10 11 12 ...
##  $ name           : chr  "Abeliuk Manasevic, Ren\xe9" "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" "3/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 ...

Hipótesis: El sexo (sex), la edad (age) y el presidente (president) influyen sobre el riesgo de salida ministerial. Para la supervivencia usa las variables enministerio (tiempo en días) y evento.

hipominis=formula(enministerio~sex+age+president)

library(survival)
chilenitos$survival=with(chilenitos,Surv(time = enministerio,event =  as.numeric(evento)))
library(magrittr) 
library(dplyr)
chilenitos%>%
    rmarkdown::paged_table()

library(ggplot2)
library(ggfortify)

KM.generico = survfit(survival ~ 1, data = chilenitos)

ejeX='DAYS\n Curva cae cuando destituyen al ministro'
ejeY='Probabilidad \n(Mantenerse en el cargo)'
titulo="Curva de Sobrevivencia"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)

chilenitos$sex <- factor(chilenitos$sex, levels = c("man", "woman"))
chilenitos$sex <- as.numeric(chilenitos$sex)
str(chilenitos)

## 'data.frame':    180 obs. of  13 variables:
##  $ id             : int  1 2 4 5 6 7 9 10 11 12 ...
##  $ name           : chr  "Abeliuk Manasevic, Ren\xe9" "Albornoz Pollman, Laura" "Alvarado Constenla, Luis" "Alvear Valenzuela, Soledad" ...
##  $ sex            : num  1 2 1 2 2 2 2 1 1 2 ...
##  $ 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" "3/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 ...
##  $ survival       : 'Surv' num [1:180, 1:2] 1461+ 1318  1461+ 1163+ 1663  2106   142+ 1005  2093   125  ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "time" "status"
##   ..- attr(*, "type")= chr "right"

KM_H1=formula(survival ~ sex)

KM.sex = survfit(KM_H1, data = chilenitos)

###
ejeX='DAYS\n Curva cae cuando destituyen al ministro'
ejeY='Probabilidad \n(Mantenerse en el cargo)'
titulo="Curva de Sobrevivencia"

autoplot(KM.sex,xlab=ejeX,ylab=ejeY, 
         main = titulo,conf.int = F)  + 
        labs(colour = "Impacto del género") + 
         scale_color_discrete(labels = c("Hombre", "Mujer"))

autoplot(KM.sex,xlab=ejeX,ylab=ejeY, 
         main = titulo,conf.int = T)+ 
        labs(colour = "Impacto del género") + 
         scale_color_discrete(labels = c("Hombre", "Mujer"))

LogRank=survdiff(KM_H1, data = chilenitos)
LogRank$pvalue

## [1] 0.8969547

Interpretación Kaplan Meier Según el análisis KM, el sexo no tiene un efecto sobre la sobrevivencia de un ministro en su cargo. Además, esto se confirma al interpretar el longrank, pues este no es menor a 0.05. Por ende, la variable sexo no tendría un efecto significativo en un potencial modelo.

Regresión COX

COX_H1= formula(survival~factor(sex)+age+president)
rcox1 <- coxph(COX_H1,data=chilenitos)
summary(rcox1)

## Call:
## coxph(formula = COX_H1, data = chilenitos)
## 
##   n= 180, number of events= 97 
## 
##                                  coef  exp(coef)   se(coef)      z Pr(>|z|)   
## factor(sex)2               -0.2887245  0.7492186  0.2673500 -1.080  0.28016   
## age                        -0.0008448  0.9991555  0.0121047 -0.070  0.94436   
## presidentMichelle Bachelet  0.2271447  1.2550115  0.2810804  0.808  0.41903   
## presidentPatricio Aylwin   -1.4360901  0.2378559  0.4482177 -3.204  0.00136 **
## presidentRicardo Lagos     -0.0644762  0.9375585  0.2456334 -0.262  0.79294   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                            exp(coef) exp(-coef) lower .95 upper .95
## factor(sex)2                  0.7492     1.3347   0.44365    1.2653
## age                           0.9992     1.0008   0.97573    1.0231
## presidentMichelle Bachelet    1.2550     0.7968   0.72342    2.1772
## presidentPatricio Aylwin      0.2379     4.2042   0.09881    0.5726
## presidentRicardo Lagos        0.9376     1.0666   0.57932    1.5173
## 
## Concordance= 0.61  (se = 0.029 )
## Likelihood ratio test= 18.39  on 5 df,   p=0.002
## Wald test            = 12.99  on 5 df,   p=0.02
## Score (logrank) test = 15.02  on 5 df,   p=0.01

modelcox=list('Riesgo - Ser destituido'=rcox1,'Riesgo- Ser destituido (exponenciado)'=rcox1)
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 - Ser destituido	Riesgo- Ser destituido (exponenciado)
factor(sex)2	-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

(presidentAylwin=abs(1-exp(coef(rcox1)[4])))

## presidentPatricio Aylwin 
##                0.7621441

En base a los resultados de la tabla, solo la variable ‘presidentPatricioAylwin’ tiene un efecto estadísticamente significativo. Por tanto, el riesgo de ser destituido siendo ministro del gobierno de Patricio Aylwin es, en promedio, 76.2% menor a comparación de otros ministros que no fueron parte de su gobierno. Asimismo, se evidencia que las demás variables no poseen un efecto estadísticamente significativo.