Profesor:Dr. José Manuel MAGALLANES REYES, Ph.D

Guía Básica

El Pipeline del Analista Político

Podemos proponer que los analistas existen para brindar explicaciones y recomendaciones sobre algún asunto de interés. Por otro lado, el decisor es responsable de eligir un curso de acción sabiendo que el analista le ha dado información incompleta. Ya el gestor se encarga de implementar las decisiones; esa implementación traerá nueva información y el analista vuelve a su trabajo.

La estadística es una ciencia matemática que guía científicamente el análisis de datos. Así el analista trata de seguir una secuencia en su labor:

  1. Apuesta por entender una variable que explicaría un problema de interés (narcotráfico, elecciones, etc). En esta etapa, hace exploración de los datos de esa variable, es decir, organiza los datos en tablas, y de ahi calcula medidas, tablas resumen y gráficos. Se entiende que la variable de interés tiene una variabilidad tal que despierta en el analista la necesidad de preguntar por qué esa variabilidad.

  2. Se plantea hipótesis respecto a qué se relaciona con la variabilidad de la variable de interés. Las hipótesis se formula no antes de haber revisado la literatura. En esta etapa hace uso de análisis bivariado o multivariado. Esta etapa se enriquece si está actualizado en las teorías que proponen cierta asociación de la variable de interés con otras variables.

  3. Aplica la técnica estadística que corresponda, tal que confirme, o no, la hipótesis que se planteó.

  4. Interpreta los resultados obtenidos.

  5. Elabora síntesis de lo actuado; propone explicaciones a lo encontrado; y elabora recomendaciones.

Hay muchas opciones para las técnicas señaladas en el punto 3. La elección de las mismas dependerá de la preparación del analista. Esta sesión te guiará para que:

  1. Recuerdes cómo y para qué hacer exploración univariada (ir)
  2. Recuerdes cómo y para qué hacer análisis bivariada (ir).
  3. Introducir el concepto (y necesidad) de la regresión multivariada (ir).

1. La tabla de datos

Supongamos que estás interesado en la “situación de los locales escolares en el Perú”. Ese simple interés te lleva a buscar datos para saber tal situación. De pronto te encuentras con estos datos (basado en CEPLAN):

Imaginemos que estos datos son lo ‘mejor’ que podrás conseguir:

  • Variable: Locales Publicos en buen estado; Tipo Medida: Porcentaje.

  • Variable: Cantidad de Contribuyentes a la SUNAT; Tipo Medida: Conteo.

  • Variable: PEA Ocupada; Tipo Medida: Conteo.

  • Variables de Población; Tipo Medida: Conteo.

2. Carga de datos:

Si los datos están en GoogleDrive, puedes leerlos 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'

dataRegiones=read.csv(linkADrive )

head(dataRegiones)

Antes de explorar los datos de manera estadística, debemos primero verificar cómo R ha intepretado el tipo de datos:

str(dataRegiones)
## '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 ...

Asumamos que la variable de interés será la situación de los locales de colegios. La columna de interés es de tipo numérico, y R también lo tiene intrepretado así.

3. Exploración Univariada de la Variable de Interés

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ígrafos es usando el comando summary():

summary(dataRegiones$buenEstado)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.70   13.90   17.30   18.61   21.80   40.50

Pero, nos faltaría algunos estadígrafos, que podemos añadir así:

library(DescTools)

allStats=c(summary(dataRegiones$buenEstado),
  sd=sd(dataRegiones$buenEstado), # variabilidad (en relacion a la media)
  skew=Skew(dataRegiones$buenEstado), # asimetria (positiva o negativa)
  kurt=Kurt(dataRegiones$buenEstado), # concentración (enpinada / aplanada)
  cv=CoefVar(dataRegiones$buenEstado)) # variabilidad (mayor o menor que uno)
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

De lo anterior queda claro que hay dispersion moderada (sd<mean / cv<1 ); ligera cola a la derecha (skew >0 / mean > median); y que está poco empinada (kurt > 0).

Veamos lo mismo en el histograma:

library(ggplot2)

base=ggplot(data=dataRegiones,
            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=dataRegiones,
            aes(y=buenEstado))
boxplot=base + geom_boxplot() 

boxplot

A esta altura sólo falta identificar a los atípicos en el boxplot, lo cual podemos recuperar con ggplot_build:

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

Los outliers son:

outliersLocales=valuesFromBox$outliers[[1]]
outliersLocales
## [1] 33.8 40.5

Nota que cuando hay outliers, se puede calcular del grafico los maximos y minimos teóricos:

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

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

dataRegiones[dataRegiones$buenEstado>31.4,]

