Práctica Laboratorio 5

1. Regresión Gaussiana

#install.packages("rio")
library(rio)
data_A = import("Ceplan.xlsx")

Filtramos la base de datos para solo tomar en cuenta las provincias

#install.packages("dplyr")
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

data_Ceplan <- filter(data_A, TIPO == "Provincia")

Comprobamos nuestras variables

str(data_Ceplan)

## 'data.frame':    195 obs. of  8 variables:
##  $ TIPO     : chr  "Provincia" "Provincia" "Provincia" "Provincia" ...
##  $ UBIGEO   : chr  "010200" "010300" "010100" "010400" ...
##  $ PROVINCIA: chr  "BAGUA" "BONGARÁ" "CHACHAPOYAS" "CONDORCANQUI" ...
##  $ IVIA     : chr  "0.51947927474975597" "0.50935518741607699" "0.39617508649826" "0.69078904390335105" ...
##  $ IDE      : chr  "0.64920285677771705" "0.71639089022042102" "0.76727899863208504" "0.40882896570212401" ...
##  $ ABASTOS  : chr  "6" "2" "4" "1" ...
##  $ IDH      : chr  "0.46103622424613799" "0.41289545257260002" "0.54266528829961902" "0.2534515691766" ...
##  $ POBREZA  : num  34.8 33.2 22 56.8 48.1 ...

Cambiamos al tipo de variable correspondiente y completamos casos

data_Ceplan$IVIA=as.numeric(data_Ceplan$IVIA)
data_Ceplan$ABASTOS=as.numeric(data_Ceplan$ABASTOS)

## Warning: NAs introducidos por coerción

data_Ceplan$IDH=as.numeric(data_Ceplan$IDH)
data_Ceplan$IDE=as.numeric(data_Ceplan$IDE)
data_Ceplan=data_Ceplan[complete.cases(data_Ceplan),]

Verificamos

str(data_Ceplan)

## 'data.frame':    183 obs. of  8 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 ...

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

reg_gauss = lm(IVIA ~ ABASTOS + IDH + IDE + POBREZA, data = data_Ceplan)
summary(reg_gauss)

## 
## Call:
## lm(formula = IVIA ~ ABASTOS + IDH + IDE + POBREZA, data = data_Ceplan)
## 
## 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

Si queremos verlo mejor:

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)

modelo_gauss=list('Vulnerabilidad' = reg_gauss)
modelsummary(modelo_gauss, title = "Regresion Gausseana",
             stars = TRUE,
             output = "kableExtra")

Regresion Gausseana
	Vulnerabilidad
