TESTEANDO LOS SUPUETOS MCO

#Quitamos notacion cientifica
options(scipen=999) 

#Cargamos datos
welfare <- read_excel("Data_set.xlsx")

#Revisamos variable
skimr::skim(welfare)
Data summary
Name welfare
Number of rows 1074
Number of columns 14
_______________________
Column type frequency:
character 1
numeric 13
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
country_id 0 1 3 3 0 25 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
year 0 1.00 1990.98 12.41 1970.00 1980.00 1991.00 2002.00 2012.00 ▇▇▇▇▇
population 129 0.88 36.36 7.01 18.99 30.79 36.57 42.26 48.73 ▂▆▇▇▇
gini 649 0.40 49.54 6.73 28.90 44.00 50.60 55.00 68.30 ▁▅▇▇▁
sector_dualism 606 0.44 -7.19 11.61 -36.66 -14.59 -9.13 -0.84 18.51 ▁▅▇▅▃
gdp 194 0.82 7345.94 4981.75 1421.71 4046.04 6264.84 8667.02 29343.90 ▇▅▁▁▁
foreign_inv 198 0.82 2.19 4.18 -55.24 0.48 1.43 3.29 39.81 ▁▁▇▆▁
ethnic_diversity 370 0.66 0.45 0.50 0.00 0.00 0.00 1.00 1.00 ▇▁▁▁▇
regime_type 0 1.00 3.03 1.22 0.00 2.00 4.00 4.00 4.00 ▁▃▁▃▇
education_budget 325 0.70 3.80 1.65 0.40 2.60 3.80 4.80 9.10 ▃▇▇▂▁
health_budget 325 0.70 2.50 1.89 0.10 1.10 2.10 3.40 15.00 ▇▃▁▁▁
socialsec_budget 467 0.57 4.10 3.57 0.00 1.50 3.00 5.70 16.80 ▇▃▂▁▁
legislative_bal 2 1.00 0.00 0.27 -0.82 -0.15 0.00 0.11 0.92 ▁▃▇▂▁
repression 0 1.00 0.10 0.30 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁ 
stat.desc(welfare)
##          country_id              year    population          gini
## nbr.val          NA    1074.000000000   945.0000000   425.0000000
## nbr.null         NA       0.000000000     0.0000000     0.0000000
## nbr.na           NA       0.000000000   129.0000000   649.0000000
## min              NA    1970.000000000    18.9925880    28.9000000
## max              NA    2012.000000000    48.7311592    68.3000000
## range            NA      42.000000000    29.7385712    39.4000000
## sum              NA 2138314.000000000 34356.0089188 21053.9933333
## median           NA    1991.000000000    36.5708504    50.6000000
## mean             NA    1990.981378026    36.3555650    49.5388078
## SE.mean          NA       0.378561808     0.2281924     0.3262479
## CI.mean          NA       0.742805393     0.4478231     0.6412646
## var              NA     153.913911986    49.2078319    45.2360207
## std.dev          NA      12.406204576     7.0148294     6.7257729
## coef.var         NA       0.006231201     0.1929506     0.1357678
##          sector_dualism             gdp  foreign_inv ethnic_diversity
## nbr.val     468.0000000      880.000000  876.0000000     704.00000000
## nbr.null      0.0000000        0.000000    0.0000000     384.00000000
## nbr.na      606.0000000      194.000000  198.0000000     370.00000000
## min         -36.6597777     1421.710000  -55.2422031       0.00000000
## max          18.5111249    29343.900000   39.8092346       1.00000000
## range        55.1709027    27922.190000   95.0514377       1.00000000
## sum       -3365.4469246  6464423.000000 1915.2881060     320.00000000
## median       -9.1254131     6264.845000    1.4272784       0.00000000
## mean         -7.1911259     7345.935227    2.1864019       0.45454545
## SE.mean       0.5367223      167.934709    0.1410747       0.01877977
## CI.mean       1.0546898      329.599822    0.2768843       0.03687115
## var         134.8171616 24817818.479232   17.4342052       0.24828656
## std.dev      11.6110793     4981.748536    4.1754287       0.49828362
## coef.var     -1.6146400        0.678164    1.9097260       1.09622396
##            regime_type education_budget health_budget socialsec_budget
## nbr.val  1074.00000000     749.00000000  749.00000000      607.0000000
## nbr.null   22.00000000       0.00000000    0.00000000       13.0000000
## nbr.na      0.00000000     325.00000000  325.00000000      467.0000000
## min         0.00000000       0.40000000    0.10000000        0.0000000
## max         4.00000000       9.10000000   15.00000000       16.8000000
## range       4.00000000       8.70000000   14.90000000       16.8000000
## sum      3259.00000000    2846.30000000 1875.38000000     2486.6000000
## median      4.00000000       3.80000000    2.10000000        3.0000000
## mean        3.03445065       3.80013351    2.50384513        4.0965404
## SE.mean     0.03723480       0.06034266    0.06897698        0.1447749
## CI.mean     0.07306129       0.11846112    0.13541151        0.2843214
## var         1.48902640       2.72728608    3.56361033       12.7225788
## std.dev     1.22025669       1.65144969    1.88775272        3.5668724
## coef.var    0.40213430       0.43457676    0.75394149        0.8707036
##          legislative_bal     repression
## nbr.val   1072.000000000 1074.000000000
## nbr.null   274.000000000  961.000000000
## nbr.na       2.000000000    0.000000000
## min         -0.816666678    0.000000000
## max          0.916149985    1.000000000
## range        1.732816663    1.000000000
## sum         -2.198186918  110.000000000
## median       0.000000000    0.000000000
## mean        -0.002050547    0.102420857
## SE.mean      0.008248672    0.009185576
## CI.mean      0.016185391    0.018023730
## var          0.072939509    0.090618552
## std.dev      0.270073155    0.301029154
## coef.var  -131.707826749    2.939139194
# after the comma we indicate the data.frame that contains the data 
model_1 <- lm(gini ~ 1 + education_budget, data = welfare)
class(model_1) # we verify that the object class is "lm" 
## [1] "lm"