Es decir, los locales en buen estados se distribuyen por lo general entre 7.7 y 31.4 por ciento, y que la mitad de los departamentos tienen a lo más un 17.3 por ciento de buenos locales. Las regiones destacadas son sólo dos, pero ninguna de ellas supera siquiera el 50 por ciento.

4. Análisis Bivariado

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

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

Calculemos ahora los indices de correlación:

f1=formula(~peaOcupada + buenEstado)

Camino parametrico:

pearsonf1=cor.test(f1,data=dataRegiones)[c('estimate','p.value')]
pearsonf1
## $estimate
##       cor 
## 0.3567442 
## 
## $p.value
## [1] 0.08002776

Camino no parametrico

spearmanf1=cor.test(f1,data=dataRegiones,method='spearman',exact=F)[c('estimate','p.value')]
spearmanf1
## $estimate
##       rho 
## 0.2785151 
## 
## $p.value
## [1] 0.1776146

Hasta aquí, dudamos que esas variables estén correlacionadas.

Otro caso importante es cuando analizamos nuestra variable versus una variable categórica. Por ejemplo, creemos una columna que indique (de manera imprecisa) costa,sierra y selva:

dataRegiones$zona=c('selva','costa','sierra', 'costa','sierra',
                      'sierra','costa','sierra','sierra','selva',
                      'costa','sierra','costa','costa','costa',
                      'selva','selva','costa','sierra','costa',
                      'sierra','selva','costa','costa','selva')

# tendriamos
dataRegiones

Pidamos un boxplot, pero por grupos:

base=ggplot(data=dataRegiones, aes(x=zona, y=buenEstado))
base + geom_boxplot() +  geom_jitter(color="black", size=0.4, alpha=0.9)

Parece que hubiese diferencia entre las categorías. Para 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 tests de normalidad. Usemos Shapiro-Wilk en cada grupo:

f2=formula(buenEstado~zona)


