EJERCICIO 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,000,000\) y otras compensaciones que superaron los \(17,000,000\). 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).

# Crear el dataframe con los datos de ejecutivos
ejecutivos <- data.frame(
  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"),
  
  Title = c("Chmn/CEO", "Chmn/Pres/CEO", "Pres/CEO", "Chmn/CEO", "Chmn/Pres/CEO",
            "Chmn/CEO", "Chmn/CEO", "Chmn/Pres/CEO", "CEO", "Chmn/Pres/CEO",
            "Chmn/Pres/CEO", "Chmn/Pres/CEO", "Chmn/CEO", "Chmn/CEO", "Exec Chmn/Founder"),
  
  Company = c("Citigroup", "McGraw-Hill Cos.", "JP Morgan Chase & Co.", "News Corp.",
              "Bank of America", "American Express Co.", "Altria Group", "Hewlett-Packard Co.",
              "WPP Group", "Home Depot", "IBM Corp.", "Yum Brands", "Cendant Corp.",
              "NBC Universal", "Viacom"),
  
  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)
)

# Estructura del dataframe
str(ejecutivos)
## 'data.frame':    15 obs. of  5 variables:
##  $ Executive: chr  "Charles Prince" "Harold McGraw III" "James Dimon" "K. Rupert Murdoch" ...
##  $ Title    : chr  "Chmn/CEO" "Chmn/Pres/CEO" "Pres/CEO" "Chmn/CEO" ...
##  $ Company  : chr  "Citigroup" "McGraw-Hill Cos." "JP Morgan Chase & Co." "News Corp." ...
##  $ Age      : num  56 57 50 75 58 54 51 49 61 57 ...
##  $ Salary   : num  1000 1172 1000 4509 1500 ...

Cambiamos la variable título a factor para categorizarlo

ejecutivos$Title <- as.factor(ejecutivos$Title)
str(ejecutivos)
## 'data.frame':    15 obs. of  5 variables:
##  $ Executive: chr  "Charles Prince" "Harold McGraw III" "James Dimon" "K. Rupert Murdoch" ...
##  $ Title    : Factor w/ 5 levels "CEO","Chmn/CEO",..: 2 3 5 2 3 2 2 3 1 3 ...
##  $ Company  : chr  "Citigroup" "McGraw-Hill Cos." "JP Morgan Chase & Co." "News Corp." ...
##  $ Age      : num  56 57 50 75 58 54 51 49 61 57 ...
##  $ Salary   : num  1000 1172 1000 4509 1500 ...
library(knitr)
## Warning: package 'knitr' was built under R version 4.4.3
kable(ejecutivos, 
      col.names = c("Ejecutivos","Cargo", "Compania", "Edad","Salario ($1,000´s)"),
      caption = "Ejercicio 1",
      align = c('l', 'c', 'c','c','c'))
Ejercicio 1
Ejecutivos Cargo Compania Edad Salario ($1,000´s)
Charles Prince Chmn/CEO Citigroup 56 1000
Harold McGraw III Chmn/Pres/CEO McGraw-Hill Cos. 57 1172
James Dimon Pres/CEO JP Morgan Chase & Co. 50 1000
K. Rupert Murdoch Chmn/CEO News Corp. 75 4509
Kenneth D. Lewis Chmn/Pres/CEO Bank of America 58 1500
Kenneth I. Chenault Chmn/CEO American Express Co. 54 1092
Louis C. Camilleri Chmn/CEO Altria Group 51 1663
Mark V. Hurd Chmn/Pres/CEO Hewlett-Packard Co. 49 817
Martin S. Sorrell CEO WPP Group 61 1562
Robert L. Nardelli Chmn/Pres/CEO Home Depot 57 2164
Samuel J. Palmisano Chmn/Pres/CEO IBM Corp. 55 1680
David C. Novak Chmn/Pres/CEO Yum Brands 53 1173
Henry R. Silverman Chmn/CEO Cendant Corp. 65 3300
Robert C. Wright Chmn/CEO NBC Universal 62 2500
Sumner Redstone Exec Chmn/Founder Viacom 82 5807

1.-Estimación del modelo de regresión

Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente la edad y Titulo o cargo, y como variable dependiente o respuesta el Salario de los ejecutivos.

Utilizando variable explicativa cualitativa y una cuantitativa

modelo_mult1 <- lm(Salary ~ Age + Title, data = ejecutivos)

modelo_mult1
## 
## Call:
## lm(formula = Salary ~ Age + Title, data = ejecutivos)
## 
## Coefficients:
##            (Intercept)                     Age           TitleChmn/CEO  
##                  -6978                     140                     852  
##     TitleChmn/Pres/CEO  TitleExec Chmn/Founder           TitlePres/CEO  
##                    719                    1305                     978

El modelo de regresión lineal estimado es: \(\hat{y}_{Salario} = -6978 + 140(Age)+852(TitleChmn/CEO)+ 719 (TitleChmn/Pres/CEO) + 1,305(TitleExec Chmn/Founder)+978(TitlePres/CEO)\)

Donde: \(\hat{\beta}_0\)=-6978, esto sería el salario esperado cuando la edad es 0 y el Title es CEO.

\(\hat{\beta}_1\)=140, es decir, por cada  año adicional de edad, el salario esperado aumenta \(140\) en miles.

\(\hat{\beta}_2\) = Si el ejecutivo es hmn/CEO, gana en promedio 852 más que un CEO, manteniendo edad constante.

\(\hat{\beta}_3\)=Si es Chmn/Pres/CEO, gana \(719\) más que un CEO, manteniendo edad constante.

\(\hat{\beta}_4\)= Si es Exec Chmn/Founder, gana \(1,305\) mil más que un CEO, manteniendo la edad constante

\(\hat{\beta}_5\) = Si es Pres/CEO, gana \(978\) más que un CEO.

2.-Summary e interpretación del modelo estimado

summary(modelo_mult1)
## 
## Call:
## lm(formula = Salary ~ Age + Title, data = ejecutivos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -714.0 -198.0    0.0  227.5  649.0 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -6977.86    1412.31  -4.941 0.000801 ***
## Age                      140.00      21.89   6.396 0.000126 ***
## TitleChmn/CEO            852.00     497.35   1.713 0.120850    
## TitleChmn/Pres/CEO       718.99     515.23   1.395 0.196343    
## TitleExec Chmn/Founder  1305.05     796.94   1.638 0.135937    
## TitlePres/CEO            977.98     694.12   1.409 0.192454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 460.3 on 9 degrees of freedom
## Multiple R-squared:  0.9341, Adjusted R-squared:  0.8975 
## F-statistic: 25.52 on 5 and 9 DF,  p-value: 4.606e-05

Modelo de regresión estimadad:

\(\hat{y}_{Salario} = -6978 + 140(Age)+852(TitleChmn/CEO)+ 719 (TitleChmn/Pres/CEO) + 1,305(TitleExec Chmn/Founder)+978(TitlePres/CEO)\)

Respecto a la variable edad, por cada año adicional a la edad, el salario aumenta \(140\), notamos que es altamente significativo (p=0.000126)

Respecto a la variable explicativa titulo (variables dummy):

La categoría de referencia es “CEO” ya que es la categoría con menos observaciones.