Linealidad en los parámetros

Para evaluar la linealidad, hacemos un gráfico de dispersión de los valores predichos contra los residuos de u.

El objetivo es evaluar si el promedio de los residuos tiende a situarse aleatoriamente por encima y por debajo de cero. Si los residuos presentan un patrón creciente o decreciente - o cualquier otro tipo - entonces la forma funcional de una de las variables en cuestión es no lineal. Para ello utilizamos ggplot2:

ggplot(mapping = aes(x = model_1$fitted.values, y = model_1$residuals)) +
  labs(x = "Predicted Values", y = "Residuals") +
  geom_point() +
  geom_hline(mapping = aes(yintercept = 0), color= "red")

También podemos hacer un gráfico de residuos parciales donde cada variable independiente del modelo se grafica contra los residuos. El objetivo es obtener un gráfico “parcial” para observar la relación entre la(s) variable(s) independiente(s) y la variable dependiente teniendo en cuenta (controlando) todas las demás variables del modelo. Una línea de puntos nos muestra la predicción del MCO, y otra línea (roja) nos muestra la relación “real.” Si observamos que una de nuestras variables no tiene una relación lineal podemos hacer transformaciones (¡a las variables!) para que la forma funcional se aproxime a la empírica. Hay que señalar que, además de la justificación empírica, esta transformación lineal debe siempre estar apoyada por un argumento teórico de por qué la relación entre las dos variables toma tal forma.

Una transformación común que se verá regularmente en los trabajos es la de las transformaciones logarítmicas de las variables. Éstas están presentes tanto en las variables dependientes como en las independientes. Por esta razón, le ofrecemos una tabla que le será útil. Esto le permitirá saber cómo cambia la interpretación de los resultados cuando una de las variables (o ambas) se transforma.

La tabla resume la interpretación de los coeficientes en presencia de transformaciones logarítmicas”.

Por ejemplo, decidimos transformar nuestra variable dependiente de tal manera que

model_1_log <- lm(log(gini) ~ 1 + education_budget, data = welfare)

screenreg(model_1_log)
## 
## ============================
##                   Model 1   
## ----------------------------
## (Intercept)         3.78 ***
##                    (0.02)   
## education_budget    0.03 ***
##                    (0.01)   
## ----------------------------
## R^2                 0.07    
## Adj. R^2            0.07    
## Num. obs.         356       
## ============================
## *** p < 0.001; ** p < 0.01; * p < 0.05

La interpretación sería la siguiente: si aumentamos el gasto en educacion en una unidad, esperaríamos que el Gini aumente un 3%, ceteris paribus. Para saber cuándo transformar nuestras variables, veremos con un ejemplo cómo podemos diagnosticar un problema de medición en nuestras variables.

crPlots(model_1)

La relación de nuestra variable de interés con las variables dependientes parece ser cada vez más cuadrática. Entonces, es razonable hacer una transformación cuadrática a la variable. Evaluemos esto gráficamente:

