1 Introduccion
2 Creamos el modelo de regresión con todas las variables preditoras
3 Nuevo modelo con selección de los mejores predictores
4 Diagnósticos del modelo lineal trabajado

1 Introduccion

En este documento presentó un ejemplo ‘de libro’ de regresión lineal multiple con su correspondiente evaluación. Dejo aquí el link al documento RMardown en GitHub.

Los datos que emplearemos provienen de la librería MASS dataset Boston. Dichos datos refieren a vivienda en EEUU y contienen 506 secciones censales de Boston del censo de 1970. La base de datos Boston Housing contiene los datos originales de Harrison y Rubinfeld (1979), el marco de datos BostonHousing 2 es la versión corregida con información espacial adicional.

Esta informaci?n esta incluida en la biblioteca mlbench o descargar el conjunto de datos. Los datos tienen las siguientes características, siendo medv la variable de objetivo o independiente:

crim - Crimen per cápita por ciudad
zn - proporción de terrenos residenciales divididos en zonas para lotes de más de 25,000 pies cuadrados
indus - proporción de acres de negocios no minoristas por ciudad
chas - variable ficticia de Charles River (= 1 si el tramo limita el río, 0 de lo contrario)
nox - concentración de óxidos nítricos (partes por 10 millones)
rm - número promedio de habitaciones por vivienda
age - proporción de unidades ocupadas por sus propietarios construidas antes de 1940
dis - Distancias desproporcionadas a cinco centros de empleo de Boston
rad - índice de accesibilidad a las autopistas radiales
tax - tasa de impuesto a la propiedad de valor completo por USD 10,000
ptratio - colegios por localidad
b 1000 (B - 0,63)^ 2, donde B es la proporción de afrodescendientes por ciudad
lstat - porcentaje de estado inferior de la población
medv - valor mediano de las viviendas ocupadas por sus propietarios en USD 1000

Repetimos que el ejemplo es de libro, aunque bien cabe considerar la dimensión ética de la información contenida con un claro sesgo racista. Retomaremos esto en otro lugar.

1.1 Cargamos script de librerias y funciones de ayuda

# tools
source("main.R")

model_equation <- function(model, ...) {
  format_args <- list(...)
  
  model_coeff <- model$coefficients
  format_args$x <- abs(model$coefficients)
  model_coeff_sign <- sign(model_coeff)
  model_coeff_prefix <- case_when(model_coeff_sign == -1 ~ " - ",
                                  model_coeff_sign == 1 ~ " + ",
                                  model_coeff_sign == 0 ~ " + ")
  model_eqn <- paste(strsplit(as.character(model$call$formula), "~")[[2]], # 'y'
                     "=",
                     paste(if_else(model_coeff[1]<0, "- ", ""),
                           do.call(format, format_args)[1],
                           paste(model_coeff_prefix[-1],
                                 do.call(format, format_args)[-1],
                                 " * ",
                                 names(model_coeff[-1]),
                                 sep = "", collapse = ""),
                           sep = ""))
  return(model_eqn)
}

# Datos
df = MASS::Boston %>% as_tibble()
df = df %>% mutate_all(as.numeric) %>% as_tibble()

1.2 Composición del dataset

A continuación veremos la composición de los datos, antes de crear el modelo de regresión lineal.

skimr::skim(df)

Data summary
Name	df
Number of rows	506
Number of columns	14
_______________________
Column type frequency:
numeric	14
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
crim	1	3.61	8.60	0.01	0.08	0.26	3.68	88.98	▇▁▁▁▁
zn	1	11.36	23.32	0.00	0.00	0.00	12.50	100.00	▇▁▁▁▁
indus	1	11.14	6.86	0.46	5.19	9.69	18.10	27.74	▇▆▁▇▁
chas	1	0.07	0.25	0.00	0.00	0.00	0.00	1.00	▇▁▁▁▁
nox	1	0.55	0.12	0.38	0.45	0.54	0.62	0.87	▇▇▆▅▁
rm	1	6.28	0.70	3.56	5.89	6.21	6.62	8.78	▁▂▇▂▁
age	1	68.57	28.15	2.90	45.02	77.50	94.07	100.00	▂▂▂▃▇
dis	1	3.80	2.11	1.13	2.10	3.21	5.19	12.13	▇▅▂▁▁
rad	1	9.55	8.71	1.00	4.00	5.00	24.00	24.00	▇▂▁▁▃
tax	1	408.24	168.54	187.00	279.00	330.00	666.00	711.00	▇▇▃▁▇
ptratio	1	18.46	2.16	12.60	17.40	19.05	20.20	22.00	▁▃▅▅▇
black	1	356.67	91.29	0.32	375.38	391.44	396.22	396.90	▁▁▁▁▇
lstat	1	12.65	7.14	1.73	6.95	11.36	16.96	37.97	▇▇▅▂▁
medv	1	22.53	9.20	5.00	17.02	21.20	25.00	50.00	▂▇▅▁▁

