library(kableExtra)library(rio) #usamos rio para abrir archivos de diversas extensionesWeb Scraping
Cuando el objeto es una tabla
library(rvest)
url = "https://en.wikipedia.org/wiki/Democracy_Index"
fullxpath = '/html/body/div[2]/div/div[3]/main/div[3]/div[3]/div[1]/table[5]'
url = read_html(url)
table = html_nodes(url, xpath= fullxpath)
table = html_table (table)
table =data.frame(table)Cuando no es una tabla
pagina_web ="https://www.gob.pe/institucion/mef/funcionarios"
pagina_web = read_html(pagina_web)
css_nombre="h3.text-2xl"
nombre_html = html_nodes(pagina_web,css_nombre)
nombre_texto = html_text(nombre_html)
css_cargo="p"
cargo_html = html_nodes(pagina_web,css_cargo)
cargo_texto = html_text(cargo_html)
df = data.frame(NOMBRE = nombre_texto, CARGO = cargo_texto)Iteración de scraping
url="https://www.gob.pe/institucion/mef/funcionarios?sheet="
css_cargo="p"
css_name="h3.text-2xl"
final_table = list() # list es una función para crear listas
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionfor(i in 1:3) { # INPUT
webpage <- read_html(paste0(url, i)) #obtenemos el código html de las 3 páginas
cargo_texto <- webpage %>%
html_nodes(css_cargo) %>% # obtener el código html del css del cargo
html_text() # lo convertimos en un vector de texto
name_texto <- webpage %>%
html_nodes(css_name) %>% # obtener el código html del css del name
html_text() # lo convertimos en un vector de texto
final_table[[i]] <- data.frame(NOMBRE=name_texto, CARGO=cargo_texto) # OUTPUT: Con esto estamos creando una lista con 3 data frames que contenga las 3 páginas scrapeadas
}
df = bind_rows(final_table) # Ahora convertimos todo en una base de datos2. Limpieza de datos
data="https://github.com/WendyAdrianzenRossi/Statistics/raw/main/midata.csv"
urban=import(data)Eliminar espacios en blanco
urban$COUNTRY =trimws(urban$COUNTRY,which=c("right"),whitespace = "[\\h\\v]") Extracción de patrones
library(stringr)
# Creamos objetos
PATRON="(\\-*\\d+\\.*\\d*)(?=\\%)"
COLSUCIA=urban$URBANIZATION
# UNA COLUMNA
urban$pop_urb=str_extract_all(string = COLSUCIA,
pattern= PATRON,
simplify = T)[,1]
# OTRA COLUMNA
urban$rate_urb=str_extract_all(string = COLSUCIA,
pattern=PATRON,
simplify = T)[,2]Uso de particiones
urban$pop_urb2=str_split(string = urban$URBANIZATION,
pattern = 'rate of urbanization:',
simplify = T)[,1]
urban$pop_urb2=str_split(urban$pop_urb2,
'% of total',
simplify = T)[,1]
urban$pop_urb2=str_split(urban$pop_urb2,
':',
simplify = T)[,2]3. Merge
data="https://github.com/PoliticayGobiernoPUCP/estadistica_anapol2/raw/master/DATA/corruption.csv"
corru=import(data)
data2="https://github.com/PoliticayGobiernoPUCP/estadistica_anapol2/raw/master/DATA/demo.csv"
demo=import(data2)
corru$Rank=NULL
demo$Rank=NULL
colnames(demo)[2]='ScoreDemo'
corru$Region=NULLInner join
corrgdp=merge(corru,demo) #Consiste en fusionar dos data frames en uno que contenga los elementos comunes de ambos.
str(corrgdp)## 'data.frame': 157 obs. of 14 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Angola" ...
## $ 2016Score : chr "15" "39" "34" "18" ...
## $ 2015Score : chr "11" "36" "36" "15" ...
## $ 2014Score : chr "12" "33" "36" "19" ...
## $ 2013Score : chr "8" "31" "36" "23" ...
## $ 2012Score : chr "8" "33" "34" "22" ...
## $ ScoreDemo : num 2.97 5.98 3.5 3.62 7.02 4.79 9.09 8.29 2.65 2.71 ...
## $ Electoralprocessandpluralism: num 2.92 7 2.58 1.75 9.17 5.67 10 9.58 0.5 0.83 ...
## $ Functioningofgovernment : num 1.14 4.71 2.21 2.86 5.36 4.64 8.93 7.86 2.14 3.21 ...
## $ Politicalparticipation : num 4.44 5.56 3.89 5.56 6.11 5.56 7.78 8.33 3.33 2.78 ...
## $ Politicalculture : num 2.5 5 5 5 6.25 2.5 8.75 6.88 3.75 4.38 ...
## $ Civilliberties : num 3.82 7.65 3.82 2.94 8.24 5.59 10 8.82 3.53 2.35 ...
## $ Regimetype : chr "Authoritarian" "Hybrid regime" "Authoritarian" "Authoritarian" ...
## $ Region : chr "Asia" "Europe" "Africa" "Africa" ...Outer join
corrgdp2=merge(corru,demo,all=TRUE) #combina todas las columnas de ambos conjuntos de datos en uno para todos los elemento
str(corrgdp2)## 'data.frame': 187 obs. of 14 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Angola" ...
## $ 2016Score : chr "15" "39" "34" "18" ...
## $ 2015Score : chr "11" "36" "36" "15" ...
## $ 2014Score : chr "12" "33" "36" "19" ...
## $ 2013Score : chr "8" "31" "36" "23" ...
## $ 2012Score : chr "8" "33" "34" "22" ...
## $ ScoreDemo : num 2.97 5.98 3.5 3.62 7.02 4.79 9.09 8.29 2.65 NA ...
## $ Electoralprocessandpluralism: num 2.92 7 2.58 1.75 9.17 5.67 10 9.58 0.5 NA ...
## $ Functioningofgovernment : num 1.14 4.71 2.21 2.86 5.36 4.64 8.93 7.86 2.14 NA ...
## $ Politicalparticipation : num 4.44 5.56 3.89 5.56 6.11 5.56 7.78 8.33 3.33 NA ...
## $ Politicalculture : num 2.5 5 5 5 6.25 2.5 8.75 6.88 3.75 NA ...
## $ Civilliberties : num 3.82 7.65 3.82 2.94 8.24 5.59 10 8.82 3.53 NA ...
## $ Regimetype : chr "Authoritarian" "Hybrid regime" "Authoritarian" "Authoritarian" ...
## $ Region : chr "Asia" "Europe" "Africa" "Africa" ...Left join
corrgdp3=merge(corru,demo,all.x = TRUE) #consiste en unir todas las filas del primer data frame con los valores correspondientes del segundo
str(corrgdp3)## 'data.frame': 177 obs. of 14 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Angola" ...
## $ 2016Score : chr "15" "39" "34" "18" ...
## $ 2015Score : chr "11" "36" "36" "15" ...
## $ 2014Score : chr "12" "33" "36" "19" ...
## $ 2013Score : chr "8" "31" "36" "23" ...
## $ 2012Score : chr "8" "33" "34" "22" ...
## $ ScoreDemo : num 2.97 5.98 3.5 3.62 7.02 4.79 9.09 8.29 2.65 NA ...
## $ Electoralprocessandpluralism: num 2.92 7 2.58 1.75 9.17 5.67 10 9.58 0.5 NA ...
## $ Functioningofgovernment : num 1.14 4.71 2.21 2.86 5.36 4.64 8.93 7.86 2.14 NA ...
## $ Politicalparticipation : num 4.44 5.56 3.89 5.56 6.11 5.56 7.78 8.33 3.33 NA ...
## $ Politicalculture : num 2.5 5 5 5 6.25 2.5 8.75 6.88 3.75 NA ...
## $ Civilliberties : num 3.82 7.65 3.82 2.94 8.24 5.59 10 8.82 3.53 NA ...
## $ Regimetype : chr "Authoritarian" "Hybrid regime" "Authoritarian" "Authoritarian" ...
## $ Region : chr "Asia" "Europe" "Africa" "Africa" ...Right join
corrgdp4=merge(corru,demo,all.y = TRUE) #En este caso, la combinación consiste en unir todas las filas del segundo data frame con las correspondientes en el primero
str(corrgdp4)## 'data.frame': 167 obs. of 14 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Angola" ...
## $ 2016Score : chr "15" "39" "34" "18" ...
## $ 2015Score : chr "11" "36" "36" "15" ...
## $ 2014Score : chr "12" "33" "36" "19" ...
## $ 2013Score : chr "8" "31" "36" "23" ...
## $ 2012Score : chr "8" "33" "34" "22" ...
## $ ScoreDemo : num 2.97 5.98 3.5 3.62 7.02 4.79 9.09 8.29 2.65 2.71 ...
## $ Electoralprocessandpluralism: num 2.92 7 2.58 1.75 9.17 5.67 10 9.58 0.5 0.83 ...
## $ Functioningofgovernment : num 1.14 4.71 2.21 2.86 5.36 4.64 8.93 7.86 2.14 3.21 ...
## $ Politicalparticipation : num 4.44 5.56 3.89 5.56 6.11 5.56 7.78 8.33 3.33 2.78 ...
## $ Politicalculture : num 2.5 5 5 5 6.25 2.5 8.75 6.88 3.75 4.38 ...
## $ Civilliberties : num 3.82 7.65 3.82 2.94 8.24 5.59 10 8.82 3.53 2.35 ...
## $ Regimetype : chr "Authoritarian" "Hybrid regime" "Authoritarian" "Authoritarian" ...
## $ Region : chr "Asia" "Europe" "Africa" "Africa" ...4. Regresiones
Regresión Lineal
library(rio)
data=import("https://github.com/gabrielramost/Lab_estadistica_2/blob/main/Sesi%C3%B3n%204/trabajadores.sav?raw=true")Cálculo del modelo
Leemos validez del modelo (anova), explicación del modelo (R2), y aporte de las variables
if (!requireNamespace("kableExtra", quietly = TRUE)) {
install.packages("kableExtra")
}
library(kableExtra)library(kableExtra)
library(modelsummary)
modelo1<-lm(salario_actual~educ+ antiguedad + experiencia, data)
modelo2<-lm(salario_actual~salario_inicial + educ + antiguedad + experiencia, data)
summary(modelo1)##
## Call:
## lm(formula = salario_actual ~ educ + antiguedad + experiencia,
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23029 -7861 -2403 5785 75981
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -27886.290 5529.479 -5.043 6.55e-07 ***
## educ 4004.576 210.628 19.013 < 2e-16 ***
## antiguedad 87.951 58.441 1.505 0.1330
## experiencia 11.936 5.803 2.057 0.0403 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12770 on 470 degrees of freedom
## Multiple R-squared: 0.4441, Adjusted R-squared: 0.4406
## F-statistic: 125.2 on 3 and 470 DF, p-value: < 2.2e-16models = list(modelo1,modelo2)
modelsummary(models, stars = TRUE, output = "kableExtra")| (1) | (2) | |
|---|---|---|
| (Intercept) | -27886.290*** | -16149.671*** |
| (5529.479) | (3255.470) | |
| educ | 4004.576*** | 669.914*** |
| (210.628) | (165.596) | |
| antiguedad | 87.951 | 161.486*** |
| (58.441) | (34.246) | |
| experiencia | 11.936* | -17.303*** |
| (5.803) | (3.528) | |
| salario_inicial | 1.768*** | |
| (0.059) | ||
| Num.Obs. | 474 | 474 |
| R2 | 0.444 | 0.810 |
| R2 Adj. | 0.441 | 0.809 |
| AIC | 10314.5 | 9806.4 |
| BIC | 10335.3 | 9831.4 |
| Log.Lik. | -5152.227 | -4897.195 |
| F | 125.176 | 501.450 |
| RMSE | 12717.55 | 7425.66 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La comparación de modelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasar de un modelo a otro). Como la comparación es significativa (vea el Pr(>F)), rechazamos igualdad de modelos: el modelo 2 sí reduce el error al incluir una variable más.
library(magrittr)
library(knitr)
tanova=anova(modelo1,modelo2)
tanova## Analysis of Variance Table
##
## Model 1: salario_actual ~ educ + antiguedad + experiencia
## Model 2: salario_actual ~ salario_inicial + educ + antiguedad + experiencia
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 470 7.6663e+10
## 2 469 2.6137e+10 1 5.0526e+10 906.66 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1kable(tanova,
caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 470 | 76662947710 | NA | NA | NA | NA |
| 469 | 26136575912 | 1 | 50526371798 | 906.6554 | 0 |
Tambien podemos indagar los efectos de los regresores con coeficientes estandarizados:
modelo2b=formula(scale(salario_actual)~scale(salario_inicial)+scale(educ)+scale(antiguedad)+scale(experiencia))
reg2b=lm(modelo2b,data)
model2b=list(reg2b)
modelsummary(model2b, title = "Regresion: modelo con \ncoeficientes estandarizados",
stars = TRUE,
output = "kableExtra")| (1) | |
|---|---|
| (Intercept) | 0.000 |
| (0.020) | |
| scale(salario_inicial) | 0.815*** |
| (0.027) | |
| scale(educ) | 0.113*** |
| (0.028) | |
| scale(antiguedad) | 0.095*** |
| (0.020) | |
| scale(experiencia) | -0.106*** |
| (0.022) | |
| Num.Obs. | 474 |
| R2 | 0.810 |
| R2 Adj. | 0.809 |
| AIC | 567.7 |
| BIC | 592.7 |
| Log.Lik. | -277.871 |
| F | 501.450 |
| RMSE | 0.43 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
summary(reg2b)##
## Call:
## lm(formula = modelo2b, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.73349 -0.24121 -0.07298 0.15475 2.69850
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.129e-16 2.008e-02 0.000 1
## scale(salario_inicial) 8.151e-01 2.707e-02 30.111 < 2e-16 ***
## scale(educ) 1.132e-01 2.798e-02 4.045 6.11e-05 ***
## scale(antiguedad) 9.515e-02 2.018e-02 4.715 3.19e-06 ***
## scale(experiencia) -1.060e-01 2.161e-02 -4.904 1.30e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4372 on 469 degrees of freedom
## Multiple R-squared: 0.8105, Adjusted R-squared: 0.8089
## F-statistic: 501.5 on 4 and 469 DF, p-value: < 2.2e-16Verificación de requisitos
1. Linealidad
Se asume una relación lineal de la Y y las Xs. Para evaluar la linealidad podemos utilizar un gráfico residual que evalúa los residuos y los valores Esperados de Y. La diferencia entre el valor observado de la variable dependiente (y) y el valor predicho (ŷ) se denomina residuo (e). Cada punto de datos tiene un residuo y se espera que el promedio de los residuos sea cercano a 0.
plot(modelo2,1) #Se espera que la línea roja tienda a la horizontal
O también:
mean(modelo2$residuals) # Se espera que el promedio de los residuos sea cercano a 0## [1] 2.600199e-13La falta de linearidad provocaría que el modelo no sirva para explicar las mismas variables con datos diferentes en otrs estudios.
2. Homocedasticidad
Se asume que la varianza de los errores del modelo de regresión se mantienen constantes. Para evaluarlo podemos utilizar un gráfico de Scale-location plot. También podemos utilizar el test de Breusch-Pagan. La homocedasticidad es lo contrario a la heterocedasticidad.
plot(modelo2,3) # Se espera que la línea roja debe tender a horizontal.
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(modelo2) #H0 = Modelo homocedástico, ojo aquí buscamos que sea mayor a 0.05 para cumplir con el requisito##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 50.265, df = 4, p-value = 3.179e-10La presencia de heterocedasticidad afecta el cálculo de los p-valores, lo que afectará la validez de todo el modelo.
3. Normalidad de residuos
Los residuos (la distancia entre el valor esperado y el observado de la variable) deben distribuirse de manera normal. Para evaluarlo podemos utilizar un gráfico de Residuos vs distribución normal (Normal Q-Q) , o podemos aplicar el test de Shapiro Wilk.
plot(modelo2,2) # En este caso la expectativa es que los puntos se acerquen lo más posible a la línea oblicua del gráfico que representa una distribución normal
shapiro.test(modelo2$residuals) #H0= Hay normalidad##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.86112, p-value < 2.2e-16La falta de normalidad limita la capacidad de hacer inferencias a partir de lo encontrado.
4. No multicolinealidad
La multicolinealidad es una situación en la que se presenta una fuerte correlación entre variables explicativas del modelo. Cuando el valor de la prueba VIF es mayor a 5 significa que esa variable tiene una fuerte correlación con otras variables incluidas en el modelo. Para evaluarlo podemos utilizar el Factor de Inflación de la Varianza
library(DescTools)##
## Attaching package: 'DescTools'## The following objects are masked from 'package:modelsummary':
##
## Format, Mean, Median, N, SD, VarVIF(modelo2) #>5 es problemático## salario_inicial educ antiguedad experiencia
## 1.813582 1.937010 1.007613 1.155814La presencia de la multicolinealidad no perjudica tanto el calculo de la dependiente, pero evita calcular bien el efecto de cada regresor.
5. Valores influyentes
Los valores influyentes con casos en nuestra data que no siguen el patrón general del resto de casos. No todos los valores atípicos influyen en el análisis de regresión. Aunque nuestros datos tengan valores extremos, es posible que no sean influyentes para determinar una línea de regresión. Eso significa que los resultados no serían muy diferentes si los incluyéramos o los excluyéramos del análisis.
plot(modelo2,5) #miramos valores fuera de las lineas punteadas
checModelo2=as.data.frame(influence.measures(modelo2)$is.inf)
checModelo2[checModelo2$cook.d & checModelo2$hat,] #visualizamos valores ## [1] dfb.1_ dfb.slr_ dfb.educ dfb.antg dfb.expr dffit cov.r cook.d
## [9] hat
## <0 rows> (or 0-length row.names)Si alguna fila aparece en el resultado, ese caso está afectando los cálculos de la regresión (sin él habría otro resultado).
Modelos lineales generalizados
A. Regresión Poisson (No lineal)
La regresión Poisson tiene sus supuestos:
Variable de respuesta es un conteo por unidad de tiempo o espacio, que puede ser descrita por la distribución Poisson.
Independencia: Las observaciones (filas) no deben tener relación entre sí.
Media=Varianza: Por definición, la media de una variable que se distribuye como Poisson debe ser igual a su varianza.
Linealidad: El logaritmo de la media de los datos, log(λ), debe ser una función lineal de los datos
library(rio)
data=import("https://github.com/gabrielramost/Lab_estadistica_2/blob/main/Sesi%C3%B3n%204/IOP_1119_01_02_D.sav?raw=true")
data=subset(data, SERIE== 2019)Variable dependiente: Indice de acaso
#Codificación
data$P12A[data$P12A == "2"] <- 0
data$P12B[data$P12B == "2"] <- 0
data$P12C[data$P12C == "2"] <- 0
data$P12D[data$P12D == "2"] <- 0
data$P12E[data$P12E == "2"] <- 0
data$P12F[data$P12F == "2"] <- 0
data$P12G[data$P12G == "2"] <- 0
data$P12H[data$P12H == "2"] <- 0
data$P12A[data$P12A == "98"] <- NA
data$P12B[data$P12B == "98"] <- NA
data$P12C[data$P12C == "98"] <- NA
data$P12D[data$P12D == "98"] <- NA
data$P12E[data$P12E == "98"] <- NA
data$P12G[data$P12G == "98"] <- NA
data$P12H[data$P12H == "98"] <- NA
# Creación del indice
data$indiceAC=data$P12A+data$P12B+data$P12C+data$P12D+data$P12E+data$P12F+data$P12G+data$P12HVariables independientes
data$SEXO=as.factor(data$SEXO)
data$SEXO=factor(data$SEXO,levels=levels(data$SEXO) ,labels=c("Hombre","Mujer"))
data$NIVEDU=as.factor(data$NIVEDU)
data$NIVEDU=factor(data$NIVEDU,levels=levels(data$NIVEDU) ,labels=c("Ninguna","Inicial","Primaria"
,"SecundariaIncompleta","SecundariaCompleta",
"SuperiorTecnicaincompleta",
"Superiortecnica completa","Superior universitaria incompleta"
,"Superior universitaria completa","Postgrado"))
data$NSEGRUP=as.factor(data$NSEGRUP)
data$NSEGRUP=factor(data$NSEGRUP,levels=levels(data$NSEGRUP) ,labels=c("AB","C","DE"))
data$EDAD = as.numeric(data$EDAD)Cálculo de modelos
Modelo 1
library(modelsummary)
modelo1<-glm(indiceAC~EDAD,
family=poisson (link = "log"),
data=data)
modelsummary(modelo1,
exponentiate = T,
statistic = 'conf.int',
title = "EXP() de la Regresión Poisson para Interpretación",
stars = TRUE,
output = "kableExtra")| (1) | |
|---|---|
| (Intercept) | 5.031*** |
| [4.260, 5.938] | |
| EDAD | 0.968*** |
| [0.963, 0.972] | |
| Num.Obs. | 679 |
| AIC | 2603.1 |
| BIC | 2612.1 |
| Log.Lik. | -1299.544 |
| F | 201.569 |
| RMSE | 1.99 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Nota que mientras en la regresión lineal no deseábamos que nuestro coeficiente esté cerca al cero, es decir, que su intervalo de confianza no incluya al cero, aquí no deseamos que el intervalo de confianza incluya al uno. Una vez exponenciado, podemos interpretar los coeficientes de manera más sencilla. Nótese que esta regresión propone un efecto multiplicativo sobre el valor medio de la respuesta (la regresión OLS o Gaussiana propone un efecto aditivo).
Interpretación: para el modelo se ve que por cada unidad que aumente la edad, la probabilidad de ser victima de acoso callejero se multiplica por 0.968. En otras palabras, disminuiría en 3.2 %
Modelo 2
modelo2<-glm(indiceAC~SEXO,
family=poisson (link = "log"),
data=data)
modelsummary(modelo2,
exponentiate = T,
statistic = 'conf.int',
title = "EXP() de la Regresión Poisson para Interpretación",
stars = TRUE,
output = "kableExtra")| (1) | |
|---|---|
| (Intercept) | 0.600*** |
| [0.520, 0.688] | |
| SEXOMujer | 3.766*** |
| [3.228, 4.416] | |
| Num.Obs. | 679 |
| AIC | 2492.1 |
| BIC | 2501.1 |
| Log.Lik. | -1244.044 |
| F | 275.730 |
| RMSE | 1.92 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Interpretación:
Si una persona es mujer, el número de eventos de acoso callejero va a aumentar en 3.7 veces. Cuando una persona es mujer, la probabilidad de sufrir acoso callejero va a aumentar en 376%.
Comparación de modelos
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(modelo1,modelo2, type ="text")##
## ==============================================
## Dependent variable:
## ----------------------------
## indiceAC
## (1) (2)
## ----------------------------------------------
## EDAD -0.033***
## (0.002)
##
## SEXOMujer 1.326***
## (0.080)
##
## Constant 1.616*** -0.511***
## (0.085) (0.072)
##
## ----------------------------------------------
## Observations 679 679
## Log Likelihood -1,299.544 -1,244.044
## Akaike Inf. Crit. 2,603.087 2,492.088
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Se puede comparar el LOG lIKELIHOOD, el mayor será el mejor, o el modelo que tenga menor akaike.
Tambien podemos comparar en función de las residual deviances, haciendo un anova. Nos quedamos con el segundo.
anova(modelo1,modelo2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: indiceAC ~ EDAD
## Model 2: indiceAC ~ SEXO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 677 1715.2
## 2 677 1604.2 0 111Equidispersión
Uno de los supuestos era que la media y la varianza sean iguales. Si la varianza es mayor que la media, hablamos de sobredispersión (subdispersion si es menor que la media). Para comprobar hacemos el test de equidispersión
overdispersion=AER::dispersiontest(modelo1,alternative='greater')$ p.value<0.05
underdispersion=AER::dispersiontest(modelo1,alternative='less')$ p.value<0.05
# tabla
testResult=as.data.frame(rbind(overdispersion,underdispersion))
names(testResult)='Es probable?'
testResult## Es probable?
## overdispersion TRUE
## underdispersion FALSELa Tabla muestra que es altamente improbable que la varianza sea igual a la media, por lo que se opta por aceptar que lo más probable es que tengamos sobredispersión.
VarY=data$indiceAC #Variable dependiente
list(media=mean(VarY,na.rm=T),varianza=var(VarY,na.rm=T),razon=var(VarY,na.rm=T)/mean(VarY,na.rm=T))## $media
## [1] 1.46539
##
## $varianza
## [1] 4.378963
##
## $razon
## [1] 2.988257Veamos de solucionarlo:
Quasipoisson
reg1quasi = glm(indiceAC~EDAD, family = quasipoisson(link = "log"), data = data)
stargazer(modelo1,reg1quasi,type = "text",intercept.bottom = FALSE,style="all2")##
## =========================================================
## Dependent variable:
## ----------------------------
## indiceAC
## Poisson glm: quasipoisson
## link = log
## (1) (2)
## ---------------------------------------------------------
## Constant 1.616*** 1.616***
## (0.085) (0.137)
## EDAD -0.033*** -0.033***
## (0.002) (0.004)
## ---------------------------------------------------------
## Observations 679 679
## Log Likelihood -1,299.544
## Akaike Inf. Crit. 2,603.087
## Residual Deviance (df = 677) 1,715.164 1,715.164
## Null Deviance (df = 678) 1,949.025 1,949.025
## =========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Vemos que los coeficientes no han cambiado, pero observemos los residuos:
library(arm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: Matrix## Loading required package: lme4##
## Attaching package: 'lme4'## The following object is masked from 'package:rio':
##
## factorize##
## arm (Version 1.13-1, built: 2022-8-25)## Working directory is C:/Users/USER/Downloadscbind(sePoi=se.coef(modelo1),seQuasiPoi=se.coef(reg1quasi))## sePoi seQuasiPoi
## (Intercept) 0.084757677 0.137359706
## EDAD 0.002304708 0.003735048Binomial negativa
Otra alternativa ante la sobredispersión es usar la regresión binomial negativa.
library(MASS)
nb <- glm.nb(indiceAC~EDAD, data = data)
stargazer(modelo1,reg1quasi,nb, type = "text",intercept.bottom = FALSE,style="all2")##
## ======================================================================================
## Dependent variable:
## ---------------------------------------------------------
## indiceAC
## Poisson glm: quasipoisson negative
## link = log binomial
## (1) (2) (3)
## --------------------------------------------------------------------------------------
## Constant 1.616*** 1.616*** 1.793***
## (0.085) (0.137) (0.168)
## EDAD -0.033*** -0.033*** -0.037***
## (0.002) (0.004) (0.004)
## --------------------------------------------------------------------------------------
## Observations 679 679 679
## Log Likelihood -1,299.544 -1,062.279
## theta 0.562*** (0.059) (p = 0.000)
## Akaike Inf. Crit. 2,603.087 2,128.559
## Residual Deviance (df = 677) 1,715.164 1,715.164 628.013
## Null Deviance (df = 678) 1,949.025 1,949.025 707.150
## ======================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Vemos que la binomial negativa responde mejor y corrige los errores:
library(dplyr)
library(knitr)
anova(modelo1,reg1quasi,nb, test = "Chisq") %>%
kable(caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 677 | 1715.1641 | NA | NA | NA |
| 677 | 1715.1641 | 0 | 0.000 | NA |
| 677 | 628.0129 | 0 | 1087.151 | NA |
La caída del Deviance es tanta para el último caso que la mejor opción es la binomial negativa. Por lo general, la binomial negativa es más utilizada que la quasipoisson, pero la binomial negativa no es apropiada para la subdispersión, mientras que la quasipoisson sí se usa para ese caso.
library(ggplot2)
dotwhisker::dwplot(list(Poisson1=modelo1,QuasiPoisson=reg1quasi,NegativeBinomial=nb),exp=T) +
scale_color_discrete(name="Modelos para: Acoso callejero") +
geom_vline(xintercept = 1,
colour = "grey60",
linetype = 2)
B. Regresión Logística
Ideas clave
La probabilidad: Es una medida que señala que tan posible es que ocurra un fenómeno o evento. Oscila entre 0 y 1. Cuanto más cerca de 0 menos probabilidad, y cuanto más cerca del 1 indica más probabilidad. Ejemplo: la probabilidad de ganar en un partido de futbol.
Odds: Es la probabilidad de que suceda un evento dividido por la probabilidad de que no suceda. Oscilan entre 0 e infinito y se pueden calcular para la ocurrencia del evento como para la no ocurrencia del evento. Ejemplo: la probabilidad de ganar una apuesta en un partido de futbol es 1.5 veces más probable que perder
ODDs Ratio: Es la razón entre dos odds. Permite comparar los odds de un evento en dos grupos Va de 0 a infinito. Los modelos de regresión logística están basados en probabilidades entre dos variables. Por esta razón, es útil conocer los ODDs ratio. Ejemplo: Ganar en una apuesta de un partido de futbol es 5 veces más probable que el ganar en una apuesta de voley
La regresión logística modela el comportamiento de la probabilidad del evento de interés. Es un tipo de análisis de regresión utilizado para predecir el resultado de una variable categórica ( dependiente) en función de las variables predictoras (independientes). Es útil para modelar la probabilidad de que ocurra un evento en función del efecto de un conjunto de variables.
La idea es que la regresión logística aproxime la probabilidad de obtener “0” (no ocurre cierto suceso) o “1” (ocurre el suceso) con el valor de la variable explicativa x.
iop=import("https://github.com/gabrielramost/Lab_estadistica_2/blob/main/Sesi%C3%B3n%204/IOP_0717_01_D.sav?raw=true")
data = subset(iop, select = c(NRO, P17, P3B, P8, P14, P30, P19, P92, EDAD, SEXO))library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:arm':
##
## logit## The following object is masked from 'package:DescTools':
##
## Recode## The following object is masked from 'package:dplyr':
##
## recodedata$no_cong = recode(data$P17, "2 = 1 ; 1 = 0 ; 8:9 = NA")
data$no_cong = as.factor(data$no_cong)
data$conf_con = factor(ifelse(data$P3B == 1,1,
ifelse(data$P3B == 2,1,
ifelse(data$P3B == 3,2,
ifelse(data$P3B == 4,2,0)))))
data$conf_con[data$conf_con == 0] <- NA
data$conf_con = factor(data$conf_con, levels = c(1:2), labels = c("Trust","NoTrust"))
data$sit_eco = ifelse(data$P8 == 1,1,
ifelse(data$P8 == 2,1,
ifelse(data$P8 == 3,2,
ifelse(data$P8 == 4,2,
ifelse(data$P8 == 5,2,0)))))
data$sit_eco[data$sit_eco == 0] <- NA
data$sit_eco = factor(data$sit_eco, levels = c(1:2), labels = c("Better","Worst"))
data$ideologia = recode(data$P14, "88:99 = NA")
data$demo = ifelse(data$P19 == 1,1,
ifelse(data$P19 == 2,1,
ifelse(data$P19 == 3,2,
ifelse(data$P19 == 4,2,0))))
data$demo[data$demo == 0] <- NA
data$demo = factor(data$demo, levels = c(1:2), labels = c("Yes","No"))
data$EDAD=as.numeric(data$EDAD)
data$sexo = as.factor(data$SEXO)
levels(data$sexo) = c("Hombre", "Mujer")Realizamos el modelo de regresión logística ¿Qué factores determinan la creencia de que la democracia puede funcionar sin el Congreso de la República? Vamos a modelar esa probabilidad
data = na.omit(data)Calculo mi modelo:
Modelando la probabilidad de considerar que la democracia puede funcionar sin congreso:
set.seed(2019)
modelo = glm(no_cong ~ EDAD, family = binomial, data = data)
modelor = glm(no_cong ~ EDAD + conf_con, family = binomial, data = data)
modelog = glm(no_cong ~ EDAD + conf_con + demo, family = binomial, data = data)
modelo1 = glm(no_cong ~ EDAD + conf_con + demo + sit_eco, family = binomial, data = data)
stargazer(modelo, modelor,modelog,modelo1, type = "text")##
## =========================================================
## Dependent variable:
## ---------------------------------------
## no_cong
## (1) (2) (3) (4)
## ---------------------------------------------------------
## EDAD -0.008* -0.008* -0.008* -0.008*
## (0.004) (0.004) (0.004) (0.004)
##
## conf_conNoTrust 0.444** 0.392** 0.397**
## (0.173) (0.175) (0.175)
##
## demoNo 0.450*** 0.460***
## (0.149) (0.150)
##
## sit_ecoWorst -0.099
## (0.155)
##
## Constant 0.007 -0.348 -0.633** -0.586**
## (0.182) (0.230) (0.250) (0.261)
##
## ---------------------------------------------------------
## Observations 971 971 971 971
## Log Likelihood -659.835 -656.462 -651.837 -651.630
## Akaike Inf. Crit. 1,323.670 1,318.923 1,311.673 1,313.260
## =========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Cuando queremos comparar modelos vemos el Akaike (AIC), se considera que un modelo es mejor si tiene un AIC menor a los demás.
¿Cuál es el efecto de las X sobre el Odds de la dependiente?
modelsummary(modelog,
exponentiate = T,
statistic = 'conf.int',
title = "Regresión Logísticas (Coeficientes Exponenciados)",
stars = TRUE,
output = "kableExtra")| (1) | |
|---|---|
| (Intercept) | 0.531* |
| [0.324, 0.865] | |
| EDAD | 0.992+ |
| [0.983, 1.000] | |
| conf_conNoTrust | 1.479* |
| [1.054, 2.094] | |
| demoNo | 1.568** |
| [1.173, 2.107] | |
| Num.Obs. | 971 |
| AIC | 1311.7 |
| BIC | 1331.2 |
| Log.Lik. | -651.837 |
| F | 6.223 |
| RMSE | 0.49 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Interpretación: + El odds de creer que la democracia puede funcionar sin el Congreso de la República, disminuye en 0.8% cada vez que la edad aumenta en uno (No significativa) + El odds de creer que la democracia puede funcionar sin el Congreso de la República, aumenta en 47.9% cuando no se tiene confianza en el congreso + El odds de creer que la democracia puede funcionar sin el Congreso de la República, aumenta en 56.8% cuando no se está satisfecho con la democracia
Calculamos efectos marginales (lenguaje probabilístico)
library(margins)
model = margins(modelog) #visualicemos esto en la consola
margins=summary(model)#Otra forma, aqui vemos los AME. visualicemos esto en la consola
margins$ef = margins$AME*100
margins## factor AME SE z p lower upper ef
## conf_conNoTrust 0.0920 0.0399 2.3081 0.0210 0.0139 0.1701 9.1999
## demoNo 0.1066 0.0345 3.0922 0.0020 0.0390 0.1742 10.6590
## EDAD -0.0020 0.0010 -1.8892 0.0589 -0.0040 0.0001 -0.1968Para no hablar en términos de odds sino en lenguaje probabilístico, calculamos los efectos marginales, que se interpretan de la siguiente manera:
- En promedio, no confiar en el congreso, aumenta en 0.0920 la probabilidad de creer que la democracia puede funcionar sin el Congreso de la República.
- En promedio, no estar satisfecho con la democracia, aumenta en 0.10 la probabilidad de creer que la democracia puede funcionar sin el Congreso de la República.
- En promedio, el incremento de un punto de le edad, disminuye en 0.0020 la probabilidad de creer que la democracia puede funcionar sin el Congreso de la República. (No significativa)
c. Regresión Cox
El objetivo de la regresión Cox es evaluar los impactos de diferentes variables sobre la “sobrevivencia”. Es decir, nos ayuda a saber si ciertos factores influecian la tasa de ocurrencia de un evento en específico en un lapso de tiempo.
library(survival)
library(survminer)## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomaVamos a trabajar con una base de datos de personas con cancer al pulmón:
head(lung)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0Estas son las variables:
inst: Institution code time: Survival time in days status: censoring status 1=censored, 2=dead age: Age in years sex: Male=1 Female=2 ph.ecog: ECOG performance score (0=good 5=dead) ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician pat.karno: Karnofsky performance score as rated by patient meal.cal: Calories consumed at meals wt.loss: Weight loss in last six months
Vamos a ver algunos descriptivos del tiempo de sobrevivencia:
summary(lung$time)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.0 166.8 255.5 305.2 396.5 1022.0Ahora de la ocurrencia del evento:
#Recodificamos para que 0 sea la ocurrencia:
lung$status1 = ifelse(lung$status == 2,1,0)
table(lung$status1)##
## 0 1
## 63 165Empezamos con el análisis:
lung$survival=with(lung,Surv(time = time,
event = status1))library(ggplot2)
library(ggfortify)
KM.generico = survfit(survival ~ 1, data = lung)
ejeX='DAYS\n Curva cae cuando ocurre un fallecimiento'
ejeY='Probabilidad \n(Sobrevivir)'
titulo="Curva de Sobrevivencia"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)
Lo que indica el gráfico es que la probabilidad de sobrevivir al cancer en el día 0 es del 100% (todas las personas del estudio están vivas). Esta diminuye conforme avanzan los días. Pasado el día número 500, la probabilidad de sobrevivir es menor al 25%.
Ahora probemos una hipótesis: los hombres tienen una menor probabilidad de sobrevivir.
KM_H1=formula(survival ~ sex)
KM.sex = survfit(KM_H1, data = lung)
###
ejeX='DAYS\n Curva cae cuando ocurre un fallecimiento'
ejeY='Probabilidad \n(Sobrevivir)'
titulo="Curva de Sobrevivencia"
autoplot(KM.sex,xlab=ejeX,ylab=ejeY,
main = titulo,conf.int = F) +
labs(colour = "Impacto del género") +
scale_color_discrete(labels = c("Hombre", "Mujer"))
Efectivamente, parece que la mejor no pinta tan bien para los hombres. Veamos los intervalos de confianza:
LogRank=survdiff(KM_H1, data = lung)
LogRank$pvalue## [1] 0.001311165autoplot(KM.sex,xlab=ejeX,ylab=ejeY,
main = titulo,conf.int = T)+
labs(colour = "Impacto del género") +
scale_color_discrete(labels = c("Hombre", "Mujer"))
Probemos ahora otra hipótesis: más joven es la persona, sus probabilidad de sobrevivir serán mayores. Para ver cómo juegan ambas variables a la vez, hay que evaluarlo con una regresión:
COX_H1= formula(survival~factor(sex)+age)
rcox1 <- coxph(COX_H1,data=lung)
summary(rcox1)## Call:
## coxph(formula = COX_H1, data = lung)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## factor(sex)2 -0.513219 0.598566 0.167458 -3.065 0.00218 **
## age 0.017045 1.017191 0.009223 1.848 0.06459 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## factor(sex)2 0.5986 1.6707 0.4311 0.8311
## age 1.0172 0.9831 0.9990 1.0357
##
## Concordance= 0.603 (se = 0.025 )
## Likelihood ratio test= 14.12 on 2 df, p=9e-04
## Wald test = 13.47 on 2 df, p=0.001
## Score (logrank) test = 13.72 on 2 df, p=0.001Vemos que el género tiene un impacto negativo sobre la probabilidad de fallecimiento. Es decir, el ser mujer (2) disminuye la probabilidad de ocurrencia del evento. Esto es estadísticamente significativo, mientras que la edad no muestra impactos. Pero, en tanto esta es una combinación lineal, no se pueden interpretar directamente:
(sex=abs(1-exp(coef(rcox1)[1])))## factor(sex)2
## 0.401434Esto ya se puede interpretar: ser mujer está asociado con un mejor pronóstico, reduciendo las probabilidades de ocurrencia del evento en 40%.
5. Análisis Factorial Exploratorio (EFA)
library(htmltab)
demoL=c("https://en.wikipedia.org/wiki/Democracy_Index",
'//*[@id="mw-content-text"]/div[1]/table[6]/tbody')
demo = htmltab(doc = demoL[1], which = demoL[2], encoding = "UTF-8")
# limpieza
demo[,]=lapply(demo[,], trimws,whitespace = "[\\h\\v]") # no blanks
library(stringr) # nombres simples
names(demo)=str_split(names(demo),">>",simplify = T)[,1]
names(demo)=trimws(names(demo),whitespace = "[\\h\\v]")
## Formateo
#guion invisible
names(demo)=gsub('\u00AD','',names(demo))
#capitalizar
names(demo)=str_to_title(names(demo))
#espacios por "_"
names(demo)=names(demo)=gsub('\\s','_',names(demo))
# Eliminemos columnas que no usaremos:
demo[,c(1,2,4,6)]=NULL
str(demo)## 'data.frame': 167 obs. of 7 variables:
## $ Country : chr "Norway" "New Zealand" "Iceland" "Sweden" ...
## $ Overall_Score : chr "9.81" "9.61" "9.52" "9.39" ...
## $ Electoral_Processand_Pluralism: chr "10.00" "10.00" "10.00" "9.58" ...
## $ Functioningof_Government : chr "9.64" "9.29" "9.64" "9.64" ...
## $ Politicalparticipation : chr "10.00" "10.00" "8.89" "8.33" ...
## $ Politicalculture : chr "10.00" "8.75" "9.38" "10.00" ...
## $ Civilliberties : chr "9.41" "10.00" "9.71" "9.41" ...# Tipo de variables:
## En demo:
# a numerica
demo[,-1]=lapply(demo[,-1],as.numeric)
# sin perdidos:
demo=na.omit(demo)Proceso del Analisis Factorial Exploratorio (EFA)
El análisis factorial exploratorio requiere que hagamos algunas observaciones previas.
- Calculemos matriz de correlación:
dontselect=c("Country","Regime_Type","Overall_Score")
select=setdiff(names(demo),dontselect)
theData=demo[,select] # sin los Scores ni nombre de país.
# esta es:
library(polycor)
corMatrix=polycor::hetcor(theData)$correlations- Explorar correlaciones:
- Sin evaluar significancia:
library(ggcorrplot)
ggcorrplot(corMatrix)
Si puedes ver bloques correlacionados, hay esperanza de un buen analisis factorial.
- Verificar si datos permiten factorizar:
Cuanto más cerca de 1 tenga el valor obtenido del test KMO, implica que la relación entres las variables es alta. Si KMO ≥ 0.9, el test es muy bueno; notable para KMO ≥ 0.8; mediano para KMO ≥ 0.7; bajo para KMO ≥ 0.6; y muy bajo para KMO < 0.5.17 nov 2019
library(psych)
psych::KMO(corMatrix) ## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = corMatrix)
## Overall MSA = 0.85
## MSA for each item =
## Electoral_Processand_Pluralism Functioningof_Government
## 0.80 0.89
## Politicalparticipation Politicalculture
## 0.94 0.82
## Civilliberties
## 0.82- Verificar si la matriz de correlaciones es adecuada
Aqui hay dos pruebas:
- Hnula: La matriz de correlacion es una matriz identidad
El test de esfericidad de Bartlett pone a prueba la hipótesis nula de que las variables analizadas no están correlacionadas en la muestra o, dicho de otro modo, que la matriz de correlación es la identidad (las intercorrelaciones entre las variables son cero). Diagonal principal igual a 1
cortest.bartlett(corMatrix,n=nrow(theData))$p.value>0.05## [1] FALSE- Hnula: La matriz de correlacion es una matriz singular. No posee inversa, su determinante es cero
library(matrixcalc)
is.singular.matrix(corMatrix)## [1] FALSE- Determinar en cuantos factores o variables latentes podríamos redimensionar la data:
Podemos hacer uso de multitud de representaciones a la hora de escoger el número óptimo de componentes principales:
fa.parallel(theData,fm = 'ML', fa = 'fa',correct = T)
## Parallel analysis suggests that the number of factors = 1 and the number of components = NASe sugiere 1, lo esperado, sigamos.
- Redimensionar a numero menor de factores
- Resultado inicial:
library(GPArotation)
resfa <- fa(theData,
nfactors = 1,
cor = 'mixed',
rotate = "varimax",
fm="minres")
print(resfa$loadings)##
## Loadings:
## MR1
## Electoral_Processand_Pluralism 0.914
## Functioningof_Government 0.911
## Politicalparticipation 0.827
## Politicalculture 0.620
## Civilliberties 0.971
##
## MR1
## SS loadings 3.675
## Proportion Var 0.735- Resultado mejorado (solo apropiado si hay más de un factor):
print(resfa$loadings,cutoff = 0.5)##
## Loadings:
## MR1
## Electoral_Processand_Pluralism 0.914
## Functioningof_Government 0.911
## Politicalparticipation 0.827
## Politicalculture 0.620
## Civilliberties 0.971
##
## MR1
## SS loadings 3.675
## Proportion Var 0.735Cuando logramos que cada variable se vaya a un factor, tenemos una estructura simple.
- Resultado visual:
fa.diagram(resfa)
- Evaluando Resultado obtenido:
- ¿Qué variables aportaron mas a los factores?
sort(resfa$communality)## Politicalculture Politicalparticipation
## 0.3843847 0.6836567
## Functioningof_Government Electoral_Processand_Pluralism
## 0.8291684 0.8345046
## Civilliberties
## 0.9435319- ¿Qué variables contribuyen a mas de un factor?
sort(resfa$complexity)## Politicalculture Functioningof_Government
## 1 1
## Electoral_Processand_Pluralism Civilliberties
## 1 1
## Politicalparticipation
## 1- Posibles valores proyectados:
¿Qué nombres les darías?
as.data.frame(resfa$scores)%>%head()## MR1
## 3 1.648552
## 4 1.737765
## 5 1.651807
## 6 1.552071
## 7 1.620391
## 8 1.525564El nombre sería democracia. Y debe estar correlacionado con el puntaje propuesto por TheEconomist:
demo$efa=resfa$scores
ggplot(data=demo,aes(x=Overall_Score,y=efa)) + geom_point() + theme_minimal()
Nota que los valores no son parecidos. Pero podríamos cambiar el rango:
library(BBmisc)##
## Attaching package: 'BBmisc'## The following object is masked from 'package:DescTools':
##
## %nin%## The following objects are masked from 'package:dplyr':
##
## coalesce, collapse, symdiff## The following object is masked from 'package:base':
##
## isFALSEdemo$efa_ok=normalize(demo$efa,
method = "range",
margin=2, # by column
range = c(0, 10))Esto te daría:
ggplot(data=demo,aes(x=Overall_Score,y=efa_ok)) + geom_point() + theme_minimal()
El problema sería que en la nueva variable aparecen casos 0 y casos 10. Es decir la data debe ser modificada:
# creo 5 valores con 10
top=rep(10,5)
# creo 5 valores con 0
bottom=rep(0,5)
# hago data frame con eso:
limitCases=as.data.frame(rbind(bottom,top))
# vuelvo a subsetear
theData=demo[,select]
names(limitCases)=names(theData)
#anexo
theData=rbind(theData,limitCases)
# nuevamente
resfa <- fa(theData,
nfactors = 1,
cor = 'mixed',
rotate = "varimax",
fm="minres")
scores=normalize(resfa$scores,
method = "range",
margin=2, # by column
range = c(0, 10))
#tenemos:
tail(scores)## MR1
## 169 1.60618819
## 170 0.89352156
## 171 0.09791259
## 172 0.18531014
## bottom 0.00000000
## top 10.00000000De ahi que:
demo$efa_ok=head(scores,-2)Ahora sí:
ggplot(data=demo,aes(x=Overall_Score,y=efa_ok)) + geom_point() + theme_minimal()
6. ANALISIS FACTORIAL CONFIRMATORIO
Si la exploración apoyaba nuestro marco teórico, podemos proponer cómo construir los indices:
model <- ' democracia =~ Electoral_Processand_Pluralism + Functioningof_Government + Politicalparticipation + Politicalculture + Civilliberties'Ahora vemos qué arroja el modelo:
# normalizar las variables:
theDataNorm=scale(theData)
library(lavaan)## This is lavaan 0.6-15
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covcfa_fit <- cfa(model, data=theDataNorm,
std.lv=TRUE,
missing="fiml")Preparo los tests:
allParamCFA=parameterEstimates(cfa_fit,standardized = T)
allFitCFA=as.list(fitMeasures(cfa_fit))Veamos resultados:
- Si cada indicador tiene una buena conexión con su latente (ver p valor):
allParamCFA[allParamCFA$op=="=~",]## lhs op rhs est se z pvalue
## 1 democracia =~ Electoral_Processand_Pluralism 0.929 0.058 15.989 0
## 2 democracia =~ Functioningof_Government 0.893 0.060 14.933 0
## 3 democracia =~ Politicalparticipation 0.830 0.063 13.259 0
## 4 democracia =~ Politicalculture 0.632 0.069 9.091 0
## 5 democracia =~ Civilliberties 0.971 0.056 17.337 0
## ci.lower ci.upper std.lv std.all std.nox
## 1 0.815 1.043 0.929 0.932 0.932
## 2 0.776 1.010 0.893 0.896 0.896
## 3 0.707 0.953 0.830 0.832 0.832
## 4 0.496 0.768 0.632 0.634 0.634
## 5 0.861 1.081 0.971 0.974 0.974Averigüemos qué tan bien salió el modelo:
- El ChiSquare es NO significativo? (p_value debe ser mayor a 0.05 para que sea bueno)
allFitCFA[c("chisq", "df", "pvalue")] # pvalue>0.05## $chisq
## [1] 36.04353
##
## $df
## [1] 5
##
## $pvalue
## [1] 9.309554e-07- El Índice Tucker Lewi es mayor a 0.9? Este índice tiende a 1 para modelos con muy buen ajuste, considerándose aceptables valores superiores a 0,90.
allFitCFA$tli # > 0.90## [1] 0.9260249- La Raíz del error cuadrático medio de aproximación es menor a 0.05? Se resta la diferencia entre el valor pronosticado y el valor real en cada punto que se pronostica y se eleva al cuadrado. Se suman todos los resultados del punto anterior.. Un valor de cero indica ajuste perfecto y un valor inferior a 0.05 se considera un buen ajuste
allFitCFA[c('rmsea.ci.lower','rmsea' ,'rmsea.ci.upper')] # 0.05 en el Int de Conf?## $rmsea.ci.lower
## [1] 0.1358434
##
## $rmsea
## [1] 0.1916714
##
## $rmsea.ci.upper
## [1] 0.2527308Ya sabemos que no hay buen augurio.
scorescfa=normalize(lavPredict(cfa_fit),
method = "range",
margin=2, # by column
range = c(0, 10))De ahi que:
demo$cfa_ok=head(scorescfa,-2)Veamos ambos scores calculados
ggplot(data=demo,aes(x=cfa_ok,y=efa_ok)) + geom_point() + theme_minimal()
7. CONGLOMERADOS
Traigamos un archivo del COVID a nivel mundial:
#limpiar memoria
rm(list = ls())
###
linkCovid="https://github.com/owid/covid-19-data/raw/master/public/data/owid-covid-data.csv"
allCovid=read.csv(linkCovid)Veamos el tipo de dato date:
str(allCovid$date)## chr [1:325443] "2020-01-03" "2020-01-04" "2020-01-05" "2020-01-06" ...La fecha está como character, por lo que no se podría sacar cálculos:
summary(allCovid$date)## Length Class Mode
## 325443 character characterSi deseas que sea fecha:
allCovid$date=as.Date(allCovid$date,format="%Y-%m-%d")Verificando:
str(allCovid$date)## Date[1:325443], format: "2020-01-03" "2020-01-04" "2020-01-05" "2020-01-06" "2020-01-07" ...Y calculando:
summary(allCovid$date)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "2020-01-01" "2020-11-22" "2021-10-09" "2021-10-08" "2022-08-25" "2023-07-13"Quedémonos con los datos del 31 de julio de 2021:
covidJulio21=allCovid[allCovid$date=="2021-07-31",]Quedemonos con algunas variables, las cuáles están definidas aqui:
Nos quedaremos con:
total_cases_per_million: Total confirmed cases of COVID-19 per 1,000,000 people.
total_deaths_per_million: Total deaths attributed to COVID-19 per 1,000,000 people.
aged_65_older: Share of the population that is 65 years and older, most recent year available.
population_density: Number of people divided by land area, measured in square kilometers, most recent year available
Seleccionemos esas variables, incluyendo iso_code y location:
selec=c("iso_code","location","total_cases_per_million","total_deaths_per_million","aged_65_older","population_density")
covidJulio21=covidJulio21[,selec]El iso_code es de tres caracteres, pero hay algunos de mayor tamaño:
table(nchar(covidJulio21$iso_code))##
## 3 8
## 236 18Nos quedaremos con los que tienen sólo tres caracteres:
covidJulio21=covidJulio21[nchar(covidJulio21$iso_code)==3,]covidJulio21 debe tener varios valores perdidos, veamos qué países los tienen:
covidJulio21[!complete.cases(covidJulio21),"location"]## [1] "American Samoa" "Andorra"
## [3] "Anguilla" "Bermuda"
## [5] "Bonaire Sint Eustatius and Saba" "British Virgin Islands"
## [7] "Cayman Islands" "Cook Islands"
## [9] "Dominica" "Faeroe Islands"
## [11] "Falkland Islands" "French Guiana"
## [13] "Gibraltar" "Greenland"
## [15] "Guadeloupe" "Guernsey"
## [17] "Hong Kong" "Isle of Man"
## [19] "Jersey" "Kiribati"
## [21] "Liechtenstein" "Macao"
## [23] "Marshall Islands" "Martinique"
## [25] "Mayotte" "Micronesia (country)"
## [27] "Monaco" "Montserrat"
## [29] "Nauru" "New Caledonia"
## [31] "Niue" "North Korea"
## [33] "Northern Mariana Islands" "Palau"
## [35] "Pitcairn" "Reunion"
## [37] "Saint Barthelemy" "Saint Helena"
## [39] "Saint Kitts and Nevis" "Saint Martin (French part)"
## [41] "Saint Pierre and Miquelon" "Samoa"
## [43] "San Marino" "Sint Maarten (Dutch part)"
## [45] "Solomon Islands" "South Sudan"
## [47] "Syria" "Taiwan"
## [49] "Tokelau" "Tonga"
## [51] "Turkmenistan" "Turks and Caicos Islands"
## [53] "Tuvalu" "Vanuatu"
## [55] "Vatican" "Wallis and Futuna"Omitiremos esos territorios:
covidJulio21=covidJulio21[complete.cases(covidJulio21),]Revisemos la data:
summary(covidJulio21[,-c(1,2)])## total_cases_per_million total_deaths_per_million aged_65_older
## Min. : 20.87 Min. : 0.698 Min. : 1.144
## 1st Qu.: 3567.62 1st Qu.: 57.708 1st Qu.: 3.450
## Median : 27345.97 Median : 361.594 Median : 6.659
## Mean : 40099.38 Mean : 716.092 Mean : 8.831
## 3rd Qu.: 68574.04 3rd Qu.:1054.208 3rd Qu.:14.342
## Max. :182769.49 Max. :5762.595 Max. :27.049
## population_density
## Min. : 1.98
## 1st Qu.: 36.16
## Median : 82.46
## Mean : 207.76
## 3rd Qu.: 208.66
## Max. :7915.73La distribución de datos tiene rangos muy diferentes. Preparemos unos boxplot:
- Sin modificación:
boxplot(covidJulio21[,-c(1,2)],horizontal=T,las=2,cex.axis=0.3)
- Tipificado:
boxplot(scale(covidJulio21[,-c(1,2)]),horizontal=T,las=2,cex.axis=0.3)
- Suavizado logarítmicamente:
boxplot(log(covidJulio21[,-c(1,2)]),horizontal=T,las=2,cex.axis=0.3)
La opción de convertir en logaritmo aleja menos las variables entre sí. Usemos eso:
covidJulio21[,-c(1,2)]=log(covidJulio21[,-c(1,2)])- Los nombres envialos como nombre de fila:
row.names(covidJulio21)=covidJulio21$location- Prepare los datos para el analisis cluster, use como semilla aleatoria al numero 123.
Hay diversas maneras de calculas matrices de distancias entre casos. Cuando las variables son todas numéricas, es comun usar la distancia Euclideana). Hay otras técnicas útiles como la Mahattan. En nuestro caso, usaremos la distancia Gower útil cuando las variables (columnas) están de diversos tipos de escalas.
library(cluster)
set.seed(123)
g.dist = daisy(covidJulio21[,c(3:6)], metric="gower")Procesos de clusterización
Hay diversas estrategías de clusterización. Veremos dos de ellas:
- La técnica de Partición: Como su nombre lo indica, la estrategia de partición busca partir los casos en grupos. El algoritmo básico establece puntos que deben atraer a los casos, tal que estos se separen. Claro está, que estos puntos atractores van moviendose conforme los grupos se van formando, hasta que al final se han partido todos los casos.
Hay diversos algoritmos que buscan una implementación de estos principios básicos. El más conocido es el de K-medias, pero para ciencias sociales tiene la desventaja que requiere que todas las variables sean numéricas, no siendo muy adecuado ante categorías. La alternativa a las necesidades en ciencias sociales es la técnica de k-medoides.
- La técnica de Jerarquización: La jerarquización busca clusterizar por etapas, hasta que todas las posibilidades de clusterizacion sean visible. Este enfoque tiene dos familias de algoritmos:
++ Jerarquización Aglomerativa: En esta estrategia se parte por considerar cada caso (fila) como un cluster, para de ahi ir creando miniclusters hasta que todos los casos sean un solo cluster. El proceso va mostrando qué tanto esfuerzo toma juntar los elementos cluster tras cluster.
++ Jerarquización Divisiva: Esta estrategia comienza con todos los casos como un gran cluster; para de ahi dividir en clusters más pequeños.
Decidir cantidad de clusters:
- Obtenga sugerencia de cantidad de clusters
## para PARTICIÓN
library(factoextra)
fviz_nbclust(covidJulio21[,c(3:6)], pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)
## PARA JERARQUICO AGLOMERATIVA
fviz_nbclust(covidJulio21[,c(3:6)], hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "agnes")
## PARA JERARQUICO DIVISIVA
fviz_nbclust(covidJulio21[,c(3:6)], hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "diana")
Usemos la cantidad sugerida de clusters:
set.seed(123)
SUGERIDOS=2
res.pam=pam(g.dist,k = SUGERIDOS,cluster.only = F) # para PAM (Estrategia de partición)
covidJulio21$pam=res.pam$cluster
res.agnes <- hcut(g.dist, k = SUGERIDOS,hc_func='agnes') #PARA AGLOMERATIVA (VIA AGNES)
covidJulio21$agnes=res.agnes$cluster
res.diana <- hcut(g.dist, k = SUGERIDOS,hc_func='diana') #PARA DIVISIVA - DIANA
covidJulio21$diana=res.diana$clusterEvaluemos el resultado usando el coeficiente de silueta:
Un caso se ha clusterizado bien si tiene valor positivo
Un caso es dificilmente clusterizable si es muy cercano a cero
Un caso está mal clusterizado si es negativo.
fviz_silhouette(res.pam)## cluster size ave.sil.width
## 1 1 70 0.45
## 2 2 110 0.48
fviz_silhouette(res.agnes)## cluster size ave.sil.width
## 1 1 73 0.42
## 2 2 107 0.49
fviz_silhouette(res.diana)## cluster size ave.sil.width
## 1 1 69 0.47
## 2 2 111 0.48
- El Dendograma puede apoyar a entender evaluación:
dend1=pltree(res.agnes,cex = 0.5,hang = -1,main = "Agnes")
dend2=pltree(res.diana,cex = 0.5,hang = -1,main = "Diana")
- Encontremos los casos MAL clusterizados:
silPAM=data.frame(res.pam$silinfo$widths)
silPAM$country=row.names(silPAM)
poorPAM=silPAM[silPAM$sil_width<0,'country']
silAGNES=data.frame(res.agnes$silinfo$widths)
silAGNES$country=row.names(silAGNES)
poorAGNES=silAGNES[silAGNES$sil_width<0,'country']
silDIANA=data.frame(res.diana$silinfo$widths)
silDIANA$country=row.names(silDIANA)
poorDIANA=silDIANA[silDIANA$sil_width<0,'country']Juntemos todos en un Dataframe:
library("qpcR")
mal_covidClus=as.data.frame(qpcR:::cbind.na(poorPAM, poorAGNES,poorDIANA))
mal_covidClus## poorPAM poorAGNES poorDIANA
## 1 United Arab Emirates Japan <NA>
## 2 Qatar Iceland <NA>
## 3 Singapore Barbados <NA>