Instrucciones

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

Problema 1

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).

1. Cargar y visualizar los datos

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

2. Estimación del modelo de regresión

modelo <- lm(Salary ~ Age, data = df)
coef(modelo)
## (Intercept)         Age 
##  -6745.4366    149.2888

3. Summary e interpretación

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

Interpretación:

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.

4. Coeficiente de determinación

summary(modelo)$r.squared
## [1] 0.903874

5. Diagrama de dispersión

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.

6. Predicción para Bill Gustin (72 años)

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

Interpretación:

Con 95% de confianza su salario estaría entre 3,544,650 y 4,462,058

7. Intervalos de confianza y predicción para varias edades

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

Interpretación:

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.

8. Propiedades de los residuales

# 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))

Interpretación:

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

9. Verificación de supuestos

Normalidad

\(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

Interpretación:

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

Homocedasticidad

$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

Interpretación:

Como el p-valor 0.4931 es mayor a 0.05 por lo tanto sí tienen varianza constante

Independencia

\(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

Interpretación:

Como el p-valor 0.846 es mayor a 0.05, por lo tanto sí hay independencia entre los residuales

10. Análisis de influencia

plot(modelo, which = 4)   # Cook's 

plot(modelo, which = 5)   # Residuales vs leverage

Interpretación:

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

11. ANOVA

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

12. Intervalos de confianza para parámetros

confint(modelo, level = 0.95)
##                  2.5 %     97.5 %
## (Intercept) -8485.7607 -5005.1126
## Age           120.1179   178.4597

Interpretación:

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

Problema 2

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.

1. Cargar y visualizar los datos

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

2. Estimación de los modelos de regresión

a) Modelo simple con HR

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

Interpretación:

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

b) Modelo simple con ERA

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

Interpretación:

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

c) Modelo múltiple con HR y ERA

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

Interpretación:

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

3. Coeficiente de determinación

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

4. Predicciones con el modelo múltiple d) Predicción para San Diego con HR=180 y ERA=4.0

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

Interpretación:

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%.

5. Intervalos de confianza para parámetros

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

Interpretación:

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))

7. Verificación de supuestos

Normalidad

\(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

Interpretación:

Como el p-valor=0.5623 es mayor que 0.05 NO rechazamos \(H_{0}: los\ residuales\ provienen\ de\ una\ distribución\ normal\)

Homocedasticidad

$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

Interpretación:

Como el p-valor=0.4956 es mayor que 0.05 NO rechazamos $H_{0}: la varianza es constante $

Independencia

\(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

Interpretación:

Como el p-valor= 0.498 es mayor que 0.05 NO rechazamos \(H_{0}: los\ errores\ son\ independientes\)

8. Análisis de puntos atípicos

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

9. ANOVA del modelo múltiple

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

Interpretación:

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.

Lea el capítulo 15 correspondiente al análisis de regresión múltiple, ¿en qué consiste la multicolinealidad?

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