Un Chmn/CEO: gana \(852,000\) más que la categoría de referencia

Un Chmn/Pres/CEO: gana \(719,000\) más que la categoria de referencia

Un Exec Chmn/Founder: gana \(1,305,000\) más que la categría de referencia

Un Pres/CEO: gana \(978,000\) más que la categoria de refrencia

\(p-valor=4.606e-05\) : Modelo estadísticamente significativo en general.

3.- Coeficiente de determinación

\(R^2=0.9341\), el modelo explica el 93.41% de la variación de los salarios.

\(R^2-ajustado\): Ajusta el \(R^2\) por el numero de pedictores en el modelo.

\(R^2-ajustado=0.8975\), el 89.75% de la variación es explicada, considerando la complejidad del modelo.

Es un excelente ajuste - las variables Edad y Títulos explican más del 80% de la variabilidad en los salarios de los ejecutivos. Sin embargo considerando el tamaño muestral del 15 observaciones podría generar un riesgo de sobreajuste.

4.- Predicción de la variable de respuesta para nuevos valores de las variables independientes.

**Valores ajustados (predichos) de la variable respuesta.

predict(modelo_mult1)
##         1         2         3         4         5         6         7         8 
## 1714.0102 1720.9951 1000.0000 4373.9672 1860.9928 1434.0147 1014.0215  601.0132 
##         9        10        11        12        13        14        15 
## 1562.0000 1720.9951 1440.9996 1161.0041 2973.9898 2553.9966 5807.0000
ejecutivos$Estimados<- predict(modelo_mult1)
ejecutivos

Ejecutivo 1: Se observo salario de \(1000\), y su predicción fue de \(1719.01\), hubo sobreestimación.

Ejecutivo 2: Se observó salario de \(1172\), y su predicción fue de \(1720.99\), hubo sobreestimación.

Ejecutivo 3: Se observó salario de \(1000\), y su predicción fue de \(1000.00\), hubo ajuste perfecto.

Ejecutivo 4: Se observó salario de \(4509\), y su predicción fue de \(4373.97\), hubo subestimación.

Ejecutivo 5: Se observó salario de \(1500\), y su predicción fue de \(1860.99\), hubo sobreestimación.

Ejecutivo 6: Se observó salario de \(1092\), y su predicción fue de \(1434.01\), hubo sobreestimación.

Ejecutivo 7: Se observó salario de \(1663\), y su predicción fue de \(1014.02\), hubo subestimación.

Ejecutivo 8: Se observó salario de \(817\), y su predicción fue de \(601.01\), hubo subestimación.

Ejecutivo 9: Se observó salario de \(1562\), y su predicción fue de \(1562.00\), hubo ajuste perfecto.

Ejecutivo 10: Se observó salario de \(2164\), y su predicción fue de \(1720.99\), hubo subestimación.

Ejecutivo 11: Se observó salario de \(1680\), y su predicción fue de \(1441.00\), hubo subestimación.

Ejecutivo 12: Se observó salario de \(1173\), y su predicción fue de \(1161.00\), hubo subestimación.

Ejecutivo 13: Se observó salario de \(3300\), y su predicción fue de \(2973.99\), hubo subestimación.

Ejecutivo 14: Se observó salario de \(2500\), y su predicción fue de \(2554.00\), hubo sobreestimación.

Ejecutivo 15: Se observó salario de \(5807\), y su predicción fue de \(5807.00\), hubo ajuste perfecto.

d) Suponga que Bill Gustin, de 72 años, es el presidente y CEO de una de las principales empresas de electrónica. Prediga su sueldo anual.

Estimación del salario anual de los ejecutivos con \(72 años\) y \(Pres/CEO\)

# Datos específicos para Bill Gustin
bill_guatin <- data.frame(
  Age = 72,
  Title = "Pres/CEO"  # Asumiendo que es Presidente y CEO.
)

# Asegurar el factor
bill_guatin$Title <- factor(bill_guatin$Title)
bill_guatin

Predicción

# Predicción
salario_predicho <- predict(modelo_mult1, newdata = bill_guatin)

salario_predicho
##       1 
## 4079.95

Esto significa que Bill Gustin, de 72 años , es el presidente y CEO tiene un sueldo anual de \(4,079.95\) en miles de pesos.

5.-Intervalos esperados para los valores esperados

ic1<-predict(object=modelo_mult1,newdata=ejecutivos[,c(2,4,5)],interval="confidence",level=0.95)
ic1
##          fit        lwr      upr
## 1  1714.0102 1234.01981 2194.001
## 2  1720.9951 1282.52805 2159.462
## 3  1000.0000  -41.37610 2041.376
## 4  4373.9672 3539.57412 5208.360
## 5  1860.9928 1407.86037 2314.125
## 6  1434.0147  900.79060 1967.239
## 7  1014.0215  379.97964 1648.063
## 8   601.0132   87.03888 1114.987
## 9  1562.0000  520.62390 2603.376
## 10 1720.9951 1282.52805 2159.462
## 11 1440.9996 1015.77952 1866.220
## 12 1161.0041  726.28075 1595.728
## 13 2973.9898 2493.99946 3453.980
## 14 2553.9966 2122.41772 2985.575
## 15 5807.0000 4765.62390 6848.376
#Meter las variables que se ocupan
#Rectificar la clase
class(ic1)
## [1] "matrix" "array"
#Convertir a data-frame
ic1<-as.data.frame(ic1)
class(ic1)
## [1] "data.frame"
datejecutivos<-data.frame(ejecutivos, ic1$lwr,ic1$upr)
datejecutivos

Con el \(95%\) de confianza sabemos que el Salario promedio verdadero de todos los ejecutivos en función de la edad y el titulo cae dentro del intervalo de confianza.

Notamos que existen valores negativos, lo cual no puede ser posible, lo que puede significar que el modelo es poco preciso para ciertos perfiles o faltan datos.

6. Intervalos de confianza para las predicciones

ic1.1<-predict(object=modelo_mult1,newdata=ejecutivos[,c(2,4,5)],interval="predict",level=0.95)
ic1.1
##          fit        lwr      upr
## 1  1714.0102  567.33893 2860.681
## 2  1720.9951  591.07584 2850.914
## 3  1000.0000 -472.72821 2472.728
## 4  4373.9672 3039.54748 5708.387
## 5  1860.9928  725.30220 2996.683
## 6  1434.0147  264.06080 2603.969
## 7  1014.0215 -205.18861 2233.232
## 8   601.0132 -560.29374 1762.320
## 9  1562.0000   89.27179 3034.728
## 10 1720.9951  591.07584 2850.914
## 11 1440.9996  316.15460 2565.845
## 12 1161.0041   32.53234 2289.476
## 13 2973.9898 1827.31858 4120.661
## 14 2553.9966 1426.73243 3681.261
## 15 5807.0000 4334.27179 7279.728
#Rectificar la clase
class(ic1.1)
## [1] "matrix" "array"
#Convertir a data-frame
ic1.1<-as.data.frame(ic1.1)
class(ic1.1)
## [1] "data.frame"
datejecutivos1<-data.frame(ejecutivos, ic1.1$lwr,ic1.1$upr)
datejecutivos1

