Facultad de Derecho y Ciencia Politica
Escuela de Ciencia Política
Guia de Clase de ESTADISTICA
Exploración de Variables Numéricas
Nuestros datos muestran los accidentes ocurridos en las empresas mineras del 2000 al 2006:
Traigamos esos datos:
link2='https://docs.google.com/spreadsheets/d/e/2PACX-1vTgDnZeJUe5Qdn-Il3Ob1u630GIYyxf5nomPu3oqrb5L31vUYCHsrtA0_tl2gEQ82mOI7-l1B9TtyZo/pub?gid=753961441&single=true&output=csv'
mineria=read.csv(link2, stringsAsFactors = F)
# que tenemos
str(mineria)## 'data.frame': 33 obs. of 4 variables:
## $ EMPRESA : chr "Ananea" "Antamina" "Ares" "Argentum" ...
## $ Responsable: chr "Contrata" "Empresa" "Empresa" "Contrata" ...
## $ Tama : chr "Pequeña" "Grande" "Mediana" "Mediana" ...
## $ accidentes : int 20 5 10 6 5 19 5 23 6 16 ...Parte 1. Exploración Gráfica
NOTA: Las variables numéricas no requieren tablas de frecuencia para ser exploradas.
- El grafico inicial a usar es el histograma:
library(ggplot2)
base1=ggplot(mineria,aes(x=accidentes))
histNum= base1 + geom_histogram(bins=7)
histNum 
La exploración numérica nos debe sugerir:
- Si la media es representativo o no.
- Si hay asimetría, y hacia donde se concentran los datos.
- Si hay valores atípicos.
Del gráfico podemos decir que, como no es simétrico, informar la media de accidentes no será muy representativa (la mediana será la mejor opción); por otro lado, como la asimetría nos muestra que los accidentes son en su mayoría pocos. Hay un numero de accidentes que se aleja hacia valores altos, pero no estamos seguros si está tan lejos del ‘centro’ para ser atípico.
Los atípicos se ven claramente en un boxplot:
base2=ggplot(mineria,aes(y=accidentes))
box=base2 + geom_boxplot() + coord_flip()
box 
El punto a la derecha confirma que hay atípicos.
Para ser más preciso en nuestra exploración, debemos calcular diversos indicadores estadísticos.
Parte 2. Exploración con Estadígrafos
Los estadigrafos aparecen rapidamente así:
summary(mineria$accidentes)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.00 5.00 8.00 10.76 13.00 42.00Podriamos graficarlos en el boxplot:
library(ggplot2)
estadigrafos=round(as.vector(summary(mineria$accidentes)),2)
box + scale_y_continuous(breaks = estadigrafos) 
¿Te das cuenta por que hay asimetría?
Si la media es mayor que la mediana la asimetría tiende a ser positiva (cola a la derecha). Hay tendencia a la asimetría negativa (cola a la izquierda) cuando la mediana es mayor que la media. Aquí se nota claramente que hay asimetría, pero podemos confirmarla calculando el coeficiente respectivo:
library(DescTools)
Skew(mineria$accidentes,conf.level = 0.05)## skew lwr.ci upr.ci
## 2.251112 2.661165 2.737899La distancia intercuartilica es importante saberla:
IQR(mineria$accidentes)## [1] 8Es decir, entre el primer y tercer cuartil hay sólo 8 valores; asi, el 50% de los valores centrales varian en 8 valores. Así, podemos proponer que un atípico es aquel que está a una distancia lejana de estos valores centrales. Este umbral, tradicionalmente se calcula así:
# cuartil tres
q3=as.numeric(summary(mineria$accidentes)[5])
# calculando umbral (distancia del q3)
umbral= q3+1.5*IQR(mineria$accidentes)
umbral## [1] 25Osea, en teoria, todo valor mayor que 25 será considerado un atípico:
mineria[mineria$accidentes>umbral,]## EMPRESA Responsable Tama accidentes
## 31 Volcan Contrata Grande 42Los accidentes pueden ser representados en el Gini:
Gini(mineria$accidentes,conf.level=0.95)## Warning in norm.inter(t, adj.alpha): extreme order statistics used as
## endpoints## gini lwr.ci upr.ci
## 0.3362676 0.2694925 0.4576164Si el Gini es 0, todas las empresas causan la misma cantidad de accidentes; si fuera 1, una sola empresa causa todos los accidentes.
Normalmente al Gini le acompaña la curva de Lorenz:
library(gglorenz) # instalar ## Registered S3 methods overwritten by 'ineq':
## method from
## plot.Lc DescTools
## lines.Lc DescToolsbase1 + gglorenz::stat_lorenz(color='red') +
geom_abline(linetype = "dashed") + coord_fixed() +
labs(x = "% Empresas ordenadas por accidentes causados",
y = "% Acumulado de Accidentes",
title = "Relación empresa / accidente",
caption = "Fuente: MINEM")
Si la curva se acerca a la diagonal, hay igualdad de distribución: cada empresa contribuye con la misma cantidad de accidentes (Gini = 0).
EJEMPLO CON MULTIPLES COLUMNAS:
LINK="https://docs.google.com/spreadsheets/d/e/2PACX-1vQw5yDzY3trmUGrF-rDJGCLxDUD51L9NwUyhrZ6QeZjzLTL8XSsQpOciXec6Bs9MWoDUKzKAiikixJl/pub?gid=1363658439&single=true&output=csv"
idh=read.csv(LINK,stringsAsFactors = F)#install:"summarytools"
summarytools::descr(idh[,-1]) ## Registered S3 method overwritten by 'pryr':
## method from
## print.bytes Rcpp## Warning in fun(libname, pkgname): couldn't connect to display ":0"## Descriptive Statistics
## idh
## N: 24
##
## X1993 X2000 X2007 X2012
## ----------------- -------- -------- -------- --------
## Mean 0.54 0.58 0.60 0.45
## Std.Dev 0.10 0.07 0.04 0.09
## Min 0.37 0.46 0.54 0.30
## Q1 0.46 0.52 0.57 0.38
## Median 0.53 0.57 0.60 0.44
## Q3 0.60 0.62 0.64 0.53
## Max 0.75 0.75 0.68 0.63
## MAD 0.11 0.08 0.04 0.09
## IQR 0.14 0.10 0.07 0.14
## CV 0.18 0.13 0.06 0.20
## Skewness 0.28 0.32 0.27 0.39
## SE.Skewness 0.47 0.47 0.47 0.47
## Kurtosis -0.63 -0.66 -1.12 -0.90
## N.Valid 24.00 24.00 24.00 24.00
## Pct.Valid 100.00 100.00 100.00 100.00#install: "psych"
psych::describe(idh[,-1])## vars n mean sd median trimmed mad min max range skew kurtosis
## X1993 1 24 0.54 0.10 0.53 0.54 0.11 0.37 0.75 0.38 0.28 -0.63
## X2000 2 24 0.58 0.07 0.57 0.57 0.08 0.46 0.75 0.29 0.32 -0.66
## X2007 3 24 0.60 0.04 0.60 0.60 0.04 0.54 0.68 0.14 0.27 -1.12
## X2012 4 24 0.45 0.09 0.44 0.45 0.09 0.30 0.63 0.33 0.39 -0.90
## se
## X1993 0.02
## X2000 0.02
## X2007 0.01
## X2012 0.02stargazer::stargazer(idh[,-1], type = "html")| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| X1993 | 24 | 0.540 | 0.096 | 0.370 | 0.465 | 0.602 | 0.750 |
| X2000 | 24 | 0.576 | 0.075 | 0.460 | 0.518 | 0.622 | 0.750 |
| X2007 | 24 | 0.603 | 0.038 | 0.540 | 0.570 | 0.635 | 0.680 |
| X2012 | 24 | 0.453 | 0.091 | 0.300 | 0.387 | 0.525 | 0.630 |
statsIDH=psych::describe(idh[,-1])
stargazer::stargazer(statsIDH, type = "html",summary = FALSE)| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
| X1993 | 1 | 24 | 0.540 | 0.096 | 0.530 | 0.536 | 0.111 | 0.370 | 0.750 | 0.380 | 0.279 | -0.629 | 0.020 |
| X2000 | 2 | 24 | 0.576 | 0.075 | 0.575 | 0.574 | 0.082 | 0.460 | 0.750 | 0.290 | 0.319 | -0.664 | 0.015 |
| X2007 | 3 | 24 | 0.603 | 0.038 | 0.600 | 0.602 | 0.044 | 0.540 | 0.680 | 0.140 | 0.273 | -1.115 | 0.008 |
| X2012 | 4 | 24 | 0.453 | 0.091 | 0.440 | 0.450 | 0.089 | 0.300 | 0.630 | 0.330 | 0.388 | -0.897 | 0.019 |
library(knitr)
library(kableExtra)
kable(statsIDH, format = "html", digits = 2)%>%
kable_styling()| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| X1993 | 1 | 24 | 0.54 | 0.10 | 0.53 | 0.54 | 0.11 | 0.37 | 0.75 | 0.38 | 0.28 | -0.63 | 0.02 |
| X2000 | 2 | 24 | 0.58 | 0.07 | 0.57 | 0.57 | 0.08 | 0.46 | 0.75 | 0.29 | 0.32 | -0.66 | 0.02 |
| X2007 | 3 | 24 | 0.60 | 0.04 | 0.60 | 0.60 | 0.04 | 0.54 | 0.68 | 0.14 | 0.27 | -1.12 | 0.01 |
| X2012 | 4 | 24 | 0.45 | 0.09 | 0.44 | 0.45 | 0.09 | 0.30 | 0.63 | 0.33 | 0.39 | -0.90 | 0.02 |
reshape::melt(idh,id.vars='X')## X variable value
## 1 Amazonas X1993 0.47
## 2 Ancash X1993 0.53
## 3 Apurímac X1993 0.40
## 4 Arequipa X1993 0.65
## 5 Ayacucho X1993 0.42
## 6 Cajamarca X1993 0.45
## 7 Cusco X1993 0.48
## 8 Huancavelica X1993 0.37
## 9 Huánuco X1993 0.45
## 10 Ica X1993 0.64
## 11 Junín X1993 0.54
## 12 La Libertad X1993 0.58
## 13 Lambayeque X1993 0.59
## 14 Lima X1993 0.75
## 15 Loreto X1993 0.53
## 16 Madre de Dios X1993 0.60
## 17 Moquegua X1993 0.62
## 18 Pasco X1993 0.52
## 19 Piura X1993 0.53
## 20 Puno X1993 0.45
## 21 San Martin X1993 0.54
## 22 Tacna X1993 0.71
## 23 Tumbes X1993 0.61
## 24 Ucayali X1993 0.52
## 25 Amazonas X2000 0.52
## 26 Ancash X2000 0.58
## 27 Apurímac X2000 0.46
## 28 Arequipa X2000 0.64
## 29 Ayacucho X2000 0.49
## 30 Cajamarca X2000 0.50
## 31 Cusco X2000 0.54
## 32 Huancavelica X2000 0.46
## 33 Huánuco X2000 0.49
## 34 Ica X2000 0.67
## 35 Junín X2000 0.58
## 36 La Libertad X2000 0.61
## 37 Lambayeque X2000 0.63
## 38 Lima X2000 0.75
## 39 Loreto X2000 0.56
## 40 Madre de Dios X2000 0.62
## 41 Moquegua X2000 0.67
## 42 Pasco X2000 0.58
## 43 Piura X2000 0.55
## 44 Puno X2000 0.51
## 45 San Martin X2000 0.55
## 46 Tacna X2000 0.68
## 47 Tumbes X2000 0.62
## 48 Ucayali X2000 0.57
## 49 Amazonas X2007 0.57
## 50 Ancash X2007 0.60
## 51 Apurímac X2007 0.56
## 52 Arequipa X2007 0.65
## 53 Ayacucho X2007 0.56
## 54 Cajamarca X2007 0.56
## 55 Cusco X2007 0.58
## 56 Huancavelica X2007 0.54
## 57 Huánuco X2007 0.57
## 58 Ica X2007 0.65
## 59 Junín X2007 0.60
## 60 La Libertad X2007 0.62
## 61 Lambayeque X2007 0.62
## 62 Lima X2007 0.68
## 63 Loreto X2007 0.59
## 64 Madre de Dios X2007 0.63
## 65 Moquegua X2007 0.65
## 66 Pasco X2007 0.59
## 67 Piura X2007 0.60
## 68 Puno X2007 0.56
## 69 San Martin X2007 0.59
## 70 Tacna X2007 0.65
## 71 Tumbes X2007 0.65
## 72 Ucayali X2007 0.60
## 73 Amazonas X2012 0.38
## 74 Ancash X2012 0.44
## 75 Apurímac X2012 0.34
## 76 Arequipa X2012 0.58
## 77 Ayacucho X2012 0.33
## 78 Cajamarca X2012 0.38
## 79 Cusco X2012 0.44
## 80 Huancavelica X2012 0.30
## 81 Huánuco X2012 0.37
## 82 Ica X2012 0.54
## 83 Junín X2012 0.45
## 84 La Libertad X2012 0.47
## 85 Lambayeque X2012 0.46
## 86 Lima X2012 0.63
## 87 Loreto X2012 0.40
## 88 Madre de Dios X2012 0.56
## 89 Moquegua X2012 0.62
## 90 Pasco X2012 0.41
## 91 Piura X2012 0.43
## 92 Puno X2012 0.39
## 93 San Martin X2012 0.44
## 94 Tacna X2012 0.56
## 95 Tumbes X2012 0.52
## 96 Ucayali X2012 0.43idhLong=reshape::melt(idh,id.vars='X')
library(ggplot2)
hist=ggplot(data=idhLong, aes(x=value)) + geom_histogram()
hist+facet_wrap(.~variable)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
EJERCICIO:
Mejore el codigo y haga un grafico adicional de su elección.
Haga el mismo analisis para cada variable del INDICE de DEMOCRACIA.