welfare_no_na <- welfare %>%
drop_na(gini, education_budget, foreign_inv, health_budget,
socialsec_budget, population, sector_dualism, ethnic_diversity,
gdp, regime_type, legislative_bal)

welfare_no_na <- welfare_no_na %>%
mutate(cseduc2 = education_budget * education_budget)

model_1_quadratic <- lm(gini ~ 1 + cseduc2 + education_budget,
data = welfare_no_na)


crPlots(model_1_quadratic)

Basándonos en un diagnóstico visual, observamos una tendencia creciente en los residuos a medida que avanza en los valores predichos. También detectamos una relación no lineal entre el gasto en educación y los niveles de desigualdad. Sospechamos que esta relación podría ser cuadrática (parábola cuadrática creciente) y, según la gráfica de residuos parciales, parece que la variable transformada está mucho más cerca de la relación lineal estimada por MCO (marcada por la línea punteada). Observa que la escala de la figura de la izquierda es de 0 a 15, mientras que la de la derecha es de 0 a 20, lo que denota una pendiente más pronunciada.

Para confirmar las observaciones visuales, a menudo se utiliza una prueba estadística para diagnosticar una especificación funcional deficiente del modelo: El test RESET de Ramsey. La idea es, precisamente, evaluar si hay un error de especificación en la ecuación de regresión. Lo que esta prueba hace es estimar de nuevo el modelo pero añadiendo los valores predichos del modelo original con algunas transformaciones no lineales de las variables. Luego, a partir de un test-F, el modelo se evalúa con las especificaciones no lineales para comprobar si tiene un mejor ajuste que el modelo original sin la transformación no lineal. La hipótesis nula establece que las nuevas variables (en este caso, csdeuc^2) no son significativas para explicar la variación de la variable dependiente; es decir, que su coeficiente es igual a cero ( beta=0 )

resettest(model_1, power = 2, type = "fitted", data = bienestar_no_na)
## 
##  RESET test
## 
## data:  model_1
## RESET = 8.5762, df1 = 1, df2 = 353, p-value = 0.003627

Según el resultado del test-F, confirmamos lo que observamos gráficamente: añadir un término cuadrático al gasto en educación mejora el ajuste de nuestra estimación. Llegamos a esta conclusión observando el valor p del test RESET de Ramsey: a un nivel de significancia estadística del 5%, se rechaza la hipótesis nula de que la incorporación del término cuadrático no mejora el ajuste del modelo.

Nota: Esta evaluación se hizo para un modelo de regresión simple (bivariante), pero el mismo procedimiento puede hacerse para modelos multivariantes.

Variación en variables independientes y colinealidad no-perfecta

Para evaluar la multicolinealidad, el primer paso es observar la matriz de correlación de las variables de nuestro modelo (tal como hicimos en la etapa de análisis estadístico descriptivo)

corr_selected <- welfare %>%
select(gini, education_budget, sector_dualism, foreign_inv, gdp,
ethnic_diversity, regime_type, health_budget, socialsec_budget,
legislative_bal, population) %>%
# calculate correlation matrix and round to 1 decimal place:
cor(use = "pairwise") %>%
round(1)
ggcorrplot(corr_selected, type = "lower", lab = T, show.legend = F)

Vemos que algunas de nuestras variables tienen fuertes correlaciones, como el gasto en seguridad social gasto_segsocial y la población (poblacion), que tiene una correlación negativa de 0,7. En cualquier caso, para detectar si la multicolinealidad es un problema, es necesario realizar un test de VIF (variance inflation factor), porque mirar los pares de correlaciones no nos ayuda a establecer si más de dos variables tienen una correlación lineal. Lo que el test VIF revela es cuánto “crecen” los errores de los coeficientes cuando el resto de las variables están presentes (cuánto aumenta la varianza del error).

model_2 <- lm(gini ~ 1 + education_budget + foreign_inv + health_budget +
socialsec_budget + population + sector_dualism + ethnic_diversity + gdp + factor(regime_type) +
legislative_bal,
data = welfare_no_na)
  
vif(model_2)
##                         GVIF Df GVIF^(1/(2*Df))
## education_budget    1.831416  1        1.353298
## foreign_inv         1.473459  1        1.213861
## health_budget       1.753385  1        1.324155
## socialsec_budget    4.842801  1        2.200637
## population          5.019055  1        2.240325
## sector_dualism      1.216700  1        1.103041
## ethnic_diversity    1.944754  1        1.394544
## gdp                 2.450633  1        1.565450
## factor(regime_type) 2.653868  3        1.176648
## legislative_bal     1.743851  1        1.320550