El intervalo de predicción para el ejecutivo pres/CEO, de 56 años con salario de 1000 anual, su intervalo de predicción, está entre 567.33 en miles sueldo minimo anual y 2860.68 como sueldo anual maximo.

7.- Propiedades de los residuales

1)Suma de los residuales muy cercano al 0 ó es 0.

suma_residuales1<- sum(modelo_mult1$residuals)
suma_residuales1
## [1] -2.842171e-14

Es muy cercano al 0 la suma de los residuales, por lo tanto cumple.

suma_cuadrado_residuales1<-sum(suma_residuales1^2)
suma_cuadrado_residuales1
## [1] 8.077936e-28

Es muy cercano al 0 la suma de los residuales, por lo tanto cumple.

2) El promedio de los residuales es muy cercano al 0.

promedio_residuales1<-mean(modelo_mult1$residuals)
promedio_residuales1
## [1] -1.901226e-15

Es muy cercano al 0, entonces, cumple.

3) La suma de los valores observados coincide con la suma de los valores estimados

suma_observados1<-sum(ejecutivos$Salary)
suma_observados1
## [1] 30939
suma_estimados1<-sum(modelo_mult1$fitted.values)
suma_estimados1
## [1] 30939

Notemos que la suma de los valores observados y de los estimados es igual.

4) La suma de los residuales ponderados por los valores de las variables explicativas es 0

suma_ponderada_x1<-sum(ejecutivos$Age*modelo_mult1$residuals)
suma_ponderada_x1
## [1] -2.532374e-11

Es cercano al 0.

5) El promedio de los valores observados es igual al promedio de los valores estimados.

Valores observados

mean_observados1<-mean(ejecutivos$Salary)
mean_observados1
## [1] 2062.6

Valores estimados

mean_estimados1<-mean(modelo_mult1$fitted.values)
mean_estimados1
## [1] 2062.6

Notemos que el promedio de los valores observados y de los estimados es igual, entonces, cumple.

6) La recta ajustada pasa por el punto \((\bar{X},\bar{Y})\)

mean_xy1<-c(mean(ejecutivos$Age),mean(ejecutivos$Salary))
mean_xy1
## [1]   59.0 2062.6

Gráfico

plot(ejecutivos$Age,ejecutivos$Salary,col=4
     ,pch=18,xlab="Edad de los ejecutivos",ylab = "Salario Anual",main = "Modelo de regresión ajustado a las variables Edad y Salario de los ejecutivos")
abline(modelo_mult1)
## Warning in abline(modelo_mult1): only using the first two of 6 regression
## coefficients
points(59.0,2062.6)

Notemos que la propiedad no cumple, es decir la recta ajustada no pasa por el punto \((\bar{X},\bar{Y})\)

7) La distribución de los residuales sea Normal

hist(modelo_mult1$residuals,col = 4, main = "Histograma de los residuales", xlab="Residuales")

Notemos que no sigue una distribucion, pico no centrado en 0, no es simetrica.

Probemos con otra grafica

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
ejecutivos$residuales<-modelo_mult1$residuals
ejecutivos$Estimados<- modelo_mult1$fitted.values
head(ejecutivos)
qqnorm(ejecutivos$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= 4, pch=19)
qqline(ejecutivos$residuales)

Vemos que los residuales NO siguen una distribución normal, los puntos se desvían por encima de la linea cola derecha, es decir, hay valores positivos mas grandes de los esperado.

Los puntos se desvian por debajo de la linea.Esto indica que los residuales negativos son menos extremos de lo que se esperarian.

Probamos con otra gráfica

ggplot(data = data.frame(Residuales = modelo_mult1$residuals), 
       aes(x = Residuales)) +
  geom_density(fill = 4, alpha = 0.6, ) +
  labs(title = "Densidad de los Residuales",
       subtitle = "Distribución de los errores del modelo lineal",
       x = "Residuales",
       y = "Densidad") +
  theme_minimal()

Vemos que hay sesgo a la izquierda hacia los valores bajos, los datos estan concentrados en el centro.

Gráfico de probabilidad normal de los residuales ajustados

library(car)
qqPlot(modelo_mult1$residuals, 
       pch = 20, 
       col = 4, 
       main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado")

## [1] 1 7

Los residuales se comportan correctamente a excepción de los punto 1 y 7.

La varianza estimada de los residuales, es el cuadrado medio del error, que obtenemos a continuación:

El MSE mide el promedio de los cuadrados de las diferencias entre los valores reales observados y los valores predichos por un modelo (residuos), en terminos simples me permite responder a la pregunta ¿Que tan lejos en promedio estan las predicciones de mi modelo de la realidad?

mse1<-(sum(modelo_mult1$residuals)^2)/15
mse1
## [1] 5.38529e-29
des_estandar1<-sqrt(mse1)
des_estandar1
## [1] 7.338454e-15

Penaliza los errores grandes, pues las diferencias grandes entre la prediccion y el valor real se vuelven mucho mas significativas. Aunque sus unidades estàn al cuadrado un MSE de cercano a 0 significa un ajuste perfecto.

La desviación estándar estimada en el error es de 7.338454e-15

8)Verificación de supuestos: Independencia, Varianza Constante y Normalidad de los Residuales

a) Validación Gráfica

plot(modelo_mult1)
## Warning: not plotting observations with leverage one:
##   3, 9, 15

La caracteristica a observar en el gráfico de residuales contra ajustados es una dispersión aleatoria de los residuales o de los punto.No debe haber patrones. Al parecer no hay patrones

La segunda gráfica es un qqplot o bien un gráfico de probabilidad normal y permite valorar la normalidad de los residuales.

La tercera gráfica de residuales estandarizados permite valorar la varianza constante, entonces debemos ver puntos dispersos sin patrón alguno. No tenemos patron alguno, pero hay puntos muy lejos.

La distancia de Cook identifica observaciones que, si se eliminaran, cambiarían significativamente los coeficientes de la regresión, ayudando a diagnosticar problemas en el ajuste del modelo. Notemos que el punto 7, es un punto influyente ya que esta arriba del 1.

Pruebas de normalidad

Para probar si los residuales proviene de una población normal se usa test de shapiro, que contrasta la hipótesis:

\(H_{0}: los\ datos\ provienen\ de\ una\ distribución\ normal\)

\(H_{1}: los\ datos\ no\ provienen\ de\ una\ distribución\ normal\)

