1.1.1.1 Cargamos los datos

Cargamos los datos, como están en GoogleDrive los podemos leer desde ahí:

rm(list = ls()) # limpiar el working environment

linkADrive='https://docs.google.com/spreadsheets/d/e/2PACX-1vS2ZSNM8BIZtoufVTO4Mw3ZmTWW1rAAtsGzFg0shHTJXX-3GmtLsgU-Nqkw5RzDgrNX31GTC9L7LnEz/pub?gid=0&single=true&output=csv'

estadoLocales=read.csv(linkADrive)

head(estadoLocales)

##   DEPARTAMENTO UBIGEO buenEstado contribuyentes_sunat peaOcupada TOTALpob
## 1     AMAZONAS  10000       18.6                75035     130019   417365
## 2       ÁNCASH  20000       13.9               302906     387976  1139115
## 3     APURÍMAC  30000        8.7               103981     140341   424259
## 4     AREQUIPA  40000       27.4               585628     645001  1460433
## 5     AYACUCHO  50000       17.0               151191     235857   650940
## 6    CAJAMARCA  60000       18.0               277457     461312  1427527
##    URBANO  RURAL
## 1  205976 211389
## 2  806065 333050
## 3  243354 180905
## 4 1383694  76739
## 5  444473 206467
## 6  567141 860386

La vista preliminar nos muestra varias columnas, pero la de nuestro interés es la última columna de la derecha. A esta altura, antes de explorar los datos de manera estadística, debemos primero verificar cómo R ha intepretado el tipo de datos:

str(estadoLocales)

## 'data.frame':    25 obs. of  8 variables:
##  $ DEPARTAMENTO        : chr  "AMAZONAS" "ÁNCASH" "APURÍMAC" "AREQUIPA" ...
##  $ UBIGEO              : int  10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 ...
##  $ buenEstado          : num  18.6 13.9 8.7 27.4 17 18 33.8 11.9 10.1 15.6 ...
##  $ contribuyentes_sunat: int  75035 302906 103981 585628 151191 277457 499257 466883 80353 185658 ...
##  $ peaOcupada          : int  130019 387976 140341 645001 235857 461312 445072 496399 111996 253200 ...
##  $ TOTALpob            : int  417365 1139115 424259 1460433 650940 1427527 1046953 1315220 367252 759962 ...
##  $ URBANO              : int  205976 806065 243354 1383694 444473 567141 1046953 865771 166163 447620 ...
##  $ RURAL               : int  211389 333050 180905 76739 206467 860386 0 449449 201089 312342 ...

Nuestra variable de interés, buenEstado, es de tipo numérico y R también lo ha interpretado así:

1.1.1.2 Tablas, Gráficos y Estadígrafos

Si los datos están en el tipo adecuado, puede iniciarse la exploración estadística.

Tenemos 25 observaciones a nivel departamental. La manera más rápida de ver los estadígragos es usando el comando summary()

summary(estadoLocales$buenEstado)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.70   13.90   17.30   18.61   21.80   40.50

Ahora, podemos nos faltarían algunos estadígrados, que podemos verlo así:

library(DescTools)

## Warning: package 'DescTools' was built under R version 4.4.2

allStats=c(summary(estadoLocales$buenEstado),
  sd=sd(estadoLocales$buenEstado),
  skew=Skew(estadoLocales$buenEstado),
  kurt=Kurt(estadoLocales$buenEstado),
  cv=CoefVar(estadoLocales$buenEstado))
allStats

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max.         sd 
##  7.7000000 13.9000000 17.3000000 18.6120000 21.8000000 40.5000000  8.2596065 
##       skew       kurt         cv 
##  0.9008641  0.2026850  0.4437786

A esta altura, los gráficos que necesitamos son el histograma:

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.2

base=ggplot(data=estadoLocales,
            aes(x=buenEstado))
histogram= base + geom_histogram(aes(y = after_stat(density)),
                 colour = 1, fill = "white",bins=10) +  
    stat_function(fun = dnorm,
                  args = list(mean = allStats['Mean'],
                              sd = allStats['sd']),col='red')
    
histogram

Y el boxplot

base=ggplot(data=estadoLocales,
            aes(y=buenEstado))
boxplot=base + geom_boxplot()

boxplot

A esta altura solo falta identificar los atípicos del boxplot, lo cual podemos recuperar con ggplot_build:

data_boxLocales=ggplot_build(boxplot)$data[[1]]
data_boxLocales

##   ymin lower middle upper ymax   outliers notchupper notchlower x flipped_aes
## 1  7.7  13.9   17.3  21.8 31.4 33.8, 40.5    19.7964    14.8036 0       FALSE
##   PANEL group ymin_final ymax_final   xmin  xmax xid newx new_width weight
## 1     1    -1        7.7       40.5 -0.375 0.375   1    0      0.75      1
##   colour  fill alpha shape linetype linewidth
## 1 grey20 white    NA    19    solid       0.5

Los outliers están en una lista, por lo que debemos escribir

data_boxLocales$outliers

## [[1]]
## [1] 33.8 40.5

Nota la utilidad de los “bigotes” del boxplot, cuando hay atípicos:

data_boxLocales[c('ymin','ymax')]

##   ymin ymax
## 1  7.7 31.4

En este caso, sabemos que los valores que exceden 31.4 serán considerados atípicos.

Hasta aquí, podrías sustentar que los valores de la variable buen estado de colegios públicos a nivel regional se distribuyen asimétricamente, con una dispersión baja, y con presencia de valores atípicos altos.

1.1.2.1 Numérica - Numérica

Consideramos que nos interesa explorar la posible relación entre nuestra variable de interés y la PEA ocupada. Traigamos esa variable. Asimismo, como son dos variables de tipo numérico, la estrategia a seguir es el análisis de correlación. Veamos este scatterplot:

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.4.2

base=ggplot(data=estadoLocales, aes(x=peaOcupada, y=buenEstado))
scatter = base + geom_point()
scatterText = scatter + geom_text_repel(aes(label=DEPARTAMENTO),size=2)
scatterText

Calculemos ahora los índices de correlación:

f1=formula(~peaOcupada + buenEstado)

pearsonf1=cor.test(f1,data=estadoLocales)[c('estimate','p.value')]
pearsonf1

## $estimate
##       cor 
## 0.3567442 
## 
## $p.value
## [1] 0.08002776

spearmanf1=cor.test(f1,data=estadoLocales,method='spearman',exact=F)[c('estimate','p.value')]
spearmanf1

## $estimate
##       rho 
## 0.2785151 
## 
## $p.value
## [1] 0.1776146

1.1.2.2 Numérica - Categórica

Creamos aquí una variable categórica (factor) que represente “populoso” como tener más de un millón de habitantes

library(magrittr)
estadoLocales$populoso=ifelse(estadoLocales$TOTALpob>1000000,'Si','No')%>%as.factor()

# veamos conteo de cada categoría
table(estadoLocales$populoso)

## 
## No Si 
## 14 11

Pidamos un boxplot, pero en este caso por grupos:

base=ggplot(data=estadoLocales, aes(x=populoso, y=buenEstado))
base + geom_boxplot(notch = T) +  geom_jitter(color="black", size=0.4, alpha=0.9)

## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?

Los boxplots tienen un notch flanqueando a la mediana, para sugerir igualdad de medianas si éstos se intersectan; de ahí que parece no haber diferencia sustantiva entre las categorías.

Verificar si hay o no igualdad entre distribuciones depende si las variables se distribuyen o no de manera normal por grupo. Como ello no es fácil de discernir visualmente, complementemos el análisis con los test de normlalidad. Usemos Shapiro-Wilk en cada grupo:

f2=formula(buenEstado~populoso)


tablag= aggregate(f2, estadoLocales,
          FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})




shapiroTest=as.data.frame(tablag[,2])
names(shapiroTest)=c("W","Prob")
shapiroTest

##           W        Prob
## 1 0.8126339 0.007182561
## 2 0.9583503 0.750951875

Para que se vea mejor en html:

library(knitr)

## Warning: package 'knitr' was built under R version 4.4.2

library(magrittr)
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.4.3

