Superliga - Modelo Estadístico
Martin Alalu
9/22/2019
library(tidyverse)
library(janitor)
library(ggpubr)
futbol <- read.csv("https://raw.githubusercontent.com/martoalalu/superliga-modelo/master/data/datos_2015.csv", encoding = 'UTF-8', sep=";")
futbol <- clean_names(futbol)Análisis Exploratorio
El dataframe tiene 49 variables y 30 observaciones, correspondientes a los equipos que jugaron la Superliga de fútbol durante el año 2015. De los datos se puede observar que el campeón al cabo de 30 fechas obtuvo 64 puntos, mientras que el último 14. La valoración máxima de un plantel fue de 43.4M de euros mientras que la mínima de 3.38M. El equipo que más invirtió lo hizo por un monto de 11.2M, y el que menos invirtó en realidad vendió más de lo que compró al inicio de la temporada, por un valor total de 2.9M. En términos relativos hubo un equipo que multiplicó por 1.08 la inversión realizada para el campeonato anterior mientras que hubo un equipo cuya inversión retrocedió UN 22% para esta temporada.
Durante esta temporada 7 equipos participaron de la Copa Sudamericana, mientras que 6 para la Libertadores y 2 equipos jugaron ambas competiciones continentales. Un total de 10 equipos ascendieron de la B Nacional para jugar la temporada 2015.
futbolPara empezar vemos cuál es la correlación entre las principales variables. Sabemos que correlación no es causalidad pero nos ayudará a ir guiando la exploración.
cor_matrix <- cor(select(futbol, pts_2012_13,pts_2013_14,pts_2014, inversion_abs, inversion_relativa,valor,pos, pts,sudamericana, ascenso,libertadores), method = c("pearson", "kendall", "spearman"))
col<- colorRampPalette(c("red", "white", "blue"))(20)
heatmap(x = cor_matrix, col = col, symm = TRUE)
Cómo era esperable hay una fuerte correlación entre los puntos obtenidos sucesivamente en los últimos 3 torneos. Se destaca la correlación entre las variables valor y puntos obtenidos en el presente torneo (0.61), e inversamente entre los puntos y el hecho de haber ascendido (o sea que haber ascendido impacta negativamente en la cantidad de puntos obtenidos al final del torneo), y por último entre valor e inversión absoluta (0.57); es decir que los equipos de mayor valor de mercado y quienes más invirtieron en términos absolutos solieron tener más puntos, y por ende estar mejor posicionados. Un poco más rezagada se encuentra la correlación entre puntos e inversión relativa (es decir medida en variación respecto al año anterior), con 0.43.
modelo_pts_valor <- lm(pts ~ valor, data = futbol)
modelo_pts_invabs <- lm(pts ~ inversion_abs, data = futbol)
modelo_pts_invrel <- lm(pts ~ inversion_relativa, data = futbol)
modelo_pts_valor##
## Call:
## lm(formula = pts ~ valor, data = futbol)
##
## Coefficients:
## (Intercept) valor
## 29.6578 0.7111modelo_pts_invabs##
## Call:
## lm(formula = pts ~ inversion_abs, data = futbol)
##
## Coefficients:
## (Intercept) inversion_abs
## 36.142 1.744modelo_pts_invrel##
## Call:
## lm(formula = pts ~ inversion_relativa, data = futbol)
##
## Coefficients:
## (Intercept) inversion_relativa
## 36.83 20.35ggpubr::ggarrange(ncol=2, nrow=2,
ggplot(data = futbol) +
geom_point(aes(x = valor, y = pts)) +
labs(title = "Valor del equipo y puntos",
subtitle = "Superliga 2015",
y = "Puntos",
x= "Valor del equipo",
caption = "con línea de regresión") +
geom_abline(aes(intercept = 29.6578, slope = 0.7111), color = "blue")+
geom_smooth(aes(x = valor, y = pts), method = "lm"),
ggplot(data = futbol) +
geom_point(aes(x = inversion_abs, y = pts)) +
labs(title = "Inversión absoluta y puntos",
subtitle = "Superliga 2015",
y = "Puntos",
x= "Inversión absoluta",
caption = "con línea de regresión") +
geom_abline(aes(intercept = 36.142, slope = 1.744), color = "blue")+
geom_smooth(aes(x = inversion_abs, y = pts), method = "lm"),
ggplot(data = futbol) +
geom_point(aes(x = inversion_relativa, y = pts)) +
labs(title = "Inversión relativa y puntos",
subtitle = "Superliga 2015",
y = "Puntos",
x= "Inversión relativa",
caption = "con línea de regresión") +
geom_abline(aes(intercept = 36.83, slope = 20.35), color = "blue")+
geom_smooth(aes(x = inversion_relativa, y = pts), method = "lm")
)
summary(modelo_pts_valor)##
## Call:
## lm(formula = pts ~ valor, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.8081 -6.8428 -0.7274 6.7124 17.7294
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.6578 3.1246 9.492 3.01e-10 ***
## valor 0.7111 0.1713 4.152 0.000279 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.725 on 28 degrees of freedom
## Multiple R-squared: 0.3811, Adjusted R-squared: 0.359
## F-statistic: 17.24 on 1 and 28 DF, p-value: 0.0002791summary(modelo_pts_invabs)##
## Call:
## lm(formula = pts ~ inversion_abs, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.142 -7.132 -2.467 4.501 25.915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.1424 2.2854 15.815 1.73e-15 ***
## inversion_abs 1.7440 0.5228 3.336 0.00241 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.46 on 28 degrees of freedom
## Multiple R-squared: 0.2844, Adjusted R-squared: 0.2589
## F-statistic: 11.13 on 1 and 28 DF, p-value: 0.002407summary(modelo_pts_invrel)##
## Call:
## lm(formula = pts ~ inversion_relativa, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.8305 -8.1469 -0.8645 4.4746 23.7623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.830 2.444 15.069 5.82e-15 ***
## inversion_relativa 20.347 7.909 2.573 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.12 on 28 degrees of freedom
## Multiple R-squared: 0.1912, Adjusted R-squared: 0.1623
## F-statistic: 6.618 on 1 and 28 DF, p-value: 0.01569Estos 3 modelos son excesivamente simples, en el sentido de que sólo consideran el efecto de una sola variable independiente por separado (valor del plantel, inversión absoluta y relativa) en la dependiente (pts), sin tener en cuenta el resto (que en este caso serían parte de los errores de estimación).
El primero intenta explicar la variación de los puntos obtenidos a partir del valor del plantel, y se los resultados podemos obtenemos que el incremento de un millón de euros en el valor del plantel trae aparejado un incremento de 0.7 puntos al final del torneo. Este hallazgo es estadísticamente significativo y permite explicar el 35% de la varianza de los puntos.
El segundo modelo sólo considera la inversión absoluta como variable independiente, la cual explica también el 25% de la varianza de los puntos y es estadísticamente significativa. Ante el incremento de un millón de euros en la inversión respecto del año pasado se esperaría un incremento de 1.7 puntos al final de la temporada.
Por último el modelo que sólo utiliza la inversión relativa al año anterior. Ante un incremento de un 10% de la inversión respecto de la temporada anterior se esperaría un incremento de 2 puntos botenidos al final del torneo. Este hallazgo es estadísticamente significativo pero sólo explica un 16% de la varianza de puntos.
Analizamos ahora con un modelo más complejo, que tenga en cuenta todas las variables disponibles y así veremos si ganamos o no en poder explicativo.
#Modelo
modelo_pts <- lm(pts ~ pts_2012_13 + pts_2013_14 + pts_2014 + inversion_abs +
inversion_relativa + valor + libertadores + sudamericana + ascenso,
data = futbol)
summary(modelo_pts)##
## Call:
## lm(formula = pts ~ pts_2012_13 + pts_2013_14 + pts_2014 + inversion_abs +
## inversion_relativa + valor + libertadores + sudamericana +
## ascenso, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.1568 -3.9097 -0.2494 3.4886 15.7965
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.058086 7.459359 2.019 0.0571 .
## pts_2012_13 -0.084927 0.103988 -0.817 0.4237
## pts_2013_14 0.008748 0.106711 0.082 0.9355
## pts_2014 0.575616 0.223596 2.574 0.0181 *
## inversion_abs -2.173272 1.268791 -1.713 0.1022
## inversion_relativa 40.110542 15.874772 2.527 0.0201 *
## valor 0.834816 0.338797 2.464 0.0229 *
## libertadores 3.896544 5.391639 0.723 0.4782
## sudamericana -1.348871 3.915157 -0.345 0.7340
## ascenso 7.400458 5.703893 1.297 0.2092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.973 on 20 degrees of freedom
## Multiple R-squared: 0.7028, Adjusted R-squared: 0.5691
## F-statistic: 5.256 on 9 and 20 DF, p-value: 0.0009804En primer lugar las variables independientes seleccionadas (puntos en las 3 temporadas anteriores, inversión absoluta y relativa, el valor del equipo, el hecho de haber jugado la Copa Libertadores o la Copa Sudamericana y haber ascendido) explican en 56% la variación en los puntos. Hay 3 variables independientes cuyo estimador es estadísticamente significativo en un 95%, los puntos obtenidos en la temporada anterior (2014), la inversión relativa y el valor del equipo. En el primer caso por cada punto obtenido en el torneo 2014 conlleva un incremento de 0.57 puntos en los puntos obtenidos durante el torneo 2015; duplicar la inversión relativa implica un incremento de 40 puntos (cabe destacar que sólo un equipo hizo tal incremento, Rosario Central, cuya inversión relativa fue de 1.08), o bien incrementar un 10% la inversión implica un incremento de 4 puntos; por último por cada incremento de un millón de euros en el valor del equipo implica un incremento de 0.83 puntos. Estos efectos señalados se dan manteniendo constantes el resto de las variables independientes del modelo. En este modelo participar o no de la Copa Libertadores o Copa Sudamericana, o bien haber ascendido no tiene ningún impacto estadísticamente significativo sobre los puntos obtenidos, algo que va en contra de la intuición futbolera (los equipos que juegan copas le dan prioridad a ellas por ser torneos más prestigiosos y se ven obligados a usar equipos alternativos en el torneo local y los equipos que ascienden tienen menos recursos y menos experiencia en primera, por ende corren con desventaja frente al resto).
Diagnóstico del modelo
En primer lugar vemos cómo se distribuyen los residuos
residuos1 <- residuals(modelo_pts)
# Agrego la columna de con los residuos al data frame y los visualizamos
futbol <- futbol %>% mutate(residuo_ml = residuos1)
ggplot(futbol) +
geom_point(aes(x = pts, y = residuo_ml)) +
geom_hline(yintercept = 0, col = "blue") +
labs(x = "Pts", y = "Residuo regresión lin. aplicada")
hist(futbol$residuo_ml)
El primer gráfico muestra la distribución de los residuos a lo largo de las observaciones. A simple vista no hay ningún patrón claro, sino que se parecerían estar distribuidos aleatoriamente. Mientras que en el segundo se ve que los mismos tienen una distribución normal, algo que resulta positivo también.
par(mfrow=c(2,2)) # Change the panel layout to 2 x 2
plot(modelo_pts)
par(mfrow=c(1,1))El primer gráfico (Residuals vs Fitted) muestra si los residuos tienen patrones no lineales. En este caso es posible descartar la no-linealidad de los residuos del modelo, es decir que los mismos se pareceb distribuirse aleatoriamente, sin presentar patrón alguno. El segundo gráfico muestra que los residuos tienen una distribución normal, algo que nuevamente es positivo. El tercero muestra la distribución de los residuos a lo largo de los rangos de los predictores, de modo tal de chequear la asunción de homocedasticidad (igual varianza de los residuos en función de X). En este caso la homocedasticidad se cumple al estar los residuos distribuidos aleatoriamente, una vez más. El último gráfico permite encontrar casos que están influenciando en el modelo de sobremanera, de modo tal que alteran los resultados del mismo (empujandolo para arriba o para abajo). Nuevamente, ningún caso tiene un valor elevado en términos de la distancia de Cook, por lo que no hay casos extremos que estén influenciando en el modelo y por ende alteren los resultados.
Ahora bien, con la función step podemos ir evaluando distintos modelos sustrayendo variables y ver cuál es el que mejor poder explicativo tiene.
step(modelo_pts)## Start: AIC=132.4
## pts ~ pts_2012_13 + pts_2013_14 + pts_2014 + inversion_abs +
## inversion_relativa + valor + libertadores + sudamericana +
## ascenso
##
## Df Sum of Sq RSS AIC
## - pts_2013_14 1 0.43 1271.8 130.41
## - sudamericana 1 7.55 1279.0 130.58
## - libertadores 1 33.20 1304.6 131.17
## - pts_2012_13 1 42.40 1313.8 131.38
## <none> 1271.4 132.40
## - ascenso 1 107.01 1378.4 132.82
## - inversion_abs 1 186.51 1457.9 134.51
## - valor 1 385.98 1657.4 138.35
## - inversion_relativa 1 405.84 1677.3 138.71
## - pts_2014 1 421.30 1692.7 138.99
##
## Step: AIC=130.41
## pts ~ pts_2012_13 + pts_2014 + inversion_abs + inversion_relativa +
## valor + libertadores + sudamericana + ascenso
##
## Df Sum of Sq RSS AIC
## - sudamericana 1 8.26 1280.1 128.60
## - libertadores 1 32.91 1304.8 129.18
## - pts_2012_13 1 61.06 1332.9 129.82
## <none> 1271.8 130.41
## - ascenso 1 107.47 1379.3 130.84
## - inversion_abs 1 186.18 1458.0 132.51
## - valor 1 387.90 1659.8 136.40
## - inversion_relativa 1 412.75 1684.6 136.84
## - pts_2014 1 458.42 1730.3 137.65
##
## Step: AIC=128.6
## pts ~ pts_2012_13 + pts_2014 + inversion_abs + inversion_relativa +
## valor + libertadores + ascenso
##
## Df Sum of Sq RSS AIC
## - libertadores 1 33.46 1313.6 127.38
## - pts_2012_13 1 57.59 1337.7 127.92
## <none> 1280.1 128.60
## - ascenso 1 109.66 1389.8 129.07
## - inversion_abs 1 210.54 1490.6 131.17
## - valor 1 404.33 1684.4 134.84
## - pts_2014 1 457.72 1737.8 135.78
## - inversion_relativa 1 467.99 1748.1 135.95
##
## Step: AIC=127.38
## pts ~ pts_2012_13 + pts_2014 + inversion_abs + inversion_relativa +
## valor + ascenso
##
## Df Sum of Sq RSS AIC
## - pts_2012_13 1 49.17 1362.7 126.48
## <none> 1313.6 127.38
## - ascenso 1 186.84 1500.4 129.37
## - inversion_abs 1 248.65 1562.2 130.58
## - inversion_relativa 1 511.20 1824.8 135.24
## - pts_2014 1 533.75 1847.3 135.61
## - valor 1 806.24 2119.8 139.74
##
## Step: AIC=126.48
## pts ~ pts_2014 + inversion_abs + inversion_relativa + valor +
## ascenso
##
## Df Sum of Sq RSS AIC
## <none> 1362.7 126.48
## - ascenso 1 151.63 1514.4 127.65
## - inversion_abs 1 235.60 1598.3 129.27
## - inversion_relativa 1 533.79 1896.5 134.40
## - pts_2014 1 645.87 2008.6 136.12
## - valor 1 762.31 2125.0 137.81##
## Call:
## lm(formula = pts ~ pts_2014 + inversion_abs + inversion_relativa +
## valor + ascenso, data = futbol)
##
## Coefficients:
## (Intercept) pts_2014 inversion_abs
## 13.5176 0.4665 -2.3508
## inversion_relativa valor ascenso
## 43.9452 0.9416 7.8035modelo_step <- lm(pts ~ pts_2014 + inversion_abs + inversion_relativa +
valor + ascenso, data = futbol)
summary(modelo_step)##
## Call:
## lm(formula = pts ~ pts_2014 + inversion_abs + inversion_relativa +
## valor + ascenso, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.303 -3.851 -0.657 3.961 14.887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.5176 6.0767 2.224 0.03578 *
## pts_2014 0.4665 0.1383 3.373 0.00252 **
## inversion_abs -2.3508 1.1540 -2.037 0.05282 .
## inversion_relativa 43.9452 14.3326 3.066 0.00530 **
## valor 0.9416 0.2570 3.664 0.00123 **
## ascenso 7.8035 4.7753 1.634 0.11528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.535 on 24 degrees of freedom
## Multiple R-squared: 0.6815, Adjusted R-squared: 0.6152
## F-statistic: 10.27 on 5 and 24 DF, p-value: 2.345e-05residuos_step <- residuals(modelo_step)
# Agrego la columna de con los residuos al data frame y los visualizamos
futbol <- futbol %>% mutate(residuo_step = residuos_step)
ggplot(futbol) +
geom_point(aes(x = pts, y = residuos_step)) +
geom_hline(yintercept = 0, col = "blue") +
labs(x = "Pts", y = "Residuo regresión lin. aplicada")
hist(futbol$residuo_step)
par(mfrow=c(2,2)) # Change the panel layout to 2 x 2
plot(modelo_step)
par(mfrow=c(1,1))El modelo seleccionado es el que tiene como variables independientes a los puntos obtenidos en 2014, la inversión absoluta, la inversión relativa, el valor del equipo al inicio de la temporada y el hecho de haber ascendido o no en el último torneo. Este tiene un poder explicativo más alto que los anteriores, llegando al 61%. Por cada incremento en una unidad de los puntos obtenidos en el torneo anterior (2014) se esperaría un incremento de 0.4 puntos en el actual torneo; por cada incremento de un 10% en la inversión relativa se esperaría un incremento de 4.3 puntos; mientras que por cada incremento de un millón de euroes en el valor del equipo se esperaría un incremento de 0.9 puntos. Al igual que en elr esto de los modelos la distribución de los errores parece tener una distribución aleatoria y normal. Lo único que se observa es la observación 19 correspondiente a River Plate, el cual en términos de distancia de Cook esta ligeramente fuera del rango, por lo que podría ser considerado un outlier. Veamos cuáles son los valores esperados por el modelo para River y los finalmente obtenidos
Transformación de variables
Si bien las escalas de las variables no difieren mucho probamos ahora el mismo modelo pero con las variables en escala logaritmica para estandarizarlas.
futbol <- futbol %>%
mutate(log_pts = log(pts),
inversion_abs_log=log(inversion_abs),
inversion_relativa_log=log(inversion_relativa),
valor_log=log(valor),
ascenso_log=log(ascenso),
pts_2014_log=log(pts_2014))
futbol$inversion_abs_log## [1] 1.8453002 0.8109302 2.0215476 1.2237754 2.3823201 0.2390169
## [7] 0.8109302 NaN 1.8484548 -1.7147984 -Inf -1.2729657
## [13] 1.0885620 -0.2876821 NaN 2.1435894 1.4586150 -0.4780358
## [19] 2.4159138 NaN NaN 1.9459101 NaN -1.0498221
## [25] -0.1863296 -0.4307829 0.0000000 -0.5978370 NaN -InfAl revisar los nuevos valores vemos que al aplicar la escala logaritmica en la variable “inversión_absoluta” tenemos NaN y “- infinito” por lo que no podremos correr el nuevo modelo. Probamos usar otro método para estandarizarlas centrado en la media.
library(robustHD)
futbol <- futbol %>%
mutate(pts_std = standardize(pts),
inversion_abs_std=standardize(inversion_abs),
inversion_relativa_std=standardize(inversion_relativa),
valor_std=standardize(valor),
ascenso_std=standardize(ascenso),
pts_2014_std=standardize(pts_2014))
modelo_std <- lm(pts_std ~ pts_2014_std + inversion_abs_std + inversion_relativa_std +
valor_std + ascenso_std, data = futbol)
summary(modelo_std)##
## Call:
## lm(formula = pts_std ~ pts_2014_std + inversion_abs_std + inversion_relativa_std +
## valor_std + ascenso_std, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.01291 -0.31706 -0.05409 0.32610 1.22562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.982e-16 1.133e-01 0.000 1.00000
## pts_2014_std 5.308e-01 1.574e-01 3.373 0.00252 **
## inversion_abs_std -7.189e-01 3.529e-01 -2.037 0.05282 .
## inversion_relativa_std 9.443e-01 3.080e-01 3.066 0.00530 **
## valor_std 8.174e-01 2.231e-01 3.664 0.00123 **
## ascenso_std 3.080e-01 1.885e-01 1.634 0.11528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6204 on 24 degrees of freedom
## Multiple R-squared: 0.6815, Adjusted R-squared: 0.6152
## F-statistic: 10.27 on 5 and 24 DF, p-value: 2.345e-05Al ver los datos vemos que el poder explicativo se mantiene en el 61%, por lo que al momento nos quedamos con el modelo anteriorimente descripto. Lo que sí podemos extraer a partir de este análisis es una medida de la importancia de las variables (ya que ahora están estandarizadas), así es posible afirmar que la inversión relativa es la variable independiente de mayor peso en el modelo, seguida por el valor del equipo y luego los puntos obtenidos en el torneo anterior (2014).
Mantiendo el modelo_step vamos a agregar una nueva variable que mide la eficiencia en la inversión en el torneo pasado para capturar en qué medida una buena inversión en el torneo 2014 impacta o no positivamente en la performance del torneo 2015. Es esperable que el rendimiento de un equipo no sea inmediato sino que tenga cierto lag, esta variable intenta capturar esta maduración esperable de un equipo.
futbol <- futbol %>%
mutate(valor_2014=valor-inversion_abs)
futbol <- futbol %>%
mutate(eficiencia_2014=valor_2014 / pts_2014)
futbol$eficiencia_2014[!is.finite(futbol$eficiencia_2014)]<- 0
modelo_step_2 <- lm(pts ~ pts_2014 + inversion_abs +
inversion_relativa + valor + ascenso + eficiencia_2014,
data = futbol)
summary(modelo_step_2)##
## Call:
## lm(formula = pts ~ pts_2014 + inversion_abs + inversion_relativa +
## valor + ascenso + eficiencia_2014, data = futbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0960 -4.2788 -0.5357 3.9409 14.7023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.9840 6.5692 1.976 0.06021 .
## pts_2014 0.4944 0.1814 2.726 0.01205 *
## inversion_abs -2.4632 1.2632 -1.950 0.06347 .
## inversion_relativa 45.2411 15.5453 2.910 0.00788 **
## valor 0.9921 0.3333 2.977 0.00675 **
## ascenso 7.8963 4.8862 1.616 0.11972
## eficiencia_2014 -2.0748 8.4507 -0.246 0.80824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.687 on 23 degrees of freedom
## Multiple R-squared: 0.6823, Adjusted R-squared: 0.5995
## F-statistic: 8.234 on 6 and 23 DF, p-value: 7.926e-05El nuevo modelo tiene un mayor poder explicativo menor que los anteriores (59%), por lo que lo descartamos.
Intervalos de confianza
Construimos los intervalos de confianza para los parámetros obtenidos en el modelo_step
reg <- function(X,Y){
boot <- sample.int(nrow(futbol),nrow(futbol),replace=T)
boot_fit = lm(pts ~ pts_2014 + inversion_abs + inversion_relativa +
valor + ascenso,data=futbol[boot,])
return(coef(boot_fit))
}
boot <- t(replicate(10000,reg(X,Y)))
head(boot)## (Intercept) pts_2014 inversion_abs inversion_relativa valor
## [1,] 14.463071 0.3735520 -2.693497 45.81814 1.0306870
## [2,] 19.445102 0.6359408 -1.253064 25.29880 0.5214799
## [3,] 16.662918 0.4415764 -2.507880 43.85518 0.9027087
## [4,] 14.478317 0.3887384 -2.849695 50.75373 1.0393630
## [5,] 6.803529 0.5850233 -1.716217 38.16449 1.1193616
## [6,] 14.472173 0.3008266 -1.852119 38.27135 1.0005829
## ascenso
## [1,] 6.933269
## [2,] 6.836471
## [3,] 8.595136
## [4,] 6.044444
## [5,] 9.851370
## [6,] 9.199738par(mfrow=c(3,2))
hist(boot[,1])
hist(boot[,2])
hist(boot[,3])
hist(boot[,4])
hist(boot[,5])
hist(boot[,6])
quantile(boot[,1],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## 1.872514 9.194253 13.106004 17.506756 29.192759quantile(boot[,2],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## 0.1934308 0.3848472 0.4722892 0.5582378 0.7438982quantile(boot[,3],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## -4.352733 -2.941781 -2.139368 -1.132739 1.515263quantile(boot[,4],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## 1.201656 30.252883 41.490028 51.442278 72.210906quantile(boot[,5],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## 0.1348448 0.7746052 0.9289176 1.0722425 1.4215708quantile(boot[,6],p=c(.025,.25,.5,.75,.975))## 2.5% 25% 50% 75% 97.5%
## -1.459279 5.248869 8.221979 10.980021 16.846245Los intervalos de confianza reflejan lo hallado anteriormente. Que con un 97,5% de certeza podemos decir que los estimadores encontrados para las variables pts_2014, inversion_relativa y valor son estadísticamente significativos, o bien que apenas hay un 2,5% de probabildiad de que el parámetro obtenido no sea estadísticamente significativo.
#Si hacemos confint del modelo dan distintos valores! ¿Es porque los otros están bootsrapeados?
confint(modelo_step)## 2.5 % 97.5 %
## (Intercept) 0.9759153 26.05923445
## pts_2014 0.1810062 0.75189869
## inversion_abs -4.7326344 0.03104736
## inversion_relativa 14.3641559 73.52627402
## valor 0.4112167 1.47198017
## ascenso -2.0521436 17.65920382Test de hipótesis vía simulación
El objetivo principal es ver si el modelo elegido performa mejor que uno solo con la media, es decir si el modelo predice mejor que el ruido mismo. Vemos si los estimadores son estadísticamente significativos a partir de la técnica de bootstrapping.
#Defino los 2 modelos, uno con la VD y todas las VI y otro solo con la VD
datreg_M1 <- select(futbol, pts, pts_2014, inversion_abs, inversion_relativa, valor, ascenso)
datreg_M0 <- select(futbol, pts)
# Modelos
# Mi modelo: M1
fit_M1 <- lm(pts ~ ., data = datreg_M1)
# El modelo M0 (anidado dentro de M1)
fit_M0 <- lm(pts ~ ., data = datreg_M0)
summary(fit_M0)##
## Call:
## lm(formula = pts ~ ., data = datreg_M0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.3333 -9.8333 -0.3333 9.4167 23.6667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.333 2.218 18.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.15 on 29 degrees of freedom# Extraigo los errores
residuos_M1 <- fit_M1$residuals
residuos_M0 <- fit_M0$residuals
# Extraigo la prediccion del modelo nulo M0, la cual voy a usar para generar un data set sintético
#Esta es la predicción basandonos en la media nada mas
B0 = fit_M0$fitted.values
# Calculo el valor del test F observado
## Calculo la suma de los errores al cuadrado (desviación respecto a la media) del modelo de interés (M1, quien contiene al otro modelo) y M0 contine solo el intercept
RSS_M1 <- sum(residuos_M1 ^ 2)
RSS_M0 <- sum(residuos_M0 ^ 2)
## Calculo los grados de libertad
p_M0 <- fit_M0$rank
p_M1 <- fit_M1$rank
N = nrow(datreg_M1)
## Calculo el numerador y el denominador
numerador = (RSS_M0 - RSS_M1) / (p_M1 - p_M0)
denominador = RSS_M1 / (N - p_M1)
## Calculo F
F_obs <- numerador / denominador
F_obs## [1] 10.27088reg <- function(datreg_M0, datreg_M1, B0, residuos_M1) {
# índice para hacer el boostrap
N <- nrow(datreg_M1)
boot_i <- sample.int(N, N, replace = T)
# y sintético (recuerden que R copia lo que ingresa a la func)
residuos_2_boot <- residuos_M1[boot_i]
datreg_M0$y <- residuos_2_boot + B0
datreg_M1$y <- residuos_2_boot + B0
# M1 con los y's sintéticos
boot_fit_M0 = lm(y ~ ., data = datreg_M0)
boot_fit_M1 = lm(y ~ ., data = datreg_M1)
# Residuos
boot_residuos_M0 <- boot_fit_M0$residuals
boot_residuos_M1 <- boot_fit_M1$residuals
# Suma residuos al cuadrado
RSS_M0 <- sum(boot_residuos_M0 ^ 2)
RSS_M1 <- sum(boot_residuos_M1 ^ 2)
# F simulacion boot
p_M0 <- boot_fit_M0$rank
p_M1 <- boot_fit_M1$rank
numerador = (RSS_M0 - RSS_M1) / (p_M1 - p_M0)
denominador = RSS_M1 / (N - p_M1)
F_boot <- numerador / denominador
# Devuelve F_boot
return(F_boot)
}
boot <-
t(replicate(2000, reg(datreg_M0, datreg_M1, B0, residuos_M1)))
# Que tan probable es?
mean(F_obs < boot)## [1] 0(mean(F_obs < boot)) < .05## [1] TRUEhist(boot,xlim = c(0,35))
abline(v=F_obs,col="red",lwd=3)
A partir de esto es posible rechazar la hipótesis nula y por ende seguir afirmar que el modelo construido es estadísticamente significativo. Es decir que el modelo que incluye a pts_2014, inversión absoluta, relativa, valor del equipo y el hecho de haber ascendido predice mejor que la media, o bien que los pts obtenidos en el torneo de Superliga 2015 no se explican por el error o ruido en este caso.