ESTADISTICA PARA EL ANALISIS POLITICO II
- Profesor: Dr. José Manuel Magallanes, Ph.D.
- Profesor del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
- Oficina 105 - Edificio CISEPA / ECONOMIA / CCSS
- Telefono: (51) 1 - 6262000 anexo 4302
- Correo Electrónico: jmagallanes@pucp.edu.pe
Sesión de cierre 2020 II
Esta sesión está pensada en tu test teórico y en tu proyecto final, por lo que la he redactado integrando diversos componentes del curso que hemos compartido.
La data
Si uds trabajan en grupo lo promero que deben haber hecho es preguntarse por como el comportamiento de una variable (la dependiente) puede explicarse, en parte, por el comportamiento de otra(s) (la(s) independientes). De ahí, revisas diversos trabajos de investigación que te permiten hipotetizar esa relación. Cuando lo tienes en claro, vas por la data que operacionalice tal hipotesis.
Si me pregunto: ¿Qué explica el nivel felicidad?, luego leo y encuentro que la democracia puede influir en la felicidad, me animo a hipotetizar que la democracia influye en la felicidad. Al buscar la data, encuentro que hay data a nivel de país (país será la unidad de análisis). Estas datas ya son archiconocidas por ti:
Traigamamos la data:
library(htmltab)
# creo una "vector" con "nombres" para los links:
happyInfo=c(url="https://en.wikipedia.org/wiki/World_Happiness_Report",
path='//*[@id="mw-content-text"]/div/table/tbody')
demoInfo=c(url="https://en.wikipedia.org/wiki/Democracy_Index",
path='//*[@id="mw-content-text"]/div/table[2]/tbody')
# carga
happy = htmltab(doc = happyInfo["url"],which = happyInfo["path"])
demo = htmltab(doc = demoInfo["url"],which = demoInfo["path"])El preprocesamiento
Nuestra data seguro necesitara algo de pre procesamiento, por lo que luego de ver los nombres (names()) y la estructura (str()) será mejor asegurarnos que nos quedemos con lo que corresponda:
# Que no haya espacios en blanco al inicio ni al fin de las celdas
happy[,]=lapply(happy[,], trimws,whitespace = "[\\h\\v]")
demo[,]=lapply(demo[,], trimws,whitespace = "[\\h\\v]")
# nombres sin espacios en blanco
library(stringr)
names(happy)=str_split(names(happy)," ",simplify = T)[,1]
names(demo)=str_split(names(demo)," ",simplify = T)[,1]
#nombres de variables solo letras (usar con cuidado)
names(demo)=gsub(pattern = "\\W",replacement = "",x = names(demo))
# Eliminemos columnas que no usaremos esta vez:
happy$Overall=NULL
demo$Rank=NULL
demo$Changes=NULL
# Asegurando tipo de variable es el adecuado:
## En demo:
demo[,2:7]=lapply(demo[,2:7],as.numeric)
# En happy:
happy[,-1]=lapply(happy[,-1],as.numeric)La integración
Tener dos indicadores con sus componentes es lo que cada uno de los miembros del grupo debe haber aportado al trabajo. Luego del paso anterior, cada miembro del grupo debe tener su data lista, pero ahora hay que integrarla:
# como ambas tiene una columna 'Score', hay que renombrar
names(happy)[names(happy)=="Score"]="ScoreHappy"
names(demo)[names(demo)=="Score"]="ScoreDemo"
#ahora sí:
HappyDemo=merge(happy,demo)Hemos perdido países en el merge, en tu trabajo final debes recuperar la mayor cantidad de casos.
Hasta aquí, ya podrías hacer una regresión, para probar tu hipotesis, usando los ‘scores’:
hipotesis1=formula(ScoreHappy~ScoreDemo)
regresion1=lm(hipotesis1,data=HappyDemo)
#encontrando:
summary(regresion1)##
## Call:
## lm(formula = hipotesis1, data = HappyDemo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63706 -0.50939 0.03414 0.62139 2.29929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.65161 0.19619 18.613 <2e-16 ***
## ScoreDemo 0.31670 0.03267 9.695 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8671 on 145 degrees of freedom
## Multiple R-squared: 0.3933, Adjusted R-squared: 0.3891
## F-statistic: 93.99 on 1 and 145 DF, p-value: < 2.2e-16De lo anterior, puedes ver que se cumple que hay relación significativa (p-value del F-statistic). Esa relación es directa: “a mayor democracia, mayor felicidad”. Puede inclusive proponer que si el nivel de democracia aumenta en un punto, el nivel de felicidad aumenta en 0.31670, con una significancia del 0.001.
Calculando tus latentes
En vez de usar los scores, usemos analisis factorial. Vayamos por el camino exploratorio: “¿Habrá conjuntos de variables que representen algun concepto?”. Nótese que las variables medidas para democracia están en una misma “escala” (del 0 al 10), pero que para el caso de felicidad usaremos un conjunto de variables que no se usan para su construcció pero que están relacionadas con la felicidad.
El primer paso es examinar la matriz de correlaciones:
theData=HappyDemo[,c(3:8,10:14)]
#matriz de correlaciones
library(polycor)
corMatrix=polycor::hetcor(theData)$correlationsUna vez calculada la matriz, verificamos si el analisis factorial será util dados los datos disponibles (por convención del curso asumiremos que podemos continuar si el MSA > 0.6):
library(psych)
KMO(corMatrix)$MSA## [1] 0.855409Para ver como cada variable contribuye al “adequacy”:
KMO(corMatrix)$MSAi## GDP Social Healthy
## 0.8273328 0.8749227 0.8878139
## Freedom Generosity Perceptions
## 0.7956738 0.5524972 0.7234691
## Electoral Functioning Politicalparticipation
## 0.8021894 0.9086594 0.9269262
## Politicalculture Civilliberties
## 0.8838357 0.8595249Hasta aquí parece que todo está bien, aunque una variable podría dar problemas (Generosity). Sólo falta confirmar que la matriz de correlación cumple los dos requisitos solicitados:
- Hnula: La matriz de correlacion es una matriz identidad
cortest.bartlett(corMatrix,n=nrow(theData))$p.value>0.05## [1] FALSE- Hnula: La matriz de correlacion es una matriz singular.
#instalaste "matrixcalc'?
matrixcalc::is.singular.matrix(corMatrix)## [1] FALSEPor teoría podríamos asumir que hay dos latentes, pero sigamos explorando cuantas sugiere esta técnica:
#instalar 'parameters' y 'n_factors'
sugerencia=parameters::n_factors(corMatrix)
#instalar 'see'
# tenemos:
plot(sugerencia)
No nos sugieren 2, sino 3 latentes. Veamos cuales son:
# instalaste "GPArotation"?
library(GPArotation)
resfa <- fa(theData,
nfactors = 3,
rotate = "varimax")
fa.diagram(resfa)
Vemos que una latente corresponde a la latente de democracia, mientras que las que hemos asumido representan la felicidad no se agrupa en una sola. Si puedes sustentar una modificación a tu hipotesis, es el momento de revisar tu teoría, aquí en particular el concepto de felicidad.
Planteemos el siguiente cambio, que felicidad es el MR3, y ahora revisas tu teoría y vez si se puede sustentar la existencia de la “tranquilidad”, que sería la MR2 (¿y que tal si descartabas las variables de MR2?). Si aplicamos las evaluaciones a este EFA, vemos que hay cierto sustento para mantenerlo.
Regresion usando el camino de las latentes:
Planteemos nuestra regresión usando SEM (structural equation model). SEM hace un trabajo doble, calcula las latentes usando CFA, y lleva a cabo a regresion con las variables latentes (y otras no latentes si hubieran). El primer paso es asegurarse que la data este en la misma escala:
HappyDemoScaled=scale(HappyDemo[,-c(1,2,9,15:16)])Luegp, en formato TEXTO (entre comillas) se describe el modelo (teoria) a probar con las latentes y la regresión. Por ejemplo, una latente como felicidad, se construye a partir de GDP + Social + Healthy usando el simbolo \(=\sim\). La regresion luego conectara variables dependiente a independientes con \(\sim\), veamos:
model <- '# describiendo las latentes:
democracia =~ Electoral + Functioning + Politicalparticipation + Politicalculture + Civilliberties
tranquilidad =~ Freedom + Generosity + Perceptions
felicidad =~ GDP + Social + Healthy
# regresion con las latentes:
felicidad~tranquilidad + democracia'Ya con el modelo definido, testeas todo:
#instalaste "lavaan"?
library(lavaan)
fit <- sem(model, data=HappyDemoScaled)Todo se ha guardado en fit. Primero recupero los parametros estimados:
# conseguir los parametros
allParamSEM=parameterEstimates(fit,standardized = T)Así, podemos ver los resultados de la regresion:
allParamSEM[allParamSEM$op=="~",]## lhs op rhs est se z pvalue ci.lower ci.upper std.lv
## 12 felicidad ~ tranquilidad 0.350 0.118 2.967 0.003 0.119 0.581 0.307
## 13 felicidad ~ democracia 0.531 0.081 6.551 0.000 0.372 0.690 0.532
## std.all std.nox
## 12 0.307 0.307
## 13 0.532 0.532Eses resultado se interpreta igual que cuando hicimos regresion con lm, aunque veremos luego que no saldrá lo mismo.
En este punto, debemos evaluar todo el SEM, usemos la función fitMeasures:
# conseguir las medidas de evaluación
allFitSEM=as.list(fitMeasures(fit))- Test 1: El ChiSquare es NO significativo? (queremos que NO lo sea)
allFitSEM[c("chisq", "pvalue")] # pvalue>0.05## $chisq
## [1] 214.1208
##
## $pvalue
## [1] 0No pasó .
- Test 2: El Índice Tucker Lewis es mayor a 0.9? (queremos que sea mayor).
allFitSEM$tli # > 0.90## [1] 0.8182477No pasó.
- La Raíz del error cuadrático medio de aproximación es menor a 0.05?
allFitSEM[c('rmsea.ci.lower','rmsea' ,'rmsea.ci.upper')] ## $rmsea.ci.lower
## [1] 0.1474601
##
## $rmsea
## [1] 0.1694822
##
## $rmsea.ci.upper
## [1] 0.1922318Tampoco pasó, pues 0.05 es menor que el rmsea.ci.lower, y debería ser mayor que el rmsea.ci.upper.
A esta altura, ya sabemos que el SEM en su totalidad nos informa que aun cuando por partes parece haber correspondencias, en su totalidad la teoría planteada tiene problemas para ser generalizada con la data disponible.
Con lavPredict() podemos calcular las latentes, y podemos luego añadirlas a la data original:
HappyDemo=as.data.frame(cbind(HappyDemo,lavPredict(fit)))Las latentes tienen valores muy diferentes a los de los scores, pero podemos transformarlos a ese rango:
library(BBmisc)
# el rango nuevo usara como minimo y maximo los valores del Score que vino con la data (eso lo hago con min y max):
HappyDemo$tranquilidad=normalize(HappyDemo$tranquilidad ,
method = "range",
margin=2, # by column
range = c(min(HappyDemo$ScoreHappy),
max(HappyDemo$ScoreHappy)))
HappyDemo$felicidad=normalize(HappyDemo$felicidad,
method = "range",
margin=2, # by column
range = c(min(HappyDemo$ScoreHappy),
max(HappyDemo$ScoreHappy)))
HappyDemo$democracia=normalize(HappyDemo$democracia,
method = "range",
margin=2, # by column
range = c(min(HappyDemo$ScoreDemo),
max(HappyDemo$ScoreDemo)))Veamos un resumen estadistico:
summary(HappyDemo)## Country ScoreHappy GDP Social
## Length:147 Min. :3.083 Min. :0.026 Min. :0.000
## Class :character 1st Qu.:4.541 1st Qu.:0.615 1st Qu.:1.057
## Mode :character Median :5.386 Median :0.983 Median :1.293
## Mean :5.423 Mean :0.913 Mean :1.216
## 3rd Qu.:6.190 3rd Qu.:1.221 3rd Qu.:1.453
## Max. :7.769 Max. :1.684 Max. :1.624
## Healthy Freedom Generosity Perceptions
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.5550 1st Qu.:0.3100 1st Qu.:0.1085 1st Qu.:0.0460
## Median :0.7980 Median :0.4170 Median :0.1770 Median :0.0850
## Mean :0.7325 Mean :0.3947 Mean :0.1855 Mean :0.1109
## 3rd Qu.:0.8875 3rd Qu.:0.5080 3rd Qu.:0.2520 3rd Qu.:0.1415
## Max. :1.1410 Max. :0.6310 Max. :0.5660 Max. :0.4530
## ScoreDemo Electoral Functioning Politicalparticipation
## Min. :1.320 Min. : 0.000 Min. :0.000 Min. : 1.110
## 1st Qu.:3.630 1st Qu.: 3.125 1st Qu.:2.930 1st Qu.: 3.890
## Median :5.810 Median : 7.000 Median :5.360 Median : 5.560
## Mean :5.593 Mean : 6.041 Mean :4.978 Mean : 5.367
## 3rd Qu.:7.335 3rd Qu.: 9.170 3rd Qu.:6.790 3rd Qu.: 6.670
## Max. :9.870 Max. :10.000 Max. :9.640 Max. :10.000
## Politicalculture Civilliberties Regimetype Region
## Min. : 1.880 Min. : 0.000 Length:147 Length:147
## 1st Qu.: 4.380 1st Qu.: 3.820 Class :character Class :character
## Median : 5.630 Median : 6.180 Mode :character Mode :character
## Mean : 5.653 Mean : 5.929
## 3rd Qu.: 6.250 3rd Qu.: 8.240
## Max. :10.000 Max. :10.000
## democracia tranquilidad felicidad
## Min. :1.320 Min. :3.083 Min. :3.083
## 1st Qu.:3.977 1st Qu.:4.925 1st Qu.:5.374
## Median :6.291 Median :5.618 Median :6.215
## Mean :5.997 Mean :5.614 Mean :6.087
## 3rd Qu.:7.824 3rd Qu.:6.196 3rd Qu.:6.868
## Max. :9.870 Max. :7.769 Max. :7.769En SEM la regresión salió bien; ahora que cuento con las latentes puedo correr la regresión con el metodo clásico:
summary(lm(felicidad~tranquilidad+democracia,
data = HappyDemo))##
## Call:
## lm(formula = felicidad ~ tranquilidad + democracia, data = HappyDemo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.54553 -0.42162 0.03224 0.43197 1.40098
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.52099 0.31646 7.966 4.52e-13 ***
## tranquilidad 0.38345 0.06326 6.061 1.13e-08 ***
## democracia 0.23563 0.02797 8.425 3.36e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6637 on 144 degrees of freedom
## Multiple R-squared: 0.5855, Adjusted R-squared: 0.5797
## F-statistic: 101.7 on 2 and 144 DF, p-value: < 2.2e-16Esto no contradice al SEM, pero no tiene la información que SEM tenía, aqui la funcion lm no sabe que está usando latentes. Correr un SEM es mejor que usar regresion lineal ante la presencia de latentes y cuando queremos probar TODO en simultaneo. Esto es muy exigente, y lo visto aquí es sólo muy poco del tema.
Veamos visualmente los scores originales y los calculados con SEM, pues debería haber correlación:
library(ggplot2)
library(ggpubr)
p1=ggscatter(HappyDemo,
x = "democracia", y = "ScoreDemo",
cor.coef = TRUE,
cor.method = "pearson") # spearman?
p2=ggscatter(HappyDemo,
x = "tranquilidad", y = "ScoreHappy",
cor.coef = TRUE,
cor.method = "pearson")
p3=ggscatter(HappyDemo,
x = "felicidad", y = "ScoreHappy",
cor.coef = TRUE,
cor.method = "pearson")
ggarrange(p1, p2, p3,
ncol = 3)
Usando clusters
Podria usar clusters en la regresión? Me puede servir para crear una ordinal:
set.seed(123)
# creemos clusters con las variable relacionadas con felicidad
library(cluster)
dist=daisy(HappyDemo[,3:8],stand = T)
res.pam=pam(x=dist,diss = T,k=3)
HappyDemo$pamhappy=res.pam$clusteringVerifiquemos cómo rankea el valor del cluster:
aggregate(ScoreHappy~pamhappy,data = HappyDemo,FUN = mean)## pamhappy ScoreHappy
## 1 1 4.322021
## 2 2 5.649377
## 3 3 6.913522Verificamos que la etiqueta no requiere recodificar, pues va como ascendente. Continuamos:
#como ordinal!
HappyDemo$pamhappy=as.ordered(HappyDemo$pamhappy)
#usemos regresión para otra hipotesis:
hipotesis2=formula(ScoreDemo~pamhappy)
# regresion
regresion2=lm(hipotesis2,data=HappyDemo)
summary(regresion2)##
## Call:
## lm(formula = hipotesis2, data = HappyDemo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.625 -1.018 0.517 1.397 2.799
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.87621 0.17269 34.028 < 2e-16 ***
## pamhappy.L 2.49887 0.33433 7.474 6.93e-12 ***
## pamhappy.Q -0.02055 0.25913 -0.079 0.937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.858 on 144 degrees of freedom
## Multiple R-squared: 0.2944, Adjusted R-squared: 0.2846
## F-statistic: 30.05 on 2 and 144 DF, p-value: 1.243e-11Observa que el efecto lineal (pamhappy.L) es significativo, pero el cuadrático no lo es (pamhappy.Q). El lineal nos importará, pues quiere decir el efecto de pasar del nivel 1 al 2, y del 2 al 3. Esto sugiere que si un país esta en el grupo 1 de felicidad, su indice de democracia en promedio sube en 2.49 si llega a subir al nivel 2 de felicidad (lo mismo puedes decir si pasa del 2 al 3).
Y si mi dependiente es categórica?
Cuando la independiente es categórica, hacemos uso de la regresión logística. Veamos sólo el caso con dependiente dicotómica.
# es feliz?
HappyDemo$esfeliz=HappyDemo$felicidad>median(HappyDemo$felicidad)La nueva variable es dicotomica:
table(HappyDemo$esfeliz)##
## FALSE TRUE
## 74 73Apliquemos la regresión logística, utilizando las latentes calculadas (menos felicidad):
hipotesis3=formula(esfeliz~democracia + tranquilidad)
regresion3=glm(hipotesis3, data=HappyDemo,family = binomial)
#resultado clásico:
summary(regresion3)##
## Call:
## glm(formula = hipotesis3, family = binomial, data = HappyDemo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9175 -0.7521 -0.1456 0.7238 2.3975
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.0487 1.6856 -4.775 1.80e-06 ***
## democracia 0.6774 0.1239 5.469 4.53e-08 ***
## tranquilidad 0.7125 0.2761 2.580 0.00987 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 203.78 on 146 degrees of freedom
## Residual deviance: 133.92 on 144 degrees of freedom
## AIC: 139.92
##
## Number of Fisher Scoring iterations: 5Vemos que la tranquilidad y la democracia afectan la probabilidad que un país sea feliz (ambas son significativas y tiene efecto positivo). A diferencia de la regresión lineal, no podemos dar una interpretación en lenguaje sencillo a esos coeficientes.
# interpracion usando marginal effects:
library(margins)
#
(model = margins(regresion3))## democracia tranquilidad
## 0.1004 0.1056Esto quiere decir en cuanto las variables afectan a la probabilidad de ser feliz, por ejemplo, cada vez que la democracia sube un punto, la probabilidad que un país sea feliz se eleva en 0.1 (10%).