shapiro.test(modelo_mult1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_mult1$residuals
## W = 0.96886, p-value = 0.8408

De acuerdo con el \(p-valor=0.8404\) se puede decir con un 95% de confianza que los residuales provienen de una poblacion con distribucion normal.

Análisis de dispersión con utilizando la edad del ejercicio como la variable independiente.

ggplot(ejecutivos,aes(x=Age,y=Salary,color= Title))+
  geom_point(sixe=4)+theme_light()+
  ggtitle("grafico de dispersion")
## Warning in geom_point(sixe = 4): Ignoring unknown parameters: `sixe`

Notemos que hay relación positiva, aunque no es contante la relación.

Independencia

Para probar Independencia un supuesto muy importante, es importante saber que una muestra de manera indempendiente ie. que no hay patrones de dependencia espacial, temporal o multinivel entre las observaciones. En este supuesto espero que ¿espero que suceda cerca del acceso a la carretera sea similar a los que estan lejos?

En este caso existen pruebas estadisticas para evalual el supuesto de independencia.El test de Durbin-Watson

contrasta la hipótesis:

\(H_{0}: los\ datos\ no\ presentan\ autocorrelacion\)

\(H_{1}: los\ datos\ presentan\ autocorrelacion\)

library(car)

durbinWatsonTest(modelo_mult1)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2516531      1.229395   0.096
##  Alternative hypothesis: rho != 0

Con un p-valor > 0.05 no se puede rechazar la hipotesis nula, por lo que se presenta autocorrelacion.

Homogeneidad de Varianzas

Grafica de homogeneidad de varianzas

plot(ejecutivos$Age, ejecutivos$Salary, main = "Diagrama de Dispersión\n(Homogeneidad Varianzas)", 
     xlab = "Edad", ylab = "Salario")
abline(modelo_mult1)
## Warning in abline(modelo_mult1): only using the first two of 6 regression
## coefficients
abline(v = median(ejecutivos$Age), col = 4, lty = 3)

var.test(residuals(modelo_mult1)[ejecutivos$Age > median(ejecutivos$Age)], 
         residuals(modelo_mult1)[ejecutivos$Age < median(ejecutivos$Age)])
## 
##  F test to compare two variances
## 
## data:  residuals(modelo_mult1)[ejecutivos$Age > median(ejecutivos$Age)] and residuals(modelo_mult1)[ejecutivos$Age < median(ejecutivos$Age)]
## F = 0.26827, num df = 5, denom df = 6, p-value = 0.1703
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.04480424 1.87189619
## sample estimates:
## ratio of variances 
##          0.2682683

$H_{0}: =1 $

$H_{1}: $

En este caso dado que muchos de los valores no se acercan a la media es posible que no genere ese error

Tomando en cuienta que \(p-valor>0.05\) por lo que son iguales las varianzas por lo que se cuenta con varianza contante.

b) Analisis de puntos atipicos

Prueba de Bonferroni para datos atipicos

Con ayuda de la funcion outlierTest() de la libreria car nos ayuda a idendificar valores atípicos de los residuos studentizados del modelo

Contraste de hipotesis:

\(H_{0}: la\ observación\ i\ no\ es\ un\ valor\ atípico\)

\(H_{1}: la\ observación\ i\ es\ un\ valor\ atípico\)

library(car)
outlierTest(modelo_mult1)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 7 2.079628           0.071156      0.85388

con un \(p-valor > 0.05\) no se rechazar la hipotesis es decir contamos con la observacion de 7 valores atipicos que son significativos.

Desde una vista mas grafica podemos explorar dos graficos

ejecutivos <- na.omit(ejecutivos)
cook <- cooks.distance(modelo_mult1)
observaciones_influyentes <- which(cook > 4/nrow(ejecutivos))
observaciones_influyentes#Valores influyentes en la predicción
## 7 
## 7

Notamos que le punto influyente es el 7.

influencePlot(modelo_mult1)

Notamos que el valor 1 y 7 son los puntos mas influyentes, pero el 7 es mas influyente.

c)Validación con pruebas de hipotesis

Prueba de Homocedasticidad

\(H_{0}: La\ varianza \ es \ constante\ en \ los\ residuales\)

\(H_{1}: La\ varianza \ no\ es \ constante\ en \ los \ residuales\)

Prueba de homocedasticidad:

\(H_{0}: Hay\ homocedasticidad \ de \ los \ residuales\)

\(H_{1}: No\ hay \ homocedasticidad \ de \ los \ residuales\)

ncvTest(modelo_mult1) #prueba de homocedosticidad

Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.8312437, Df = 1, p = 0.36191 Como p-valor es mayor a \(0.05\), entonces no se rechaza la hipotesis Nula, es decir la varianza es constante, es decir hay homocedasticidad de los residuales

Test de Breusch-Pagan (homocedasticidad  los residuos)

\(H_{0}: La\ varianza \ es \ constante\ ( Hay Homocedasticidad de los residuales)\)

\(H_{1}: La\ varianza \ no\ es \ constante\ ( Hay Heterodesaticidad de los residuales)\)

library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Warning: package 'zoo' was built under R version 4.4.3
#studentized Breusch-Pagan test modelo
bptest(modelo_mult1)
studentized Breusch-Pagan test

data: modelo_mult1 BP = 5.4604, df = 5, p-value = 0.3623 Como el \(p-valor\) > \(0.05\) entonces no se rechaza la hipotesis nula, es decir, la varianza es constante ( Hay homocedasticidad de los residuales), es decir se cumple el supuesto de la homocedasticidad.

Prueba de correlación

Comprobar la independencia para los residuos estudentizados del modelo ajustado. Se realiza a través del Test de Durbin-Watson (asume bajo la hipótesis nula que no existe correlación)

\(H_{0}: No\ existe \ correlación \ entre \ los \ residuales\)

\(H_{1}: Existe\ correlación\ entre \ los \ residuales\)

dwtest(Salary~Age, alternative = "two.sided", data=ejecutivos)
Durbin-Watson test

data: Salary ~ Age DW = 1.9257, p-value = 0.8077 alternative hypothesis: true autocorrelation is not 0

Como \(p-valor=0.8077\), entonces no se tiene evidencia para rechazar hipotesis nula, es decir, no existe correlación, entonces se cumple el supuesto de independencia.

9.- Anova del Modelo

El análisis de la varianza, para las pruebas de hipotesis

Contraste de hipotesis:

\(H_{0}: \beta_{1}=0\)

\(H_{1}: \beta_{1} \neq0\)

\(H_{0}: \beta_{2}=0\)

\(H_{1}: \beta_{2} \neq0\)

anova(modelo_mult1)

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{1}\neq0\), es decir, la edad y el salario están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor >0.05\) por lo que no se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{2}\neq0\), es decir, el Title y el salario no están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

10.- Pruebas de hipotesis para los parametros

El anova también permite realizar la prueba de hipotesis para la significancia de los coeficientes de las variables explicativas, en este caso, al tener solo una variable independiente, sólo hay un contraste:

\(H_{0}: \beta_{1}=0\)

\(H_{1}: \beta_{1} \neq0\)

anova(modelo_mult1)

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{1}\neq0\), es decir, el rating y el precio están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor >0.05\) por lo que no se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{2}\neq0\), es decir, el Title y el salario no están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

La gráfica con los intervalos de confianza

library(ggplot2)
ggplot(data = ejecutivos, mapping = aes(x = Age, y = Salary)) +
  geom_point(color = "firebrick", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  labs(title = "Salary ~ Age", x = "Age", y = "Salary") +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

Hay una tendencia , ya que a medida que aumenta la edad el salario tambien. Es decir, hay casos en los que tienen menos edad y tienen maás sueldo a comparación de otros. Tambien vemos que la mayoria de los puntos estan fuera de la banda de confianza.

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

# Crear el dataframe
equipos <- data.frame(
  Team = c("Arizona", "Atlanta", "Chicago", "Cincinnati", "Colorado", 
           "Florida", "Houston", "Los Ángeles", "Milwaukee", "Montreal",
           "New York", "Philadelphia", "Pittsburgh", "San Diego",
           "San Francisco", "St. Louis"),
  Proportion_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)
)

