Realizar un html en R con el título de ejercicios de análisis de regresión lineal simple, en cada ejercicio, además de lo que se solicita, integrar los siguientes puntos a revisar en el proceso de análisis de regresión lineal simple: 1. Estimación del modelo de regresión 2. Summary e Interpretación del modelo estimado 3. Coeficiente de determinación 4. Predicción (Estimación) de la variable de respuesta para nuevos valores de la variable independiente. 5. Intervalos de confianza para los valores esperados 6. Intervalos de confianza para las predicciones 7. Propiedades de los residuales 8. Verificación de supuestos a. Validación grafica b. Análisis de puntos atípicos c. Validación con pruebas de hipótesis 9. Anova del modelo 10. Prueba de hipótesis para los parámetros
Con base en la revisión anual de sueldos de Advertising Age, Mark Hurd, de 49 años, presidente (Chairman) y presidente ejecutivo (CEO) de Hewlett-Packard Co., recibió un sueldo anual de $817 000, un bono de más de 5 millones y otras compensaciones que superaron los $17 millones. Su compensación total fue ligeramente mejor que el pago total promedio de unCEO, $12.4 millones. La tabla siguiente muestra la edad (Age) y el sueldo anual (Salary) en miles de dólares de Mark Hurd y otros 14 ejecutivos (Executive) con su respectivo cargo (Title), quienes dirigen empresas que cotizan en la bolsa (Advertising Age, 5 de diciembre de 2006).
library(tidyverse)
df <- tibble(
Executive = c("Charles Prince","Harold McGraw III","James Dimon","K. Rupert Murdoch",
"Kenneth D. Lewis","Kenneth I. Chenault","Louis C. Camilleri","Mark V. Hurd",
"Martin S. Sorrell","Robert L. Nardelli","Samuel J. Palmisano","David C. Novak",
"Henry R. Silverman","Robert C. Wright","Sumner Redstone"),
Age = c(56,57,50,75,58,54,51,49,61,57,55,53,65,62,82),
Salary = c(1000,1172,1000,4509,1500,1092,1663,817,1562,2164,1680,1173,3300,2500,5807)
)
df
## # A tibble: 15 × 3
## Executive Age Salary
## <chr> <dbl> <dbl>
## 1 Charles Prince 56 1000
## 2 Harold McGraw III 57 1172
## 3 James Dimon 50 1000
## 4 K. Rupert Murdoch 75 4509
## 5 Kenneth D. Lewis 58 1500
## 6 Kenneth I. Chenault 54 1092
## 7 Louis C. Camilleri 51 1663
## 8 Mark V. Hurd 49 817
## 9 Martin S. Sorrell 61 1562
## 10 Robert L. Nardelli 57 2164
## 11 Samuel J. Palmisano 55 1680
## 12 David C. Novak 53 1173
## 13 Henry R. Silverman 65 3300
## 14 Robert C. Wright 62 2500
## 15 Sumner Redstone 82 5807
modelo <- lm(Salary ~ Age, data = df)
coef(modelo)
## (Intercept) Age
## -6745.4366 149.2888
summary(modelo)
##
## Call:
## lm(formula = Salary ~ Age, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -799.18 -318.73 57.78 295.88 794.71
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6745.4 805.6 -8.374 1.35e-06 ***
## Age 149.3 13.5 11.056 5.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 462.7 on 13 degrees of freedom
## Multiple R-squared: 0.9039, Adjusted R-squared: 0.8965
## F-statistic: 122.2 on 1 and 13 DF, p-value: 5.546e-08
El modelo de regresión muestra que existe una relación muy fuerte entre la edad y el salario de los ejecutivos. Por cada año adicional de edad, el salario aumenta aproximadamente en 149,300 dólares y aunque el ajuste es excelente, todavía hay algunas diferencias notables entre los valores reales y los predichos por el modelo, con errores que llegan a ser de casi 800,000 dólares en algunos casos particulares.
summary(modelo)$r.squared
## [1] 0.903874
ggplot(df, aes(x = Age, y = Salary)) +
geom_point(size = 2) +
geom_smooth(method = "lm", se = FALSE, color = "pink") +
labs(title = "Diagrama de dispersión: Salary vs Age",
x = "Edad (años)", y = "Sueldo (miles de USD)")
## `geom_smooth()` using formula = 'y ~ x'
## Interpretación: Tu gráfico muestra una relación positiva muy fuerte
entre la edad y el salario de los ejecutivos. A medida que aumenta la
edad los salarios también aumentan de manera consistente, la línea de
tendencia que atraviesa los puntos tiene una pendiente pronunciada hacia
arriba confirmando que esta relación es positiva.
Bill <- data.frame(Age = 72)
# Predicción
predict(modelo, Bill)
## 1
## 4003.354
# Intervalo de confianza
predict(modelo, Bill, interval = "confidence", level = 0.95)
## fit lwr upr
## 1 4003.354 3544.65 4462.058
# Intervalo de predicción
predict(modelo, Bill, interval = "prediction", level = 0.95)
## fit lwr upr
## 1 4003.354 2903.62 5103.088
Con 95% de confianza su salario estaría entre 3,544,650 y 4,462,058
edades <- data.frame(Age = seq(45, 85, by = 5))
# IC para valor esperado
ic <- predict(modelo, edades, interval = "confidence")
# Intervalo de predicción
ip <- predict(modelo, edades, interval = "prediction")
cbind(edades, ic, ip)
## Age fit lwr upr fit lwr upr
## 1 45 -27.44259 -510.5418 455.6567 -27.44259 -1137.5733 1082.688
## 2 50 719.00119 350.8618 1087.1406 719.00119 -346.1429 1784.145
## 3 55 1465.44497 1182.2217 1748.6682 1465.44497 426.5895 2504.300
## 4 60 2211.88876 1952.1749 2471.6026 2211.88876 1179.1949 3244.583
## 5 65 2958.33254 2646.5085 3270.1566 2958.33254 1911.3179 4005.347
## 6 70 3704.77632 3292.9945 4116.5582 3704.77632 2623.7723 4785.780
## 7 75 4451.22010 3917.8896 4984.5506 4451.22010 3318.3270 5584.113
## 8 80 5197.66388 4532.9342 5862.3936 5197.66388 3997.3011 6398.027
## 9 85 5944.10767 5142.9605 6745.2548 5944.10767 4663.1541 7225.061
El modelo predice valores negativos para 45 años, lo cual no tiene sentido para salarios, esto sugiere que el modelo podría no ser adecuado para edades muy jóvenes en este contexto.
Mientras más edad mayor el salario esperado, pero también mayor la incertidumbre en la predicción para casos individuales.
# Calcular residuales y valores ajustados
residuales <- resid(modelo)
ajustados <- fitted(modelo)
#gráfico
plot(ajustados, residuales,
xlab = "Valores ajustados",
ylab = "Residuales",
main = "Residuales vs Valores Ajustados",
pch = 19, col = "purple",
ylim = c(-1500, 1500)) # Ajustar límites del eje Y
# línea horizontal en cero
abline(h = 0, col = "red", lwd = 2)
# Añadir línea de tendencia suavizada
lines(lowess(ajustados, residuales), col = "lightblue", lwd = 2)
# leyenda
legend("topright",
legend = c("Residuales", "Línea cero", "Tendencia"),
col = c("purple", "red", "lightblue"),
pch = c(19, NA, NA),
lwd = c(NA, 2, 2))
Se observa un patrón donde la variabilidad de los residuales no es constante a lo largo de los valores ajustados sino que presenta una dispersión desigual que se extienden desde aproximadamente -1500 esto puede indicar la existencia de valores anormales potencialmente influyentes en el modelo
\(H_{0}: los\ residuales\ provienen\ de\ una\ distribución\ normal\)
\(H_{1}: los\ residuales\ no\ provienen\ de\ una\ distribución\ normal\)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(carData)
shapiro.test(residuales)
##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.94932, p-value = 0.5139
Como EL valor p-valor es 0.5139 es MAYOR que 0.05 por lo tanto los residuales del modelo si siguen una distribución normal
$H_{0}: la varianza es constante $ $H_{1}: la varianza no es constante $
library(lmtest)
##
## Attaching package: 'lmtest'
## The following object is masked _by_ '.GlobalEnv':
##
## ip
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 0.46971, df = 1, p-value = 0.4931
Como el p-valor 0.4931 es mayor a 0.05 por lo tanto sí tienen varianza constante
\(H_{0}: los\ errores\ son\ independientes\)
\(H_{1}: los\ errores\ NO\ son\ independientes\)
library(car)
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
durbinWatsonTest(modelo)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.04812618 1.925742 0.788
## Alternative hypothesis: rho != 0
Como el p-valor 0.846 es mayor a 0.05, por lo tanto sí hay independencia entre los residuales
plot(modelo, which = 4) # Cook's
plot(modelo, which = 5) # Residuales vs leverage
La mayoría de observaciones 1-6, 8, 11, 13 tienen valores bajos de distancia de Cook menores a 0.1 estas observaciones NO son influyentes y no afectan significativamente el modelo mientras que las observaciones 7 y 15 son las más influyentes
anova(modelo)
## Analysis of Variance Table
##
## Response: Salary
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 26165094 26165094 122.24 5.546e-08 ***
## Residuals 13 2782630 214048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mi F value = 122.24 valor alto y p-value = 0.00000005546 esto quiere decir que la edad si es un predictor significativo del salario de ejecutivos
confint(modelo, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) -8485.7607 -5005.1126
## Age 120.1179 178.4597
Con 95% de confianza podemos afirmar que por cada año adicional de edad el salario aumenta entre 120,120 y 178,460 dólares anuales
En el beisbol, el éxito de un equipo suele valorarse en función del desempeño en bateo y en lanzamiento. Una medida del desempeño en el bateo es la cantidad de jonrones que anota el equipo mientras que en lanzamiento es el promedio de carreras permitidas por el equipo que lanza. En general, se cree que los equipos que anotan más jonrones y tienen un promedio menor de carreras permitidas ganan un mayor porcentaje de juegos. Los datos siguientes muestran la proporción de juegos ganados (Proportion Won), la cantidad de jonrones (HR, home runs) del equipo (Team) y el promedio de carreras permitidas (ERA, earned run average) de 16 equipos de la Liga Nacional que participaron en la temporada de las Grandes Ligas de Beisbol de 2003 (sitio web de USA Today, 7 de enero de 2004).
Además de los 10 puntos solicitados, valora las siguientes: a) Obtenga la ecuación de regresión estimada para predecir la proporción de juegos ganados en función de la cantidad de jonrones. b) Desarrolle la ecuación de regresión estimada para predecir la proporción de juegos ganados dado el promedio de carreras permitidas por los miembros del equipo que lanza. c) Obtenga la ecuación de regresión estimada para predecir la proporción de juegos ganados en función de la cantidad de jonrones y del promedio de carreras permitidas por los miembros del equipo que lanza. d) En la temporada 2003, San Diego ganó sólo 39.5% de sus juegos, el más bajo de la Liga Nacional. Para mejorar el récord del año siguiente, el equipo buscó nuevos jugadores que incrementaran la cantidad de jonrones a 180 y disminuyera el promedio de carreras permitidas por el equipo que lanza a 4.0. Use la ecuación de regresión estimada obtenida en el inciso c) para estimar el porcentaje de juegos que ganaría San Diego si tuviera 180 jonrones y su promedio de carreras permitidas fuera de 4.0.
library(tidyverse)
beisbol <- tibble(
Team = c("Arizona","Atlanta","Chicago","Cincinnati","Colorado","Florida","Houston","Los Angeles",
"Milwaukee","Montreal","New York","Philadelphia","Pittsburgh","San Diego","San Francisco","St. Louis"),
Won = c(0.519,0.623,0.543,0.426,0.457,0.562,0.537,0.525,
0.420,0.512,0.410,0.531,0.463,0.395,0.621,0.525),
HR = c(152,235,172,182,198,157,191,124,
196,144,124,166,163,128,180,196),
ERA = c(3.857,4.106,3.842,5.127,5.269,4.059,3.880,3.162,
5.058,4.027,4.517,4.072,4.664,4.904,3.734,4.642)
)
beisbol
## # A tibble: 16 × 4
## Team Won HR ERA
## <chr> <dbl> <dbl> <dbl>
## 1 Arizona 0.519 152 3.86
## 2 Atlanta 0.623 235 4.11
## 3 Chicago 0.543 172 3.84
## 4 Cincinnati 0.426 182 5.13
## 5 Colorado 0.457 198 5.27
## 6 Florida 0.562 157 4.06
## 7 Houston 0.537 191 3.88
## 8 Los Angeles 0.525 124 3.16
## 9 Milwaukee 0.42 196 5.06
## 10 Montreal 0.512 144 4.03
## 11 New York 0.41 124 4.52
## 12 Philadelphia 0.531 166 4.07
## 13 Pittsburgh 0.463 163 4.66
## 14 San Diego 0.395 128 4.90
## 15 San Francisco 0.621 180 3.73
## 16 St. Louis 0.525 196 4.64
modelo_a <- lm(Won ~ HR, data = beisbol)
summary(modelo_a)
##
## Call:
## lm(formula = Won ~ HR, data = beisbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10807 -0.05877 0.02147 0.04226 0.10714
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.354017 0.095906 3.691 0.00242 **
## HR 0.000888 0.000558 1.591 0.13386
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06666 on 14 degrees of freedom
## Multiple R-squared: 0.1532, Adjusted R-squared: 0.09268
## F-statistic: 2.532 on 1 and 14 DF, p-value: 0.1339
Este modelo no funciona bien para predecir victorias, es decir que batear muchos jonrones no garantiza que un equipo gane más juegos ya que los jonrones por sí solos explican solo el 15.3% de las victorias y no son estadísticamente significativos p-valor = 0.134
modelo_b <- lm(Won ~ ERA, data = beisbol)
summary(modelo_b)
##
## Call:
## lm(formula = Won ~ ERA, data = beisbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.076783 -0.021889 -0.006412 0.034078 0.101827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.86474 0.09661 8.951 3.62e-07 ***
## ERA -0.08367 0.02223 -3.764 0.0021 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05107 on 14 degrees of freedom
## Multiple R-squared: 0.503, Adjusted R-squared: 0.4675
## F-statistic: 14.17 on 1 and 14 DF, p-value: 0.002095
Este modelo es mucho mejor el ERA explica el 50.3% de las victorias y es muy significativo p-valor = 0.002 por cada punto que aumenta el ERA, las victorias disminuyen en 0.0837. Esto confirma que los equipos con mejores pitcheos ganan más juegos
modelo_c <- lm(Won ~ HR + ERA, data = beisbol)
summary(modelo_c)
##
## Call:
## lm(formula = Won ~ HR + ERA, data = beisbol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04478 -0.01230 0.00497 0.01187 0.04935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7091884 0.0600608 11.808 2.54e-08 ***
## HR 0.0014006 0.0002453 5.710 7.17e-05 ***
## ERA -0.1025967 0.0127556 -8.043 2.11e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0283 on 13 degrees of freedom
## Multiple R-squared: 0.8583, Adjusted R-squared: 0.8365
## F-statistic: 39.37 on 2 and 13 DF, p-value: 3.046e-06
Este es el mejor modelo explica el 85.8% de las victorias. Ambas variables son altamente significativas: Por cada jonrón adicional las victorias aumentan 0.0014 Por cada punto que aumenta el ERA las victorias disminuyen 0.1026
r2_a <- summary(modelo_a)$r.squared
r2_b <- summary(modelo_b)$r.squared
r2_c <- summary(modelo_c)$r.squared
c(R2_HR = r2_a, R2_ERA = r2_b, R2_HR_ERA = r2_c)
## R2_HR R2_ERA R2_HR_ERA
## 0.1531724 0.5029622 0.8583070
nuevo <- data.frame(HR = 180, ERA = 4.0)
# Predicción
predict(modelo_c, nuevo)
## 1
## 0.5509179
# Intervalo de confianza
predict(modelo_c, nuevo, interval = "confidence", level = 0.95)
## fit lwr upr
## 1 0.5509179 0.5318672 0.5699686
# Intervalo de predicción
predict(modelo_c, nuevo, interval = "prediction", level = 0.95)
## fit lwr upr
## 1 0.5509179 0.4868842 0.6149516
Con 180 jonrones y ERA de 4.0 se predice que ganarán el 55.09% de sus juegos. El intervalo de confianza sugiere que podrían ganar entre 53.2% y 57.0%.
confint(modelo_c, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 0.5794350215 0.838941813
## HR 0.0008706971 0.001930594
## ERA -0.1301533941 -0.075039945
Con 95% de confianza podemos afirmar que: El intercepto está entre 0.579 y 0.839 Por cada jonrón adicional las victorias aumentan entre 0.00087 y 0.00193 Por cada punto que aumenta el ERA las victorias disminuyen entre 0.130 y 0.075
##6. Diagnóstico de residuales
par(mfrow=c(2,2))
plot(modelo_c)
par(mfrow=c(1,1))
\(H_{0}: los\ residuales\ provienen\ de\ una\ distribución\ normal\)
\(H_{1}: los\ residuales\ no\ provienen\ de\ una\ distribución\ normal\)
shapiro.test(resid(modelo_c))
##
## Shapiro-Wilk normality test
##
## data: resid(modelo_c)
## W = 0.95439, p-value = 0.5623
Como el p-valor=0.5623 es mayor que 0.05 NO rechazamos \(H_{0}: los\ residuales\ provienen\ de\ una\ distribución\ normal\)
$H_{0}: la varianza es constante $ $H_{1}: la varianza no es constante $
bptest(modelo_c)
##
## studentized Breusch-Pagan test
##
## data: modelo_c
## BP = 1.4039, df = 2, p-value = 0.4956
Como el p-valor=0.4956 es mayor que 0.05 NO rechazamos $H_{0}: la varianza es constante $
\(H_{0}: los\ errores\ son\ independientes\)
\(H_{1}: los\ errores\ NO\ son\ independientes\)
durbinWatsonTest(modelo_c)
## lag Autocorrelation D-W Statistic p-value
## 1 0.1609368 1.643358 0.466
## Alternative hypothesis: rho != 0
Como el p-valor= 0.498 es mayor que 0.05 NO rechazamos \(H_{0}: los\ errores\ son\ independientes\)
influencePlot(modelo_c, id.method="identify", main="Puntos influyentes")
## Warning in plot.window(...): "id.method" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter
## Warning in box(...): "id.method" is not a graphical parameter
## Warning in title(...): "id.method" is not a graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter
## StudRes Hat CookD
## 2 0.2653483 0.42250676 0.01849352
## 6 2.0231050 0.08014246 0.09602080
## 8 -1.5833117 0.37772799 0.45454822
## 9 -1.9349405 0.18995051 0.24164139
anova(modelo_c)
## Analysis of Variance Table
##
## Response: Won
## Df Sum Sq Mean Sq F value Pr(>F)
## HR 1 0.011253 0.011253 14.053 0.002433 **
## ERA 1 0.051806 0.051806 64.695 2.108e-06 ***
## Residuals 13 0.010410 0.000801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El ANOVA confirma que tanto los jonrones como el promedio de carreras permitidas son predictores estadísticamente significativos de la proporción de juegos ganados siendo el ERA el factor más determinante para el éxito de un equipo.
La multicolinealidad en el análisis de regresión múltiple ocurre cuando dos o más variables independientes están muy correlacionadas entre sí, esto significa que no son independientes en sentido estadístico y comparten información similar
El problema principal es que aunque el modelo completo pueda ser significativo en una prueba F la alta correlación entre las variables hace que sea difícil medir el efecto individual de cada una por eso las pruebas t pueden indicar que ninguna de las variables por separado es significativa aun cuando en conjunto sí lo sean, esto pasa porque las variables aportan básicamente la misma información y el modelo no puede distinguir cuál es la contribución única de cada una.
Una regla práctica para detectar este problema es fijarse en el coeficiente de correlación entre las variables independientes, si este valor es mayor a 0.7 o menor a -0.7 entre cualquier par de variables es una señal de alerta de multicolinealidad