kable(cbind(tablag[1],shapiroTest))%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F, position = "left")

Habría normalidad en un grupo (populoso sí) y no en otro (populoso no). Usemos entonces tanto la prueba de Mann-Whitney (no paramétrica) como la prueba t (paramétrica) para analizar diferencias

(student_T=t.test(f2,data=estadoLocales)[c('estimate','p.value')])

## $estimate
## mean in group No mean in group Si 
##             17.6             19.9 
## 
## $p.value
## [1] 0.5025783

(Mann_Whitnery=wilcox.test(f2,data=estadoLocales,exact=F)['p.value'])

## $p.value
## [1] 0.3960432

1.2.1.1 Preparando la data

#carga de data
linkToData='https://docs.google.com/spreadsheets/d/e/2PACX-1vSRpCC8gKIxMxpK0wjgLcl-GQWdw6sAeB16Sixkq6kZXeTfMS8_n70twbEbQ2tMcGp8tm-6x8qf8ieo/pub?gid=234740779&single=true&output=csv'
pavi=read.csv(linkToData)
str(pavi)

## 'data.frame':    1096 obs. of  9 variables:
##  $ municipio       : chr  "m1" "m2" "m3" "m4" ...
##  $ apropiaciondolar: num  102.17 4.19 1.59 0 0 ...
##  $ consejocomunal  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ejecucion       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ uribista        : int  0 1 1 1 1 1 1 1 1 NA ...
##  $ pctopo          : num  14.82 14.51 15.08 6.15 47.31 ...
##  $ priorizado      : int  0 1 0 0 0 0 0 0 1 0 ...
##  $ poblacioncienmil: num  20.9155 0.2406 0.0398 0.0558 0.2723 ...
##  $ nbi             : num  12.2 33.8 28.5 33.1 27.1 ...

Recordar que es mejor usar datos sin valores perdidos en los análisis multivariados

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

Los datos debemos darle el formato adecuado. En este caso hay que cambiar a categóricos algunas variables:

seleccion=c("consejocomunal","ejecucion","uribista","priorizado")
pavi[,seleccion]=lapply(pavi[,seleccion],as.factor)

1.2.1.2 Explorando las variables

paviStats=summary(pavi[,-1])
paviStats

##  apropiaciondolar  consejocomunal ejecucion uribista     pctopo      
##  Min.   :  0.000   0:827          0:844     0:325    Min.   : 0.000  
##  1st Qu.:  0.000   1: 49          1: 32     1:551    1st Qu.: 6.285  
##  Median :  0.000                                     Median :19.607  
##  Mean   :  8.926                                     Mean   :27.306  
##  3rd Qu.: 11.436                                     3rd Qu.:44.455  
##  Max.   :132.643                                     Max.   :99.419  
##  priorizado poblacioncienmil        nbi       
##  0:639      Min.   : 0.00629   Min.   : 6.84  
##  1:237      1st Qu.: 0.07460   1st Qu.:27.33  
##             Median : 0.14593   Median :40.26  
##             Mean   : 0.45639   Mean   :41.72  
##             3rd Qu.: 0.27433   3rd Qu.:53.78  
##             Max.   :69.26836   Max.   :97.32

¿Qué puedes decir de estos resultados?

library(ggcorrplot)

## Warning: package 'ggcorrplot' was built under R version 4.4.3

colNums=names(pavi)[c(2,6,8,9)]  #las columnas numéricas
numXs=pavi[,colNums]
ggcorrplot(cor(numXs),lab = T,show.diag = F)

library(magrittr)
library(kableExtra)

colCats=setdiff(names(pavi), colNums)[-1]
diffPara=c()
diffNoPara=c()

for (col in colCats){
    diffPara=c(diffPara,t.test(pavi[,"apropiaciondolar"]~pavi[,col])['p.value']<=0.05)
    diffNoPara=c(diffNoPara,wilcox.test(pavi[,"apropiaciondolar"]~pavi[,col])['p.value']<=0.05)
}
data.frame(cbind(colCats,diffPara,diffNoPara),
           row.names = 1:length(colCats))%>%
           kable(caption = "Diferencia de 'VD:Apropiacion' por Grupo")%>%
            kableExtra::kable_styling(full_width = FALSE)

par(mfrow = c(2, 2))  

for (col in colCats) { 
  boxplot(pavi$apropiaciondolar~pavi[,col],
       main = col,xlab="",ylab = "")
}

1.2.2.1 Hipótesis 1

Una hipótesis: - El Beneficio recibido en un municipio ha sido afectado por el porcentaje de votos recibidos por los candidatos de la oposición a Uribe a la cámara de representantes, controlado por tamaño de población.

# hipotesis en R
modelo1=formula(apropiaciondolar~pctopo+poblacioncienmil)

Probamos la hipótesis

reg1=lm(modelo1,data=pavi)
summary(reg1)

## 
## Call:
## lm(formula = modelo1, data = pavi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.809  -8.671  -7.448   2.602  90.426 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       8.80319    0.79976  11.007   <2e-16 ***
## pctopo           -0.03146    0.02150  -1.463    0.144    
## poblacioncienmil  2.15125    0.20073  10.717   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.85 on 873 degrees of freedom
## Multiple R-squared:  0.1181, Adjusted R-squared:  0.1161 
## F-statistic: 58.45 on 2 and 873 DF,  p-value: < 2.2e-16

El resultado anterior puede ser presentado de mejor manera usando modelsummary package, como se puede ver en la Tabla 2.1 (el error típico está entre parentesis).

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')
##   options(modelsummary_factory_latex = 'kableExtra')
##   options(modelsummary_factory_html = 'kableExtra')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)

## 
## Adjuntando el paquete: 'modelsummary'

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

model1=list('apropiacion (I)'=reg1)
modelsummary(model1, title = "Regresion: modelo 1",
             stars = TRUE,
             output = "kableExtra")

Al probar esta hipótesis vemos…

1. que pctopo tiene signo negativo (relación inversa con la VD)
2. que la magnitud de ese efecto es -0.031, lo que indica cuanto varía apropiaciondolar en promedio cuando pctopo se incrementa en una unidad, controlando por poblacioncienmil
3. que pctopo NO tiene efecto significativo

Justamente el R cuadrado ajustado (0.1180881) nos brinda un porcentaje (multiplicalo por 100) que da una puesta de nuestra cercanía a una situación perfecta (cuando vale 1).

1.2.2.2 Hipótesis 2

¿Y si queremos ver el efecto de consejo comunal (consecocomunal)?

modelo2=formula(apropiaciondolar~pctopo+consejocomunal+poblacioncienmil)

Esta nueva hipótesis desea evaluar si la vista de Uribe a un Consejo Comunal influye en la asignación de presupuesto. El resultado de este proceso lo vemos a continuación:

reg2=lm(modelo2,data=pavi)
summary(reg2)

## 
## Call:
## lm(formula = modelo2, data = pavi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -62.378  -7.933  -6.834   1.573  91.247 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       8.04119    0.79129  10.162  < 2e-16 ***
## pctopo           -0.02987    0.02103  -1.420    0.156    
## consejocomunal1  14.79601    2.32127   6.374 2.98e-10 ***
## poblacioncienmil  1.91236    0.19987   9.568  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.51 on 872 degrees of freedom
## Multiple R-squared:  0.1573, Adjusted R-squared:  0.1545 
## F-statistic: 54.28 on 3 and 872 DF,  p-value: < 2.2e-16

model2=list('apropiacion (II)'=reg2)
modelsummary(model2, title = "Regresion: modelo 2",
             stars = TRUE,
             output = "kableExtra")

Al probar esta hipótesis vemos que… - pctopo tiene signo negativo; NO tiene efecto significativo; y la magnitud de ese efecto es -0.03, lo que indica cuando variaría apropiaciondolar en promedio cuando pctopo se incremente en una unidad, controlando por las demás variables. - consejocomunal SÍ tiene efecto significativo al 0.001: ese efecto es directo, pues el coeficiente calculado es positivo; y la magnitud de ese efecto es 14.796, lo que indica cuando varía apropiaciondolar en promedio cuando consejocomunal es 1 y no 0; también controlando por las demas variables.