De los datos precedentes podemos advertir distribución de los datos.

1.3 Si hubiera datos faltantes eliminaríamos dichas observaciones para generar el modelo

df = na.omit(df)

1.4 Correlación

Inspeccionamos la correlación entre la variable y otros datos estadísticos.

Esta información es muy importante para el estudio de regresión pues variables muy correlacionadas pueden indicar la necesidad de optimizar el modelo.

2 Creamos el modelo de regresión con todas las variables preditoras

Creamos el modelo de regresión lineal multiple, partiendo de considerar como predictores a todas las varibles disponibles y analizaremos sus resultados.

Observations	506
Dependent variable	medv
Type	OLS linear regression

F(13,492)	108.077
R²	0.741
Adj. R²	0.734

	Est.	S.E.	t val.	p
(Intercept)	36.459	5.103	7.144	0.000
crim	-0.108	0.033	-3.287	0.001
zn	0.046	0.014	3.382	0.001
indus	0.021	0.061	0.334	0.738
chas	2.687	0.862	3.118	0.002
nox	-17.767	3.820	-4.651	0.000
rm	3.810	0.418	9.116	0.000
age	0.001	0.013	0.052	0.958
dis	-1.476	0.199	-7.398	0.000
rad	0.306	0.066	4.613	0.000
tax	-0.012	0.004	-3.280	0.001
ptratio	-0.953	0.131	-7.283	0.000
black	0.009	0.003	3.467	0.001
lstat	-0.525	0.051	-10.347	0.000
Standard errors: OLS

2.1 Intervalos de confianza para los coeficientes β0 y β1

A continuación veremos los coefcientes y entre parentesis los intervalos de confianza respectivos. La estimación de todo coeficiente de regresión tiene asociada un error estándar, por lo tanto todo coeficiente de regresión tiene su correspondiente intervalo de confianza. Si dicho intervalo contiene 0, la variable carece de significancia para el modelo.

summ(model, confint = TRUE, digits = 3)

Observations	506
Dependent variable	medv
Type	OLS linear regression

F(13,492)	108.077
R²	0.741
Adj. R²	0.734

	Est.	2.5%	97.5%	t val.	p
(Intercept)	36.459	26.432	46.487	7.144	0.000
crim	-0.108	-0.173	-0.043	-3.287	0.001
zn	0.046	0.019	0.073	3.382	0.001
indus	0.021	-0.100	0.141	0.334	0.738
chas	2.687	0.994	4.380	3.118	0.002
nox	-17.767	-25.272	-10.262	-4.651	0.000
rm	3.810	2.989	4.631	9.116	0.000
age	0.001	-0.025	0.027	0.052	0.958
dis	-1.476	-1.867	-1.084	-7.398	0.000
rad	0.306	0.176	0.436	4.613	0.000
tax	-0.012	-0.020	-0.005	-3.280	0.001
ptratio	-0.953	-1.210	-0.696	-7.283	0.000
black	0.009	0.004	0.015	3.467	0.001
lstat	-0.525	-0.624	-0.425	-10.347	0.000
Standard errors: OLS

2.2 Calculamos R parcial

rsq::rsq.partial(model) %>% knitr::knit_print()

## $adjustment
## [1] FALSE
## 
## $variable
##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"  
## 
## $partial.rsq
##  [1] 0.021482035621 0.022714068014 0.000227109394 0.019381758450 0.042119850461
##  [6] 0.144502577602 0.000005581299 0.100105008922 0.041456694550 0.021398864103
## [11] 0.097305611557 0.023845646636 0.178718014696

2.3 Conclusiones

m_output = tidy(s, conf.int = TRUE)
var_low_pvalue = str_c(m_output$term[m_output$p.value < 0.05], collapse = ", ")
var_high_pvalue = str_c(m_output$term[m_output$p.value > 0.05], collapse = ", ")

