REGRESIÓN
SESIÓN 8 - REGRESIÓN
Apuntes de Hugo: http://rpubs.com/HUGO-ALONSO-VALVERDE/551842
PASO 1: PREPARANDO DATOS
- Extrayendo la data:
link='https://docs.google.com/spreadsheets/d/e/2PACX-1vQLagf4DBS-oed_VP26Rq8H1yaaw0CgDkUwR1QwZ8_f3MSk2vfXLHKQezHRfGF2gw/pub?gid=1800569641&single=true&output=csv'
PaviVotos=read.csv(link, stringsAsFactors = F)- Variables:
str(PaviVotos)## 'data.frame': 1096 obs. of 8 variables:
## $ apropiaciondolar : chr "102,17" "62,80" "0,00" "0,00" ...
## $ priorizado : chr "no" "no" "si" "no" ...
## $ pctopo : chr "14,82" "33,57" "55,42" "63,40" ...
## $ uribista : chr "no" "no" "no" "no" ...
## $ con_consejocomunal: chr "no" "no" "no" "no" ...
## $ ejecucion : chr "no" "no" "no" "no" ...
## $ poblacioncienmil : chr "20,92" "0,39" "0,11" "0,17" ...
## $ nbi : chr "12,17" "21,64" "31,25" "45,33" ...Los nombres de la variables NO deben tener espacios, tildes ni Ñs.
Las números no deben aparecer entre comillas, si eso sucede, quiere decir que el numero tiene simbolos como “comas”, elimina eso desde el original. Tambien puede ser porque hay alguna celda con algun valor inapropiado (quizas una X, un “-”, etc.), esos se deben eliminar y dejar la celda vacía.
PaviVotos$apropiaciondolar= gsub(',','',PaviVotos$apropiaciondolar)
PaviVotos$pctopo= gsub(',','',PaviVotos$pctopo)
PaviVotos$poblacioncienmil= gsub(',','',PaviVotos$poblacioncienmil)
PaviVotos$nbi=gsub(',','',PaviVotos$nbi)Las variables categoricas que se muestran, han sido leídas como texto. Hay que volverlas categóricas usando la función “as.factor”:
PaviVotos[,c(2,4:6)]=lapply(PaviVotos[,c(2,4:6)],as.factor)PaviVotos[,c(1,3,7,8)]=lapply(PaviVotos[,c(1,3,7,8)],as.numeric)- Eliminar valores perdidos
PaviVotos=PaviVotos[complete.cases(PaviVotos),]PASO 2: ESTADÍSTICOS DE LAS VARIABLES
summary(PaviVotos)## apropiaciondolar priorizado pctopo uribista con_consejocomunal
## Min. : 0.0 no:790 Min. : 0 :184 no:1001
## 1st Qu.: 0.0 si:270 1st Qu.: 624 no:325 si: 59
## Median : 0.0 Median :2039 si:551
## Mean : 844.8 Mean :2797
## 3rd Qu.: 1056.0 3rd Qu.:4570
## Max. :13264.0 Max. :9942
## ejecucion poblacioncienmil nbi
## no:1021 Min. : 0.00 Min. : 536
## si: 39 1st Qu.: 7.00 1st Qu.:2830
## Median : 14.00 Median :4127
## Mean : 41.42 Mean :4294
## 3rd Qu.: 27.00 3rd Qu.:5550
## Max. :6927.00 Max. :9881URBISTA tiene un valor extraño
levels(PaviVotos$uribista)## [1] "" "no" "si"De ahí que
levels(PaviVotos$uribista)[levels(PaviVotos$uribista)==''] = NAOtra alternativa
library(Hmisc)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, units# NO MUESTRO RESULTADO
Hmisc::describe(PaviVotos)## PaviVotos
##
## 8 Variables 1060 Observations
## ---------------------------------------------------------------------------
## apropiaciondolar
## n missing distinct Info Mean Gmd .05 .10
## 1060 0 212 0.842 844.8 1338 0 0
## .25 .50 .75 .90 .95
## 0 0 1056 3322 4232
##
## lowest : 0 3 4 15 20, highest: 9540 9664 9972 10217 13264
## ---------------------------------------------------------------------------
## priorizado
## n missing distinct
## 1060 0 2
##
## Value no si
## Frequency 790 270
## Proportion 0.745 0.255
## ---------------------------------------------------------------------------
## pctopo
## n missing distinct Info Mean Gmd .05 .10
## 1060 0 947 1 2797 2810 30.95 99.90
## .25 .50 .75 .90 .95
## 624.00 2039.00 4570.25 6780.70 8027.15
##
## lowest : 0 2 3 7 10, highest: 9397 9408 9436 9604 9942
## ---------------------------------------------------------------------------
## uribista
## n missing distinct
## 876 184 2
##
## Value no si
## Frequency 325 551
## Proportion 0.371 0.629
## ---------------------------------------------------------------------------
## con_consejocomunal
## n missing distinct
## 1060 0 2
##
## Value no si
## Frequency 1001 59
## Proportion 0.944 0.056
## ---------------------------------------------------------------------------
## ejecucion
## n missing distinct
## 1060 0 2
##
## Value no si
## Frequency 1021 39
## Proportion 0.963 0.037
## ---------------------------------------------------------------------------
## poblacioncienmil
## n missing distinct Info Mean Gmd .05 .10
## 1060 0 138 0.999 41.43 59.72 3.00 4.00
## .25 .50 .75 .90 .95
## 7.00 14.00 27.00 55.00 99.05
##
## lowest : 0 1 2 3 4, highest: 971 1286 2092 2299 6927
##
## Value 0 50 100 150 200 250 300 350 400 450
## Frequency 775 211 35 10 3 4 3 6 3 2
## Proportion 0.731 0.199 0.033 0.009 0.003 0.004 0.003 0.006 0.003 0.002
##
## Value 500 550 700 950 1300 2100 2300 6950
## Frequency 1 1 1 1 1 1 1 1
## Proportion 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##
## For the frequency table, variable is rounded to the nearest 50
## ---------------------------------------------------------------------------
## nbi
## n missing distinct Info Mean Gmd .05 .10
## 1060 0 981 1 4294 2126 1560 1974
## .25 .50 .75 .90 .95
## 2830 4127 5550 6829 7734
##
## lowest : 536 684 712 869 894, highest: 9609 9698 9732 9779 9881
## ---------------------------------------------------------------------------PASO 3: PRUEBA DE NORMALIDAD
Quiero saber:
Es la variable dependiente normal: para elegir pruebas no parametricas bivariadas. Si hay asimetría, habrán atípicos: los atípicos complican todo! 1. Verificar normalidad:
Usemos prueba Shapiro-Wilk:
shapiro.test(PaviVotos$apropiaciondolar)##
## Shapiro-Wilk normality test
##
## data: PaviVotos$apropiaciondolar
## W = 0.595, p-value < 2.2e-16Como el p-value es menor que 0.05, se confirma que la variable no es similar a la curva normal.
- Verificar asimetría y si hay atípicos:
# Variable dependiente:
# Apropiacion presupuestal en millones US$: 'apropiaciondolar'
library(DescTools) # instalar antes##
## Attaching package: 'DescTools'## The following objects are masked from 'package:Hmisc':
##
## %nin%, Label, Mean, QuantileSkew(PaviVotos$apropiaciondolar)## [1] 2.708773Ahora sabemos que hay cola a la derecha (asimetría positiva). La medida de resumen es la mediana cuando hay asimetría.
Podemos verlo en el los gráficos:
- Histograma
library(ggplot2)
base=ggplot(data=PaviVotos, aes(x=apropiaciondolar))
base+geom_histogram(bins=20) #usar 20 siempre
- Diagrama de cajas
Con el diagrama de cajas, identificamos si hay atípicos
base=ggplot(data=PaviVotos, aes(y=apropiaciondolar))
base+geom_boxplot() + coord_flip()
PASO 4. ANÁLISIS BIVARIADO
El analisis bivariado es la Variable dependiente versus todas las independientes:
# HIPOTESIS:
#### a) criterios técnicos: 'priorizado'
#### Ha: "appropiacion" tiene relacion directa con "priorizado".
#### b) clientelismo:
####### b.1 porcentaje de la oposicion:'pctopo'
#### Hb1: "appropiacion" tiene relacion directa con "pctopo"
####### b.2 municipio lo gano o no partido afin al ejecutivo: 'uribista'
#### Hb2: "appropiacion" tiene relacion directa con "uribista"
#### c) populismo:
####### c.1 los municipios en los que se han realizado Consejos Comunales: 'consejocomunal'
#### Hc1: "appropiacion" tiene relacion directa con 'consejocomunal'
# c.2 los municipios donde se habían prometido 'ejecucion'
#### Hc2: "appropiacion" tiene relacion directa con 'ejecucion'Entonces:
- Ha: “appropiacion” tiene relacion directa con “priorizado”:
Veamos la apropiación promedio por grupo:
Ha=formula(apropiaciondolar ~ priorizado)
aggregate(Ha, data=PaviVotos,median) ## priorizado apropiaciondolar
## 1 no 0.0
## 2 si 65.5Veamos gráfico de cada grupo:
base=ggplot(data=PaviVotos, aes(x=priorizado,y=apropiaciondolar))
base+geom_boxplot() + coord_flip()
Prueba no paramétrica que informa si un grupo tiene distribución que se diferencia de la otra:
wilcox.test(Ha, data=PaviVotos)##
## Wilcoxon rank sum test with continuity correction
##
## data: apropiaciondolar by priorizado
## W = 96581, p-value = 0.01152
## alternative hypothesis: true location shift is not equal to 0Por el p-valor rechazamos con una significancia de 0.01 que las medianas sean las mismas: osea, habría efecto de la variable priorizado en apropiación.
Hb1: “appropiacion” tiene relacion directa con “pctopo”
Hb1=formula(~ apropiaciondolar + pctopo)
cor.test(Hb1,data=PaviVotos,method = "spearm",exact=F)##
## Spearman's rank correlation rho
##
## data: apropiaciondolar and pctopo
## S = 205100315, p-value = 0.2796
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.03323799El p.valor nos sugiere rechazar la existencia de correlación con una significancia de 0.279, de ahi que el rho (-0.033) puede interpretarse como 0 (el Spearman va de -1 a 1).
Veamos la gráfica:
library(ggpubr)## Loading required package: magrittrggscatter(PaviVotos,
x = "pctopo", y = "apropiaciondolar",
cor.coef = TRUE,
cor.method = "spearman",
add = "reg.line",
add.params = list(color = "blue", fill = "lightgray"),
conf.int = TRUE) 
Hb2: “appropiacion” tiene relacion directa con “uribista”
Veamos estadísticos:
Hb2=formula(apropiaciondolar ~ uribista)
aggregate(Hb2, data=PaviVotos,median)## uribista apropiaciondolar
## 1 no 0
## 2 si 0Serían distribuciones iguales?
wilcox.test(Hb2,data=PaviVotos)##
## Wilcoxon rank sum test with continuity correction
##
## data: apropiaciondolar by uribista
## W = 95048, p-value = 0.09883
## alternative hypothesis: true location shift is not equal to 0Por el p-valor rechazaríamos la igualdad de medianas con una significancia de 0.09. Osea, hay efecto de la variable uribista.
Veamoslo gráficamente:
base=ggplot(data=PaviVotos, aes(x=uribista,y=apropiaciondolar))
base+geom_boxplot() + coord_flip()
Hc1: “appropiacion” tiene relacion directa con ‘consejocomunal’
Hc1=formula(apropiaciondolar ~ con_consejocomunal)
aggregate(Hc1, data=PaviVotos,median) ## con_consejocomunal apropiaciondolar
## 1 no 0
## 2 si 1692Parece que hay diferencias:
wilcox.test(Hc1,data=PaviVotos)##
## Wilcoxon rank sum test with continuity correction
##
## data: apropiaciondolar by con_consejocomunal
## W = 14004, p-value = 1.32e-13
## alternative hypothesis: true location shift is not equal to 0Por el p-valor rechazamos la igualdad de medianas con una significancia de 0.13e-14. Osea, hay efecto de la variable con_consejocomunal.
Gráficamente:
base=ggplot(data=PaviVotos, aes(x=con_consejocomunal,y=apropiaciondolar))
base+geom_boxplot() + coord_flip()
Hc2: “appropiacion” tiene relacion directa con ‘ejecucion’
Hc2=formula(apropiaciondolar ~ ejecucion)
aggregate(Hc2, data=PaviVotos,median)## ejecucion apropiaciondolar
## 1 no 0
## 2 si 99Gráficamente:
base=ggplot(data=PaviVotos, aes(x=ejecucion,y=apropiaciondolar))
base+geom_boxplot() + coord_flip()
wilcox.test(Hc2,data=PaviVotos)##
## Wilcoxon rank sum test with continuity correction
##
## data: apropiaciondolar by ejecucion
## W = 17497, p-value = 0.1612
## alternative hypothesis: true location shift is not equal to 0Por el p-valor aceptamos que las medianas no serían diferentes. Osea, NO hay efecto de la variable ejecucion.
Paso 5: Análisis Regresión
Dado que data tiene todas las variables:
names(PaviVotos)## [1] "apropiaciondolar" "priorizado" "pctopo"
## [4] "uribista" "con_consejocomunal" "ejecucion"
## [7] "poblacioncienmil" "nbi"Podemos usar:
paviModeloOrd=lm(apropiaciondolar~.,data=PaviVotos)
summary(paviModeloOrd)##
## Call:
## lm(formula = apropiaciondolar ~ ., data = PaviVotos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6020.0 -820.8 -587.7 251.2 9249.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1413.10657 167.64728 8.429 < 2e-16 ***
## priorizadosi -221.44016 118.71322 -1.865 0.06247 .
## pctopo -0.03096 0.02126 -1.456 0.14562
## uribistasi -257.30427 109.18533 -2.357 0.01867 *
## con_consejocomunalsi 1405.08489 232.44536 6.045 2.22e-09 ***
## ejecucionsi 295.64548 280.75147 1.053 0.29261
## poblacioncienmil 1.83856 0.20023 9.182 < 2e-16 ***
## nbi -0.09290 0.02941 -3.159 0.00164 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1538 on 868 degrees of freedom
## (184 observations deleted due to missingness)
## Multiple R-squared: 0.1751, Adjusted R-squared: 0.1684
## F-statistic: 26.32 on 7 and 868 DF, p-value: < 2.2e-16Por lo tanto:
Cuando un proyecto sí está priorizado la asignación presupuestal promedio disminuye (relación inversa). El efecto que un proyecto priorizado disminuye la asignación promedio en -2.21440 millones de dolares. Este efecto es significativo al 0.1.
El porcentaje de voto de la oposición no tiene efecto significativo (ni al 0.05, ni a 0.1).
Cuando un municipio tiene alcalde que sí es Uribista la asignación presupuestal promedio disminuye (relación inversa). El efecto que haya alcalde Uribista disminuye la asignación promedio en -2.57304 millones de dolares. Este efecto es significativo al 0.05.
Cuando un municipio sí tiene Consejo Comunal la asignación presupuestal promedio aumenta (relación directa). El efecto que haya Consejo Comunal aumenta la asignación promedio en 14.05085 millones de dolares. Este efecto es significativo al 0.001.
Cuando un Consejo Comunal sí recibió promesa de ejecución no tiene efecto sobre la asignación presupuetal.
PRÁCTICA DIRIGIDA - REGRESIÓN
PREGUNTAS PLANTEADAS
- Contamos con las siguientes variables:
- pobla: poblacion de la provincia
- esperanza: esperanza de vida promedio de las personas de la provincia
- accesoedu: % de poblacion de la provincia con acceso a la educacion
- percapitaf: Ingreso percapita promedio familia
- PPK: Votos PPK 2a Vuelta 2016
- FP: Votos FP 2a Vuelta 2016
- Lo han contratado para una investigación. La institucion que lo ha contratado supone que:
- H1: PPK tuvo más apoyo que FP mientras mejor acceso a la educación.
- H2: PPK tuvo más apoyo que FP mientras mayor era a la esperanza de vida.
- H3: PPK tuvo menos apoyo que FP mientras mayores ingresos.
- H4: ¿Será cierto que el promedio de ingresos es más alto donde ganó PPK y perdió Keiko?
No olvide controlar sus modelos inferenciales multivariados.
CARGANDO Y PRE PROCESANDO LA DATA
linkCSV="https://docs.google.com/spreadsheets/d/e/2PACX-1vT-XeYAxuZWm0-TXFRNosUKdat8N86GYeT-s8RcqyRU5B6DxNxe_FP0gRHN6ZWqDcHop1kjH_9_Ekh_/pub?gid=492868216&single=true&output=csv"
TABLA=read.csv(linkCSV, stringsAsFactors = F)names(TABLA)## [1] "ubiReg" "ubiProv" "ubiDis" "depa" "prov"
## [6] "dist" "pobla" "esperanza" "accesoedu" "percapitaf"
## [11] "PPK" "FP"- Revisamos nuestras variables
str(TABLA)## 'data.frame': 1834 obs. of 12 variables:
## $ ubiReg : int 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 ...
## $ ubiProv : int 10200 10200 10200 10200 10200 10200 10300 10300 10300 10300 ...
## $ ubiDis : int 10202 10201 10203 10204 10205 10206 10302 10303 10304 10305 ...
## $ depa : chr "AMAZONAS" "AMAZONAS" "AMAZONAS" "AMAZONAS" ...
## $ prov : chr "BAGUA" "BAGUA" "BAGUA" "BAGUA" ...
## $ dist : chr "ARAMANGO" "BAGUA" "COPALLIN" "EL PARCO" ...
## $ pobla : int 11587 26067 6501 1443 23820 8020 349 282 922 883 ...
## $ esperanza : chr "76,75" "74,71" "77,97" "77,42" ...
## $ accesoedu : chr "5,38" "8,33" "5,77" "6,24" ...
## $ percapitaf: chr "404,6" "662" "452,5" "551,4" ...
## $ PPK : int 1823 4949 1490 604 6282 2342 135 92 234 283 ...
## $ FP : int 3072 5809 1321 400 2059 2765 118 162 189 155 ...- antes esperanza, accesoedu y percapitaf aparecían como texto (chr) y presentablan decimales (que R no reconoce), para cambiarlo se hace lo siguiente:
- Eliminamos los decimales
TABLA$percapitaf= gsub(',','',TABLA$percapitaf)
TABLA$esperanza= gsub(',','',TABLA$esperanza)
TABLA$accesoedu= gsub(',','',TABLA$accesoedu)- Volvemos las variables de categóricas a numerales
library(readr)
TABLA[,c(8:10)]=lapply(TABLA[,c(8:10)],parse_number)IMPORTANTE SEGUIR PROCEDIMIENTO PARA EVITAR ERRORES CON GRÁFICOS
RESPONDIENDO A LAS HIPÓTESIS
Debemos analizar antes lo siguiente:
1.1) ¿PPK le ganó a FP en todas las provincias?
- ¿Cómo saber dónde ganó? Debo saber dónde PPK sacó más votos que FP, eso se calcula así:
# Debemos de saber la diferencia de votos entre PPK y FP:
TABLA$difvotos=TABLA$PPK-TABLA$FP
head(aggregate(difvotos ~ prov, data = TABLA, sum),20)## prov difvotos
## 1 ABANCAY 137
## 2 ACOBAMBA 577
## 3 ACOMAYO 4646
## 4 AIJA -397
## 5 ALTO AMAZONAS 13871
## 6 AMBO -2903
## 7 ANDAHUAYLAS -513
## 8 ANGARAES 2866
## 9 ANTA 9633
## 10 ANTABAMBA 135
## 11 ANTONIO RAYMONDI 1681
## 12 AREQUIPA 274149
## 13 ASCOPE -18314
## 14 ASUNCION 330
## 15 ATALAYA -7635
## 16 AYABACA 1627
## 17 AYMARAES -1337
## 18 AZANGARO 29420
## 19 BAGUA 2064
## 20 BARRANCA -20134Entonces, si pido un summary y veo negativos, sabre que no ganó en todas las provincias:
dataProv=aggregate(difvotos ~ prov, data = TABLA, sum)
summary(dataProv$difvotos)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -79776.0 -6174.5 -364.0 309.7 3609.5 274149.0Aquí sé en cuales:
dataProv[dataProv$difvotos>0,'prov']## [1] "ABANCAY" "ACOBAMBA"
## [3] "ACOMAYO" "ALTO AMAZONAS"
## [5] "ANGARAES" "ANTA"
## [7] "ANTABAMBA" "ANTONIO RAYMONDI"
## [9] "AREQUIPA" "ASUNCION"
## [11] "AYABACA" "AZANGARO"
## [13] "BAGUA" "CALCA"
## [15] "CALLAO" "CAMANA"
## [17] "CANAS" "CANCHIS"
## [19] "CANDARAVE" "CANGALLO"
## [21] "CARABAYA" "CARLOS FERMIN FITZCARRALD"
## [23] "CASTILLA" "CAYLLOMA"
## [25] "CELENDIN" "CHACHAPOYAS"
## [27] "CHINCHA" "CHINCHEROS"
## [29] "CHOTA" "CHUCUITO"
## [31] "CHUMBIVILCAS" "CHUPACA"
## [33] "CONDESUYOS" "CONDORCANQUI"
## [35] "COTABAMBAS" "CUSCO"
## [37] "CUTERVO" "DANIEL A. CARRION"
## [39] "DATEM DEL MARANON" "DOS DE MAYO"
## [41] "EL COLLAO" "ESPINAR"
## [43] "GENERAL SANCHEZ CERRO" "GRAU"
## [45] "HUACAYBAMBA" "HUALGAYOC"
## [47] "HUAMALIES" "HUAMANGA"
## [49] "HUANCA SANCOS" "HUANCANE"
## [51] "HUANCAVELICA" "HUANCAYO"
## [53] "HUARAZ" "HUARI"
## [55] "ILO" "ISLAY"
## [57] "JAEN" "JORGE BASADRE"
## [59] "JUNIN" "LA UNION"
## [61] "LAMPA" "LAURICOCHA"
## [63] "LIMA" "MARANON"
## [65] "MARISCAL LUZURIAGA" "MARISCAL NIETO"
## [67] "MAYNAS" "MELGAR"
## [69] "MOHO" "PALLASCA"
## [71] "PARURO" "PASCO"
## [73] "PAUCARTAMBO" "POMABAMBA"
## [75] "PUNO" "QUISPICANCHI"
## [77] "RECUAY" "RODRIGUEZ DE MENDOZA"
## [79] "SAN IGNACIO" "SAN MARTIN"
## [81] "SAN PABLO" "SAN ROMAN"
## [83] "SANTA" "SIHUAS"
## [85] "TACNA" "TALARA"
## [87] "TARATA" "TAYACAJA"
## [89] "URUBAMBA" "VICTOR FAJARDO"
## [91] "YAROWILCA" "YAULI"
## [93] "YUNGUYO"Aquí sé cuantas:
length(dataProv[dataProv$difvotos>=0,'prov'])## [1] 931.2) ¿PPK le ganó a FP en todos los departamentos?
dataDepa=aggregate(difvotos ~ depa, data = TABLA, sum)
summary(dataDepa$difvotos)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -217134 -27094 -6609 2415 23266 294889¿Cuántos ganó?
length(dataDepa[dataDepa$difvotos>0,"depa"])## [1] 122) ¿Existe correlación entre la ventaja que le sacó PPK a FP y nivel de ingresos? ¿Qué coeficiente uso?
Correlación puede ser via tecnicas paramétricas o no paramétricas. De ahí que debo decidir segun si difvotos es normal o no:
shapiro.test(TABLA$difvotos)##
## Shapiro-Wilk normality test
##
## data: TABLA$difvotos
## W = 0.46253, p-value < 2.2e-16Luego se hacer Spearman:
cor.test(TABLA$difvotos,TABLA$percapitaf, method = "spearman")## Warning in cor.test.default(TABLA$difvotos, TABLA$percapitaf, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: TABLA$difvotos and TABLA$percapitaf
## S = 1107425273, p-value = 0.0009472
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.07712924Graficamente:
library(ggpubr)
preg2=ggscatter(data=TABLA,
x = "difvotos", y = "percapitaf",
cor.coef = TRUE,
cor.method = "spearman")
preg2
Para las siguientes preguntas, hacemos una regresión:
todasLasHipo=lm(difvotos~accesoedu+esperanza+percapitaf, data=TABLA)
summary(todasLasHipo)##
## Call:
## lm(formula = difvotos ~ accesoedu + esperanza + percapitaf, data = TABLA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56071 -791 147 833 69659
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.196e+03 4.965e+02 -2.410 0.0161 *
## accesoedu 1.266e+00 5.307e-01 2.385 0.0172 *
## esperanza -7.039e-02 6.080e-02 -1.158 0.2472
## percapitaf 2.549e-01 5.409e-02 4.713 2.62e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5376 on 1830 degrees of freedom
## Multiple R-squared: 0.02853, Adjusted R-squared: 0.02693
## F-statistic: 17.91 on 3 and 1830 DF, p-value: 1.854e-113) ¿Se cumple la hipotesis 1? ¿Qué nos puedes informar sobre esa hipotesis?
H1:PPK tuvo más apoyo que FP mientras mejor acceso a la educación.
Se confirma la H1. La relación es directa y el coeficiente (180.47) es significativo al 0.1.
4) ¿Se cumple la hipotesis 2? ¿Qué nos puedes informar sobre esa hipotesis?
H2:PPK tuvo más apoyo que FP mientras mayor era a la esperanza de vida.
No se confirma la H2. La relación es inversa y el coeficiente (-105.16) es significativo al 0.001.
5) ¿Se cumple la hipotesis 3? ¿Qué nos puedes informar sobre esa hipotesis?
H3:PPK tuvo menos apoyo que FP mientras mayores ingresos.
No se confirma la H3. La relación es directa y el coeficiente (3.77) es significativo al 0.001.
Con variable de CONTROL:
todasLasHipoyControl=lm(difvotos~accesoedu+esperanza+percapitaf+pobla, data=TABLA)
summary(todasLasHipoyControl)##
## Call:
## lm(formula = difvotos ~ accesoedu + esperanza + percapitaf +
## pobla, data = TABLA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37211 -1153 -16 867 75095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.678e+03 4.900e+02 -3.425 0.000629 ***
## accesoedu 1.765e+00 5.236e-01 3.372 0.000763 ***
## esperanza -5.445e-02 5.965e-02 -0.913 0.361450
## percapitaf 3.776e-01 5.491e-02 6.877 8.34e-12 ***
## pobla -2.256e-02 2.617e-03 -8.621 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5271 on 1829 degrees of freedom
## Multiple R-squared: 0.06646, Adjusted R-squared: 0.06442
## F-statistic: 32.55 on 4 and 1829 DF, p-value: < 2.2e-16- Revisando:
- H1:PPK tuvo más apoyo que FP mientras mejor acceso a la educación.
Se confirma la H1. La relación es directa y el coeficiente (264) es significativo al 0.01.
- H2:PPK tuvo más apoyo que FP mientras mayor era a la esperanza de vida.
No se confirma la H2. La relación es inversa y el coeficiente (-94.62) es significativo al 0.001.
- H3:PPK tuvo menos apoyo que FP mientras mayores ingresos.
No se confirma la H3. La relación es directa y el coeficiente (5.003) es significativo al 0.001.
Nota que:
Sin control:Adjusted R-squared: 0.05078 Con control: Adjusted R-squared: 0.09589 Con control, el Adjusted R-squared se eleva. Ese un mejor modelo de regresión.
Lo confirmas asi:
anova(todasLasHipo,todasLasHipoyControl)## Analysis of Variance Table
##
## Model 1: difvotos ~ accesoedu + esperanza + percapitaf
## Model 2: difvotos ~ accesoedu + esperanza + percapitaf + pobla
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1830 5.2887e+10
## 2 1829 5.0822e+10 1 2.065e+09 74.317 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como el test de Analisis de Varianza es significativo (Pr(>F) = 2.2e-16), Modelo con control es mejor.
6) ¿Se cumple la hipótesis 4? ¿Será cierto que el promedio de ingresos es mas alto donde gano PPK y perdio Keiko?
Columna donde 1 es gana y 0 es no gana:
TABLA$ganaPPK=as.numeric(TABLA$difvotos>0)Calculo el promedio por grupo:
aggregate(percapitaf~ganaPPK, data = TABLA, mean)## ganaPPK percapitaf
## 1 0 3534.292
## 2 1 3628.659El promedio parece mayor, pero veamos errores:
library(ggpubr)
ggerrorplot(data=TABLA, x = "ganaPPK",
y = "percapitaf",
desc_stat = "mean_ci"
)
Como los errores se intersectan, hago prueba de contrastes:
Primero, verifico si la numerica (ingresos) es normal:
shapiro.test(TABLA$percapitaf)##
## Shapiro-Wilk normality test
##
## data: TABLA$percapitaf
## W = 0.88933, p-value < 2.2e-16Por lo anterior, corresponde la U de Mann Whitney
kruskal.test(percapitaf~ganaPPK,data=TABLA)##
## Kruskal-Wallis rank sum test
##
## data: percapitaf by ganaPPK
## Kruskal-Wallis chi-squared = 4.6157, df = 1, p-value = 0.03168Con ese p-valor (0.02989) debo rechazar que las medias sean iguales. Eso contradice el gráfico anterior. Veamos el boxplot:
ggplot(data=TABLA, aes(x=as.factor(ganaPPK), y=percapitaf)) + geom_boxplot(notch = T)
La presencia de atipicos confunde a la U de Mann Whitney, de hecho la mediana de Keiko es mas alta. En conclusión, por ambos gráficos podemos decir que la diferencia de promedios no es signigicativa.