Nótese la lectura del efecto cuando la variable independiente es categórica (o factor). Primero, nota que R indica el valor 1 con el nombre de la variable: eso indica que el valor 0 (cuando el consejocomunal de ese municipio NO fue visitado) es la categoría de referencia; es decir, el coeficiente nos indica cuanto se modifica el valor promedio de la dependiente cuando se pasa de 0 a 1. Si la variable independiente es politómica (no ordinal), aparecerá cada categoría menos la de referencia, y el efecto siempre se debe interpretar como el efecto de las variables mostrada versus la de referencia.

1.2.2.3 Hipótesis 3

# le agregamos uribista
modelo3=formula(apropiaciondolar~pctopo+consejocomunal+uribista+poblacioncienmil)

reg3=lm(modelo3,data=pavi)
summary(reg3)

## 
## Call:
## lm(formula = modelo3, data = pavi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.004  -7.546  -6.482   1.847  91.958 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       9.80542    1.11705   8.778  < 2e-16 ***
## pctopo           -0.03753    0.02126  -1.765   0.0779 .  
## consejocomunal1  14.73795    2.31613   6.363 3.19e-10 ***
## uribista1        -2.45126    1.09801  -2.232   0.0258 *  
## poblacioncienmil  1.89052    0.19965   9.469  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.47 on 871 degrees of freedom
## Multiple R-squared:  0.1621, Adjusted R-squared:  0.1583 
## F-statistic: 42.14 on 4 and 871 DF,  p-value: < 2.2e-16

model3=list('apropiacion (III)'=reg3)
modelsummary(model3, title = "Regresion: modelo 3",
             stars = TRUE,
             output = "kableExtra")

Al probar esta hipótesis vemos que…

• pctopo tiene signo negativo; y ahora SÍ tiene efecto significativo; y la magnitud de ese efecto es -0.038, lo que indica cuanto variaría apropiaciondolar en promedio cuando pctopo se incrementa en una unidad, controlando por las demás variables.

• consejocomunal MANTIENE efecto significativo al 0.001; ese efecto es directo, pues el coeficiente calculado es positivo; y la magnitud de ese efecto es 14.738, lo que indica cuanto varía apropiaciondolar en promedio cuando consejocomunal es 1 y no 0, también controlando por las demás variables.

• uribista tiene efecto significativo al 0.05; ese efecto es inverso, pues el coeficiente calculado es negativo; y la magnitud de ese efecto es -2.451, lo que indica cuanto varía apropiaciondolar en promedio cuando uribista es 1 y no 0, también controlando por las demás variables

1.2.2.4 Estandarización de Coeficientes

Del resultado de la última tabla NO podemos directamente decir que consejo comunal tiene más efecto que los demás por el solo hecho de que el valor estimado sea mayor al resto. Para saber cuál tiene más efecto, cuando los predictores tienen, como en este caso, unidades diferentes, estandarizamos los datos y volvemos a correr la regresión.

modelo3_st=formula(scale(apropiaciondolar)~scale(pctopo)+scale(as.numeric(consejocomunal))+scale(as.numeric(uribista))+scale(poblacioncienmil))

modelo3_st=lm(modelo3_st,data=pavi)

modelo3_st=list('apropiacion (III_st)'=modelo3_st)
modelsummary(modelo3_st, title = "Regresion: modelo 3 con \ncoeficientes estandarizados",
             stars = TRUE,
             output = "kableExtra")

Interpretación: De los resultados de este tabla se despeja que, dejando de lado la variable de control (población por cien mil), la que tiene mayor efecto es consejo comunal, algo que no era evidente en la anterior tabla. Esto nos indica cuantas desviaciones estándar varía la variable dependiente (apropiacion dolar) cuando la dependiente varía en una (1) desviación estándar. Nota además que hemos tenido que estandarizar variables categóricas, lo cuál NO es posible. Ello no debe preocupar pues esta etapa es sólo para comparar tamaño de efecto entre variable.

Podemos simplificar los pasos si utlizamos funciones del paquete lb.beta. La función que lleva el mismo nombre nos da el resultados rapidamente

library(lm.beta)

## Warning: package 'lm.beta' was built under R version 4.4.2

model3beta=list('apropiacion (III)'=lm.beta(reg3))
modelsummary(model3beta, title = "Regresion: modelo 3 con \ncoeficientes estandarizados usando lm.beta()",
             stars = TRUE,
             output = "kableExtra")

1.2.3.1 Tabla ANOVA para comparar modelos

En todos los modelos, el error no tiene coeficiente, representamos su variación usando el error típico de los residuos o residual standard error (RSE). Nótese que este ha variado de un modelo ha otro, y en el último modelo se tiene un R2Adj mayor (el menor RSE). Aquí vale la pena preguntarse si esta disminución del error es significativa:

library(magrittr)
library(knitr)
tanova=anova(reg1,reg2,reg3)

kable(tanova,
      caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)

La comparación dem odelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasado de un modelo a otro). Cuando la comparación es significativa (observar el Pr(<F)), rechazamos igualdad entre modelos: el modelo 2 presenta menos error que el modelo 1, pero el modelo 2, sí presenta mayot error que el modelo 3, por lo que el modelo 3 debe ser elegido (hasta aquí).

1.2.4.1 Linealidad

Se asume relación lineal entre la variable dependiente (Y) y las independientes Xs. Para ello analizamos la relación entre los residuos y los valores que predice el modelo de regresión.

# linea roja debe tender a horizontal
plot(reg2, 1)

La falta de linealidad provocaría que el modelo no sirva para explicar las mismas variables con datos diferentes en otros estudios.

1.2.4.2 Homocedasticidad

Se asume que la dispersión de los errores de la estimación (apropiaciondolar) mantiene una variación homogénea. Si analizamos las raíces de los errores estandarizados versus los valores estimados, la distancia de los puntos a la linea de referencia debe ser similar.

# linea roja debe tender a horizontal
plot(reg2, 3)

library(lmtest)

## Cargando paquete requerido: zoo

## 
## Adjuntando el paquete: 'zoo'

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

library(kableExtra)
# null: modelo homocedastico
resBP=bptest(reg2)
data.frame(list('BP'=resBP$statistic,
             'df'=resBP$parameter,
             "p-value"=resBP$p.value))%>%
    kable(caption = resBP$method)%>%kable_styling(full_width = F)

1.2.4.3 Normalidad de residuos

Los residuos deben distribuirse de manera normal. Un qq-plot nos permite relevar si hay normalidad

# puntos cerca a la diagonal?
plot(reg2, 2)

#NULL: Datos se distribuyen de manera normal
resSW=shapiro.test(reg2$residuals)
data.frame(list('SW'=resSW$statistic,
             "p-value"=resSW$p.value))%>%
    kable(caption = resSW$method)%>%kable_styling(full_width = F)

1.2.4.4 No Multicolinealidad

Si los predictores tienen una correlación muy alta entre sí, hay multicolinealidad:

• Mayor a 5 es problemático

library(DescTools)

# > 5 es problematico
VIF(reg2) %>%kable(col.names = "VIF",caption ="Evaluando Multicolinealidad usando VIF (Variance Inflation Factors)" )%>%kable_styling(full_width = F)

La presencia de la multicolinealidad no perjudica tanto al calculo, pero evitar calcular bien el efecto de cada regresor.

1.2.4.5 Valores Influyentes

Hay casos particulares, que tienen la capacidad de trastocar lo que el modelo representa. A veces detectándolos y suprimiéndolos podemos ver un mejor modelo. La siguiente figura muestra la posible existencia de valores influyentes.

plot(reg2, 5)

Si queremos ver más de cerca lo posible los casos influyentes, le prestamos atención al índice de Cook y a los valores predeciros (hat values) de la siguiente tabla:

checkReg2=as.data.frame(influence.measures(reg2)$is.inf)
checkReg2[checkReg2$cook.d & checkReg2$hat,c('cook.d','hat')]%>%kable(caption = "Valores Influyentes criticos")%>%kable_styling(full_width = F)