# Ver la estructura del dataframe
str(equipos)
## 'data.frame':    16 obs. of  4 variables:
##  $ Team          : chr  "Arizona" "Atlanta" "Chicago" "Cincinnati" ...
##  $ Proportion_Won: num  0.519 0.623 0.543 0.426 0.457 0.562 0.537 0.525 0.42 0.512 ...
##  $ HR            : num  152 235 172 182 198 157 191 124 196 144 ...
##  $ ERA           : num  3.86 4.11 3.84 5.13 5.27 ...
library(knitr)

kable(equipos, 
      col.names = c("Equipo","Proporción de Juegos ganados", "HR (Home rones)", "ERA (carreas permitidas)"),
      caption = "Ejercicio 2",
      align = c('l', 'c', 'c','c'))
Ejercicio 2
Equipo Proporción de Juegos ganados HR (Home rones) ERA (carreas permitidas)
Arizona 0.519 152 3.857
Atlanta 0.623 235 4.106
Chicago 0.543 172 3.842
Cincinnati 0.426 182 5.127
Colorado 0.457 198 5.269
Florida 0.562 157 4.059
Houston 0.537 191 3.880
Los Ángeles 0.525 124 3.162
Milwaukee 0.420 196 5.058
Montreal 0.512 144 4.027
New York 0.410 124 4.517
Philadelphia 0.531 166 4.072
Pittsburgh 0.463 163 4.664
San Diego 0.395 128 4.904
San Francisco 0.621 180 3.734
St. Louis 0.525 196 4.642

1.-Estimación del modelo de regresión

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.

Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente la Cantidad de Home runs, y como variable dependiente o respuesta el Proporción de juegos ganados de los equipos

Tomando la varible explicativa Home runs

modelo_mult2 <- lm(Proportion_Won ~ HR , data = equipos)

modelo_mult2
## 
## Call:
## lm(formula = Proportion_Won ~ HR, data = equipos)
## 
## Coefficients:
## (Intercept)           HR  
##    0.354017     0.000888

El modelo de regresión lineal estimado es: $_{Proporción Won} = 0.354017 + 0.000888(HR)

Donde: \(\hat{\beta}_0\)=0.354017, esto sería la proporción de juegos ganados esperados cuando no hay Home runs 0

\(\hat{\beta}_1\)=0.000888, es decir, por cada home runs adicional, la proporción de juegos ganados incrementa $ 0.354017$.

 ggplot(equipos,aes(x=HR,y=Proportion_Won))+
  geom_point(sixe=4)+theme_light()+
  ggtitle("grafico de dispersion")
## Warning in geom_point(sixe = 4): Ignoring unknown parameters: `sixe`

Notemos que No la relación entre Proporción de juegos ganados en función de la cantidad de Home runs es nula.

Gráfico de Modelo de regresión lineal ajustada

plot(equipos$HR, equipos$Proportion_Won, col=5, pch=18,
     xlab = "Home runs", ylab="Proporción de juegos ganados", main = "Modelo de regresión ajustado  a las \nvariables Proporción y Home runs")
abline(modelo_mult2)

Notemos que los puntos estan alejados del modelo de regresión lineal ajustado, es decir, este no es el mejor modelo, podría haber mejor modelo con respecto a las demaás variables explicativas.

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.

Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente la Carreras permitidas , y como variable dependiente o respuesta el Proporción de juegos ganados de los equipos.

Utilizando variables explicativas cuantitativas

modelo_mult3 <- lm(Proportion_Won ~ ERA, data = equipos)

modelo_mult3
## 
## Call:
## lm(formula = Proportion_Won ~ ERA, data = equipos)
## 
## Coefficients:
## (Intercept)          ERA  
##     0.86474     -0.08367

El modelo de regresión lineal estimado es: $_{Proporción Won} = 0.86474 -0.08367 (ERA)

Donde: \(\hat{\beta}_0\)=0.86474 , esto sería la proporción de juegos ganados esperados cuando no hay ninguna carrera permitida.

\(\hat{\beta}_1\) =-0.08367 es decir, por cada ERA aumente, la proporción de juegos ganados $0.86474 $. Si un equipo sube de 4 a 5 carreras permitidas, se esperaría que su proporción de victorias disminuya \(-0.08367\).

 ggplot(equipos,aes(x=ERA,y=Proportion_Won))+
  geom_point(sixe=4)+theme_light()+
  ggtitle("grafico de dispersion")
## Warning in geom_point(sixe = 4): Ignoring unknown parameters: `sixe`

Notemos que la relación podría ser negativa entre Proporción de juegos ganados en función de las carreras permitidas. Es decir, a menor carreras permitidas mayor proporción de carreras ganadas.

Gráfico de modelo de regresión lineal ajustada

plot(equipos$ERA, equipos$Proportion_Won, col=6, pch=18,
     xlab = "Partidas permitidas", ylab="Proporción de juegos ganados", main = "Modelo de regresión ajustado  a las \nvariables Proporción y Promedio de partidas permitidas")
abline(modelo_mult3)

Notemos que a más partidas permitidas menor es la proporcionde juegos ganados.

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.

Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente la Carreras permitidas y Cantidad de Home runs, y como variable dependiente o respuesta el Proporción de juegos ganados de los equipos.

Utilizando variables explicativas cuantitativas

modelo_mult4 <- lm(Proportion_Won ~ HR + ERA, data = equipos)

modelo_mult4
## 
## Call:
## lm(formula = Proportion_Won ~ HR + ERA, data = equipos)
## 
## Coefficients:
## (Intercept)           HR          ERA  
##    0.709188     0.001401    -0.102597

El modelo de regresión lineal estimado es: $_{Proporción Won} = 0.709188 + 0.001401(HR)-0.102597(ERA)

Donde: \(\hat{\beta}_0\)=0.709188, esto sería la proporción de juegos ganados esperados cuando no hay Home runs 0 y no hay carreras permitidas

\(\hat{\beta}_1\)=0.001401, es decir, por cada home runs adicional, la proporción de juegos ganados incrementa \(0.001401\). Si un equipo no logra ningún home run, y su ERA es 0, se esperaría que gane el \(70.91%\) de sus juegos

\(\hat{\beta}_2\) =-0.102597 es decir, por cada ERA aumente, la proporción de juegos ganados \(-0.102597\). Si un equipo sube de 4 a 5 carreras permitidas, se esperaría que su proporción de victorias disminuya \(-0.102597\).

 ggplot(equipos,aes(x=HR,y=Proportion_Won,color=ERA))+
  geom_point(sixe=4)+theme_light()+
  ggtitle("grafico de dispersion")
## Warning in geom_point(sixe = 4): Ignoring unknown parameters: `sixe`

Notemos que la proporción de juegos ganados,se ve que es donde los equipos tienen más honrones y un promedio menor de carreras permitidas.

Tomamos en cuenta otro gráfico

ggplot(data= equipos,
       aes(HR, Proportion_Won,color=ERA))+
  geom_point(size=4)+
  labs(title ="Modelo del Proporción de juegos en función de la cantidad de honrones ",
       subtitle = "De acuerdo con el tipo de el promedio de carreras permitidas")+
  stat_smooth(method = "lm",formula=y~x)+ #lineas y bandas de suavizado (smooth)
  facet_wrap(~ERA) #graficos distintos 