(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

Los resultados del análisis nos muestran lo siguiente:

En el presente modelo, todas las variables poseen un efecto estadísticamente significativo . En primer lugar, la variable ABASTOS sostiene una relación inversa con nuestra variable dependiente (IVIA). Su magnitud resulta ser de -0.001, por lo que cuando ABASTOS aumenta en 1 punto, la variable dependiente disminuyte en 0.001. En segundo lugar, la variable IDH también sostiene una relación inversa con la variable dependiente. Su magnitud resultar ser de -0.996, lo que significa que cuando IDH aumenta en 1 punto, la variable IVA disminuye en 0.996. En tercer lugar, la variable IDE también sostiene una relación inversa con nuestra variable dependiente. Su magnitud es de -0.236, lo que señala que cuando IDE aumenta en 1 punto, IVIA disminuye en 0.236. Finalmente, la variable POBREZA, a diferencia del resto, sí sostiene una relación positiva con la variable dependiente. Cuando la POBREZA aumenta en 1 punto, IVIA aumenta en 0.002. Asimismo, el R2 ajustado nos señala que el modelo explica el 90% de la variabilidad de IVIA

Procedemos con los diagnósticos del modelo

Linealidad:

plot(reg_gauss,1)

Sacamos el promedio de los residuos (la diferencia entre valor esperado y observado)

mean(reg_gauss$residuals)

## [1] 1.911203e-18

Homocedasticidad:

plot(reg_gauss, 3)

Buscamos un p-valor mayor a .05. Esto nos muestra que nuestro modelo es homocedástico. La varianza de los errores es constante y no muestra un patrón.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

bptest(reg_gauss)

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

Normalidad de residuos

plot(reg_gauss, 2)

Buscamos un p-valor menor a .05. Esto indica que los residuos (la distancia entre el valor esperado y el valor observado) se distribuyen de manera normal.

shapiro.test(reg_gauss$residuals)

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

No multicolinealidad

Que no haya una correlación muy alta entre predictores. Idealmente menores a 3, aceptable hasta menor a 5.

library(DescTools)

## 
## Attaching package: 'DescTools'

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

VIF(reg_gauss)

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

Valores influyentes

plot(reg_gauss, 5)

check = as.data.frame(influence.measures(reg_gauss)$is.inf)
check[check$cook.d & check$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)

Análisis:

En nuestro modelo comprobamos que sí existe linealidad, ya que el promedio que se obtiene (1.911203e-18) resulta ser menor a 0.05.
Por otro lado, nuestro modelo no cuenta como homocedasticidad, ya que el valor obtenido (0.001176) es menor a 0.05. Nuestro modelo es, pues, heterocedástico
Asimismo, concluimos que los residuos no se distribuyen de una manera normal, ya que el valor 0.3492 resulta ser mayor a 0.05.
Además,los resultados nos señalan que no existe correlación entre nuestros predictores, ya que todos son menores a 5.
Por último, nuestra regresión no presenta valores influeyentes que puedan afectar el resultado final.

2. Regresión Poisson

data_B = import("Ceplan_2.xlsx")

Filtramos la base de datos

data_Ceplan2 <- filter(data_B, TIPO == "Provincia")

Comprobamos el tipo de datos

str(data_Ceplan2)

## 'data.frame':    195 obs. of  9 variables:
##  $ TIPO        : chr  "Provincia" "Provincia" "Provincia" "Provincia" ...
##  $ Ubigeo      : chr  "010200" "010300" "010100" "010400" ...
##  $ PROVINCIA   : chr  "BAGUA" "BONGARÁ" "CHACHAPOYAS" "CONDORCANQUI" ...
##  $ POBLACIÓN   : num  82193 27085 60419 49800 48140 ...
##  $ DESNUTRICION: num  24.9 20.1 13.7 41.5 19.4 ...
##  $ IDE         : chr  "0.64920285677771705" "0.71639089022042102" "0.76727899863208504" "0.40882896570212401" ...
##  $ SENATI      : num  0 0 1 0 0 0 1 0 2 0 ...
##  $ ABASTOS     : chr  "6" "2" "4" "1" ...
##  $ POBREZA     : num  34184 9751 12943 43128 24472 ...

Cambiamos las variables respectivas y completamos los casos

data_Ceplan2$IDE=as.numeric(data_Ceplan2$IDE)
data_Ceplan2$ABASTOS=as.numeric(data_Ceplan2$ABASTOS)

## Warning: NAs introducidos por coerción

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

Ponemos la variable pobreza como número entero, ya que se trata de una regresión de conteos

data_Ceplan2$POBREZA <- as.integer(data_Ceplan2$POBREZA)

Le cambiamos el nombre a la variable “Población”, ya que la tilde puede generar problemas

names(data_Ceplan2)=gsub(pattern = "POBLACIÓN",
                           replacement = "POBLACION",
                           x = names(data_Ceplan2))

Volvemos a comprobar

str(data_Ceplan2)

## '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" ...
##  $ POBLACION   : 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 ...

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

hp1 = formula(POBREZA ~ ABASTOS + IDE + SENATI + DESNUTRICION)
reg_poisson=glm(hp1, data = data_Ceplan2, 
        offset=log(POBLACION),
        family = poisson(link = "log"))


modelP=list('POBREZA'=reg_poisson)
modelsummary(modelP, 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

¿Qué nos dice la tabla?: La tabla nos señala que todos los predictores resultan significativos a un nivel de 0.0001. Sin embargo, tenemos que exponenciar los resultados para determinar el efecto real.

Exponenciamos los coeficientes

formatoNum <- function(x) format(x, digits = 4, scientific = FALSE) #fórmula para limitar a 4 decimales

modelsummary(modelP, fmt=formatoNum, title = "Regresion Poisson Exponenciada",
             stars = TRUE,
             statistic = 'conf.int',
             exponentiate = T, # exponenciar!!!!!
             output = "kableExtra")

Regresion Poisson Exponenciada
	POBREZA
(Intercept)	0.7369***
	[0.7294, 0.7445]
ABASTOS	0.9995***
	[0.9995, 0.9995]
IDE	0.1210***
	[0.1194, 0.1226]
SENATI	0.8517***
	[0.8507, 0.8527]
DESNUTRICION	1.0365***
	[1.0363, 1.0366]
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

Justificación de offset(Poblacion): El offset (en este caso la variable Población) resulta necesario en esta regresión, ya que el propio número de habitantes en pobreza guarda completa relación con la población en total, ya que solo así se puede determinar la proporción real.

Vemos el impacto de cada variable sobre la probabilidad del modelo

ABASTOS

exp(coef(reg_poisson)[['ABASTOS']])

## [1] 0.9994862

((1-0.9994862)*-1)*100

## [1] -0.05138

Interpretación: Por cada unidad en aumento de ABASTOS, la cantidad de personas en pobreza disminuye en 0.05%.

IDE

exp(coef(reg_poisson)[['IDE']])

## [1] 0.1210089

((1-0.1210089)*-1)*100

## [1] -87.89911

Interpretación: Por cada unidad en aumento de IDE, la cantidad de personas en pobreza disminuye en 87.9%

SENATI

exp(coef(reg_poisson)[['SENATI']])

## [1] 0.851743

((1-0.851743)*-1)*100

## [1] -14.8257

Interpretación: Por cada unidad en aumento de SENATI, la cantidad de personas en pobreza disminuye en 14.82%

DESNUTRICIÓN

exp(coef(reg_poisson)[['DESNUTRICION']])

## [1] 1.036464

((1-1.036464)*-1)*100

## [1] 3.6464

Interpretación: Por cada unidad en aumento de DESNUTRICIÓN, la cantida de personas en pobreza aumenta en 3.64%.

Comprobamos uno de los supuestos en la Regresión Poisson: la media y la varianza son iguales

library(magrittr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

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

overdispersion=AER::dispersiontest(reg_poisson,alternative='greater')$ p.value<0.05
underdispersion=AER::dispersiontest(reg_poisson,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

El presente modelo presenta sobredispersión, por lo que que tenemos que realizar la Quasipoisson y la Binomial Negativa

2.1. Quasipoisson

reg_quasipoisson = glm(hp1, data = data_Ceplan2,
          offset=log(POBLACION),
          family = quasipoisson(link = "log"))

modelQP=list('QUASIPOISSON'=reg_quasipoisson)

modelsummary(modelQP, 
             title = "Regresión QuasiPoisson",
             stars = TRUE,
             output = "kableExtra")

Regresión QuasiPoisson
	QUASIPOISSON
(Intercept)	-0.305
	(0.323)
ABASTOS	-0.001
	(0.001)
IDE	-2.112***
	(0.413)
SENATI	-0.160***
	(0.037)
DESNUTRICION	0.036***
	(0.004)
Num.Obs.	183
F	99.766
RMSE	14554.43
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Revisamos coeficientes (comparación Poisson y Quasipoisson)

cbind(coef_Poi=coef(reg_poisson),coef_QuasiPoi=coef(reg_quasipoisson))

##                   coef_Poi coef_QuasiPoi
## (Intercept)  -0.3052428465 -0.3052428465
## ABASTOS      -0.0005139486 -0.0005139486
## IDE          -2.1118909897 -2.1118909897
## SENATI       -0.1604704541 -0.1604704541
## DESNUTRICION  0.0358151232  0.0358151232

Revisamos errores típicos (comparación Poisson y Quasipoisson)

library(arm)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

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

## Loading required package: Matrix

## Loading required package: lme4

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:rio':
## 
##     factorize

## 
## arm (Version 1.14-4, built: 2024-4-1)

## Working directory is F:/PUCP/PUCP Sociales/Ciclo 6/Estadistica 2/Prueba Laboratorio 5 - Oficial

cbind(se_Poi=se.coef(reg_poisson),se_QuasiPoi=se.coef(reg_quasipoisson))

##                    se_Poi  se_QuasiPoi
## (Intercept)  5.232535e-03 0.3231692710
## ABASTOS      1.527966e-05 0.0009436949
## IDE          6.692732e-03 0.4133532617
## SENATI       6.001680e-04 0.0370672818
## DESNUTRICION 6.479055e-05 0.0040015625

Exponenciamos los coeficientes de la regresión Quasipoisson

formatoNum <- function(x) format(x, digits = 4, scientific = FALSE) #fórmula para limitar a 4 decimales

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

EXP() de la Regresión Quasi Poisson para Interpretación
	QUASIPOISSON
(Intercept)	0.7369
	[0.38919, 1.3817]
ABASTOS	0.9995
	[0.99760, 1.0013]
IDE	0.1210***
	[0.05418, 0.2739]
SENATI	0.8517***
	[0.79147, 0.9153]
DESNUTRICION	1.0365***
	[1.02835, 1.0446]
Num.Obs.	183
F	99.766
RMSE	14554.43
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

2.2. Binomial Negativa

hp2=formula(POBREZA ~ ABASTOS + IDE + SENATI + DESNUTRICION + offset(log(POBLACION)))

reg_binomialnegativa=glm.nb(hp2,data=data_Ceplan2)

modelBN=list('Binomial Negativa'=reg_binomialnegativa)

modelsummary(modelBN,fmt=formatoNum,
             exponentiate = F, 
             statistic = 'conf.int',
             title = "EXP() de la Regresion Binomial Negativa",
             stars = TRUE,
             output = "kableExtra")

EXP() de la Regresion Binomial Negativa
	Binomial Negativa
(Intercept)	-0.5025895
	[-1.204778, 0.215726]
ABASTOS	-0.0007181
	[-0.004052, 0.003042]
IDE	-1.9490077***
	[-2.896151, -1.022766]
SENATI	-0.1460524**
	[-0.243406, -0.044755]
DESNUTRICION	0.0417278***
	[ 0.032401, 0.051099]
Num.Obs.	183
AIC	3824.2
BIC	3843.5
Log.Lik.	-1906.115
F	58.083
RMSE	15019.97
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Exponenciamos los coeficientes de la regresión Binomial Negativa

formatoNum <- function(x) format(x, digits = 4, scientific = FALSE) #fórmula para limitar a 4 decimales

modelsummary(modelBN,fmt=formatoNum,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "EXP() de la Regresión BN para Interpretación",
             stars = TRUE,
             output = "kableExtra")

EXP() de la Regresión BN para Interpretación
	Binomial Negativa
(Intercept)	0.6050
	[0.29976, 1.2408]
ABASTOS	0.9993
	[0.99596, 1.0030]
IDE	0.1424***
	[0.05524, 0.3596]
SENATI	0.8641**
	[0.78395, 0.9562]
DESNUTRICION	1.0426***
	[1.03293, 1.0524]
Num.Obs.	183
AIC	3824.2
BIC	3843.5
Log.Lik.	-1906.115
F	58.083
RMSE	15019.97
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Juntamos los 3 modelos: Poisson, QuasiPoisson y Binomial Negativa

models_total=list('Poisson asegurados'=reg_poisson,
                 'QuasiPoisson asegurados'=reg_quasipoisson,
                 'Binomial Negativa asegurados'=reg_binomialnegativa)


modelsummary(models_total,fmt=formatoNum,
             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	QuasiPoisson asegurados	Binomial Negativa asegurados
(Intercept)	0.7369***	0.7369	0.6050
	[0.7294, 0.7445]	[0.38919, 1.3817]	[0.29976, 1.2408]
ABASTOS	0.9995***	0.9995	0.9993
	[0.9995, 0.9995]	[0.99760, 1.0013]	[0.99596, 1.0030]
IDE	0.1210***	0.1210***	0.1424***
	[0.1194, 0.1226]	[0.05418, 0.2739]	[0.05524, 0.3596]
SENATI	0.8517***	0.8517***	0.8641**
	[0.8507, 0.8527]	[0.79147, 0.9153]	[0.78395, 0.9562]
DESNUTRICION	1.0365***	1.0365***	1.0426***
	[1.0363, 1.0366]	[1.02835, 1.0446]	[1.03293, 1.0524]
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

Verificamos si aún contamos con sobredispersión en los nuevos modelos

# poisson case
performance::check_overdispersion(reg_poisson)

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

## Overdispersion detected.

# quasipoisson case
performance::check_overdispersion(reg_quasipoisson)

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

## Overdispersion detected.

# negative binomial case
performance::check_overdispersion(reg_binomialnegativa)

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

## No overdispersion detected.

Comparamos los tres modelos con ANOVA

anova(reg_poisson,reg_quasipoisson,reg_binomialnegativa,test = "Chisq") %>%
kable(caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)

Tabla ANOVA para comparar modelos
Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
178	736313.3944	NA	NA	NA
178	736313.3944	0	0.0	NA
178	188.0349	0	736125.4	NA

#Interpretación: Como señalamos anteriormente, nuestro modelo presenta sobredispersión y, tras el análisis de quasipoisson y la binomial negativa, llegamos a la conclusión que el mejor modelo para este caso resulta ser el de la regresión Binomial Negativa.

3. Regresión Logística

En este caso empleamos la misma base de datos que el ejercicio 1; sin embargo, debemos crear la variable dicotómica

data_Ceplan$VULNERABLE <- ifelse(data_Ceplan$IVIA >= 0.5, 1, 0)

str(data_Ceplan)

## '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 ...

Hipótesis: El ser vulnerable a la inseguridad alimentaria 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).

hp3 = formula (VULNERABLE ~ ABASTOS + IDH + IDE + POBREZA)
reg_log = glm(hp3, data = data_Ceplan, family = binomial)
summary(reg_log)

## 
## Call:
## glm(formula = hp3, family = binomial, data = data_Ceplan)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.99122  -0.05640   0.01882   0.15116   2.86668  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  23.95841    7.70805   3.108 0.001882 ** 
## ABASTOS       0.02974    0.07392   0.402 0.687477    
## IDH         -43.27762   11.13137  -3.888 0.000101 ***
## IDE         -10.07738    7.00832  -1.438 0.150458    
## POBREZA       0.07304    0.05250   1.391 0.164153    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 251.716  on 182  degrees of freedom
## Residual deviance:  56.965  on 178  degrees of freedom
## AIC: 66.965
## 
## Number of Fisher Scoring iterations: 8

Lo vemos mejor

#install.packages("modelsummary")
library(modelsummary)
modelrl= list("Ser vulnerable"= reg_log)
modelsummary(modelrl, title = "Regresión Logística", stars = TRUE, output = "kableExtra")

Regresión Logística
	Ser vulnerable
(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

Intepretación de los coeficientes: En la tabla se señala que solo una de las cuatro variables empleadas en la regresión resulta significativa: IDH.

Evaluamos el modelo

data_Ceplan$predictedProbs <- predict(reg_log, data_Ceplan,type = "response")
head(data_Ceplan$predictedProbs)

## [1] 0.546590568 0.795312840 0.003943837 0.999997851 0.996658123 0.977952660

data_Ceplan$predicted=ifelse(data_Ceplan$predictedProbs>0.5,"1","0")
head(data_Ceplan[,c("predicted","VULNERABLE")],10)

##    predicted VULNERABLE
## 1          1          1
## 2          1          1
## 3          0          0
## 4          1          1
## 5          1          1
## 6          1          1
## 7          0          1
## 8          1          1
## 9          0          0
## 10         1          1

Asimismo, creamos una matriz de confusión para saber qué tanto acertamos

library(cvms)
ConfMatrix=confusion_matrix(predictions =  data_Ceplan$predicted,
                      targets= data_Ceplan$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

KAPPA

ConfMatrix$Kappa

## [1] 0.8783318

El índice Kappa nos demuestra que la predicción está más cerca del 1 (0.87) que del 0, por lo que encontramos que la mayoría de predicciones no fallan.

library(margins)
library(dplyr)
library(kableExtra)
margins(reg_log)

## Average marginal effects

## glm(formula = hp3, family = binomial, data = data_Ceplan)

##   ABASTOS    IDH     IDE  POBREZA
##  0.001414 -2.058 -0.4792 0.003473

marginalsData=summary(margins(reg_log))
marginalsData%>% kable(caption = "Efectos Marginales Promedio (AME)") %>%kableExtra::kable_styling(full_width = T)

Efectos Marginales Promedio (AME)
factor	AME	SE	z	p	lower	upper
ABASTOS	0.0014139	0.0035122	0.402578	0.6872587	-0.0054699	0.0082977
IDE	-0.4791767	0.3241919	-1.478065	0.1393904	-1.1145812	0.1562278
IDH	-2.0578376	0.3592354	-5.728382	0.0000000	-2.7619260	-1.3537492
POBREZA	0.0034730	0.0024163	1.437293	0.1506349	-0.0012629	0.0082089

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

Interpretación margins: Como se puede apreciar, solo el intervalod econfianza del IDH no incluye al 0, por lo que comprendemos que resulta ser el único covariado significativo. Asimismo, como identificamos en los margins, el IDH posee un efecto significativo y negativo en torno a la probabilidad de encontrarse en vulnerable a la inseguridad ailmentaria, por cada punto que aumenta el IDH, la probabilidad de ser vulnerable disminuye en 2.05.

4. Regresión Cox

data_Chile = import("MinistrosChile.csv")

str(data_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\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.

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

Creamos survival

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

library(ggplot2)
library(ggfortify)

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

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

Dinámica de sobrevivencia según sexo

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

## '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 = data_Chile)

###
ejeX='DAYS\n Curva cae cuando sale un ministro'
ejeY='Probabilidad \n(Sobrevivir 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"))

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

## [1] 0.8969547

Interpretación: como se observa en el gráfico del análisis KM y como señala el logrank (ya que 0.89 supera el 0.05 requerido), encontramos que la variable sexo no posee un efecto significativo en la dinámica de sobrevivencia.

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

Regresión Cox

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

## Call:
## coxph(formula = COX_H1, data = data_Chile)
## 
##   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

Modelos sin exponenciar y exponenciado

modelcox=list('Riesgo - Ser destituido del cargo'=reg_cox,'Riesgo- Ser destituido del cargo (exponenciado)'=reg_cox)
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 del cargo	Riesgo- Ser destituido del cargo (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(reg_cox)[4])))

## presidentPatricio Aylwin 
##                0.7621441

Interpretación: Según la tabla, identificamos que solo la variable “presidentePatricioAylwin” resulta significativa, por lo que el riesgo de ser despedido del cargo de minisitro durante su gobierno es en promedio 76% menor en torno a los minisitros que no fueron parte de su gobierno.

Práctica Laboratorio 5 - Alfredo

Alfredo Ludmir Aro Terleira

2024-05-01