Si alguna fila aparece en el resultado, en ese caso está afectando los cálculos de la regresión (sin él habría otro resultado)

2.1.3.1 Hipótesis replanteada

Reimplementemos nuestra hipótesis original usando la regresión Poisson. Los resultados lso vemos en la siguiente tabla, comparándolos con el resultado de la regresión lineal controlada por la población:

rp1=glm(h1, data = censo2007_wide_Totales, 
        offset=log(poblacion_Total), #exposure 
        family = poisson(link = "log"))

# displaying results
modelslmpoi=list('OLS asegurados (II)'=rl2,
                 'POISSON asegurados'=rp1)

modelsummary(modelslmpoi, title = "Regresiones OLS y Poisson",
             stars = TRUE,
             output = "kableExtra")

Ahora que tenemos ambas regresiones, vemos que el modelo Poisson le devuelve efecto a la independiente. Nótese que la Poisson está modelando los conteos, teniendo en cuenta la exposición (exposure), añadida usando offset. Esto no siempre es necesario, pero en este caso sí lo era, pues necesitamos controlar la exposure de manera explícita (población). Sin embargo, el AIC de la Poisson es aún mayor que el de la Gaussiana.

2.1.3.2 Interpretación

Ya que sabemos cuándo usar una regresión Poisson, alteremos la hipótesis anterior:

A nivel distrital, la cantidad de personas con algún seguro de salud está afectada por el nivel de analfabetismo y por la presencia de trabajadores independientes

h2=formula(conSeguro_Total~analfa_Total + indep_Total)
    
rp2=glm(h2, data = censo2007_wide_Totales, 
        offset=log(poblacion_Total),
        family = poisson(link = "log"))


modelsPois=list('POISSON asegurados (I)'=rp1,
                'POISSON asegurados (II)'=rp2)
modelsummary(modelsPois, 
             title = "Regresiones Poisson anidadas",
             stars = TRUE,
             output = "kableExtra")

La interpretación de la anterior tabla NO es tan directa como lo era en la regresión lineal. Demos una interpretación simple, sin hacer ningún cálculo, de los resultados del modelo II:

• En el modelo II, ambos predictores son significativos.
• En el modelo II, a mayor analfabetismo, mayor cantidad de asegurados.
• En el modelo II, a mayor cantidad de trabjadores independientes, menor cantidad de asegurados.

Sin embargo, la interpretación precisa de los coeficientes de la tabla requiere más cálculos; tengamos primero en cuenta lo que la regresión Poisson ha calculado usando la ecuación. Cuando el exposure es constante modelamos conteos (Y); cuando no lo es modelso ratios. Pero, como vemos en la Ecuación, los coeficientes necesitan ser exponenciados para saber el efecto sobre Y

Veamos la siguiente tabla:

# formula para limitar a 4 digitos decimales, 
# sin que se muestre notación científica:

formatoNum <- function(x) format(x, digits = 4, scientific = FALSE)

modelsummary(modelsPois,
             fmt=formatoNum, # uso mi formula
             exponentiate = T, # exponenciar!!!!!
             statistic = 'conf.int',
             title = "Regresión Poisson - coeficientes exponenciados",
             stars = TRUE,
             output = "kableExtra")

#más simple:
#cbind(exp(coef(rp2)),exp(confint(rp2)))

La tabla anterior tiene los coeficientes exponenciados, así mismo, muestra los intervalos de confianza (exponenciados) en vez de los errores típicos. Nota que mientras en la regresión lineal no deseábamos que nuestro coeficiente esté cerca al cero, es decir, que su intervalo de confianza no incluya al cero, aquí no deseamos que el intervalo de confianza incluya al uno

Una vez exponenciado, podemos interpretar los coeficientes de manera más sencilla. Así, para el modelo II se ve que:

• por cada unidad que aumente el analfabetismo, la cantidad esperada de asegurados se multiplicar por 1.016, es decir, aumentaría en un 1.6% (100 x |1-1.016|) en promedio.

• por cada unidad que aumente los trabajadores independientes, la cantidad esperada de asegurados se multiplica por 0.99, es decir, disminuiría en 1% (100 x |1-0.99|) en promedio.

Nótese que la regresión propone un efecto multiplicativo sobre el promedio de la respuesta (la regresión OLS o Gaussiana propone un efecto aditivo).

2.1.4.1 Quasi Poisson

Útil para subdispersión y sobredispersión.

rqp = glm(h2, data = censo2007_wide_Totales,
          offset=log(poblacion_Total),
          family = quasipoisson(link = "log"))

modelsPQP=list('POISSON asegurados (II)'=rp2,'QUASIPOISSON asegurados (II)'=rqp)

modelsummary(modelsPQP, 
             title = "Regresiones Poisson y QuasiPoisson",
             stars = TRUE,
             output = "kableExtra")

La tabla nos muestra cosas interesantes:

• Los coeficientes son los mismos para ambos modelos:

cbind(coef_Poi=coef(rp2),coef_QuasiPoi=coef(rqp))

##                 coef_Poi coef_QuasiPoi
## (Intercept)  -0.61523114   -0.61523114
## analfa_Total  0.01649913    0.01649913
## indep_Total  -0.01073999   -0.01073999

• Pero no los errores típicos

library(arm)

## Warning: package 'arm' was built under R version 4.4.2

## Cargando paquete requerido: MASS

## Cargando paquete requerido: Matrix

## Warning: package 'Matrix' was built under R version 4.4.2

## Cargando paquete requerido: lme4

## 
## Adjuntando el paquete: 'lme4'

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

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

## Working directory is C:/Users/USUARIO/Documents/UNIVERSIDAD/Manuales/Estadística2

cbind(se_Poi=se.coef(rp2),se_QuasiPoi=se.coef(rqp))

##                    se_Poi  se_QuasiPoi
## (Intercept)  8.910449e-04 0.0182152135
## analfa_Total 4.609011e-05 0.0009421985
## indep_Total  3.037210e-05 0.0006208826

La regresión quasipoisson lidia con la sobredispersión al recalcular los errores típicos, lo que afectaría la significancia de los predictores; de ahí que calcula nuevos intervalos de confianza:

#Exponenciamos!!!

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


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

2.1.4.2 Binomial Negativa

Solo resulta útil para la sobredispersión

h2off=formula(conSeguro_Total ~ analfa_Total + indep_Total + offset(log(poblacion_Total)))

rbn=glm.nb(h2off,data=censo2007_wide_Totales)

modelsQP_BN=list('Poisson asegurados (II)'=rp2,
                 'QuasiPoisson asegurados (II)'=rqp,
                 'Binomial Negativa asegurados (II)'=rbn)


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

Notamos en la tabla que los coeficientes obtenidos en la regresión binomial negativa son diferentes a los demás. Además, los AIC son mucho mejores para caso de la Binomial Negativa; los AIC/BIC podemos verlos así de manera alternativa

data.frame(Model = c("poisson", "quasipoisson", "negative-binomial"),
           AIC = c(AIC(rp2), AIC(rqp), AIC(rbn)),
           BIC = c(BIC(rp2), BIC(rqp), BIC(rbn)),stringsAsFactors = FALSE
)%>%
kable(caption = "AIC / BIC para los modelos de conteos")%>%kableExtra::kable_styling(full_width = FALSE)

2.1.5.1 Estandarización de coeficientes

Finalmente, podemos calcular los coeficientes estandarizados para saber cuál de los predictores tiene mayor efecto

sdVD=sd(censo2007_wide_Totales$conSeguro_Total)
sdVIs=apply(censo2007_wide_Totales[,c("analfa_Total","indep_Total")],2,sd)
DF=list(Poisson=sdVIs*coef(rp2)[c(2,3)]/sdVD,
     CuasiPoisson=sdVIs*coef(rqp)[c(2,3)]/sdVD,
     BinomNegativa=sdVIs*coef(rbn)[c(2,3)]/sdVD)%>%
       data.frame()

DF%>% kable(caption = "Coeficientes Standarizados (ordenar vía valores absolutos)")%>%
          kableExtra::kable_styling(full_width = F)

2.2.2.1 Probabilidad y Odds para caso de la mujer:

La probabilidad (entre 0 y 1) que una mujer sea voluntaria es la división del ser voluntaria entre el total de mujeres:

ProbMujVol=volsexTable[2,1]/sum(volsexTable[,1])
MASS::fractions(ProbMujVol)

## [1] 349/780

Es decir:

ProbMujVol

## [1] 0.4474359

Representemos el odds que sea voluntaria: la división entre ser o no ser voluntaria

OddsMujVol=volsexTable[2,1]/volsexTable[1,1]
MASS::fractions(OddsMujVol)

## [1] 349/431

Es decir:

OddsMujVol

## [1] 0.8097448

El odds suele representarse además como la razón entre dos probabilidades: la probabilidad que SÍ ocurra un evento dividido por la probabilidad que NO ocurra ese evento:

ProbMujVol/(1-ProbMujVol)

## [1] 0.8097448

2.2.2.2 Probabilidad y Odds para caso del hombre

La probabilidad que un hombre sea voluntario

ProbHomVol=volsexTable[2,2]/sum(volsexTable[,2])
MASS::fractions(ProbHomVol)

## [1] 248/641

Es decir:

ProbHomVol

## [1] 0.3868955

Ahora, calculamos los odds:

OddsHomVol=ProbHomVol/(1-ProbHomVol)
OddsHomVol

## [1] 0.6310433

2.2.2.3 Comparando Mujer y Hombre

Hasta aquí sabemos cuál odds ha salido mayor. Eso es información clara. Para precisar esa información, podemos dividir los odds, lo que nos da el odds ratio, que en este caso nos dirá qué diferencia produce el sexo en el voluntario:

(OR_MujHom=OddsMujVol/OddsHomVol)

## [1] 1.283184

Con ese valor, ya sabemos que el odds de ser mujer supera al odds del hombre por un factor de 0.28. El odds ratio (OR) puede ir de 0 a infinito. Un OR de 1 implica que no hay diferencias.

Veámoslo como fracción:

MASS::fractions(OddsMujVol/OddsHomVol)

## [1] 137157/106888

Si encuentras la cantidad de voluntarias del numerador, se espera la cantidad de voluntarios del denominador.

Veamos la tabla de contingencia de manera gráfica:

mosaicplot( t(volsexTable),col = c("orange", "green"),main = "")

La diferencia es clara, pero ¿será estadísticamente significativa?

2.2.3.1 Consideraciones Previas

La regresión es sencilla de calcular, pero la presencia de variables categóricas puede traer complicaciones. Veamos como R está utilizando la variable dependiente:

# la representación numérica
as.numeric(as.factor(c("no","yes")))

## [1] 1 2

Qué pasaría si tuvieramos otras columnas:

vol$IsVolunteer=ifelse(vol$volunteer=="yes",T,F)
vol$IsVolunteer_1=as.integer(ifelse(vol$volunteer=="yes",1,0))
vol$IsVolunteer_claro=as.factor(ifelse(vol$volunteer=="yes",
                             "claro","nunca"))

# nos queda:
library(magrittr)
vol[,c("volunteer","IsVolunteer","IsVolunteer_1","IsVolunteer_claro")]%>%
    rmarkdown::paged_table(options = )

Reescribamos la hipótesis inicial con las nuevas variables:

h1b=formula(IsVolunteer~sex)
h1c=formula(IsVolunteer_1~sex)
h1d=formula(IsVolunteer_claro~sex)

Veamos:

• Dicotómica como Booleana (lógica)

coef(glm(h1b, data=vol,family = binomial))

## (Intercept)     sexmale 
##  -0.2110362  -0.2493447

• Dicotómica como Entero:

coef(glm(h1c, data=vol,family = binomial))

## (Intercept)     sexmale 
##  -0.2110362  -0.2493447

• Dicotómica como otras palabras

coef(glm(h1d, data=vol,family = binomial))

## (Intercept)     sexmale 
##   0.2110362   0.2493447

Como se ve, R usa bien los casos Booleanos y enteros, pero ante la presencia de otros textos en la respuesta sólo usará el orden alfabético.

Otro tema es el uso de categóricas en las variables independientes. El resultado anterior ha usado a mujer como referencia: nos informa cuánto afecta ser hombre al logaritmo del odd de ser voluntario, en comparación con ser mujer (la referencia). La primera tabla nos confirma que ser hombre disminuye (aquí el signo ayuda) el logaritmo del odds (y por ende la probabilidad) de ser voluntario. Veamos alternativamente la tabla, con la misma regresión pero con la categoría hombre como referencia.

# cambio de referencia
vol$sex=relevel(vol$sex,ref = "male")

#regresion
rlog1=glm(h1, data=vol,family = binomial)
modelrl=list('Ser Voluntario (I)'=rlog1)
modelsummary(modelrl,
             title = "Regresión Logística (con cambio de referencia en sexo)",
             stars = TRUE,
             output = "kableExtra")

La tabla nos muestra el efecto (significativo al 0.05) positivo (nuevamente el signo ayuda) de ser mujer en ser voluntario. Recuerda que este tema de la categoría de referencia se aplica a todo tipo de regresión.

2.2.3.2 Interpretación

Los coeficientes de la anterior tabla, como modelan al logaritmo del odds, no son fáciles de interpretar. Pero, si le aplicamos exponencial obtendrás un odds ratio: veamos el caso del coeficiente 0.249 de la anterior tabla.

sexF=coef(rlog1)["sexfemale"]
exp(sexF)

## sexfemale 
##  1.283184

Este valor ya lo conocíamos, pero ahora sabemos que el efecto de sexo es significativo gracias a la regresión logística, alque no sabíamos con las tabla de contingencia. Lo que es más, sabemos que el hecho de ser mujer eleva en promedio en un factor de 0.28 al odds ratio de ser voluntario que si el entrevistado fuera hombre.

Algo importante ahora es que con la regresión logística ya podemos tener predictores numéricos, lo que escapa de la utilidad de las tablas de contingencia. Así, podemos en el modelo II añadir una numérica los predictores: El ofrecerse como voluntario está relacionado con el sexo de la persona, y su nivel de neuroticismo.

La siguiente tabla nos muestra los resultados para el modelo I y II, ya con los coeficientes exponenciados:

### second hypothesis
h2=formula(volunteer~sex + neuroticism)
rlog2=glm(h2, data=vol,family = binomial)

modelsrl=list('Ser Voluntario (I)'=rlog1,
             'Ser Voluntario (II)'=rlog2)

# formato creado para modelsummary
formatoNumero = function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
             fmt=formatoNumero, # usa función que creé antes
             exponentiate = T, # coeficientes sin logaritmo
             statistic = 'conf.int', # mostrar ICs
             title = "Regresión Logísticas (Coeficientes Exponenciados)",
             stars = TRUE,
             output = "kableExtra")

En el modelo II notamos que el efecto de ser mujer disminuye, pero no mucho. Probemos ahora el modelo III, el que completa la hipótesis de los investigadores: El ofrecerse como voluntario está relacionado con el sexo de la persona, su nivel de neuroticismo y su nivel de extraversión.

h3=formula(volunteer~sex + neuroticism + extraversion)
rlog3=glm(h3, data=vol,family = binomial)

summary(rlog3)

## 
## Call:
## glm(formula = h3, family = binomial, data = vol)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.351658   0.239818  -5.636 1.74e-08 ***
## sexfemale     0.235161   0.111185   2.115   0.0344 *  
## neuroticism   0.006362   0.011357   0.560   0.5754    
## extraversion  0.066325   0.014260   4.651 3.30e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1933.5  on 1420  degrees of freedom
## Residual deviance: 1906.1  on 1417  degrees of freedom
## AIC: 1914.1
## 
## Number of Fisher Scoring iterations: 4

modelsrl=list('Ser Voluntario (I)'=rlog1,
             'Ser Voluntario (II)'=rlog2,
             'Ser Voluntario (III)'=rlog3)

f <- function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
             fmt=formatoNumero,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "Comparando Regresión Logísticas (Coeficientes Exponenciados)",
             stars = TRUE,
             gof_map = c("nobs","aic","bic","rmse","logLik"), #comparar
             gof_omit = c("F"),
             output = "kableExtra")