Notemos que a más honrones en las proporciones de carreras ganadas, menor es el promedio de carreras permitidas.

2.-Summary e interpretación del modelo estimado

summary(modelo_mult4)
## 
## Call:
## lm(formula = Proportion_Won ~ HR + ERA, data = equipos)
## 
## 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

Modelo de regresión estimada:

$_{Proporción Won} = 0.709188 + 0.001401(HR)-0.102597(ERA)

Donde: \(\hat{\beta}_0\)=0.709188, esto sería la proporción de juegos ganados esperados cuando no hay Home runs 0 y no hay carreras permitidas

\(\hat{\beta}_1\)=0.001401, es decir, por home runs adicional, la proporción de juegos ganados incrementa \(0.001401\). Si un equipo no logra ningún home run, y su ERA es 0, se esperaría que gane el \(70.91%\) de sus juegos

\(\hat{\beta}_2\) =-0.102597 es decir, por cada ERA aumente, la proporción de juegos ganados \(-0.102597\). Si un equipo sube de 4 a 5 carreras permitidas, se esperaría que su proporción de victorias disminuya \(-0.102597\).

Toamndo en cuenta el p-valor de la variables explicativas cuantitativas, la variable más influyente es ERA con un \(p-valor=2.11e-06\), a comparación de HR que tiene un \(p-valor=7.11e-05\) .

\(p-valor=3.046e-06\) : Modelo estadísticamente significativo en general.

3.- Coeficiente de determinación

\(R^2=0.8583\), el modelo explica el 85.83% de la variación de la proporción de juegos ganados.

\(R^2-ajustado\): Ajusta el \(R^2\) por el numero de pedictores en el modelo.

\(R^2-ajustado=0.8365\), el 83.65% de la variación es explicada por la cantidad de home runs y por el promedio de las carreras permitidas, hay un \(16.35%\) de variabilidad en la proporción de juegos ganados que no se explican por las dos variables. Pero aún así es un ajuste bueno.

4.- Predicción de la variable de respuesta para nuevos valores de las variables independientes.

**Valores ajustados (predichos) de la variable respuesta.

predict(modelo_mult4)
##         1         2         3         4         5         6         7         8 
## 0.5263712 0.6170782 0.5559230 0.4380928 0.4459344 0.5126499 0.5786366 0.5584578 
##         9        10        11        12        13        14        15        16 
## 0.4647810 0.4977246 0.4194393 0.5239219 0.4589828 0.3853370 0.5782086 0.5074612
equipos$Estimados<- predict(modelo_mult4)
equipos

Notemos que los valores estimados estan relativamente cerca de los valores reales, algunos estan subestimados y otros sobrestiman, pero no tanto como para preocuparse.

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

Tomando en cuenta que San diego tuviera \(180\) jonrones y su promedio de carreras permitidas \(4\)

# Datos específicos para Bill Gustin
San_diego<- data.frame(
  HR = 180,
  ERA = 4  # Asumiendo que es Presidente y CEO.
)

San_diego

Predicción

# Predicción
Proporciondejuegos_predicha<- predict(modelo_mult4, newdata = San_diego)

Proporciondejuegos_predicha
##         1 
## 0.5509179

Esto significa que el equipo de San Diego tuviera \(180\) jonrones y su promedio de carreras permitidas \(4\) se esperaría un proporción del 55% de juegos ganados.

5.-Intervalos esperados para los valores esperados

ic4<-predict(object=modelo_mult4,newdata=equipos[,c(2,3,4)],interval="confidence",level=0.95)
ic4
##          fit       lwr       upr
## 1  0.5263712 0.5060660 0.5466763
## 2  0.6170782 0.5773407 0.6568157
## 3  0.5559230 0.5356752 0.5761708
## 4  0.4380928 0.4114482 0.4647373
## 5  0.4459344 0.4149832 0.4768855
## 6  0.5126499 0.4953431 0.5299566
## 7  0.5786366 0.5546411 0.6026321
## 8  0.5584578 0.5208850 0.5960306
## 9  0.4647810 0.4381367 0.4914252
## 10 0.4977246 0.4772638 0.5181853
## 11 0.4194393 0.3892085 0.4496701
## 12 0.5239219 0.5074035 0.5404404
## 13 0.4589828 0.4400625 0.4779030
## 14 0.3853370 0.3511554 0.4195185
## 15 0.5782086 0.5544897 0.6019275
## 16 0.5074612 0.4862100 0.5287124
#Meter las variables que se ocupan
#Rectificar la clase
class(ic4)
## [1] "matrix" "array"
#Convertir a data-frame
ic4<-as.data.frame(ic4)
class(ic4)
## [1] "data.frame"
datequipos<-data.frame(equipos, ic4$lwr,ic4$upr)
datequipos

Con el \(95%\) de confianza sabemos que la Proporción de juegos ganados verdaderolos home runs y las carreen función de las carreras permitidas cae dentro del intervalo de confianza.

6. Intervalos de confianza para las predicciones

ic4.1<-predict(object=modelo_mult4,newdata=equipos[,c(2,3,4)],interval="predict",level=0.95)
ic4.1
##          fit       lwr       upr
## 1  0.5263712 0.4619532 0.5907892
## 2  0.6170782 0.5441642 0.6899922
## 3  0.5559230 0.4915231 0.6203230
## 4  0.4380928 0.3714046 0.5047809
## 5  0.4459344 0.3774117 0.5144570
## 6  0.5126499 0.4491132 0.5761865
## 7  0.5786366 0.5129619 0.6443113
## 8  0.5584578 0.4867006 0.6302150
## 9  0.4647810 0.3980929 0.5314690
## 10 0.4977246 0.4332573 0.5621918
## 11 0.4194393 0.3512390 0.4876396
## 12 0.5239219 0.4605955 0.5872484
## 13 0.4589828 0.3949878 0.5229777
## 14 0.3853370 0.3152958 0.4553781
## 15 0.5782086 0.5126345 0.6437828
## 16 0.5074612 0.4427387 0.5721836
#Rectificar la clase
class(ic4.1)
## [1] "matrix" "array"
#Convertir a data-frame
ic4.1<-as.data.frame(ic4.1)
class(ic4.1)
## [1] "data.frame"
datequipos1<-data.frame(equipos, ic4.1$lwr,ic4.1$upr)
datequipos1

El intervalo de predicción para el equipo de arizona, con 152 home runs,y cantidad de carrera permitidas con proporcion de juegos ganados de \(51.9%\), su intervalo de predicción, está entre 46.19% de proporción de juegos ganados minimo y 59% de proporción de partidas ganadas.

7.- Propiedades de los residuales

1)Suma de los residuales muy cercano al 0 ó es 0.

suma_residuales4<- sum(modelo_mult4$residuals)
suma_residuales4
## [1] -1.301043e-17

Es muy cercano al 0 la suma de los residuales, por lo tanto cumple.

suma_cuadrado_residuales4<-sum(suma_residuales4^2)
suma_cuadrado_residuales4
## [1] 1.692712e-34

Es muy cercano al 0 la suma de los residuales, por lo tanto cumple.

2) El promedio de los residuales es muy cercano al 0.

promedio_residuales4<-mean(modelo_mult4$residuals)
promedio_residuales4
## [1] -8.131516e-19