A continuación, hacemos una consulta sobre si la raíz cuadrada de VIF para cada variable es inferior a 2 (la raíz cuadrada porque nos interesa el error estándar y no la varianza). Como regla general, la puntuación debe ser inferior a 2. Si es superior a 2, significa que la varianza es alta. Por lo tanto, hay problemas de multicolinealidad.

sqrt(vif(model_2)) > 2
##                      GVIF    Df GVIF^(1/(2*Df))
## education_budget    FALSE FALSE           FALSE
## foreign_inv         FALSE FALSE           FALSE
## health_budget       FALSE FALSE           FALSE
## socialsec_budget     TRUE FALSE           FALSE
## population           TRUE FALSE           FALSE
## sector_dualism      FALSE FALSE           FALSE
## ethnic_diversity    FALSE FALSE           FALSE
## gdp                 FALSE FALSE           FALSE
## factor(regime_type) FALSE FALSE           FALSE
## legislative_bal     FALSE FALSE           FALSE

Según la consulta, parece que no tenemos problemas de multicolinealidad.

Homoscedasticidad

Para evaluar este supuesto, se suelen seguir dos pasos:

  1. Diagnóstico visual: Queremos observar si los residuos (la distancia entre los puntos y la línea de regresión) son constantes para diferentes valores de x. Primero, hacemos un simple gráfico de dispersión de las variables independientes que nos interesan y la variable dependiente:
ggplot(welfare_no_na, aes(education_budget, gini)) +
geom_point() +
geom_smooth(method = "lm", color = "black")
## `geom_smooth()` using formula = 'y ~ x'

Al parecer, en los niveles bajos de gasto en educación, la variabilidad de los niveles de desigualdad es significativamente mayor que en los niveles más altos de gasto en educación. Podemos hacer un mejor diagnóstico visual si utilizamos el modelo estimado (y no sólo la relación entre las dos variables) y graficamos los residuos. En primer lugar, lo hacemos para el modelo bivariante, simplemente usando plot() con el argumento which = 1 (hay algunos otros diagnósticos disponibles, que no discutiremos aquí):

plot(model_1, which = 1)

Despues, lo hacemos para el modelo mutivariado:

plot(model_2, which = 1)

Estos gráficos marcan los valores ajustados contra los residuos que se grafican.Sin embargo, observamos claramente que los residuos no son constantes para los diferentes valores de la variable gasto en educación. Estamos ante un caso de heteroscedasticidad.

También podemos evaluar cada variable del modelo y así identificar en qué variables específicas está presente la heteroscedasticidad. De nuevo, lo que esperamos es que la línea roja coincida con la línea punteada (en cero). Usaremos la función residualPlots en el paquete car. Una forma de usar una función sin subir toda la librería es con :: como puedes ver aquí:

car::residualPlots(model_2)

##                     Test stat Pr(>|Test stat|)    
## education_budget      -1.7963          0.07442 .  
## foreign_inv           -0.0911          0.92750    
## health_budget         -0.8524          0.39534    
## socialsec_budget       0.9866          0.32542    
## population            -2.2980          0.02292 *  
## sector_dualism        -0.1908          0.84893    
## ethnic_diversity      -0.2501          0.80282    
## gdp                    1.2251          0.22240    
## factor(regime_type)                               
## legislative_bal        1.4071          0.16144    
## Tukey test            -4.2589       0.00002054 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. En el segundo paso hacemos un diagnóstico estadístico. Hay diferentes maneras de evaluar la homoscedasticidad, pero el test de Breusch-Pagan es la más utilizada. La lógica de este test es la siguiente: se hace una regresión, donde la variable dependiente consiste en los residuos cuadrados para evaluar si las variables independientes del modelo tienen alguna relación con u o no. Lo que esperamos es que el efecto sea 0, porque si la varianza del error es constante, el error (residuos) no debería variar en relación con los valores de x′s. En resumen, ¡no queremos rechazar la hipótesis nula!
bptest(model_2, studentize = T)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_2
## BP = 32.369, df = 12, p-value = 0.001213

Como el valor p es menor de 0,05, la hipótesis nula es rechazada. Por lo tanto, estamos en un escenario de heteroscedasticidad.

Soluciones a la heteroscedasticidad