El p-value del estadístico F es significativo para un α=0.05, luego al menos un predictor tiene asociación significativa con la variable respuesta, el modelo es útil.El p-value de las variables predictoras: indus, age no es significativo para un α=0.05, dado ello podrían eliminarse para mejorar el modelo.El p-value de las variables predictoras: (Intercept), crim, zn, chas, nox, rm, dis, rad, tax, ptratio, black, lstat son significativos para un α=0.05, luego dichos predictores tiene asociación significativa con la variable respuesta y pueden considerarse para integrar el modelo.Sobre el ‘residual standar error (RSE)’: Cualquier predicción del modelo se aleja en promedio 4.75 unidades del verdadero valor de regresión.R2: En los modelos lineales múltiples, cuantos más predictores se incluyan en el modelo mayor es el valor de R2, ya que, por poco que sea, cada predictor va a explicar una parte de la variabilidad observada en Y. Es por esto que R2 no puede utilizarse para comparar modelos con distinto número de predictores. En cambio R2_ajustado introduce una penalización al valor de R2 por cada predictor que se introduce en el modelo. El valor de la penalización depende del número de predictores utilizados y del tamaño de la muestra, es decir, del número de grados de libertad. En el caso bajo estudio el modelo empleado es capaz de explicar el 73.38% de la variabilidad observada en los datos.

Ecuación del modelo:

cat(paste0("La ecuación lineal de regresión es: ", model_equation(model, digits = 3, trim = TRUE)))

La ecuación lineal de regresión es: medv = 36.459488 - 0.108011 * crim + 0.046420 * zn + 0.020559 * indus + 2.686734 * chas - 17.766611 * nox + 3.809865 * rm + 0.000692 * age - 1.475567 * dis + 0.306049 * rad - 0.012335 * tax - 0.952747 * ptratio + 0.009312 * black - 0.524758 * lstat

3 Nuevo modelo con selección de los mejores predictores

La evaluación de un modelo de regresión múltiple así como la elección de qué predictores se deben de incluir en el modelo es uno de los pasos más importantes en la modelización estadística. A la hora de seleccionar los predictores que deben formar parte del modelo se pueden seguir varios métodos, desde emplear el criterio del analista (se introducen predictores determinados en base a conocimiento del campo) hasta métodos más sofisticados como método paso a paso (stepwise) que emplean criterios matemáticos para decidir qué predictores contribuyen significativamente al modelo y en qué orden se introducen. Dentro de este método se diferencias tres estrategias: Dirección forward (el modelo inicial no contiene ningún predictor, solo el parámetro β0, a partir de este se generan todos los posibles modelos introduciendo una sola variable de entre las disponibles por iteracion, aquella variable que mejore en mayor medida el modelo se selecciona) dirección backward (el modelo se inicia con todas las variables disponibles incluidas como predictores, y se prueba a eliminar una a una cada variable, si se mejora el modelo, queda excluida) y Doble o mixto (se trata de una combinación de la selección forward y backward. Se inicia igual que el forward pero tras cada nueva incorporación se realiza un test de extracción de predictores no útiles como en el backward).

El método paso a paso requiere de algún criterio matemático para determinar si el modelo mejora o empeora con cada incorporación o extracción. Existen varios parámetros empelados, de entre los que destacan el Cp, AIC, BIC y R2ajustado, cada uno de ellos con ventajas e inconvenientes. El método Akaike(AIC) tiende a ser más restrictivo e introducir menos predictores que el R2-ajustado.

model_subset = step(
              object    = lm(medv ~ ., data = df),
              direction = "backward",
              scope     = list(upper = ~., lower = ~1),
              trace     = T)