Es muy cercano al 0, entonces, cumple.

3) La suma de los valores observados coincide con la suma de los valores estimados

suma_observados4<-sum(equipos$Proportion_Won)
suma_observados4
## [1] 8.069
suma_estimados4<-sum(modelo_mult4$fitted.values)
suma_estimados4
## [1] 8.069

Notemos que la suma de los valores observados y de los estimados es igual.

4) La suma de los residuales ponderados por los valores de las variables explicativas es 0

suma_ponderada_x4<-sum((equipos$HR + equipos$ERA)*modelo_mult4$residuals)
suma_ponderada_x4
## [1] -1.554312e-15

Es cercano al 0.

5) El promedio de los valores observados es igual al promedio de los valores estimados.

Valores observados

mean_observados4<-mean(equipos$Proportion_Won)
mean_observados4
## [1] 0.5043125

Valores estimados

mean_estimados4<-mean(modelo_mult4$fitted.values)
mean_estimados4
## [1] 0.5043125

Notemos que el promedio de los valores observados y de los estimados es igual, entonces, cumple.

6) La recta ajustada pasa por el punto \((\bar{X},\bar{Y})\)

mean_xy4<-c(mean(equipos$HR),mean(equipos$ERA),mean(equipos$Proportion_Won))
mean_xy4
## [1] 169.2500000   4.3075000   0.5043125

Gráfico

plot(ejecutivos$Age,ejecutivos$Salary,col=4
     ,pch=18,xlab="Edad de los ejecutivos",ylab = "Salario Anual",main = "Modelo de regresión ajustado a las variables Edad y Salario de los ejecutivos")
abline(modelo_mult1)
## Warning in abline(modelo_mult1): only using the first two of 6 regression
## coefficients
points(59.0,2062.6)

Notemos que la propiedad no cumple, es decir la recta ajustada no pasa por el punto \((\bar{X},\bar{Y})\)

7) La distribución de los residuales sea Normal

hist(modelo_mult4$residuals,col = 4, main = "Histograma de los residuales", xlab="Residuales")

Notemos que no sigue una distribucion, pico no centrado en 0, no es simetrica.

Probemos con otra grafica

library(ggplot2)
equipos$residuales<-modelo_mult4$residuals
equipos$Estimados<- modelo_mult4$fitted.values
head(equipos)
qqnorm(equipos$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= 4, pch=19)
qqline(equipos$residuales)

Vemos que los residuales siguen una distribución normal, los puntos se desvían por encima de la linea cola izquierda, es decir, hay valores negativos mas pequeños de lo esperado.

Los puntos se desvian por debajo de la linea.Esto indica que los residuales negativos son menos extremos de lo que se esperarian y tambien se desvian hacia valores mas altos.

Probamos con otra gráfica

ggplot(data = data.frame(Residuales = modelo_mult4$residuals), 
       aes(x = Residuales)) +
  geom_density(fill = 4, alpha = 0.6, ) +
  labs(title = "Densidad de los Residuales",
       subtitle = "Distribución de los errores del modelo lineal",
       x = "Residuales",
       y = "Densidad") +
  theme_minimal()

Notemos que la grafica tiene sesgo en sus estremos, sin embargo la mayoria de sus residuales se mantienen entre -0.3,0.20 mas o menos.

Gráfico de probabilidad normal de los residuales ajustados

library(car)
qqPlot(modelo_mult4$residuals, 
       pch = 20, 
       col = 4, 
       main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado")

## [1] 6 9

Los residuales se comportan correctamente a excepción de los punto 6 y 9, sin embargo hay varios puntos fuera de la banda de confianza.

La varianza estimada de los residuales, es el cuadrado medio del error, que obtenemos a continuación:

El MSE mide el promedio de los cuadrados de las diferencias entre los valores reales observados y los valores predichos por un modelo (residuos), en terminos simples me permite responder a la pregunta ¿Que tan lejos en promedio estan las predicciones de mi modelo de la realidad?

mse4<-(sum(modelo_mult4$residuals)^2)/16
mse4
## [1] 1.057945e-35
des_estandar4<-sqrt(mse4)
des_estandar4
## [1] 3.252607e-18

Penaliza los errores grandes, pues las diferencias grandes entre la prediccion y el valor real se vuelven mucho mas significativas. Aunque sus unidades estàn al cuadrado un MSE de cercano a 0 significa un ajuste perfecto.

Entonces, el ajuste es perfecto.

La desviación estándar estimada en el error es de 3.252607e-18

8)Verificación de supuestos: Independencia, Varianza Constante y Normalidad de los Residuales

a) Validación Gráfica

plot(modelo_mult4)

La caracteristica a observar en el gráfico de residuales contra ajustados es una dispersión aleatoria de los residuales o de los punto.No debe haber patrones. Al parecer no hay patrones

La segunda gráfica es un qqplot o bien un gráfico de probabilidad normal y permite valorar la normalidad de los residuales.

La tercera gráfica de residuales estandarizados permite valorar la varianza constante, entonces debemos ver puntos dispersos sin patrón alguno. No tenemos patron alguno, pero hay puntos lejos.

La distancia de Cook identifica observaciones que, si se eliminaran, cambiarían significativamente los coeficientes de la regresión, ayudando a diagnosticar problemas en el ajuste del modelo. Notemos que el punto 6,9 y 15, son puntos influyente ya que esta arriba del 1.

Pruebas de normalidad

Para probar si los residuales proviene de una población normal se usa test de shapiro, que contrasta la hipótesis:

\(H_{0}: los\ datos\ provienen\ de\ una\ distribución\ normal\)

\(H_{1}: los\ datos\ no\ provienen\ de\ una\ distribución\ normal\)

shapiro.test(modelo_mult4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_mult4$residuals
## W = 0.95439, p-value = 0.5623

De acuerdo con el \(p-valor=0.5623\) se puede decir con un 95% de confianza que los residuales provienen de una poblacion con distribucion normal.

Análisis de dispersión con utilizando la edad del ejercicio como la variable independiente.

ggplot(equipos,aes(x=HR,y=Proportion_Won,color= ERA))+
  geom_point(sixe=4)+theme_light()+
  ggtitle("grafico de dispersion")
## Warning in geom_point(sixe = 4): Ignoring unknown parameters: `sixe`

Notemos que hay menor home runs cuando las carreras permitidas son mayores, y la proporción disminuye. La relación no es contanet.

Independencia

Para probar Independencia un supuesto muy importante, es importante saber que una muestra de manera indempendiente ie. que no hay patrones de dependencia espacial, temporal o multinivel entre las observaciones. En este supuesto espero que ¿espero que suceda cerca del acceso a la carretera sea similar a los que estan lejos?

En este caso existen pruebas estadisticas para evalual el supuesto de independencia.El test de Durbin-Watson

contrasta la hipótesis:

\(H_{0}: los\ datos\ no\ presentan\ autocorrelacion\)

\(H_{1}: los\ datos\ presentan\ autocorrelacion\)

library(car)

durbinWatsonTest(modelo_mult4)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1609368      1.643358   0.548
##  Alternative hypothesis: rho != 0

Con un p-valor > 0.05 no se puede rechazar la hipotesis nula, por lo que se presenta autocorrelacion.

Homogeneidad de Varianzas

Grafica de homogeneidad de varianzas

plot(equipos$HR, equipos$Proportion_Won, main = "Diagrama de Dispersión\n(Homogeneidad Varianzas)", 
     xlab = "Home runs", ylab = "Proporcion de juegos ganados")
abline(modelo_mult4)
## Warning in abline(modelo_mult4): only using the first two of 3 regression
## coefficients
abline(v = median(equipos$HR), col = 4, lty = 3)

library(lmtest)

# Test de Breusch-Pagan para homogeneidad de varianzas
# Este test es mejor para regresión múltiple
bp_test <- bptest(modelo_mult4)
print(bp_test)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_mult4
## BP = 1.4039, df = 2, p-value = 0.4956
# Interpretación:
# H0: Homocedasticidad (varianzas constantes)
# H1: Heterocedasticidad (varianzas no constantes)
# p-value < 0.05 sugiere problemas de heterocedasticidad

$H_{0}: =1 $

$H_{1}: $

En este caso dado que muchos de los valores no se acercan a la media es posible que no genere ese error

Tomando en cuienta que \(p-valor>0.05\) por lo que son iguales las varianzas por lo que se cuenta con varianza contante.

b) Analisis de puntos atipicos