Una vez que identificamos que tenemos heteroscedasticidad, es necesario resolverla. La primera alternativa es corregir la forma funcional. Puede estar frente a un caso en el que la no constancia del término de error es el resultado de una relación no lineal entre las variables, un problema que ya hemos aprendido a resolver, por ejemplo, exponiendo las variables. La segunda alternativa se produce a menudo cuando la naturaleza empírica de la relación hace que el error no sea constante. Sabemos que no podemos calcular los errores estándar de los estimadores como siempre lo hacemos en MCO: como la varianza del error es no constante, es necesario modificar la forma en que calculamos los errores.

Mientras que hay varias maneras de reforzar los errores (incluso podemos hacerlo a mano), R nos permite calcularlos fácilmente con el comando coeftest del paquete lmtest. Además, el paquete sandwich, con su función vcovHC, nos permite incorporar las especificaciones de la robusta matriz de varianza-covarianza.

HC0 = es el original de White (Wooldridge 2016) HC1= Es el que usa el software de Stata HC3 =Es el más conservador, por lo tanto, es muy recomendable

model_2_robust_3 <- coeftest(model_2, vcov = vcovHC(model_2, "HC3"))
## Warning in meatHC(x, type = type, omega = omega): HC3 covariances become
## numerically unstable if hat values are close to 1 as for observations 123
model_2_robust_1 <- coeftest(model_2, vcov = vcovHC(model_2, "HC1"))
model_2_robust_0 <- coeftest(model_2, vcov = vcovHC(model_2, "HC0"))
models_robust <- list(model_2, model_2_robust_0,
model_2_robust_1, model_2_robust_3)
screenreg(models_robust,
custom.model.names = c("w/o robust SE",
"robust HC0", "robust HC1", "robust HC3"))
## 
## =======================================================================
##                       w/o robust SE  robust HC0  robust HC1  robust HC3
## -----------------------------------------------------------------------
## (Intercept)            85.94 ***      85.94 ***   85.94 ***   85.94    
##                        (8.73)         (8.77)      (9.14)               
## education_budget        1.59 ***       1.59 **     1.59 **     1.59    
##                        (0.45)         (0.50)      (0.52)               
## foreign_inv             0.24           0.24        0.24        0.24    
##                        (0.18)         (0.14)      (0.14)               
## health_budget          -0.83 **       -0.83 ***   -0.83 ***   -0.83    
##                        (0.26)         (0.22)      (0.23)               
## socialsec_budget       -0.83 ***      -0.83 **    -0.83 **    -0.83    
##                        (0.20)         (0.25)      (0.26)               
## population             -0.93 ***      -0.93 ***   -0.93 ***   -0.93    
##                        (0.17)         (0.20)      (0.21)               
## sector_dualism         -0.17 ***      -0.17 ***   -0.17 ***   -0.17    
##                        (0.03)         (0.03)      (0.03)               
## ethnic_diversity        3.68 ***       3.68 ***    3.68 ***    3.68    
##                        (1.04)         (0.92)      (0.96)               
## gdp                    -0.00 **       -0.00 *     -0.00 *     -0.00    
##                        (0.00)         (0.00)      (0.00)               
## factor(regime_type)2   -2.29          -2.29       -2.29       -2.29    
##                        (4.75)         (1.36)      (1.41)               
## factor(regime_type)3   -2.90          -2.90 *     -2.90 *     -2.90    
##                        (4.70)         (1.16)      (1.20)               
## factor(regime_type)4   -5.14          -5.14 ***   -5.14 ***   -5.14    
##                        (4.62)         (0.84)      (0.88)               
## legislative_bal       -10.40 ***     -10.40 ***  -10.40 ***  -10.40    
##                        (2.22)         (2.29)      (2.38)               
## -----------------------------------------------------------------------
## R^2                     0.59                                           
## Adj. R^2                0.56                                           
## Num. obs.             167                                              
## =======================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Todas las alternativas dan lugar a errores sólidos similares. La diferencia viene dada por las diferentes especificaciones sobre la matriz robusta de varianza-covarianza (HC).

Un caso especial de Heteroscedasticidad: varianza de error asociada a los clusters

Por ejemplo, los países de América Latina podrían relacionarse por pertenecer a regiones similares (América del Sur frente a América Central o el Caribe, regiones andinas frente a no andinas, etc.). Así pues, sus errores estándar podrían correlacionarse en función de la región a la que pertenecen. Entonces, tenemos que la varianza del error condicionada por la región no es constante.