## Start:  AIC=1589.64
## medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + 
##     tax + ptratio + black + lstat
## 
##           Df Sum of Sq   RSS    AIC
## - age      1      0.06 11079 1587.7
## - indus    1      2.52 11081 1587.8
## <none>                 11079 1589.6
## - chas     1    218.97 11298 1597.5
## - tax      1    242.26 11321 1598.6
## - crim     1    243.22 11322 1598.6
## - zn       1    257.49 11336 1599.3
## - black    1    270.63 11349 1599.8
## - rad      1    479.15 11558 1609.1
## - nox      1    487.16 11566 1609.4
## - ptratio  1   1194.23 12273 1639.4
## - dis      1   1232.41 12311 1641.0
## - rm       1   1871.32 12950 1666.6
## - lstat    1   2410.84 13490 1687.3
## 
## Step:  AIC=1587.65
## medv ~ crim + zn + indus + chas + nox + rm + dis + rad + tax + 
##     ptratio + black + lstat
## 
##           Df Sum of Sq   RSS    AIC
## - indus    1      2.52 11081 1585.8
## <none>                 11079 1587.7
## - chas     1    219.91 11299 1595.6
## - tax      1    242.24 11321 1596.6
## - crim     1    243.20 11322 1596.6
## - zn       1    260.32 11339 1597.4
## - black    1    272.26 11351 1597.9
## - rad      1    481.09 11560 1607.2
## - nox      1    520.87 11600 1608.9
## - ptratio  1   1200.23 12279 1637.7
## - dis      1   1352.26 12431 1643.9
## - rm       1   1959.55 13038 1668.0
## - lstat    1   2718.88 13798 1696.7
## 
## Step:  AIC=1585.76
## medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq   RSS    AIC
## <none>                 11081 1585.8
## - chas     1    227.21 11309 1594.0
## - crim     1    245.37 11327 1594.8
## - zn       1    257.82 11339 1595.4
## - black    1    270.82 11352 1596.0
## - tax      1    273.62 11355 1596.1
## - rad      1    500.92 11582 1606.1
## - nox      1    541.91 11623 1607.9
## - ptratio  1   1206.45 12288 1636.0
## - dis      1   1448.94 12530 1645.9
## - rm       1   1963.66 13045 1666.3
## - lstat    1   2723.48 13805 1695.0

summary(model_subset)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171  0.00000000000272727 ***
## crim         -0.108413   0.032779  -3.307             0.001010 ** 
## zn            0.045845   0.013523   3.390             0.000754 ***
## chas          2.718716   0.854240   3.183             0.001551 ** 
## nox         -17.376023   3.535243  -4.915  0.00000120941304009 ***
## rm            3.801579   0.406316   9.356 < 0.0000000000000002 ***
## dis          -1.492711   0.185731  -8.037  0.00000000000000684 ***
## rad           0.299608   0.063402   4.726  0.00000299679930933 ***
## tax          -0.011778   0.003372  -3.493             0.000521 ***
## ptratio      -0.946525   0.129066  -7.334  0.00000000000092351 ***
## black         0.009291   0.002674   3.475             0.000557 ***
## lstat        -0.522553   0.047424 -11.019 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 0.00000000000000022

Ecuación del nuevo modelo seleccionado en base al criterio AIC:

cat(paste0("La ecuación lineal de regresión es: ", model_equation(model_subset, digits = 3, trim = TRUE)))

La ecuación lineal de regresión es: medv = 36.34115 - 0.10841 * crim + 0.04584 * zn + 2.71872 * chas - 17.37602 * nox + 3.80158 * rm - 1.49271 * dis + 0.29961 * rad - 0.01178 * tax - 0.94652 * ptratio + 0.00929 * black - 0.52255 * lstat

library(leaps)
regsubsets.out <-
  regsubsets(medv ~ .,
         data = df,
         nbest = 1,       # 1 best model for each number of predictors
         nvmax = NULL,    # NULL for no limit on number of variables
         force.in = NULL, force.out = NULL,
         method = "exhaustive")
summary(regsubsets.out)

