Herramientas Computacionales
para la Investigación Interdisciplinaria Reproducible
- Profesor: Dr. José Manuel Magallanes, PhD (jmagallanes@pucp.edu.pe)
Profesor del Departamento de Ciencias Sociales, Pontificia Universidad Católica del Peru.
Senior Data Scientist del eScience Institute and Visiting Professor at Evans School of Public Policy and Governance, University of Washington.
Fellow Catalyst, Berkeley Initiative for Transparency in Social Sciences, UC Berkeley.
Sesión 3: Integración
En esta sesión final, vamos a utilizar el archivo que se creó en Python, donde hay diversos indicadores polÃticos a nivel de paÃs. El archivo en cuestión será analizado en R, usando RStudio.
RStudio es un entorno para usar R. Pero lo peculiar es que no sólo permite usar un archivo clásico de sintaxis (para lo cual RStudio es menos relevante), sino que también puede producir documentos en formato HTML y PDF. En este caso, veremos como convertir un archivo de códigos simples de R en un archivo más ‘sofisticado’. Este archivo sofisticado combinará código de R y Latex, para producir un PDF.
Dado que usaremos Latex desde RStudio, podemos tomar el control detallado de lo que queremos mostrar en el PDF final. En nuestro caso, no sólo veremos detalles de cómo producir lo que R genera, sino que además lo integraremos con un archivo de referencias bibliográficas producido desde Zotero; almacenaremos todo en GitHub; y veremos cómo conectar este repositorio con Overleaf para compartir versiones intermedias con usuarios menos especializados.
Nótese que cada vez que llamo a un ‘paquete’ de R usaré el comando library; este comando sólo activa un paquete (o biblioteca) si es que esta está previamente instalada. Para instalar debes usar el menu Tools y elige la opción install.packages.
Comandos básicos en R para exploración Univariada:
# carga de datos
filename="indexes.csv"
dataidx=read.csv(filename,
stringsAsFactors = FALSE)
# ver primeras filas los datos:
head(dataidx)## Country WorldFreedom EconomicFreedom PressFreedom Democracy
## 1 Abkhazia 3 2 3 2
## 2 Afghanistan 1 2 2 1
## 3 Albania 3 3 3 2
## 4 Algeria 1 1 2 1
## 5 Andorra 5 3 4 4
## 6 Angola 1 1 2 1Los datos ya llegaron ‘limpios’ desde Python.
Como tenemos números en vez de factores, podemos preparar etiquetas para las gráficas:
# creando etiquetas:
level5=c("muy malo","malo","medio","bueno","muy bueno")
level4=c("muy malo","malo","bueno","muy bueno")
level3=c("muy malo","medio","muy bueno")Aquà podemos ver la distribución de una variable:
demoTable=table(dataidx[,5])
demoTable##
## 1 2 4 5
## 60 45 82 19De ahà que:
demoTable=table(dataidx[,5])
names(demoTable)=level4
demoTable## muy malo malo bueno muy bueno
## 60 45 82 19Ahora las frecuencias relativas:
demoTableRel=round(prop.table(demoTable)*100,1)
demoTableRel## muy malo malo bueno muy bueno
## 29.1 21.8 39.8 9.2Y aquà el plot que representa esta distribución:
title='Distribución de la Democracia'
paleta='red'
barplot(demoTableRel,main=title,
col=paleta,ylim = c(0,100),
ylab = "%")
Los mismo hacemos para ver la Libertad económica en el mundo en una tabla:
ecoTable=table(dataidx[,3])
names(ecoTable)=level5
ecoTable## muy malo malo medio bueno muy bueno
## 21 78 74 28 6Sus frecuencias relativas:
ecoTableRel=round(prop.table(ecoTable)*100,1)
ecoTableRel## muy malo malo medio bueno muy bueno
## 10.1 37.7 35.7 13.5 2.9Y aquà el plot:
title='Distribución de la Libertad Económica'
paleta='red'
barplot(ecoTableRel,main=title,
col=paleta,ylim = c(0,100),
ylab = "%")
La Libertad general en el mundo en una tabla:
worldTable=table(dataidx[,2])
names(worldTable)=level3
worldTable## muy malo medio muy bueno
## 55 62 89Ahora las frecuencias relativas:
worldTableRel=round(prop.table(worldTable)*100,1)
worldTableRel## muy malo medio muy bueno
## 26.7 30.1 43.2Y aquà el plot:
title='Distribución de la Libertad en el Mundo'
paleta='red'
barplot(worldTableRel,main=title,
col=paleta,ylim = c(0,100),
ylab = "%")
La Libertad de prensa en el mundo en una tabla:
pressTable=table(dataidx[,4])
names(pressTable)=level5
pressTable## muy malo malo medio bueno muy bueno
## 22 53 66 48 17Sus frecuencias relativas:
pressTableRel=round(prop.table(pressTable)*100,1)
pressTableRel## muy malo malo medio bueno muy bueno
## 10.7 25.7 32.0 23.3 8.3Y aquà el plot que representa esta distribución:
title='Distribución de la Libertad de Prensa'
paleta='red'
barplot(pressTableRel,main=title,
col=paleta,ylim = c(0,100),
ylab = "%")
Podemos mostrar los estadÃsticos de cada variable:
summary(dataidx[,-1])## WorldFreedom EconomicFreedom PressFreedom Democracy
## Min. :1.00 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.00 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:1.000
## Median :3.00 Median :3.000 Median :3.000 Median :2.000
## Mean :3.33 Mean :2.614 Mean :2.927 Mean :2.782
## 3rd Qu.:5.00 3rd Qu.:3.000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :5.00 Max. :5.000 Max. :5.000 Max. :5.000
## NA's :1 NA's :1 NA's :1Comandos básicos en R para exploración Bivariada:
Si asumimos que estamos interesados en el impacto de los otros indices en el nivel de Democracia, podemos primero ver la relación que tiene esta variable con todas las demás:
explanans=names(dataidx)[c(2:4)]
corrDem=cor(x=dataidx[,5],
y=dataidx[,explanans],
use = "na.or.complete",
method = "spearman")
corrDem## WorldFreedom EconomicFreedom PressFreedom
## [1,] 0.8964068 0.6021479 0.7821048Veamos la correlación entre las variables independientes:
corrTable=round(cor(dataidx[explanans],
use = "na.or.complete"),2)
# ocultar media matriz
corrTable[upper.tri(corrTable)]<-""
as.data.frame(corrTable)## WorldFreedom EconomicFreedom PressFreedom
## WorldFreedom 1
## EconomicFreedom 0.49 1
## PressFreedom 0.83 0.53 1Finalmente, vemos los modelos propuestos. Primero sin la libertad mundial como independiente:
LinRegA = lm(Democracy ~ ., data = dataidx[,c(3:5)])
summary(LinRegA)##
## Call:
## lm(formula = Democracy ~ ., data = dataidx[, c(3:5)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.99066 -0.61319 0.05363 0.43110 2.22022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.64197 0.19912 -3.224 0.00147 **
## EconomicFreedom 0.37747 0.07736 4.879 2.15e-06 ***
## PressFreedom 0.83341 0.06509 12.804 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.88 on 203 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6371, Adjusted R-squared: 0.6335
## F-statistic: 178.2 on 2 and 203 DF, p-value: < 2.2e-16Luego con la libertad mundial
LinRegB = lm(Democracy ~ ., data = dataidx[,c(2:5)])
summary(LinRegB)##
## Call:
## lm(formula = Democracy ~ ., data = dataidx[, c(2:5)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.78162 -0.36268 -0.07215 0.30011 1.91679
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.35412 0.13782 -2.569 0.0109 *
## WorldFreedom 0.70394 0.04642 15.164 < 2e-16 ***
## EconomicFreedom 0.29053 0.05335 5.446 1.49e-07 ***
## PressFreedom 0.01166 0.07020 0.166 0.8683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6033 on 202 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8303, Adjusted R-squared: 0.8278
## F-statistic: 329.4 on 3 and 202 DF, p-value: < 2.2e-16Otros plots importantes
Podemos pintar un mapa con la información que tenemos. Veamos primero el mapa:
library(rgdal)## Loading required package: sp## rgdal: version: 1.2-8, (SVN revision 663)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4folder='world_map'
file='world_map.shp'
mapaFile=file.path(folder,file)
mapaMundo = rgdal::readOGR(mapaFile,stringsAsFactors=F) ## OGR data source with driver: ESRI Shapefile
## Source: "world_map/world_map.shp", layer: "world_map"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings: POP2005plot(mapaMundo)
Para pintarlo, podemos crear conglomerados. Tenemos para ello que juntar añadir la información de los Ãndices al mapa (no al revés). Veamos cuales son las claves:
head(mapaMundo@data)## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGION
## 0 AC AG ATG 28 Antigua and Barbuda 44 83039 19 29
## 1 AG DZ DZA 12 Algeria 238174 32854159 2 15
## 2 AJ AZ AZE 31 Azerbaijan 8260 8352021 142 145
## 3 AL AL ALB 8 Albania 2740 3153731 150 39
## 4 AM AM ARM 51 Armenia 2820 3017661 142 145
## 5 AO AO AGO 24 Angola 124670 16095214 2 17
## LON LAT
## 0 -61.783 17.078
## 1 2.632 28.163
## 2 47.395 40.430
## 3 20.068 41.143
## 4 44.563 40.534
## 5 17.544 -12.296Usamos el campo respectivo para juntar los datos:
# añadiendo información de indices al mapa:
mapaMundoAll=merge(mapaMundo,
dataidx,
by.x='NAME',
by.y='Country',all.x=F)De ahà hagamos varios pasos:
# nombres de las variables a utilizar
dimensions=c("EconomicFreedom","PressFreedom","Democracy")#creando subconjunto
dataCluster=mapaMundoAll@data[,c(dimensions)]# indicando que la data numerica es ordinal:
dataCluster=as.data.frame(lapply(dataCluster,as.ordered))# llamando librerÃa:
library(cluster)
# creando matriz de distancias
dist=daisy(dataCluster,metric = "gower")
# aplicando algoritmo
pam_fit <- pam(dist, diss = TRUE, k = 3)
# añadiendo los clusters al mapa:
mapaMundoAll$cluster=pam_fit$clusteringYa tenemos los clusters, pero hay que recordar que los números asignados no son necesariamente ‘reveladores’:
aggregate(Democracy~cluster,
data=mapaMundoAll,
FUN=mean,
na.rm=T)## cluster Democracy
## 1 1 3.805195
## 2 2 1.465753
## 3 3 4.266667Luego, debemos recodificar la columna cluster:
library(car)
mapaMundoAll$cluster<-recode(mapaMundoAll$cluster,
"1=2;2=3;3=1")Con los clusters calculados, podemos pintar el mundo:
# que se pintara:
varToPlot=mapaMundoAll$cluster
#cuantos colores:
numberOfColors = length(unique(varToPlot))
#qué colores:
library(RColorBrewer)
colorForScale='Set2'
paleta = brewer.pal(numberOfColors, colorForScale)
# a dibujar:
plot(mapaMundo,col='grey',border=0)
plot(mapaMundoAll, col = paleta[varToPlot],border=F,add=T)
legend('left', legend = c("TOP","MEDIUM","LOW"),
fill = paleta,
cex = 0.6,
bty = "n",
title="Conglomerado")
AUSPICIO:
El desarrollo de estos contenidos ha sido posible gracias al grant del Berkeley Initiative for Transparency in the Social Sciences (BITSS) at the Center for Effective Global Action (CEGA) at the University of California, Berkeley
RECONOCIMIENTO
El autor reconoce el apoyo que el eScience Institute de la Universidad de Washington le ha brindado desde el 2015 para desarrollar su investigación en Ciencia de Datos.