tablag= aggregate(f2, dataRegiones,
          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

Habría normalidad en todos los grupo. Usemos entonces tanto la prueba de Kruskall (no paramétrica) como la prueba anova para analizar las diferencias:

aovZona=aov(f2,data=dataRegiones)
summary(aovZona)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## zona         2  872.2   436.1   12.54 0.000232 ***
## Residuals   22  765.1    34.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aovZona)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = f2, data = dataRegiones)
## 
## $zona
##                    diff       lwr       upr     p adj
## selva-costa   -8.698485 -16.21700 -1.179969 0.0214473
## sierra-costa -13.381818 -20.26541 -6.498226 0.0001991
## sierra-selva  -4.683333 -12.68394  3.317276 0.3239459
kruskal.test(f2,data=dataRegiones)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  buenEstado by zona
## Kruskal-Wallis chi-squared = 15.004, df = 2, p-value = 0.000552
pairwise.wilcox.test(dataRegiones$buenEstado, dataRegiones$zona, p.adjust.method = "bonferroni")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  dataRegiones$buenEstado and dataRegiones$zona 
## 
##        costa   selva  
## selva  0.03580 -      
## sierra 0.00056 0.17782
## 
## P value adjustment method: bonferroni

Con toda esta información, iriamos concluyendo también que los locales en buen estado de la costa son más numeros que las otras zonas del país, y que la sierra con la selva podrían estar al mismo nivel.

5. Regresión Lineal

La regresión es una técnica donde hay que definir una variable dependiente y una o más independientes. Las independientes pueden tener rol predictor, dependiendo del diseño de investigación, aunque por defecto tiene un rol asociativo.

La regresión sí quiere informar cuánto una variable (independiente) está asociada a la variación de otra (dependiente), de ahí que es una técnica para probar hipótesis direccionales o asimétricas (las correlaciones tiene hipótesis simétricas).

La regresión devuelve un modelo relacional entre variables, es decir una ecuación, que recoge cómo una o más variables explicarían a otra. Para nuestro caso la variable dependiente es el estado de los locales:

modelo1=formula(buenEstado~peaOcupada + zona)

Por ejemplo, para la hipótesis ‘el estado de los locales escolares públicos en una región depende de la PEA ocupada regional y de la ubicación de la región’, la regresión arrojaría este resultado:


Hasta aquí, vemos que lo que nos informaba el análisis bivariado se mantiene en la regresión: la PEA no tiene efecto, pero la zona sí (no estar en la costa lo tiene).

Sin embargo, es aquí donde reflexionamos si los datos crudos que tenemos podrían necesitar alguna transformación:

# la pea como porcentaje
dataRegiones$peaOcu_pct=dataRegiones$peaOcupada/dataRegiones$TOTALpob


modelo2=formula(buenEstado~peaOcu_pct + zona)
regre2=lm(modelo2,data = dataRegiones)
summary(regre2)
## 
## Call:
## lm(formula = modelo2, data = dataRegiones)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.6580 -2.8489 -0.4126  2.6089 10.9903 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -5.978     11.794  -0.507  0.61753   
## peaOcu_pct    76.510     28.886   2.649  0.01502 * 
## zonaselva     -4.106      3.169  -1.296  0.20906   
## zonasierra    -9.359      2.864  -3.267  0.00368 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.226 on 21 degrees of freedom
## Multiple R-squared:  0.6497, Adjusted R-squared:  0.5997 
## F-statistic: 12.98 on 3 and 21 DF,  p-value: 5.136e-05

Aquí aparace algo muy interesante, primero que peaOcu_pct tiene efecto, pues es significativo (p-valor es menor que 0.05); segundo, que ese efecto es directo, pues el coeficiente calculado es positivo (signo de columna Estimate); y tercero que la magnitud de ese efecto es 76.51, lo que indica cuánto aumenta, en promedio, la variables dependiente, cuando la variable independiente se incremente en una unidad. Ahora bien, la zona selva ya no difiere de la costa ne esta nueva regresión.

Esto es información suficiente para representar esa relación con una ecuación:

\[ BuenEstado\_pct= -5.9781324 + 76.5102039 \cdot PEA\_pct + -4.1063326 \cdot Populoso + \epsilon\]

El Y verdadero es BuenEstado_pct, pero la regresión produce un \(\hat{BuenEstado_pct}\) estimado, de ahi la presencia del \(\epsilon\). Justamente el R cuadrado ajustado (0.6497263) nos brinda un porcentaje (multiplicalo por 100), que nos da una pista de nuestra cercanía a una situación perfecta (cuando vale 1).

Regresion con indicadores

Aqui hay dos indices conocido, democracya y desarrollo humano:

dataPreparada="https://github.com/Estadistica-AnalisisPolitico/Sesion6/raw/main/idhdemo.csv"
idhdemo=read.csv(dataPreparada)
head(idhdemo)

Vemos tanto el indice como sus indicadores:

names(idhdemo)
##  [1] "country"          "hdiRanking"       "hdi"              "hdiLife"         
##  [5] "hdiSchoolExpec"   "hdiMeanEduc"      "hdiGni"           "ideRanking"      
##  [9] "ideRegime"        "ide"              "ideElectoral"     "ideFunctioning"  
## [13] "ideParticipation" "ideCulture"       "ideLiberties"

Quedemonos con los indicadores de los indices:

dontselect=c("country","hdiRanking","hdi",
             "ideRanking","ideRegime","ide")
select=setdiff(names(idhdemo),dontselect) 
theData=idhdemo[,select]

corMatrix=polycor::hetcor(theData)$correlations

Veamos si tenemos suficiente data:

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

Pidamos 2 factores:

library(GPArotation)
resfa <- psych::fa(theData,
            nfactors = 2,
            cor = 'mixed',
            rotate = "varimax", #oblimin?
            fm="minres")
print(resfa$loadings,cutoff = 0.49)
## 
## 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       0.498      
## ideLiberties     0.901      
## 
##                  MR1   MR2
## SS loadings    3.522 3.280
## Proportion Var 0.391 0.364
## Cumulative Var 0.391 0.756

Los indices recuperan los indicadores:

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

Si lo anterior es válido, estos son los indices:

idhdemo$ide_efa=resfa$scores[,1]
idhdemo$idh_efa=resfa$scores[,2]

head(idhdemo)

Verificando una regresion con los indices con/sin EFA:

idh_ide_efa=formula(ide_efa~idh_efa)
r_EFA=lm(idh_ide_efa,data=idhdemo)
summary(r_EFA)
## 
## Call:
## lm(formula = idh_ide_efa, data = idhdemo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9808 -0.8325  0.3173  0.8105  1.4671 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.290e-17  7.533e-02   0.000     1.00
## idh_efa      6.478e-02  8.006e-02   0.809     0.42
## 
## Residual standard error: 0.9676 on 163 degrees of freedom
## Multiple R-squared:  0.004001,   Adjusted R-squared:  -0.00211 
## F-statistic: 0.6547 on 1 and 163 DF,  p-value: 0.4196
idh_ide=formula(ide~hdi)
r1_noEFA=lm(idh_ide,data=idhdemo)
summary(r1_noEFA)
## 
## Call:
## lm(formula = idh_ide, data = idhdemo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6195 -1.0091  0.4006  1.2527  3.3699 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.8224     0.6613  -2.756  0.00652 ** 
## hdi           9.7393     0.8923  10.915  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.817 on 163 degrees of freedom
## Multiple R-squared:  0.4223, Adjusted R-squared:  0.4187 
## F-statistic: 119.1 on 1 and 163 DF,  p-value: < 2.2e-16

al INICIO