round(exp(coef(rlog3)["sexfemale"]),4)

## sexfemale 
##    1.2651

100*round(abs(1-exp(coef(rlog3)["sexfemale"])),4)

## sexfemale 
##     26.51

round(exp(coef(rlog3)["extraversion"]),4)

## extraversion 
##       1.0686

100*round(abs(1-exp(coef(rlog3)["extraversion"])),4)

## extraversion 
##         6.86

vol$sexFem_num=ifelse(vol$sex=='male',0,1) # 1 es para mujer
sdVIs=apply(vol[,c("sexFem_num","neuroticism", "extraversion")],2,sd)
DF=list(LogitSt=sdVIs*coef(rlog3)[c(2,3,4)])%>%
       data.frame()

# DF tiene los coeficientes estandarizados

DF%>% kable(caption = "Coeficientes Estandarizados (ordenar vía valores absolutos)")%>%
          kableExtra::kable_styling(full_width = F)

2.2.3.3 Decidiendo mejor modelo

R entrega algunos índice muy buenos para comparar modelos logísticos, el BIC, AIC, RMSE, y Log.Lik. Los primeros tres indican que el modelo es mejor si tienen valores bajos; mientras que el último será mejor si su valor aumenta.

Sin embargo, las diferencias pueden ser mínimas entre estos índices, por lo que es apropiado usar el test de razón de verosimilitud (likelihood ratio test-LRT), cuyo resultado es el siguiente:

library(lmtest)