Recordemos que nuestro objetivo es estimar el efecto del gasto en educación en el Gini de los países latinoamericanos. Observemos esta relación en la región para evaluar si, a primera vista, existen clusters:

ggplot(welfare_no_na, aes(education_budget, gini)) +
geom_point() +
facet_wrap(~country_id)

Parece que hay clusters por países. Es decir, el gasto en educación por país suele mantenerse dentro de un rango que varía ligeramente. Cuando lo observamos así no es obvio, ya que hay muchos países. Sin embargo, parece que todavía hay ciertos clusters por país (las observaciones se agrupan por país; no parecen ser independientes).

ggplot(welfare_no_na, aes(education_budget, gini, color = country_id)) + 
  geom_point() + 
  theme_bw()

Para hacer la estimación de MCO con el error de cluster, usaremos el comando lm.cluser del paquete miceadds. Este comando agrupa los errores estándar según la variable de agrupación indicada. En resumen, lo que estamos haciendo es permitir la presencia de una correlación de errores dentro de los clusters, en este caso, países (aflojando la suposición de homoscedasticidad).

Los errores estándar robustos por clusters pueden aumentar o disminuir los errores estándar. Es decir, los errores estándar por conglomerados pueden ser mayores o menores que los errores estándar convencionales. La dirección en la que cambian los errores estándar depende del signo de la correlación de los errores intragrupo. Para este paso, es necesario instalar los paquetes miceadds y multiwayvcoc. Con la opción “cluster,” indicamos qué variable agrupará los errores:

library(miceadds)
model_2_cluster <- miceadds::lm.cluster(
data = welfare_no_na,
formula = gini ~ 1 + education_budget + sector_dualism + foreign_inv +
gdp + ethnic_diversity + regime_type + health_budget +
socialsec_budget + legislative_bal,
cluster = "country_id"
)
summary(model_2_cluster)
## R^2= 0.50995 
## 
##                        Estimate   Std. Error     t value
## (Intercept)      49.17707060409 3.5693644857 13.77754242
## education_budget  1.13865155652 0.6215053531  1.83208648
## sector_dualism   -0.14825914200 0.0520718528 -2.84720312
## foreign_inv       0.60984859152 0.1100757057  5.54026511
## gdp               0.00006759012 0.0001477041  0.45760475
## ethnic_diversity  3.90789809112 1.6972660370  2.30246644
## regime_type      -1.55665524244 0.7972430496 -1.95254790
## health_budget    -0.94795222851 0.3551198348 -2.66938688
## socialsec_budget  0.00158985995 0.1807967677  0.00879363
## legislative_bal  -9.69327011503 4.3862149316 -2.20993961
##                                                             Pr(>|t|)
## (Intercept)      0.0000000000000000000000000000000000000000003479156
## education_budget 0.0669385368961974425161542967543937265872955322266
## sector_dualism   0.0044105210644507586784102137755780859151855111122
## foreign_inv      0.0000000302014089774586019700361157278717394447654
## gdp              0.6472364305467457334941627777880057692527770996094
## ethnic_diversity 0.0213088818579273905906568131740641547366976737976
## regime_type      0.0508731897971326854634988023917685495689511299133
## health_budget    0.0075989864630499892594883704077801667153835296631
## socialsec_budget 0.9929837885031662647605799065786413848400115966797
## legislative_bal  0.0271093537577103520230803468393787625245749950409

Al utilizar los clusters, el coeficiente de nuestra variable independiente gasto_educ disminuyó de 1,56 a 1,19, pero mantuvo una alta significancia estadística (valor t de >12).

La normalidad en la distribución del error

El comando qqPlot viene por defecto en R, y genera una gráfica de probabilidad normal que muestra la distribución de los datos frente a una distribución normal esperada teórica. Por lo tanto, lo que necesitamos examinar en el gráfico es que las observaciones (que son los residuos) no deben salir de las líneas punteadas (que delimitan la distribución normal):

qqPlot(model_2$residuals, col.lines = "black")

## [1] 160 118

El comando ggdensity del paquete ggpubr nos permite construir diagramas de densidad. Así, podemos trazar los residuos para evaluar visualmente si siguen una distribución aproximadamente normal.

ggdensity(model_2$residuals, main = "Density plot of the residuals")

Después de evaluar los supuestos y encontrar las soluciones (cuando sea necesario), podemos tener una mayor certeza en nuestra estimación y, como resultado, en la relación encontrada entre las variables. No obstante, una explicación completa de nuestro descubrimiento debe profundizar en el por qué y cómo se relacionan las dos variables entre sí.