## Subset selection object
## Call: regsubsets.formula(medv ~ ., data = df, nbest = 1, nvmax = NULL, 
##     force.in = NULL, force.out = NULL, method = "exhaustive")
## 13 Variables  (and intercept)
##         Forced in Forced out
## crim        FALSE      FALSE
## zn          FALSE      FALSE
## indus       FALSE      FALSE
## chas        FALSE      FALSE
## nox         FALSE      FALSE
## rm          FALSE      FALSE
## age         FALSE      FALSE
## dis         FALSE      FALSE
## rad         FALSE      FALSE
## tax         FALSE      FALSE
## ptratio     FALSE      FALSE
## black       FALSE      FALSE
## lstat       FALSE      FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: exhaustive
##           crim zn  indus chas nox rm  age dis rad tax ptratio black lstat
## 1  ( 1 )  " "  " " " "   " "  " " " " " " " " " " " " " "     " "   "*"  
## 2  ( 1 )  " "  " " " "   " "  " " "*" " " " " " " " " " "     " "   "*"  
## 3  ( 1 )  " "  " " " "   " "  " " "*" " " " " " " " " "*"     " "   "*"  
## 4  ( 1 )  " "  " " " "   " "  " " "*" " " "*" " " " " "*"     " "   "*"  
## 5  ( 1 )  " "  " " " "   " "  "*" "*" " " "*" " " " " "*"     " "   "*"  
## 6  ( 1 )  " "  " " " "   "*"  "*" "*" " " "*" " " " " "*"     " "   "*"  
## 7  ( 1 )  " "  " " " "   "*"  "*" "*" " " "*" " " " " "*"     "*"   "*"  
## 8  ( 1 )  " "  "*" " "   "*"  "*" "*" " " "*" " " " " "*"     "*"   "*"  
## 9  ( 1 )  "*"  " " " "   "*"  "*" "*" " " "*" "*" " " "*"     "*"   "*"  
## 10  ( 1 ) "*"  "*" " "   " "  "*" "*" " " "*" "*" "*" "*"     "*"   "*"  
## 11  ( 1 ) "*"  "*" " "   "*"  "*" "*" " " "*" "*" "*" "*"     "*"   "*"  
## 12  ( 1 ) "*"  "*" "*"   "*"  "*" "*" " " "*" "*" "*" "*"     "*"   "*"  
## 13  ( 1 ) "*"  "*" "*"   "*"  "*" "*" "*" "*" "*" "*" "*"     "*"   "*"

3.1 Estructura del mejor modelo según BIC

mejor_modelo_bic  = which(summary(regsubsets.out)$bic == min(summary(regsubsets.out)$bic))
# Para determinar la estructura del modelo 3 identificado en la salida anterior usamos:
summary(regsubsets.out)$which[mejor_modelo_bic, ]

## (Intercept)        crim          zn       indus        chas         nox 
##        TRUE        TRUE        TRUE       FALSE        TRUE        TRUE 
##          rm         age         dis         rad         tax     ptratio 
##        TRUE       FALSE        TRUE        TRUE        TRUE        TRUE 
##       black       lstat 
##        TRUE        TRUE

model_full = model

3.2 Comparacion modelo con subselección de variables con criterio: BIC

BIC(model_full,model_subset)

3.3 Comparacion modelo con subselección de variables con criterio: AIC

AIC(model_full,model_subset)

3.4 Anova para comparar la relevancia del modelo considerando los predictores que se emplean en cada caso

if(anova_result$`Pr(>F)`[2]<0.05){'H0 debe rechazarse: hay una diferencia significativa en el modelo agregando todas las variables predictoras'} else {'No hay evidencia para rechazar H0, luego el modelo no mejora considerando todas las variables predictoras'}

[1] “No hay evidencia para rechazar H0, luego el modelo no mejora considerando todas las variables predictoras”

3.5 Se puede modelar todo el espacio de posibilidad?

Considerando las posiblidad de cómputo disponibles podemos ajustar todos los modelos posibles, y seleccionar aquél con mejor performance según el criterio de referencia (i.e. Rajustado, AIC, etc.). En este caso, presentamos los 10 mejores modelos ordenados por valor de AIC.

library(olsrr)
all_models <- lm(medv ~ ., data = df)
k = ols_step_all_possible(all_models)
k %>% arrange(aic) %>% select(-predictors) %>% head(10) %>% #glimpse() %>% tidy()
  glimpse()

## Rows: 10
## Columns: 13
## $ mindex  <int> 8100, 8178, 8179, 8191, 7814, 7815, 7816, 8101, 7817, 8102
## $ n       <int> 11, 12, 12, 13, 10, 10, 10, 11, 10, 11
## $ rsquare <dbl> 0.7405823, 0.7406412, 0.7405837, 0.7406427, 0.7352631, 0.73483…
## $ adjr    <dbl> 0.7348058, 0.7343282, 0.7342694, 0.7337897, 0.7299149, 0.72948…
## $ predrsq <dbl> 0.7214716, 0.7207760, 0.7196558, 0.7189544, 0.7187463, 0.71648…
## $ cp      <dbl> 10.11455, 12.00275, 12.11176, 14.00000, 18.20493, 19.01141, 19…
## $ aic     <dbl> 3023.726, 3025.611, 3025.724, 3027.609, 3031.997, 3032.808, 30…
## $ sbic    <dbl> 1588.436, 1590.384, 1590.490, 1592.438, 1596.197, 1596.973, 15…
## $ sbc     <dbl> 3078.671, 3084.783, 3084.895, 3091.007, 3082.715, 3083.527, 30…
## $ msep    <dbl> 11351.00, 11371.49, 11374.01, 11394.59, 11560.30, 11578.86, 11…
## $ fpe     <dbl> 22.96389, 23.04966, 23.05476, 23.14088, 23.34226, 23.37974, 23…
## $ apc     <dbl> 0.2720210, 0.2730369, 0.2730974, 0.2741175, 0.2765029, 0.27694…
## $ hsp     <dbl> 0.04550083, 0.04567542, 0.04568554, 0.04586121, 0.04624618, 0.…