lrtest(rlog1,rlog2, rlog3) %>%
kable(caption = "Tabla LRT para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)

Como vemos en la anterior tabla, pasar del modelo 1 al modelo 2 no resulta significativo (probabilidad de similitud entre los modelos es mayor que 90%), por lo que no se rechaza la hipótesis nula de igualdad de modelos. Pero sí resultado significativo pasar del modelo 2 al modelo 3.

Podemos comparar los modelo I, II y III de manera gráfica:

library(ggplot2)
dotwhisker::dwplot(list(Logit_I=rlog1,Logit_II=rlog2,Logit_III=rlog3),
                   exp=T) + #exponenciados!
            scale_y_discrete(labels=c("extraversion",
                                      "neuroticism",
                                      "sex (ref: male)")) +
            scale_color_discrete(name="Modelos para:\nSer Voluntario") +
            geom_vline(xintercept = 1,
                       colour = "grey60",
                       linetype = 2)

Acá se puede apreciar claramente por qué el neuroticismo no es significativo.

2.2.3.4 Evaluando modelo seleccionado

Cuando usamos una regresión logística queremos predecir probabilidades, en este caso, ante una combinación de los predictores X, qué probabilidad se obtiene con esta regresión para ser voluntario. Esos son los valores ajustados (fitted) que entrega el modelo, y nos servirá para evaluar si nos hemos acercado al Y original. Calculemos tales probabilidades:

vol$predictedProbs <- predict(rlog3, vol,type = "response")
# veamos algunos de ellos
head(vol$predictedProbs)

## [1] 0.4619546 0.4080078 0.4356976 0.5648592 0.4914463 0.4210156

Aquí podemos usar estos valores para crear el 0 o 1:

vol$predicted=ifelse(vol$predictedProbs>0.5,"yes","no")

Veamos algunas filas de estas dos columnas:

head(vol[,c("predicted","volunteer")],10)

##    predicted volunteer
## 1         no        no
## 2         no        no
## 3         no        no
## 4        yes        no
## 5         no        no
## 6         no        no
## 7         no        no
## 8         no        no
## 9         no        no
## 10        no        no

Las columnas predicted busca ser la misma que la columna volunteer, pero en algunos casos fallará.

Para saber qué tanto se acertó creamos una Matriz de Confusión usando el paquete cvms en la siguiente figura:

library(cvms)

## Warning: package 'cvms' was built under R version 4.4.3

ConfMatrix=confusion_matrix(predictions =  vol$predicted,
                      targets= vol$volunteer)
library(cvms)
plot_confusion_matrix(ConfMatrix,
                      class_order=c("yes","no"))

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

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

La diagonal secundaria o antidiagonal debería tener solo 0s si tuviésemos una predicción perfecta, y los demás valores estar en la diagonal principal. El modelo, debería servir para reducir sustantivamente ese error.

Note que en la matriz de confusión se puede calcular valores clave: Verdaderos Positivos (VP), Verdaderos Negativos (VN), Falsos Positivos (FP), Falsos Negativos (FN).

Y a partir de ellos: - La Sensitividad (VP/(VP+FN)): Utilidad del modelo para determinar el “voluntario” (modelo que predice a un voluntario que realmente era voluntario en la data orignial)

ConfMatrix$Sensitivity

## [1] 0.1474037

• La Especificidad (VN/(VN+FP)): Utilidad delm odelo para determinar el “no voluntario” (modelo predice a un no voluntario que realmente no era voluntario en la data original)

ConfMatrix$Specificity

## [1] 0.8932039

Una medida general de resumen es el índice de Kappa: Este índice vale 1 si la predicción es perfecta, y 0 si faltan todas las predicciones:

ConfMatrix$Kappa

## [1] 0.04497635

2.2.3.5 Prediciendo valores concretos

Si nos quedamos con el modelo III, ya que es el mejor, podemos utilizarlo para predicciones. Veamos:

• Ejemplo 1: Predecir probabilidad de voluntariado de una mujer con nivel de neuroticismo=13 y extraversion=8.

ToInput1 <- data.frame(sex=factor(c("female")), 
                      neuroticism=13, extraversion=8)
predict(object = rlog3, newdata = ToInput1, type = "response")

##         1 
## 0.3767916

• Ejemplo 2: Predecir probabilidad de voluntariado de una mujer y hombre con los mismos niveles anteriores de neuroticismo y extraversión.

ToInput2 <- data.frame(sex=factor(c("female","male")), 
                      neuroticism=13, extraversion=8)
predict(object = rlog3, newdata = ToInput2, type = "response")

##         1         2 
## 0.3767916 0.3233650

• Ejemplo 3: Predecir probabilidad de voluntariado de dos mujeres, con nivel de neuroticismo de 13 y 10, y extraversión de 8 y 21, respectivamente:

ToInput3 <- data.frame(sex=factor(c("female","female")), 
                      neuroticism=c(13,10), extraversion=c(8,21))
predict(object = rlog3, newdata = ToInput3, type = "response")

##         1         2 
## 0.3767916 0.5841798

2.2.3.6 Efectos marginales

Recuérdese qué tenemos en volunteer, predicted y predictedProbs:

head(vol[,c('volunteer','predicted','predictedProbs')])

##   volunteer predicted predictedProbs
## 1        no        no      0.4619546
## 2        no        no      0.4080078
## 3        no        no      0.4356976
## 4        no       yes      0.5648592
## 5        no        no      0.4914463
## 6        no        no      0.4210156

Cuando analizábamos la regresión Gaussiana decíamos que si el covariado aumentaba en una unidad, el coeficiente del covariado representaba el cambio promedio en la dependiente. Ese es el efecto marginal del covariado.

Esa interpretación no se puede saber con ninguno de los valores calculados hasta ahora para la regresión logística. Para ello, necesitamos calcular de manera expresa esos efectos marginales

# interpracion usando marginal effects:
library(margins)

## Warning: package 'margins' was built under R version 4.4.2

# 
margins(rlog3)

## Average marginal effects

## glm(formula = h3, family = binomial, data = vol)

##  neuroticism extraversion sexfemale
##      0.00152      0.01585   0.05623

Estos coeficientes sí indican cuanto aumentaría en promedio la probabilidad de ser voluntario, cuando el predictor aumenta en una unidad. Sin embargo, recordemos que sabemos de antemano que no todos estos covariados son significativos. Veamos los intervalos de confianza:

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

Podemos verlo gráficamente también:

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

Ya sabemos que si el intervalo incluye al cero, ese covariado no será significativo.

2.3.2.1 Análisis Kaplan-Meier (KM)

KM es el procedimiento descriptivo básico que se utiliza para ver la dinámica de sobrevivencia. Así, vemos que la siguiente figura nos muestra el comportamiento genérico de permanecer libre.

library(ggplot2)
library(ggfortify)

## Warning: package 'ggfortify' was built under R version 4.4.3

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

Esta figura nos da una idea de cómo se comporta esta población; por ejemplo, la gráfica nos dice que si pasan 40 semanas, la probabilidad de seguir libre está cerca al 80%.

El análisis KM es más interesante para ver una comparación:

KM_H1=formula(survival ~ tuvoApoyoDinero)

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 apoyo financiero?"

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

La posición de las curvas nos hace pensar que le cuesta más a los que no tuvieron fondos permanecer libres. Para una mayor confianza en tal hipótesis, podemos hacer la prueba de Mantel-Cox (LogRank), obteniendo este nivel de significancia:

LogRank=survdiff(KM_H1, data = carcel)
# ver p-valor
LogRank$pvalue

## [1] 0.05011612

Según la H0 de KM: no hay diferencias entre grupos, con el p-valor obtenido la diferencia no es significativa al 0.05. (pero sí al 0.1) Esta figura lo aclara más:

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

De nuevo, como sólo hay dos variables, es difícil saber qué más interviene. De ahí que necesitamos un modelo que permita análisis multivariado, o sea la Regresión Cox.

2.3.2.2 Regresión Cox

Como toda regresión, esta técnica permite utilizar regresores o predictores o covariados, pero no modela la duración sino el riesgo de que el evento suceda (ser re arrestado)

COX_H1= formula(survival~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)

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

La tabla muestra los coeficientes originales y exponenciados (las razones de riesgo o hazard ratios). Veamos cada uno:

• Se puede inferir que dar financiamiento disminuye el riesgo de ser re arrestado con una significancia del 0.05 La columna a la derecha (los HRs) nos da el promedio de ese efecto en el riesgo como la distancia absoluta al uno (1); en este caso:

(apoyoDinero=abs(1-exp(coef(rcox1)[1])))

## tuvoApoyoDinero1 
##        0.3399735

De ahí que usando el hazard ratio, podemos sostener que el riesgo de volver a la cárcel para alguien con apoyo financiero es en promedio 34% menor a la de los que no lo tienen.

• Se puede inferir que es muy poco probable que el nivel educativo tenga efecto en el riesgo de ser re arrestado.

• Se puede inferir que haber estado previamente en la cárcel aumenta el riesgo de ser re arrestado con una significancia del 0.01. El promedio de ese riesgo es la distancia absoluta al uno (1); en este caso:

(carcelantes=abs(1-exp(coef(rcox1)[6])))

## vecesEnCarcel 
##     0.0907188

De ahí que podemos sostener que el riesgo de volver a la cárcel para alguien con apoyo financiero es en promedio 9.07% mayor cada vez que aumenta en uno los años de encarcelamientos previos.

Los hazard ratios de la tabla anterior (columna a la derecha, exponenciados) podemos verlos gráficamente en la siguiente figura:

library(survminer)

## Warning: package 'survminer' was built under R version 4.4.3

## Cargando paquete requerido: ggpubr

## Warning: package 'ggpubr' was built under R version 4.4.2

## 
## Adjuntando el paquete: 'ggpubr'

## The following object is masked from 'package:cvms':
## 
##     font

## 
## Adjuntando el paquete: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

ggforest(rcox1, data = carcel,main = "¿Quiénes tienen mayor riesgo de volver a ser encarcelados?")

2.3.2.3 Comparación

Podemos plantear un modelo sin educación

COX_H1= formula(survival~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)
COX_H2= formula(survival~tuvoApoyoDinero+vecesEnCarcel)

#regression
rcox2 <- coxph(COX_H2,data=carcel)
modelcox=list('Riesgo de Re arrestado (I)'=rcox2,'Riesgo de Re arrestado (II)'=rcox1)

#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
modelsummary(modelcox,
             #fmt=f,
             exponentiate = T, 
             statistic = 'conf.int',
             title = "Regresión Cox (sólo Hazard Ratios)",
             stars = TRUE,
             output = "kableExtra")

Vemos que las medidas globales de ajuste no son tan diferentes, por lo que es mejor evaluar la significancia de esas diferencias.

anova(rcox2,rcox1)%>%
knitr::kable(caption = "Tabla anova para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)

Así, añadir nivel educativo no es significativo al 0.05, pero sí lo es al 0.1. Esto daría espacio para sostener la inclusión de nivel educativo en el modelo.

3.2.1.1 Verificar si los datos permiten factorizar

Si el Overall MSA resulta cercano al 1 indica que es muy adecuado hacer un EFA. Valores bajos indican que las correlaciones pueden ser débiles o espurias.

• 0.90 Excelente

• 0.80-0.89 Muy bueno
• 0.70-0.79 Bueno
• 0.60-0.69 Aceptable
• 0.50-0.59 Mediocre
• < 0.50 No aceptable

library(psych)

## Warning: package 'psych' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'psych'

## The following object is masked from 'package:polycor':
## 
##     polyserial

## The following objects are masked from 'package:arm':
## 
##     logit, rescale, sim

## The following object is masked from 'package:modelsummary':
## 
##     SD

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

## The following objects are masked from 'package:DescTools':
## 
##     AUC, ICC, SD

psych::KMO(corMatrix) 

## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = corMatrix)
## Overall MSA =  0.9
## MSA for each item = 
##          hdiLife   hdiSchoolExpec      hdiMeanEduc           hdiGni 
##             0.90             0.92             0.91             0.88 
##     ideElectoral   ideFunctioning ideParticipation       ideCulture 
##             0.84             0.93             0.96             0.87 
##     ideLiberties 
##             0.88

3.2.1.2 Verificar si la matriz de correlaciones es adecuada

Aquí hay dos pruebas:

• Hnula: La matriz de correlación es una matriz identidad

cortest.bartlett(corMatrix,n=nrow(theData))$p.value>0.05

## [1] FALSE

Este test evalúa la hipótesis nula de que la matriz de correlaciones es una matriz identidad, es decir, que no hay correlación entre variables (lo cual sería malo para el EFA). Nosotros queremos rechazar esta hipótesis (y para ello necesitamos que el p-valor < 0.05). El FALSE sinifica que el p-valor es menor a 0.05, por lo tanto rechazamos la hipótesis nula.

• Hnula: La matriz de correlación es una matriz singular

library(matrixcalc)

is.singular.matrix(corMatrix)

## [1] FALSE

Esta prueba evalúa si la matriz es singular, es decir, no invertible, lo cual causaría problemas matemáticos al hacer el EFA. Para que el análisis factorial funciones correctamente, la matriz debe NO ser singular (es decir, debe ser invertible). Si sale FALSE nos indica que la matriz NO es singular

3.2.1.3 Número de factores o variables latentes

Necesitamos determinar en cuántos factores o variables latentes podríamos redimensionar la data: En este caso, la función fa.parallel nos dará la sugerencia:

fa.parallel(theData, fa = 'fa',correct = T,plot = F)

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect.  Try a
## different factor score estimation method.
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect.  Try a
## different factor score estimation method.

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected.  Examine the results carefully

## Parallel analysis suggests that the number of factors =  2  and the number of components =  NA

La prueba nos sugiere 2, lo esperado en teoría.

3.2.1.4 Redimensionar a número menor de factores

La rotación en el análisis factorial se utiliza para mejorar la interpretabilidad de los factores extraídos. No cambia los resultados fundamentales (como la varianza total explicada), pero sí cambia cómo se distribuyen las cargas (loadings) entre los factores.

Varimax: Rotación de tipo ortogonal (factores no correlacionados), supone que los factores son independientes entre sí, se emplea cuando quieres mantener independencia entre dimensiones, el resultado es más claro (cada variable tiende a cargar alto en un solo factor)

Oblimin: Rotación de tipo oblicua (factores pueden correlacionarse), supone que los factores pueden estar relacionados, se emplea cuando sospechas que los factores tienen relación conceptual, el resultado puede mostrar una variable cargando en más de un factor si están relacionados.

• Resultado inicial:

library(GPArotation)

## 
## Adjuntando el paquete: 'GPArotation'

## The following objects are masked from 'package:psych':
## 
##     equamax, varimin

resfa <- fa(theData,
            nfactors = 2,
            cor = 'mixed',
            rotate = "varimax", #oblimin?
            fm="minres")
print(resfa$loadings)

## 
## Loadings:
##                  MR1   MR2  
## hdiLife          0.307 0.836
## hdiSchoolExpec   0.368 0.811
## hdiMeanEduc      0.316 0.782
## hdiGni           0.283 0.813
## ideElectoral     0.906 0.249
## ideFunctioning   0.801 0.470
## ideParticipation 0.768 0.331
## ideCulture       0.498 0.382
## ideLiberties     0.901 0.334
## 
##                  MR1   MR2
## SS loadings    3.522 3.280
## Proportion Var 0.391 0.364
## Cumulative Var 0.391 0.756

• Resultados mejorado (solo apropiado si hay más de un factor):

print(resfa$loadings,cutoff = 0.5)

## 
## Loadings:
##                  MR1   MR2  
## hdiLife                0.836
## hdiSchoolExpec         0.811
## hdiMeanEduc            0.782
## hdiGni                 0.813
## ideElectoral     0.906      
## ideFunctioning   0.801      
## ideParticipation 0.768      
## ideCulture                  
## ideLiberties     0.901      
## 
##                  MR1   MR2
## SS loadings    3.522 3.280
## Proportion Var 0.391 0.364
## Cumulative Var 0.391 0.756

Cuando logramos que cada variable se vaya a un factor, tenemos una estructura simple.

• Resultado visual: El resultado lo podemos ver de forma gráfica

fa.diagram(resfa,main = "Resultados del EFA")

#### Evaluando Resultado obtenido:

• ¿Qué variables aportaron más a los factores?

  • Proporción de la varianza de cada variable que es explicada por los factores extraídos.
  • Va de 0 a 1: mientras más alta, mejor representa esa variable el modelo factorial.
  • Variables con comunalidad cercana a 1: se explican bien con los factores.
  • Variables con comunalidad <0.4: quizás no encajan bien en la estructura.

sort(resfa$communality)

##       ideCulture ideParticipation      hdiMeanEduc           hdiGni 
##        0.3941682        0.6993302        0.7119301        0.7405067 
##          hdiLife   hdiSchoolExpec   ideFunctioning     ideElectoral 
##        0.7933223        0.7936103        0.8629055        0.8835971 
##     ideLiberties 
##        0.9226321

• ¿Qué variables contribuyen a la construcción de más de un factor?

  • Mide cuántos factores distintos “usan” una variable para construir su carga factorial.
  • Medida numérica del grado en que una variable está asociada con más de un factor.
  • Complejidad = 1 –> la variable carga en solo un factor, lo deseable para una interpretación clara
  • Complejidad > 2 –> la variable carga en más de un factor, lo cual puede dificultar la interpretación

sort(resfa$complexity)

##     ideElectoral           hdiGni          hdiLife     ideLiberties 
##         1.150198         1.238563         1.265339         1.270410 
##      hdiMeanEduc ideParticipation   hdiSchoolExpec   ideFunctioning 
##         1.317113         1.358907         1.395523         1.614353 
##       ideCulture 
##         1.874466

• ¿Tucker Lewis > 0.9?

  • Es un índice de ajuste comparativo que evalúa si el modelo factorial mejora sustancialmente con respecto a un modelo nulo (sin relaciones entre variables)
  • Se interpreta similar al R cuadrado en regresión.

  • 0.95 Excelente ajuste

  • 0.90-0.95 Buena ajuste
  • 0.80-0.89 Aceptable, pero mejorable
  • < 0.80 Mal ajuste

resfa$TLI

## [1] 0.9263416

• ¿RMS cerca a cero?

  • Medida cuadrática de las diferencias entre las correlaciones observadas y las reproducidas por el modelo.
  • Evalúa qué tan bien el modelo reproduce las correlaciones reales.
  • =< 0.05 Excelente ajuste
  • 0.05-0.08 Aceptable

  • 0.08 Mal ajuste

resfa$rms

## [1] 0.0346878

• ¿RMSEA cerca a cero?

  • Mide la discrepancia por grado de libertad.
  • Evalúa el ajuste del modelo a la estructura poblacional (no solo muestral)
  • =< 0.05 Muy buen ajuste
  • 0.05-0.08 Buen ajuste
  • 0.08-0.10 Aceptable/regular

  • 0.10 Mal ajuste

resfa$RMSEA

##      RMSEA      lower      upper confidence 
##  0.1319641  0.1014027  0.1651144  0.9000000

• ¿BIC?
  • Bayseian Information Criterion
  • No se interpreta en valores absolutos
  • Sirve para comparar modelos: cuanto más bajo, mejor
  • No hay valor bueno o malo
  • Se usa para comparar el BIC entre modelos con diferentes número de factores (si empleamos uno con 3,4 o 5)

resfa$BIC

## [1] -23.30257

4.5.1.1 Decidir cantidad de clusters

La siguiente figura nos sirve para determinar la cantidad de cllusters a solicitar (usando el estadístico gap):

## para PAM

library(factoextra)

## Warning: package 'factoextra' was built under R version 4.4.3

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

fviz_nbclust(dataClus, pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)

4.5.1.2 Clusterizar via PAM:

La técnica de k-medoids se implementa en la función PAM. Esta función retorna diversos valores, en este caso crearemos una columna con la etiqueta del cluster. Usemos la sugerencia de la anterior figura y hallamos:

library(kableExtra)
set.seed(123)
res.pam=pam(g.dist,3,cluster.only = F)

#nueva columna
dataClus$pam=res.pam$cluster

# ver

head(dataClus,15)%>%kbl()%>%kable_styling()

4.5.1.3 Evaluando el uso de PAM

Una manera práctica de ver el desempeño del algoritmo es calcular las silhouettes. Para el caso reciente, vemamos la siguiente figura:

fviz_silhouette(res.pam,print.summary = F)

La figura muestra barras, donde cada una es un país (caso). Mientras más alta la barra, la pertenencia a ese cluster es clara. La barra negativa indica un país mal clusterizado. Para este caso, estos serían los mal clusterizados:

silPAM=data.frame(res.pam$silinfo$widths)
silPAM$country=row.names(silPAM)
poorPAM=silPAM[silPAM$sil_width<0,'country']%>%sort()
poorPAM

## [1] "Chile"  "Latvia"

4.5.1.4 Verificando etiqueta de clusters

Exploremos el promedio de cada cluster:

aggregate(.~ pam, data=dataClus,mean)

##   pam    hdiLife hdiSchoolExpec hdiMeanEduc     hdiGni
## 1   1 -1.0811743     -1.0681695  -1.1822831 -0.8111138
## 2   2  0.1340138      0.1819924   0.3393855 -0.2301917
## 3   3  1.2520905      1.1502463   1.0317146  1.5181170

El número asignado al cluster no tiene significado necesariamente, por lo que recomiendo hacer cálculos para dárselo. En este caso, las etiquetas ascienden al igual que el promedio, por lo que no es necesario recodificar la etiqueta.

Antes de continuar, guardemos la columna de PAM en la data integrada, y eliminemos la de dataClaus

idhdemo$pamIDHpoor=idhdemo$country%in%poorPAM
idhdemo$pamIDH=as.ordered(dataClus$pam)
dataClus$pam=NULL

4.5.2.1 Estrategia Aglomerativa

En esta estrategia se parte por considerar cada caso (fila) como un cluster, para de ahí ir creando miniclusters hasta que todos los casos sean un solo cluster. El proceso va mostrando qué tanto esfuerzo toma juntar los elementos cluster tras cluster.

## PARA JERARQUICO

fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "agnes")

set.seed(123)
library(factoextra)

res.agnes<- hcut(g.dist, k = 3,hc_func='agnes',hc_method = "ward.D")

dataClus$agnes=res.agnes$cluster

# ver

head(dataClus,15)%>%kbl()%>%kable_styling()