Prueba de Bonferroni para datos atipicos

Con ayuda de la funcion outlierTest() de la libreria car nos ayuda a idendificar valores atípicos de los residuos studentizados del modelo

Contraste de hipotesis:

\(H_{0}: la\ observación\ i\ no\ es\ un\ valor\ atípico\)

\(H_{1}: la\ observación\ i\ es\ un\ valor\ atípico\)

library(car)
outlierTest(modelo_mult4)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 6 2.023105           0.065924           NA

con un \(p-valor >0.05\) no se rechazar la hipotesis es decir contamos con la observacion de 6 , valor atipicos que son significativos.

Desde una vista mas grafica podemos explorar dos graficos

cook <- cooks.distance(modelo_mult4)
observaciones_influyentes <- which(cook > 4/nrow(equipos))
observaciones_influyentes#Valores influyentes en la predicción
## 8 
## 8

Notamos que le punto influyente es el 8.

influencePlot(modelo_mult4)

Notamos que el valor 6 , 9 y 8 son los puntos mas influyentes, pero el 8 es mas influyente.

c)Validación con pruebas de hipotesis

Prueba de Homocedasticidad

\(H_{0}: La\ varianza \ es \ constante\ en \ los\ residuales\)

\(H_{1}: La\ varianza \ no\ es \ constante\ en \ los \ residuales\)

Prueba de homocedasticidad:

\(H_{0}: Hay\ homocedasticidad \ de \ los \ residuales\)

\(H_{1}: No\ hay \ homocedasticidad \ de \ los \ residuales\)

ncvTest(modelo_mult4) #prueba de homocedosticidad

Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 1.114671, Df = 1, p = 0.29107 Como p-valor es mayor a \(0.05\), entonces no se rechaza la hipotesis Nula, es decir la varianza es constante, es decir hay homocedasticidad de los residuales

Test de Breusch-Pagan (homocedasticidad  los residuos)

\(H_{0}: La\ varianza \ es \ constante\ ( Hay Homocedasticidad de los residuales)\)

\(H_{1}: La\ varianza \ no\ es \ constante\ ( Hay Heterodesaticidad de los residuales)\)

library(lmtest)
#studentized Breusch-Pagan test modelo
bptest(modelo_mult4)
studentized Breusch-Pagan test

data: modelo_mult4 BP = 1.4039, df = 2, p-value = 0.4956

Como el \(p-valor\) > \(0.05\) entonces no se rechaza la hipotesis nula, es decir, la varianza es constante ( Hay homocedasticidad de los residuales), es decir se cumple el supuesto de la homocedasticidad.

Prueba de correlación

Comprobar la independencia para los residuos estudentizados del modelo ajustado. Se realiza a través del Test de Durbin-Watson (asume bajo la hipótesis nula que no existe correlación)

\(H_{0}: No\ existe \ correlación \ entre \ los \ residuales\)

\(H_{1}: Existe\ correlación\ entre \ los \ residuales\)

dwtest(Proportion_Won~HR+ERA, alternative = "two.sided", data=equipos)
Durbin-Watson test

data: Proportion_Won ~ HR + ERA DW = 1.6434, p-value = 0.488 alternative hypothesis: true autocorrelation is not 0

Como \(p-valor=0.488\), entonces no se tiene evidencia para rechazar hipotesis nula, es decir, no existe correlación, entonces se cumple el supuesto de independencia.

9.- Anova del Modelo

El análisis de la varianza, para las pruebas de hipotesis

Contraste de hipotesis:

\(H_{0}: \beta_{1}=0\)

\(H_{1}: \beta_{1} \neq0\)

\(H_{0}: \beta_{2}=0\)

\(H_{1}: \beta_{2} \neq0\)

anova(modelo_mult4)

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{1}\neq0\), es decir, los home runs estan relacionados con la Proporcion de juegos ganados están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{2}\neq0\), es decir, el promedio de carreras permitidas esta muy relacionada con la proporción de juegos ganados , están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

10.- Pruebas de hipotesis para los parametros

El anova también permite realizar la prueba de hipotesis para la significancia de los coeficientes de las variables explicativas, en este caso, al tener solo una variable independiente, sólo hay un contraste:

\(H_{0}: \beta_{1}=0\)

\(H_{1}: \beta_{1} \neq0\)

anova(modelo_mult4)

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{1}\neq0\), es decir, los home runs estan relacionados con la Proporcion de juegos ganados están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

Este contraste se realiza con un estadistico con distribución F, se puede ver en la tabla anova el \(p-valor <0.05\) por lo que se rechaza la hipotesis nula, indicando que el parámetro \(\beta_{2}\neq0\), es decir, el promedio de carreras permitidas esta muy relacionada con la proporción de juegos ganados , están en una relación lineal.(Prueba de liealidad en el rango de valores observados de la variable independiente).

La gráfica con los intervalos de confianza

library(ggplot2)
ggplot(data = equipos, mapping = aes(x = HR, y = Proportion_Won, color= ERA)) +
  geom_point(color = "firebrick", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  labs(title = "Propotion_Won ~ HR+ERA", x = "Home runs", y = "Proporción de juegos ganados") +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

Al parecer no hay tendencia, notemos que los puntos estan fuera de la banda de confianza, lo que significa que los residuales estan alejados del proporcion real de partidas ganas

Tambien vemos que la mayoria de los puntos estan fuera de la banda de confianza, lo que podría estar subestimando la variabilidad.

¿En qué consiste la multicolinealidad?

La multicolinealidad ocurre cuando dos o más variables predictoras en un modelo de regresión están altamente correlacionadas entre sí.

Un ejemplo:

Imagina que quieres predecir el precio de una casa y usas estas variables:

Número de habitaciones

Metros cuadrados de construcción

Problema: Estas dos variables suelen estar muy correlacionadas (a más habitaciones, más metros cuadrados). El modelo no puede separar si el precio aumenta por las habitaciones o por los metros cuadrados, ya que ambas “dicen lo mismo”.

Resultado: Los coeficientes se vuelven inestables y su interpretación pierde sentido.