El mejor modelo de la tabla precedente tiene las siguientes variables prdictoras

mejor = k %>% arrange(aic) %>% slice(1)
mejor$predictors

[1] “crim zn chas nox rm dis rad tax ptratio black lstat”

model_final = model_subset

3.6 Gráfico de importancia de las variables predictoras

library(relaimpo) 
crlm <- calc.relimp(model_final)
plot(crlm)

4 Diagnósticos del modelo lineal trabajado

Una de las mejores formas de confirmar que las condiciones necesarias para un modelo de regresión lineal simple por mínimos cuadrados se cumplen es mediante el estudio de los residuos del modelo.

4.1 Relación de las variables del modelo consideradas particularmente y los residuos

En los siguientes gráficos de las variables consideradas particularmente.

car::crPlots(model_final)

4.2 Normalidad de Residuos

A continuación realizamores pruebas sobre los residuos.

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.33179, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

shapirotest = shapiro.test(model_final$residuals)
if(shapirotest$p.value<0.05){'H0 debe rechazarse: no hay normalidad en los residuos'} else {'No hay evidencia para rechazar H0, luego los residuos son normales'}

[1] “H0 debe rechazarse: no hay normalidad en los residuos”

ncv_test =  car::ncvTest(model_final)
print(ncv_test)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 15.52101, Df = 1, p = 0.000081593

if(ncv_test$p <0.05){'H0 debe rechazarse: hay evidencia de falta de homosedasticidad'} else {'No hay evidencias de falta de homocedasticidad'}

[1] “H0 debe rechazarse: hay evidencia de falta de homosedasticidad”

4.3 Independencia

#Independencia
#H0:rho=0, no hay correlaci?n entre los residuos
durbinwat_test = car::durbinWatsonTest(model_final)
durbinwat_test

##  lag Autocorrelation D-W Statistic p-value
##    1       0.4544515      1.077875       0
##  Alternative hypothesis: rho != 0

if(durbinwat_test$p<0.05){'H0 debe rechazarse: hay correlacion entre los residuos'} else {'No hay evidencia para rechazar H0, luego los residuos no presentan correlación'}

[1] “H0 debe rechazarse: hay correlacion entre los residuos”

4.4 Factor de Inflación de la Varianza

Veremos a continuación que existen valores de VIF mayores a 5, que se señalan colinealidad moderada.

vif = car::vif(model_final) %>% tidy() %>% mutate(x = round(x, digit=2), colinealidad_moderada_entre5y10 = x>5 & x<10,
                                            colienalidad_severa_mayor10 = x > 10)
vif

if(any(vif$colinealidad_moderada_entre5y10 | vif$colienalidad_severa_mayor10)){'Hay evidencia de colinealidad entre ciertas variables'} else {'No hay evidencias de colinealidad'}

[1] “Hay evidencia de colinealidad entre ciertas variables”

4.5 Datos outliers

A continuación graficaremos los datos empleando la distancia de Cook para resaltar aquellas observaciones influyentes en el modelo.

4.6 Distribución de los residuos studentizedos

library(dplyr)
df$studentized_residual <- rstudent(model_final)
ggplot(data = df, aes(x = predict(model_final), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución de los residuos studentizedos",
     x = "predicción modelo") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))

4.7 Grafico de síntesis de valores influyentes

car::influencePlot(model_final)

Los análisis muestran observaciones influyentes (StudenRes mayores abs(3))

Regresión Lineal Múltiple

Claudio Sebastian Castillo

14 de octubre de 2022