Las ejercitadoras elípticas se están convirtiendo en una de las máquinas de ejercicio más populares. Su movimiento de bajo impacto es suave y estable, lo que las vuelve la elección preferida por las persona con problemas en rodillas y tobillos. Sin embargo, elegir la ejercitadora adecuada puede resultar un proceso difícil. El precio y la calidad son factores importantes en cualquier decisión de compra. ¿Están asociados los precios altos con las ejercitadoras elípticas de alta calidad? Consumer Reports realizó amplias pruebas para desarrollar una clasificación general basada en facilidad de uso, ergonomía, construcción y rango de ejercicio. A continuación, se muestran los datos de precio (Price) y calificación (Rating) de ocho ejercitadoras elípticas probadas, de las cuales se detallan marca y modelo (Brand and Model) (Consumer Reports,febrero de 2008).
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.4.3
library(ggplot2)
datos1<- data.frame(
Brand_Model = c("Precor 5.31", "Keys Fitness CG2", "Octane Fitness Q37e",
"LifeFitness X1 Basic", "NordicTrack AudioStrider 990",
"Schwinn 430", "Vision Fitness X6100", "ProForm XP 520 Razor"),
Price = c(3700, 2500, 2800, 1900, 1000, 800, 1700, 600),
Rating = c(87, 84, 82, 74, 73, 69, 68, 55)
)
str(datos1)
## 'data.frame': 8 obs. of 3 variables:
## $ Brand_Model: chr "Precor 5.31" "Keys Fitness CG2" "Octane Fitness Q37e" "LifeFitness X1 Basic" ...
## $ Price : num 3700 2500 2800 1900 1000 800 1700 600
## $ Rating : num 87 84 82 74 73 69 68 55
kable(datos1,
col.names = c("Marca y Modelo", "Precio ($)", "Calificación"),
caption = "Datos de Ejercitadoras Elípticas - Precio vs Calificación",
align = c('l', 'c', 'c'))
Marca y Modelo | Precio ($) | Calificación |
---|---|---|
Precor 5.31 | 3700 | 87 |
Keys Fitness CG2 | 2500 | 84 |
Octane Fitness Q37e | 2800 | 82 |
LifeFitness X1 Basic | 1900 | 74 |
NordicTrack AudioStrider 990 | 1000 | 73 |
Schwinn 430 | 800 | 69 |
Vision Fitness X6100 | 1700 | 68 |
ProForm XP 520 Razor | 600 | 55 |
Tomando en cuenta el precio como variable independiente:
plot(datos1$Price, datos1$Rating,
main = " Precio y Calificación ",
xlab = "Precio ($)", ylab = "Calificación (Rating)",
pch = 19, col = "blue", cex = 1.5)
grid()
La gráfica muestra una tendencia positiva entre el precio y la calificación: a medida que el precio aumenta, la calificación también tiende a aumentar. Esto indica que, en general, las ejercitadoras más costosas obtienen mejores calificaciones.
Sin embargo, no todos los puntos se alinean perfectamente con esta tendencia. Esta dispersión sugiere que, aunque el precio es un buen predictor de la calificación, no es el único factor que influye en la calificación.
Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente el precio de las ejercitadoras elipticas y como variable explicativa la calificación de las ejercitadoras elipticas.
El modelo de regresión lineal simple, estimado esta dado con la función lm()
modelo1<-lm(Rating~Price, data=datos1)
modelo1
Call: lm(formula = Rating ~ Price, data = datos1)
Coefficients: (Intercept) Price
58.158492 0.008449
De esta forma el modelo de regresión lineal ajustado es:
\(Rating=58.158492+0.008449*Price\)
\(\hat{\beta}_{0}\) = 58.158492 , esto sería la calificación esperada cuando el precio es 0.
\(\hat{\beta}_{1}\)= 0.008449 es la pendiente. Lo que significa que por cada $1 de aumento en X (precio), Y (Rating) se Incrementa 0.008449 puntos
plot(datos1$Price,datos1$Rating, col=5, pch=18, xlab = "Precio", ylab="Rating", main="Modelo de regresión ajustado a las variables\n Precio y Rating")
abline(modelo1)
Vemos que los puntos están ligeramente alejados de la linea de ajuste, pero aunasi se podria considerar un buen modelo ajustado.
summary(modelo1)
##
## Call:
## lm(formula = Rating ~ Price, data = datos1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2278 -2.9447 -0.0132 4.2417 6.3927
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.158492 4.014418 14.487 6.78e-06 ***
## Price 0.008449 0.001885 4.482 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.383 on 6 degrees of freedom
## Multiple R-squared: 0.77, Adjusted R-squared: 0.7317
## F-statistic: 20.09 on 1 and 6 DF, p-value: 0.004184
Tenemos residuales minimos de -8.22 y como residuales maximos de 6.39.
El modelo de regresión estimada:
\(Rating=58.158492+0.008449*Price\)
\(\hat{\beta}_{0}\) = 58.158492 , esto sería la calificación esperada cuando el precio es 0.
\(\hat{\beta}_{1}\)= 0.008449 es la pendiente.
La interpretación del modelo de regresion estimada es que:
Por cada $100 de aumento en X (precio), Y (Rating) se Incrementa 0.8449 puntos
En promedio, las predicciones se desvían ±5.383 puntos en el rating.
Bondad de ajuste del modelo de regresión lineal simple ajustado
La interpretación del valor de \(R^2-ajustada\) es 73.17% de la variabilidad en el rating se explica por el precio de las ejercitadoras, hay un 26.83% de variabilidad en el rating que no se explican por el precio de las ejercitadoras. Lo que indica un ajuste bastante bueno.
Estimación de rating de las ejercitadoras con un precio de \(1500\)
predict(modelo1, data.frame(Price= 1500))
1
70.8317
Con el 95% de confianza, una ejercitadora de $1,500 tendrá un rating con una estimacion puntual de 70.83 puntos.
Estimación de un intervalo de confianza del rating de las ejercitadoras con un precio de \(1,500\)
predict(modelo1, data.frame(Price= 1500), interval = "confidence")
fit lwr upr
1 70.8317 65.8637 75.79969
Con el 95% de confianza, estimamos que una ejercitadora de \(1,500\) pesos tendrá un rating entre 65.86 y 75.79 puntos, con una estimación puntual de 70.83 puntos.
Predicción de un intervalo de confianza del rating de las ejercitadoras con un precio de \(1,500\)
predict(modelo1, data.frame(Price= 1500), interval = "prediction")
fit lwr upr
1 70.8317 56.75361 84.90979
Con 95% de confianza, se predice que una ejercitadora de \(1,500\) pesos tendrá \(56.75\) puntos como limite inferior del intervalo de predicción y \(56.75\) puntos como limite superior del intervalo de predicción.
**1) Suma de los residuales =0**
suma_residuales<-sum(modelo1$residuals)
suma_residuales
[1] 0
suma_cuadrado_residuales <- sum(suma_residuales^2)
suma_cuadrado_residuales
[1] 0 Cumple la primera propiedad de que la suma de los residuales sean igual a 0
**2) Promedio de los residuales =0**
promedio_residuales<-mean(modelo1$residuals)
promedio_residuales
[1] 0 Cumple tambien la propiedad de que el promedio de los residuales sea igual a 0.
**3) La suma de los valores observados coincide con la suma de los valores estimados
suma_observados<- sum(datos1$Rating)
suma_observados
[1] 592
suma_estimados<- sum(modelo1$fitted.values)
suma_estimados
[1] 592
De esta forma se tiene que, el promedio de los valores observados= al promedio de los valores estimados.
**4) La suma de los residuales ponderados por los valores de la variable
explicativa es 0.
suma_ponderada_x<-sum(datos1$Price*modelo1$residuals)
suma_ponderada_x
[1] -2.842171e-12 Cumple esta propiedad, ya que la suma de los residuales ponderados por los valores de la variable explicativa (Precio) es 0.
**5) La suma de los residuales ponderados por el valor ajustado de la
variable respuesta es 0:
suma_ponderada_y<- sum(modelo1$fitted.values*modelo1$residuals)
suma_ponderada_y
[1] -7.460699e-14
Esta propiedad de igual forma se cumple, ya que la suma de los residuales ponderados es muy cercano al 0
**6) La recta ajustada pasa por el punto $(\bar{X},\bar{Y})$.
El vector de promedios:
prom_xy<-c(mean(datos1$Price),mean(datos1$Rating))
prom_xy
[1] 1875 74
Gráfico con las observaciones y la recta ajustada que pasa por el vector de promedios de la variable explicativa y la variable de respuesta:
plot(datos1$Price, datos1$Rating,
col = 7, pch=18,
xlab = "Precio",
ylab = "Calificación",
main = "Modelo de regresión ajustado a las variables\n Precio y Calificación
de las ejercitadoras")
abline(modelo1)
points(1875, 74,pch=19,col=5)
Notemos que en promedio se tiene una ejecitadora de $1875 con un promedio de calificación de 74 puntos.
**7) La distribución de los residuales sea un Normal**
Recordemos en los residuales, tienen una distribución, veamos como es la forma de sus distribución:
hist(modelo1$residuals, col=4, main="Histograma de los residuales", xlab = "Residuales")
Se ven los datos de los residuales concentrados y siguiendo una distribución normal.
Intentemos obtener la gráfica de densidad para los residuales
library(ggplot2)
datos1$residuales<-modelo1$residuals
datos1$Estimados<- modelo1$fitted.values
head(datos1)
Brand_Model Price Rating residuales Estimados
1 Precor 5.31 3700 87 -2.4190681 89.41907 2 Keys Fitness CG2 2500 84 4.7194972 79.28050 3 Octane Fitness Q37e 2800 82 0.1848559 81.81514 4 LifeFitness X1 Basic 1900 74 -0.2112201 74.21122 5 NordicTrack AudioStrider 990 1000 73 6.3927039 66.60730 6 Schwinn 430 800 69 4.0824647 64.91754
ggplot(datos1, aes(x = residuales)) +
geom_density(alpha = 0.5) +
xlab("Residuals") +
ylab("Densidad") +
ggtitle("Densidad aproximada de los residuales\n del modelo de regresión lineal ajustado")
En este gráfico de ve que los residuales estan concentrados entre -4 y 4, podría considerarse una distibución simetrica.
Otro gráfico que evidencia si la distribución de los datos es normal o no, es el “qqplot”,
qqnorm(datos1$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= 4, pch=19)
qqline(datos1$residuales)
Los puntos siguen bastante bien la línea diagonal, lo que indica que los errores están distribuidos de forma normal.
library(car)
qqPlot(datos1$residuales, pch=20, col=5, main="Gráfico de probabilidad normal\n de los residuales del modelo ajustado")
[1]
8 5
Los residuales se comportan correctamente a excepcion de los puntos 5 y 8 , los extremos estan dentro de la banda de confianza.
Para probar si los residuales provienen de una población normal se usa test de shapiro, que contrasta la hipotesis
\(H_{0}: Los\ residuales \ provienen \ de \ una \ distribución \ normal\) \(H_{0}: Los\ residuales \ no\ provienen \ de \ una \ distribución \ normal\)
shapiro.test(datos1$residuales)
Shapiro-Wilk normality test
data: datos1$residuales W = 0.96152, p-value = 0.8245 De acuerdo con el \(p-valor=0.8245\) se puede decir con un 95% de confianza que los residuales Provienen de una población con ditribución normal.
Es decir que cumple la propiedad de la distribucion sea normal.
Todo modelo de regresión lineal tiene supuestos que son: independencia varianza constante normalidad de los residuales relación lineal Para valorar estos supuestos hay un resumen gráfico que nos da un modelo:
plot(modelo1)
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.
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.
Las gráficas por separado:*
Gráfica de Valores Estimados vs Valores observados
plot(datos1$Rating, modelo1$fitted.values, col=5, main="Observados vs Estimados", xlab="Valores observados", ylab ="Valores estimados", pch=20); abline(0,1)
Los puntos están relativamente cerca de la línea diagonal, lo que indica que el modelo está haciendo buenas predicciones. Tambien que los valores observados son similares a los estimados
Prueba de Homogeneidad de Varianzas de los residuales:
\(H_{0}: \dfrac{\sigma^2_{1}}{\sigma^{2}}=1\)
\(H_{1}: \dfrac{\sigma^2_{1}}{\sigma^{2}}\neq1\)
plot(datos1$Price, datos1$Rating,
main = "Diagrama de Dispersión\n (Homogeneidad de Varianzas)",
xlab = "Price", ylab = "Rating");abline(modelo1)
abline(v = median(datos1$Price), col = 4, lty = 3)
var.test(residuals(modelo1)[datos1$Price > median(datos1$Price)],
residuals(modelo1)[datos1$Price < median(datos1$Price)])
##
## F test to compare two variances
##
## data: residuals(modelo1)[datos1$Price > median(datos1$Price)] and residuals(modelo1)[datos1$Price < median(datos1$Price)]
## F = 0.18639, num df = 3, denom df = 3, p-value = 0.2012
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.01207274 2.87775852
## sample estimates:
## ratio of variances
## 0.1863932
Como \(p-valor > 0.05\), no se rechaza la hipotesis nula, tenemos varianza constante.
Prueba de bonferroni para datos atipicos
Valorando los puntos influyentes
library(car)
outlierTest(modelo1)
No Studentized residuals with Bonferroni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferroni p 8 -2.607448 0.047817 0.38254
influenceIndexPlot(modelo1) #Valores influyentes en la predicción
Estadisticos para valorar los supuestos
influencePlot(modelo1) #puntos influyentes
StudRes Hat CookD 1 -0.6234518 0.5334151 0.2473878 5 1.4669597 0.2188841
0.2529477 8 -2.6074481 0.3243409 0.8298305
Las burbujas 1,5 y 8 son mas grandes lo que indica que son los puntos mas influyentes en el modelo.
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(modelo1) #prueba de homocedosticidad
Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 1.835843, Df = 1, p = 0.17544 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 de 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(modelo1)
studentized Breusch-Pagan test
data: modelo1 BP = 3.7505, df = 1, p-value = 0.05279 Como el \(p-valor\) es igual a 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(Rating~Price, alternative = "two.sided", data=datos1)
Durbin-Watson test
data: Rating ~ Price DW = 1.1985, p-value = 0.06942 alternative hypothesis: true autocorrelation is not 0
Como \(p-valor=0.06\), entonces no se tiene evidencia para rechazar hipotesis nula, es decir, no existe correlación, entonces se cumple el supuesto de independencia.
El análisis de la varianza, para las pruebas de hipotesis
Contraste de hipotesis:
\(H_{0}: \beta_{1}=0\)
\(H_{1}: \beta_{1} \neq0\)
anova(modelo1)
Analysis of Variance Table
Response: Rating Df Sum Sq Mean Sq F value Pr(>F)
Price 1 582.12 582.12 20.087 0.004184 ** Residuals 6 173.88 28.98
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Vemos que la variabilidad explicada por el Price de las ejercitadoras es mayor a la variabilidad de los residuales, es decir, el precio es significativo en el modelo ajustado.
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(modelo1)
Analysis of Variance Table
Response: Rating Df Sum Sq Mean Sq F value Pr(>F)
Price 1 582.12 582.12 20.087 0.004184 ** Residuals 6 173.88 28.98
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
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).
La gráfica con los intervalos de confianza
library(ggplot2)
ggplot(data = datos1, mapping = aes(x = Price, y = Rating)) +
geom_point(color = "firebrick", size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
labs(title = "Rating ~ Price", x = "Price", y = "Rating") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Hay una tendencia positiva, a medida que aumenta el precio, el rating también aumenta. Tambien vemos que la mayoria de los puntos estan dentro de la banda de confianza.
Una tienda de equipo para ejercitarse que vende principalmente equipo caro puso un letrero sobre el área de exhibición que dice: “Calidad: usted obtiene lo que paga.” Con base en su análisis de los datos, ¿considera usted que el letrero refleja de manera justa la relación precio-calidad de las ejercitadoras elípticas?
Es parcialmente cierto, ya que existe una tendencia positiva, los equipos más caros suelen tener mejor calificación, aunque no se relaciona de manera perfecta.
Un gerente de ventas obtuvo los siguientes datos sobre ventas anuales (Annual Sales) y años de experiencia (Years of Experience) de 10 vendedores (Salesperson).
library(kableExtra)
library(ggplot2)
datos2<- data.frame(
Sales_person = c("1", "2", "3", "4", "5", "6", "7", "8","9","10"),
Years_of_Experience = c(1, 3, 4, 4, 6, 8, 10 , 10, 11, 13),
Anual_Sales = c(80, 97, 92, 102, 103, 111, 119, 123,117,136)
)
str(datos2)
## 'data.frame': 10 obs. of 3 variables:
## $ Sales_person : chr "1" "2" "3" "4" ...
## $ Years_of_Experience: num 1 3 4 4 6 8 10 10 11 13
## $ Anual_Sales : num 80 97 92 102 103 111 119 123 117 136
kable(datos2,
col.names = c("Vendedores", "Años de Experiencia", "Ventas Anuales ($1000s)"),
caption = "Datos de Ventas Anuales - Años de Experiencia vs Ventas anuales",
align = c('l', 'c', 'c'))
Vendedores | Años de Experiencia | Ventas Anuales ($1000s) |
---|---|---|
1 | 1 | 80 |
2 | 3 | 97 |
3 | 4 | 92 |
4 | 4 | 102 |
5 | 6 | 103 |
6 | 8 | 111 |
7 | 10 | 119 |
8 | 10 | 123 |
9 | 11 | 117 |
10 | 13 | 136 |
Tomando en cuenta el Los años de experiencia como variable independiente:
plot(datos2$Years_of_Experience, datos2$Anual_Sales,
main = " Años de experiencia y Ventas anuales ",
xlab = "Años de experiencia", ylab = "Ventas anuales",
pch = 19, col = "green", cex = 1.5)
grid()
Esta gráfica indica una tendencia positiva entre los años de experiencia y la calificacion, ya que a medida de que aumentan los años de experiencia, tambien tienden a aumentar las ventas anuales. Aunque hay cierta dispersion, es decir, aunque los años de experiencia es un buen predictos de las ventas anuales, no es el unico factor que influye en las ventas anuales.
Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente el Años de experiencia de los vendedores y como variable explicativa las ventas anuales de los vendedores.
El modelo de regresión lineal simple, estimado esta dado con la función lm()
modelo2<-lm(Anual_Sales~ Years_of_Experience, data = datos2)
modelo2
Call: lm(formula = Anual_Sales ~ Years_of_Experience, data = datos2)
Coefficients: (Intercept) Years_of_Experience
80 4
De esta forma el modelo de regresión lineal ajustado es:
\(Ventas\ anuales=80 + 4*Años\ de \ experiencia\)
\(\hat{\beta}_{0}\) = 80 , esto sería la venta esperada cuando tienen 0 años de experiencia.
\(\hat{\beta}_{1}\)= 4 es la pendiente.
Lo que significa que por cada 1 de aumento en X (años de experiencia), Y (Ventas anuales) se Incrementa 4 unidades
plot(datos2$Years_of_Experience,datos2$Anual_Sales, col="green", pch=18,
xlab = "Años de experiencia", ylab="Ventas anuales", main="Modelo de regresión ajustado \n a las variables\n Años de experiencia y Ventas anuales")
abline(modelo2)
El modelo sugiere que la experiencia tiene un impacto positivo en las ventas. La dispersión de los puntos alrededor de la línea estan relativamente cerca de esta, es decir, el modelo tiene buen ajuste, aunque tambien puede influir otras variables sobre las ventas.
summary(modelo2)
##
## Call:
## lm(formula = Anual_Sales ~ Years_of_Experience, data = datos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.00 -3.25 -1.00 3.75 6.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.0000 3.0753 26.01 5.12e-09 ***
## Years_of_Experience 4.0000 0.3868 10.34 6.61e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.61 on 8 degrees of freedom
## Multiple R-squared: 0.9304, Adjusted R-squared: 0.9217
## F-statistic: 106.9 on 1 and 8 DF, p-value: 6.609e-06
Tenemos residuales minimos de -7 y como residuales maximos de 6
El modelo de regresión estimada:
\(Ventas\ anuales=80 + 4*Años\ de \ experiencia\)
\(\hat{\beta}_{0}\) = 80 , esto sería la venta esperada cuando tienen 0 años de experiencia.
\(\hat{\beta}_{1}\)= 4 es la pendiente.
La interpretación del modelo de regresion estimada es que:
Lo que significa que por cada 1 de aumento en X (años de experiencia), Y (Ventas anuales) se Incrementa 4 (en miles) unidades en las ventas
En promedio, las predicciones se desvían ±4.61 unidades de las ventas reales.
Bondad de ajuste del modelo de regresión lineal simple ajustado
La interpretación del valor de \(R^2-ajustada\) es 92.17% de la variabilidad en la ventas anuales se explica por los años de experiencia , hay un 7.83% de variabilidad en las ventas reales anuales que no se explican por los años de experiencia. Lo que indica un ajuste bastante bueno.
Estimación de rating de las ejercitadoras con un precio de \(1500\)
predict(modelo2, data.frame(Years_of_Experience= 9))
1 116
Con el 95% de confianza, el modelo predice que un vendedor con \(9\) años de experiencia tendría unas ventas anuales estimadas de \(116\) en miles de pesos.
Estimación de un intervalo de confianza de las ventas anuales de un vendedor con \(9\) años de experiencia
predict(modelo2, data.frame(Years_of_Experience= 9), interval = "confidence")
fit lwr upr 1 116 112.1943 119.8057
Con el 95% de confianza, estimamos que un vendedor con \(9\) años de experiencia tendrá ventas anuales entre 112.19 y 119.80 en miles de pesos, con una estimación puntual de ventas anuales de 116 miles de pesos.
predict(modelo2, data.frame(Years_of_Experience= 9), interval = "prediction")
fit lwr upr 1 116 104.7092 127.2908
Con 95% de confianza, se predice que un vendedor con \(9\) años de experiencia tendrá ventas anuales minimas de \(104.70\) en miles de pesos del intervalo de predicción y ventas anuales maximas de \(127.29\) en miles de pesos del intervalo de predicción.
**1) Suma de los residuales =0**
suma_residuales<-sum(modelo2$residuals)
suma_residuales
[1] 1.221245e-15
Cumple la propiedad de que sus residuales sean 0, ya que esta muy cercano al 0
suma_cuadrado_residuales <- sum(suma_residuales^2)
suma_cuadrado_residuales
[1] 1.49144e-30
Cumple la primera propiedad de que la suma de los residuales sean igual a 0
**2) Promedio de los residuales =0**
promedio_residuales<-mean(modelo2$residuals)
promedio_residuales
[1] 1.220812e-16
Cumple tambien la propiedad de que el promedio de los residuales sea a 0.
**3) La suma de los valores observados coincide con la suma de los valores estimados**
suma_observados<- sum(datos2$Anual_Sales)
suma_observados
[1] 1080
suma_estimados<- sum(modelo2$fitted.values)
suma_estimados
[1] 1080
De esta forma se tiene que, el promedio de los valores observados= al promedio de los valores estimados.
**4) La suma de los residuales ponderados por los valores de la variable
explicativa es 0.**
suma_ponderada_x<-sum(datos2$Years_of_Experience*modelo2$residuals)
suma_ponderada_x
[1] 7.105427e-15
Cumple esta propiedad, ya que la suma de los residuales ponderados por los valores de la variable explicativa (Años de experiencia) es 0.
**5) La suma de los residuales ponderados por el valor ajustado de la
variable respuesta es 0:**
suma_ponderada_y<- sum(modelo2$fitted.values*modelo2$residuals)
suma_ponderada_y
[1] -1.278977e-13
Esta propiedad de igual forma se cumple, ya que la suma de los residuales ponderados es muy cercano al 0
**6) La recta ajustada pasa por el punto $(\bar{X},\bar{Y})$.**
El vector de promedios:
prom_xy<-c(mean(datos2$Years_of_Experience),mean(datos2$Anual_Sales))
prom_xy
[1] 7 108
Gráfico con las observaciones y la recta ajustada que pasa por el vector de promedios de la variable explicativa y la variable de respuesta:
plot(datos2$Years_of_Experience, datos2$Anual_Sales,
col = 3, pch=18,
xlab = "Años de experiencia",
ylab = "Ventas anuales",
main = "Modelo de regresión ajustado a las variables\n Años de experiencia
y ventas anuales ")
abline(modelo2)
points(7, 108,pch=19,col=2)
Vemos que en promedio se necesitan \(7\) años de experiencia para vender anualmente un promedio de \(108\) en miles de pesos.
**7) La distribución de los residuales sea un Normal**
Recordemos en los residuales, tienen una distribución, veamos como es la forma de sus distribución:
hist(modelo2$residuals, col="green", main="Histograma de los residuales", xlab = "Residuales")
Se ve un distribución con datos un poco dispersos en los valores extremos del histograma tanto en valores negativos como en positivos.
Intentemos obtener la gráfica de densidad para los residuales
library(ggplot2)
datos2$residuales<-modelo2$residuals
datos2$Estimados<- modelo2$fitted.values
head(datos2)
Sales_person Years_of_Experience Anual_Sales residuales Estimados 1 1 1 80 -4 84 2 2 3 97 5 92 3 3 4 92 -4 96 4 4 4 102 6 96 5 5 6 103 -1 104 6 6 8 111 -1 112
ggplot(datos2, aes(x = residuales)) +
geom_density(alpha = 0.5) +
xlab("Residuals") +
ylab("Densidad") +
ggtitle("Densidad aproximada de los residuales\n del modelo de regresión lineal ajustado")
Notamos un sesgo al lado izquierdo, hacia los valores negativos, es decir que las predicciones tienden a estar por debajo de loa valores reales, ademas de que no se distribuyen de manera simetrica.
Otro gráfico que evidencia si la distribución de los datos es normal o no, es el “qqplot”,
qqnorm(datos2$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= "green", pch=19)
qqline(datos2$residuales)
Indica un buen modelo ajustado, ya que la mayoria de los puntos se acercan a la linea de regresión.
library(car)
qqPlot(datos2$residuales, pch=20, col="green", main="Gráfico de probabilidad normal\n de los residuales del modelo ajustado")
[1]
9 4 Indica que los puntos estan cerca de la linea, lo cual es bueno, la
mayoria de los puntos estan cerca de esta, solo hay dos puntos que estan
alejados un poco, pero no para alertarse.
Para probar si los residuales provienen de una población normal se usa test de shapiro, que contrasta la hipotesis
\(H_{0}: Los\ residuales \ provienen \ de \ una \ distribución \ normal\) \(H_{0}: Los\ residuales \ no\ provienen \ de \ una \ distribución \ normal\)
shapiro.test(datos2$residuales)
Shapiro-Wilk normality test
data: datos2$residuales W = 0.93759, p-value = 0.5265
De acuerdo con el \(p-valor=0.5265\) se puede decir con un 95% de confianza no hay suficiente evidencia para rechazar la hipotesis nula, lo que significa que los residuales Provienen de una población con ditribución normal.
Es decir que cumple la propiedad de la distribucion sea normal.
Todo modelo de regresión lineal tiene supuestos que son: independencia varianza constante normalidad de los residuales relación lineal Para valorar estos supuestos hay un resumen gráfico que nos da un modelo:
plot(modelo2)
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. En la primera gráfica no indica patrones.
La segunda gráfica es un qqplot o bien un gráfico de probabilidad normal y permite valorar la normalidad de los residuales. Solo hay dos valores alejados de la linea de regresión.
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, lo que indica que podría haber varianza constante
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.
Las gráficas por separado:
Gráfica de Valores Estimados vs Valores observados
plot(datos2$Anual_Sales, modelo2$fitted.values, col="green", main="Observados vs Estimados", xlab="Valores observados", ylab ="Valores estimados", pch=20); abline(0,1)
Los puntos están cerca de la línea diagonal, lo que indica que el modelo está haciendo buenas predicciones. Tambien que los valores observados son similares a los estimados
En la primera gráfica tiene distribución aleatoria es decir no indica patrones.
Prueba de Homogeneidad de Varianzas de los residuales:
\(H_{0}: \dfrac{\sigma^2_{1}}{\sigma^{2}}=1\)
\(H_{1}: \dfrac{\sigma^2_{1}}{\sigma^{2}}\neq1\)
plot(datos1$Price, datos1$Rating,
main = "Diagrama de Dispersión\n (Homogeneidad de Varianzas)",
xlab = "Años de experiencia", ylab = "Ventas anuales");abline(modelo2)
abline(v = median(datos2$Years_of_Experience), col = 4, lty = 3)
var.test(residuals(modelo2)[datos2$Years_of_Experience > median(datos2$Years_of_Experience)],
residuals(modelo2)[datos2$Years_of_Experience< median(datos2$Years_of_Experience)])
##
## F test to compare two variances
##
## data: residuals(modelo2)[datos2$Years_of_Experience > median(datos2$Years_of_Experience)] and residuals(modelo2)[datos2$Years_of_Experience < median(datos2$Years_of_Experience)]
## F = 0.80687, num df = 4, denom df = 4, p-value = 0.8403
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.084009 7.749578
## sample estimates:
## ratio of variances
## 0.806867
Como \(p-valor > 0.05\), no se tiene evidencia suficiente para rechaza la hipotesis nula, es decir, tenemos varianza constante.
Prueba de bonferroni para datos atipicos
Valorando los puntos influyentes
library(car)
outlierTest(modelo2)
No Studentized residuals with Bonferroni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferroni p 9 -2.010637 0.084286 0.84286
influenceIndexPlot(modelo2) #Valores influyentes en la predicción
Estadisticos para valorar los supuestos
influencePlot(modelo2) #puntos influyentes
StudRes Hat CookD 1 -1.092127 0.3535211 0.3184474 4 1.540232 0.1633803
0.1977229 9 -2.010637 0.2126761 0.3955659 10 1.092127 0.3535211
0.3184474
Las burbujas 1,4,9 y 10 son mas grandes lo que indica que son los puntos mas influyentes en el modelo.
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(modelo2) #prueba de homocedosticidad
Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.01096545, Df = 1, p = 0.9166
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 de 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(modelo2)
studentized Breusch-Pagan test
data: modelo2 BP = 0.02775, df = 1, p-value = 0.8677
Como el \(p-valor\) es mayor a 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(Anual_Sales~Years_of_Experience, alternative = "two.sided", data=datos2)
Durbin-Watson test
data: Anual_Sales ~ Years_of_Experience DW = 3.2235, p-value = 0.05844 alternative hypothesis: true autocorrelation is not 0
Como \(p-valor=0.05844\), entonces no se tiene evidencia para rechazar hipotesis nula, es decir, no existe correlación, entonces se cumple el supuesto de independencia.
El análisis de la varianza, para las pruebas de hipotesis
Contraste de hipotesis:
\(H_{0}: \beta_{1}=0\)
\(H_{1}: \beta_{1} \neq0\)
anova(modelo2)
Analysis of Variance Table
Response: Anual_Sales Df Sum Sq Mean Sq F value Pr(>F)
Years_of_Experience 1 2272 2272.00 106.92 6.609e-06 *** Residuals 8 170
21.25
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Vemos que la variabilidad explicada por el Años de experiencia de de los vendedores es mayor a la variabilidad de los residuales, es decir, los años de esperiencia es significativo en el modelo ajustado.
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(modelo2)
Analysis of Variance Table
Response: Anual_Sales Df Sum Sq Mean Sq F value Pr(>F)
Years_of_Experience 1 2272 2272.00 106.92 6.609e-06 *** Residuals 8 170
21.25
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
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, las ventas anuales y el los años de experiencia 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 = datos2, mapping = aes(x = Years_of_Experience, y = Anual_Sales)) +
geom_point(color = "green", size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
labs(title = "Anual_Sales ~ Years_of_Experience", x = "Years_of_ Experience", y = "Anual_Sales") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Vemos que algunos puntos salen del banda de confianzas.Hay una tendencia positiva, a medida que aumenta los años de experiencia, las ventas también aumentan. Aunque no solo dependen de la variable explicativa años de servicio. El modelo puede usarse para hacer predicciones aproximadas, pero no garantiza exactitud en cada caso.
Con frecuencia, los viajeros de las aerolíneas empacan tanto como pueden en su maleta para evitar las tarifas por sobreequipaje. Encontrar una maleta rodante durable, con gran capacidad y fácil de mover puede ser difícil. La tabla siguiente muestra los resultados de pruebas realizadas por Consumer Reports en 10 maletas rodantes; las puntuaciones (Score) más altas indican mejores resultados en las pruebas en general (sitio web Consumer Reports, octubre de 2008). La tabla incluye marca (Brand) y precio (Price) de las maletas.
library(kableExtra)
library(ggplot2)
datos3 <- data.frame(
Brand = c("Briggs & Riley", "Hartman", "Heys", "Kenneth Cole Reaction",
"Liz Claiborne", "Samsonite", "Titan", "TravelPro", "Tumi", "Victorinox"),
Price = c(325, 350, 67, 120, 85, 180, 360, 156, 595, 400),
Score = c(72, 74, 54, 54, 64, 57, 66, 67, 87, 77)
)
kable(datos3,
col.names = c("Marca", "Precio ($)", "Puntuación"),
caption = "Comparativa de marcas de equipaje: precio vs puntuación",
align = c('l', 'c', 'c'))
Marca | Precio ($) | Puntuación |
---|---|---|
Briggs & Riley | 325 | 72 |
Hartman | 350 | 74 |
Heys | 67 | 54 |
Kenneth Cole Reaction | 120 | 54 |
Liz Claiborne | 85 | 64 |
Samsonite | 180 | 57 |
Titan | 360 | 66 |
TravelPro | 156 | 67 |
Tumi | 595 | 87 |
Victorinox | 400 | 77 |
Tomando en cuenta el El precio de las maleta como variable independiente:
plot(datos3$Price, datos3$Score,
main = " Precio y Puntuaciones de las maletas ",
xlab = "Precio de maleta", ylab = "Puntuaciones de las maletas",
pch = 19, col = "purple3", cex = 1.5)
grid()
Notamos en la gráfica de dispersión que no hay una tendencia clara, ya que algunas maletas caras tienen baja puntuacion, al igual que las maletas baratas que algunos tienen mas puntuación. Vemos que los puntos estan bastante dispersos.
Indica que Hay maletas baratas con puntuaciones altas y maletas caras con puntuaciones bajas, lo que rompe la idea de que “más caro es mejor”.De igual forma no se observa una tendencia ascendente ni descendente, por lo tanto, no hay correlación evidente entre las dos variables.
Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente el precio de las maletas y como variable explicativa la puntuación de las maletas.
El modelo de regresión lineal simple, estimado esta dado con la función lm()
modelo3<-lm(Score~Price, data=datos3)
modelo3
Call: lm(formula = Score ~ Price, data = datos3)
Coefficients: (Intercept) Price
52.31050 0.05644
De esta forma el modelo de regresión lineal ajustado es:
\(Puntuación= 52.31050 + 0.05644*Precio\)
\(\hat{\beta}_{0}\) = 52.31050 , esto sería la puntuación esperada cuando cuando la variable x vale 0.
\(\hat{\beta}_{1}\)= 0.05644 es la pendiente.
Lo que significa que por cada 1 peso de aumento en X (precio), Y (la puntuación ) se Incrementa 0.05644 puntos
plot(datos3$Price,datos3$Score, col="purple", pch=18,
xlab = "Precio", ylab="Puntuación", main="Modelo de regresión ajustado \n a las variables\n precio y puntuacion de las maletas")
abline(modelo3)
Notemos que el modelo de regresión ajustado, logra que los puntos se encuentren cerca de la diagonal y sigan una tendencia positiva, de igual forma podría indicar que a mayor precio mejor puntuación pero no siempre , es decir el precio no explica toda la variabilidad de las puntuaciones, lo cual indica que podría haber otras variables explicativas relacionadas con las puntuaciones.
El modelo puede ser útil para predecir puntuaciones aproximadas según el precio, pero no garantiza precisión caso por caso.
summary(modelo3)
##
## Call:
## lm(formula = Score ~ Price, data = datos3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.630 -4.336 1.226 2.068 6.892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.310495 3.017956 17.333 1.25e-07 ***
## Price 0.056442 0.009768 5.778 0.000415 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.968 on 8 degrees of freedom
## Multiple R-squared: 0.8067, Adjusted R-squared: 0.7826
## F-statistic: 33.39 on 1 and 8 DF, p-value: 0.0004153
Tenemos residuales minimos de -8.630 y como residuales maximos de 6.89.
El modelo de regresión estimada:
\(Score=52.3104 + 0.056442 *Price\)
\(\hat{\beta}_{0}\) = 52.3104 , esto sería la calificación esperada cuando el precio es 0.
\(\hat{\beta}_{1}\)= 0.056442 es la pendiente.
Es decir, por cada cada unidad monetaria que aumenta el precio de una maleta, su puntuación promedio aumenta en aproximadamente 0.056 puntos.
La interpretación del modelo de regresion estimada es que:
Por cada $100 de aumento en X (precio), Y (Puntuación) se Incrementa 5.64 puntos
En promedio, las predicciones se desvían ±4.968 puntos en la puntuación.
Bondad de ajuste del modelo de regresión lineal simple ajustado
La interpretación del valor de \(R^2-ajustada\) es 78.26% de la variabilidad en el puntaje se explica por el precio de las maletas , hay un 21.74% de variabilidad en el puntaje que no se explican por el precio de las maletas. Lo que indica un ajuste bueno.
Estimación de la puntuación de la maleta de la marca Eagle Creek Hovercraft con un precio de \(225\)
predict(modelo3, data.frame(Price= 225))
1
65.01003
Con el 95% de confianza, el modelo predice que la maleta con precio de \(225\) tendría una puntuación e estimadas de \(65.01\) puntos.
Estimación de un intervalo de confianza del puntaje de una maleta con precio de \(225\)
predict(modelo3, data.frame(Price= 225), interval = "confidence")
fit lwr upr
1 65.01003 61.28321 68.73686
Con el 95% de confianza, estimamos que una maleta con precio de \(225\) tendrá una puntuación entre 61.28 y 68.73 puntos, con una estimación puntual de 65.01 puntos.
Predicción de un intervalo de confianza del puntaje de una maleta con precio de \(225\)
predict(modelo3, data.frame(Price= 225), interval = "prediction")
fit lwr upr
1 65.01003 52.96248 77.05759
Con 95% de confianza, se predice que una maleta con precio \(225\) tendrá puntaje minimo de \(52.96\) puntos del intervalo de predicción y puntaje maximo de \(77.0575\) puntos del intervalo de predicción.
**1) Suma de los residuales =0**
suma_residuales<-sum(modelo3$residuals)
suma_residuales
[1] -4.440892e-16
Cumple la propiedad de que sus residuales sean 0, ya que esta muy cercano al 0
suma_cuadrado_residuales <- sum(suma_residuales^2)
suma_cuadrado_residuales
[1] 1.972152e-31
Cumple la primera propiedad de que la suma de los residuales sean igual a 0
**2) Promedio de los residuales =0**
promedio_residuales<-mean(modelo2$residuals)
promedio_residuales
[1] 1.220812e-16
Cumple tambien la propiedad de que el promedio de los residuales sea a 0.
**3) La suma de los valores observados coincide con la suma de los valores estimados**
suma_observados<- sum(datos3$Score)
suma_observados
[1] 672
suma_estimados<- sum(modelo3$fitted.values)
suma_estimados
[1] 672
De esta forma se tiene que, el promedio de los valores observados= al promedio de los valores estimados.
**4) La suma de los residuales ponderados por los valores de la variable
explicativa es 0.**
suma_ponderada_x<-sum(datos3$Price*modelo3$residuals)
suma_ponderada_x
[1] -1.136868e-13
Cumple esta propiedad, ya que la suma de los residuales ponderados por los valores de la variable explicativa (Precio) es 0.
**5) La suma de los residuales ponderados por el valor ajustado de la
variable respuesta es 0:**
suma_ponderada_y<- sum(modelo3$fitted.values*modelo3$residuals)
suma_ponderada_y
[1] -1.421085e-14
Esta propiedad de igual forma se cumple, ya que la suma de los residuales ponderados es muy cercano al 0
**6) La recta ajustada pasa por el punto $(\bar{X},\bar{Y})$.**
El vector de promedios:
prom_xy<-c(mean(datos3$Price),mean(datos3$Score))
prom_xy
[1] 263.8 67.2
Gráfico con las observaciones y la recta ajustada que pasa por el vector de promedios de la variable explicativa y la variable de respuesta:
plot(datos3$Price, datos3$Score,
col = "purple", pch=18,
xlab = "Precio",
ylab = "Puntaje",
main = "Modelo de regresión ajustado a las variables\n Precio y puntaje de
las maletas")
abline(modelo3)
points(263.8, 67.2,pch=19,col=2)
Vemos que en promedio las maletas con precio de \(263.8\) tiene un puntaje promedio de \(67.2\) puntos
**7) La distribución de los residuales sea un Normal**
Recordemos en los residuales, tienen una distribución, veamos como es la forma de sus distribución:
hist(modelo3$residuals, col="purple", main="Histograma de los residuales", xlab = "Residuales")
Se ve una distribución de los datos de los residuales un poco dispersos, hacia valores positivos y negativos, sin embargo el intervalo con más frecuencia es aproximadamente de \(0\) a \(2\).
Intentemos obtener la gráfica de densidad para los residuales
library(ggplot2)
datos3$residuales<-modelo3$residuals
datos3$Estimados<- modelo3$fitted.values
head(datos3)
Brand Price Score residuales Estimados
1 Briggs & Riley 325 72 1.345725 70.65427 2 Hartman 350 74 1.934665 72.06533 3 Heys 67 54 -2.092136 56.09214 4 Kenneth Cole Reaction 120 54 -5.083583 59.08358 5 Liz Claiborne 85 64 6.891901 57.10810 6 Samsonite 180 57 -5.470127 62.47013
ggplot(datos3, aes(x = residuales)) +
geom_density(alpha = 0.5) +
xlab("Residuals") +
ylab("Densidad") +
ggtitle("Densidad aproximada de los residuales\n del modelo de regresión lineal ajustado")
Notamos un sesgo al lado derecho, hacia los valores positivos, es decir que las predicciones tienden a estar por arriba de los valores reales, ademas de que no se distribuyen de manera simetrica.
Otro gráfico que evidencia si la distribución de los datos es normal o no, es el “qqplot”,
qqnorm(datos3$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= "purple", pch=19)
qqline(datos3$residuales)
Los puntos siguen bastante bien la línea diagonal, lo que indica que los errores están aproximadamente distribuidos de forma normal.
library(car)
qqPlot(datos3$residuales, pch=20, col="purple", main="Gráfico de probabilidad normal\n de los residuales del modelo ajustado")
[1]
5 7
Indica que los residuales estan cerca de la linea, lo cual es bueno, la mayoria de los puntos estan cerca de esta, a excepción de los punt 5 y 7 que se alejan un de la linea diagona, sin embargo estan dentro de las bandas de confianza,
Para probar si los residuales provienen de una población normal se usa test de shapiro, que contrasta la hipotesis:
\(H_{0}: Los\ residuales \ provienen \ de \ una \ distribución \ normal\) \(H_{0}: Los\ residuales \ no\ provienen \ de \ una \ distribución \ normal\)
shapiro.test(datos3$residuales)
Shapiro-Wilk normality test
data: datos3$residuales W = 0.92836, p-value = 0.432
De acuerdo con el \(p-valor=0.432\) se puede decir con un 95% de confianza no hay suficiente evidencia para rechazar la hipotesis nula, lo que significa que los residuales Provienen de una población con ditribución normal.
Es decir que cumple la propiedad de la distribucion sea normal.
Todo modelo de regresión lineal tiene supuestos que son: independencia varianza constante normalidad de los residuales relación lineal Para valorar estos supuestos hay un resumen gráfico que nos da un modelo:
plot(modelo3)
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. En la primera gráfica no indica patrones y lo putos estan dispersos.
La segunda gráfica es un qqplot o bien un gráfico de probabilidad normal y permite valorar la normalidad de los residuales. Solo hay 3 valores alejados de la linea de regresión.
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, lo que indica que hay varianza constante
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.
Las gráficas por separado:*
Gráfica de Valores Estimados vs Valores observados
plot(datos3$Score, modelo3$fitted.values, col="purple", main="Observados vs Estimados", xlab="Valores observados", ylab ="Valores estimados", pch=20); abline(0,1)
Los puntos están aproximados a la línea diagonal, lo que indica que el modelo es bueno, pero las variables explicativas podrían ser más.
En la primera gráfica tiene distribución aleatoria es decir no indica patrones.
Prueba de Homogeneidad de Varianzas de los residuales:
\(H_{0}: \dfrac{\sigma^2_{1}}{\sigma^{2}}=1\)
\(H_{1}: \dfrac{\sigma^2_{1}}{\sigma^{2}}\neq1\)
plot(datos3$Price, datos3$Score,
main = "Diagrama de Dispersión\n (Homogeneidad de Varianzas)",
xlab = "Precio de las maletas", ylab = "Puntaje de las maletas");abline(modelo3)
abline(v = median(datos3$Price), col = "purple", lty = 3)
var.test(residuals(modelo3)[datos3$Price> median(datos3$Price)],
residuals(modelo3)[datos3$Price< median(datos3$Price)])
##
## F test to compare two variances
##
## data: residuals(modelo3)[datos3$Price > median(datos3$Price)] and residuals(modelo3)[datos3$Price < median(datos3$Price)]
## F = 0.38795, num df = 4, denom df = 4, p-value = 0.3814
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.04039202 3.72604213
## sample estimates:
## ratio of variances
## 0.3879463
Como \(p-valor > 0.05\), no se tiene evidencia suficiente para rechaza la hipotesis nula, es decir, tenemos varianza constante.
Prueba de bonferroni para datos atipicos
Valorando los puntos influyentes
library(car)
outlierTest(modelo3)
No Studentized residuals with Bonferroni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferroni p 5 1.772613 0.11958 NA
influenceIndexPlot(modelo3) #Valores influyentes en la predicción
Estadisticos para valorar los supuestos
influencePlot(modelo3) #puntos influyentes
StudRes Hat CookD 3 -0.4616318 0.2497136 0.03933185 5 1.7726129
0.2235794 0.35685531 7 -1.5583377 0.1357735 0.16185708 9 0.3038965
0.5240251 0.05734414
Las burbujas 5,7 son mas grandes lo que indica que son los puntos mas influyentes en el modelo.
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(modelo3) #prueba de homocedosticidad
Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.9216424, Df = 1, p = 0.33704 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 de 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(modelo3)
studentized Breusch-Pagan test
data: modelo3 BP = 2.3117, df = 1, p-value = 0.1284
Como el \(p-valor\) es mayor a 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(Score~Price, alternative = "two.sided", data=datos3)
Durbin-Watson test
data: Score ~ Price DW = 2.55, p-value = 0.3513 alternative hypothesis: true autocorrelation is not 0
Como \(p-valor=0.3513\), entonces no se tiene evidencia para rechazar hipotesis nula, es decir, no existe correlación, entonces se cumple el supuesto de independencia.
El análisis de la varianza, para las pruebas de hipotesis
Contraste de hipotesis:
\(H_{0}: \beta_{1}=0\)
\(H_{1}: \beta_{1} \neq0\)
anova(modelo3)
Analysis of Variance Table
Response: Score Df Sum Sq Mean Sq F value Pr(>F)
Price 1 824.14 824.14 33.389 0.0004153 *** Residuals 8 197.46
24.68
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Vemos que la variabilidad explicada por el precio de las maletas es mayor a la variabilidad de los residuales, es decir, el precio de las maletas es significativo en el modelo ajustado.
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(modelo3)
Analysis of Variance Table
Response: Score Df Sum Sq Mean Sq F value Pr(>F)
Price 1 824.14 824.14 33.389 0.0004153 *** Residuals 8 197.46
24.68
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
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 puntaje y el precio de las maletas 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 = datos3, mapping = aes(x = Price, y = Score)) +
geom_point(color = "purple", size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
labs(title = "Puntaje ~ Precio", x = "Precio", y = "Puntaje") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Existe una relación positiva entre el precio y la puntuación,en general, las maletas más caras tienden a recibir mejores valoraciones. La dispersión indica que el precio no explica todo; hay otros factores que también afectan la puntuación (como calidad, marca, diseño). El modelo puede usarse para hacer predicciones aproximadas, pero no garantiza exactitud en cada caso.
Para el Internal Revenue Service (Servicio de Administración Tributaria de Estados Unidos), el carácter razonable de las deducciones declaradas por un contribuyente depende de su ingreso bruto ajustado.Deducciones grandes que comprenden donaciones de caridad o por atención médica son más apropiadas para contribuyentes que tengan un ingreso bruto ajustado grande. Si las deducciones de una persona son mayores que las deducciones declaradas promedio correspondientes a un determinado nivel de ingresos, aumentan las posibilidades de que se le realice una auditoría. Los datos (en miles de dólares) sobre ingreso bruto ajustado y el monto promedio o razonable de deducciones declaradas se listan a continuación.
library(kableExtra)
datos4 <- data.frame(
Ingreso_bruto_ajustado = c(22, 27, 32, 48, 65, 85, 120),
Deducciones_razonables = c(9.6, 9.6, 10.1, 11.1, 13.5, 17.7, 25.5)
)
# Mostrar la estructura
str(datos4)
## 'data.frame': 7 obs. of 2 variables:
## $ Ingreso_bruto_ajustado: num 22 27 32 48 65 85 120
## $ Deducciones_razonables: num 9.6 9.6 10.1 11.1 13.5 17.7 25.5
# Crear la tabla con formato
kable(datos4,
col.names = c("Ingreso Bruto Ajustado (miles de $)",
"Monto Razonable de Deducciones (miles de $)"),
caption = "Relación entre Ingreso Bruto Ajustado y Deducciones Declaradas",
align = c('c', 'c'))
Ingreso Bruto Ajustado (miles de $) | Monto Razonable de Deducciones (miles de $) |
---|---|
22 | 9.6 |
27 | 9.6 |
32 | 10.1 |
48 | 11.1 |
65 | 13.5 |
85 | 17.7 |
120 | 25.5 |
Tomando en cuenta El ingreso bruto ajustado como variable independiente:
plot(datos4$Ingreso_bruto_ajustado, datos4$Deducciones_razonables,
main = " Ingreso bruto ajustado y Monto razonable de las deducciones declaradas ",
xlab = "Ingreso bruto ajustado", ylab = "Monto razonable de las deducciones declaradas",
pch = 19, col = "yellow", cex = 1.5)
grid()
Esta gráfica indica una tendencia positiva entre el ingreso bruto ajustado y el monto razonable de las deducciones declaradas, ya que a medida de que aumenta el ingreso bruto ajustado, tambien tienden a aumentar el monto razonable de las deducciones declaradas.
Ajustamos o estimamos la recta de regresión lineal simple considerando como variable independiente el Ingreso bruto ajustado y como variable explicativa Monto razonable de las deducciones declaradas .
El modelo de regresión lineal simple, estimado esta dado con la función lm()
modelo4<-lm(Deducciones_razonables~Ingreso_bruto_ajustado, data = datos4)
modelo4
Call: lm(formula = Deducciones_razonables ~ Ingreso_bruto_ajustado, data = datos4)
Coefficients: (Intercept) Ingreso_bruto_ajustado
4.6768 0.1613
De esta forma el modelo de regresión lineal ajustado es:
\(Monto\ de \ deducción\ razonable = 4.6768 + 0.1613* Ingreso\ bruto \ ajustado\)
\(\hat{\beta}_{0}\) = 4.6768 , esto sería el monto de deducción razonable si el ingreso bruto ajustado es 0.
\(\hat{\beta}_{1}\)= 0.1613 es la pendiente.
Lo que significa que por cada $100 de aumento en X ( ingreso bruto ajustado), Y (Monto razonable de las deducciones declaradas ) se Incrementa 16.13 unidades
plot(datos4$Ingreso_bruto_ajustado,datos4$Deducciones_razonables, col="yellow", pch=18,
xlab = "Ingreso bruto ajustado", ylab="Monto razonable de las deducciones declaradas ", main="Modelo de regresión ajustado \n a las variables\n Ingreso bruto ajustado y Monto razonable de las deducciones declaradas ")
abline(modelo4)
Notemos que el modelo de regresión ajustado, logra que los puntos se encuentren cerca de la diagonal y sigan una tendencia positiva, de igual forma indica que a mayor Ingreso bruto ajustado mayor Monto razonable de las deducciones declaradas, es decir el Ingreso explica la variabilidad de las Monto razonable de las deducciones declaradas. El modelo puede ser útil para predecir puntuaciones aproximadas según el Ingreso bruto ajustado.
summary(modelo4)
##
## Call:
## lm(formula = Deducciones_razonables ~ Ingreso_bruto_ajustado,
## data = datos4)
##
## Residuals:
## 1 2 3 4 5 6 7
## 1.3744 0.5679 0.2613 -1.3196 -1.6619 -0.6881 1.4660
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.67675 1.03339 4.526 0.006251 **
## Ingreso_bruto_ajustado 0.16131 0.01568 10.285 0.000149 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.372 on 5 degrees of freedom
## Multiple R-squared: 0.9549, Adjusted R-squared: 0.9458
## F-statistic: 105.8 on 1 and 5 DF, p-value: 0.0001493
El modelo de regresión estimada:
\(Monto\ de \ deducción\ razonable = 4.6768 + 0.1613* Ingreso\ bruto \ ajustado\)
\(\hat{\beta}_{0}\) = 4.6768 , esto sería el monto de deducción razonable si el ingreso bruto ajustado es 0.
\(\hat{\beta}_{1}\)= 0.1613 es la pendiente.
La interpretación del modelo de regresion estimada es que:
Lo que significa que por cada $100 de aumento en X ( ingreso bruto ajustado), Y (Monto razonable de las deducciones declaradas ) se Incrementa 16.13 unidades
En promedio, las predicciones se desvían ±1.37 unidades del Monto razonable de las deducciones declaradas .
Bondad de ajuste del modelo de regresión lineal simple ajustado
La interpretación del valor de \(R^2-ajustada\) es 94.58% de la variabilidad en el Monto razonable de las deducciones declaradas se explica por el ingreso bruto ajustado, hay un 5.42% de variabilidad en el Monto razonable de las deducciones declaradas que no se explican por el ingreso bruto ajustado. Lo que indica un ajuste bastante bueno.
Estimación de Monto razonable de las deducciones declaradas de un contribuyente cuyo ingreso bruto ajustado es de $52,500
predict(modelo4, data.frame(Ingreso_bruto_ajustado= 52.5))
1
13.14553
Con el 95% de confianza, el modelo predice que un contribuyente con uningreso bruto ajustado de \(52,500\) pesos tendría un Monto razonable de las deducciones declaradas estimadas de \(13,145.53\) .
Según la tendencia de los datos, el monto razonable de deducciones para un ingreso de $52,500 sería aproximadamente $13,145.53.. Las deducciones declaradas de $20,000 excede a la deduccion estimada.Lo cuál Si, estaría justificada una auditoria, ya que la diferencia entre lo declarado y lo estimado por el modelo nos dice que el contribuyente está reportando deducciones excesivas.Aunque el modelo no es una regla legal, sirve como herramienta de referencia para detectar posibles incongruencias.
Estimación de un intervalo de confianza del Monto razonable de las deducciones declaradas de un contribuyente cuyo ingreso bruto ajustado es de $52,500
predict(modelo4, data.frame(Ingreso_bruto_ajustado= 52.5), interval = "confidence")
fit lwr upr
1 13.14553 11.80064 14.49042
Con el 95% de confianza, estimamos que un contribuyente cuyo ingreso bruto ajustado es de $52,500 tendrá un Monto razonable de las deducciones declaradas entre $11,800 y $14,490 , con una estimación puntual de Monto razonable de las deducciones declaradas de $13,140
predict(modelo4, data.frame(Ingreso_bruto_ajustado= 52.5), interval = "prediction")
fit lwr upr
1 13.14553 9.372015 16.91905
Con 95% de confianza, se predice que que un contribuyente cuyo ingreso bruto ajustado es de $52,500 tendrá un Monto razonable de las deducciones declaradas minimas de \(9,372.015\) del intervalo de predicción y Monto razonable de las deducciones declaradas maximas de \(16,919.05\) del intervalo de predicción.
**1) Suma de los residuales =0**
suma_residuales<-sum(modelo4$residuals)
suma_residuales
[1] -7.771561e-16
Cumple la propiedad de que sus residuales sean 0, ya que esta muy cercano al 0
suma_cuadrado_residuales <- sum(suma_residuales^2)
suma_cuadrado_residuales
[1] 6.039716e-31
Cumple la primera propiedad de que la suma de los residuales sean igual a 0
**2) Promedio de los residuales =0**
promedio_residuales<-mean(modelo4$residuals)
promedio_residuales
[1] -1.110223e-16
Cumple tambien la propiedad de que el promedio de los residuales sea a 0.
**3) La suma de los valores observados coincide con la suma de los valores estimados**
suma_observados<- sum(datos4$Deducciones_razonables)
suma_observados
[1] 97.1
suma_estimados<- sum(modelo4$fitted.values)
suma_estimados
[1] 97.1
De esta forma se tiene que, el promedio de los valores observados= al promedio de los valores estimados.
**4) La suma de los residuales ponderados por los valores de la variable
explicativa es 0.**
suma_ponderada_x<-sum(datos4$Ingreso_bruto_ajustado*modelo4$residuals)
suma_ponderada_x
[1] -1.953993e-14
Cumple esta propiedad, ya que la suma de los residuales ponderados por los valores de la variable explicativa (Ingreso bruto ajustado) es 0.
**5) La suma de los residuales ponderados por el valor ajustado de la
variable respuesta es 0:**
suma_ponderada_y<- sum(modelo4$fitted.values*modelo4$residuals)
suma_ponderada_y
[1] -8.437695e-15
Esta propiedad de igual forma se cumple, ya que la suma de los residuales ponderados es muy cercano al 0
**6) La recta ajustada pasa por el punto $(\bar{X},\bar{Y})$.**
El vector de promedios:
prom_xy<-c(mean(datos4$Ingreso_bruto_ajustado),mean(datos4$Deducciones_razonables))
prom_xy
[1] 57.00000 13.87143
Gráfico con las observaciones y la recta ajustada que pasa por el vector de promedios de la variable explicativa y la variable de respuesta:
plot(datos4$Ingreso_bruto_ajustado, datos4$Deducciones_razonables,
col = 3, pch=18,
xlab = "Ingreso bruto ajustado",
ylab = "Monto razonable de las deducciones declaradas ",
main = "Modelo de regresión ajustado a las variables\n Ingreso bruto ajustado
y Monto razonable de las deducciones declaradas ")
abline(modelo4)
points(57, 13.87143,pch=19,col="yellow")
Vemos que en promedio para un ingreso bruto ajustado de \(57,000\) pesos se da en promedio un Monto razonable de las deducciones declaradas de \(13,871.43\) pesos.
**7) La distribución de los residuales sea un Normal**
Recordemos en los residuales, tienen una distribución, veamos como es la forma de sus distribución:
hist(modelo4$residuals, col="green", main="Histograma de los residuales", xlab = "Residuales")
Se ve una distribución con datos concentrados en los valores de -2 a 2, tomando en cuenta que la frecuencia de los residuales es menor en el rango de -1 a 0.
Intentemos obtener la gráfica de densidad para los residuales
library(ggplot2)
datos4$residuales<-modelo4$residuals
datos4$Estimados<- modelo4$fitted.values
head(datos4)
Ingreso_bruto_ajustado Deducciones_razonables residuales Estimados 1 22 9.6 1.3744266 8.225573 2 27 9.6 0.5678758 9.032124 3 32 10.1 0.2613251 9.838675 4 48 11.1 -1.3196373 12.419637 5 65 13.5 -1.6619097 15.161910 6 85 17.7 -0.6881127 18.388113
ggplot(datos4, aes(x = residuales)) +
geom_density(alpha = 0.5) +
xlab("Residuals") +
ylab("Densidad") +
ggtitle("Densidad aproximada de los residuales\n del modelo de regresión lineal ajustado")
Los residuales no distribuyen de manera simetrica.
Otro gráfico que evidencia si la distribución de los datos es normal o no, es el “qqplot”,
qqnorm(datos4$residuales, main = "Gráfico de probabilidad normal\n de los residuales del modelo ajustado", col= "yellow", pch=19)
qqline(datos4$residuales)
Los puntos siguen bastante bien la línea diagonal, lo que indica que los errores están aproximadamente distribuidos de forma normal.
library(car)
qqPlot(datos4$residuales, pch=20, col="yellow", main="Gráfico de probabilidad normal\n de los residuales del modelo ajustado")
[1]
5 7
Indica que los puntos estan cerca de la linea, lo cual es bueno, la mayoria de los puntos estan cerca de esta, solo hay dos puntos que estan alejados un poco,pero un siguen dentro de la banda de confianza.
Para probar si los residuales provienen de una población normal se usa test de shapiro, que contrasta la hipotesis
\(H_{0}: Los\ residuales \ provienen \ de \ una \ distribución \ normal\) \(H_{0}: Los\ residuales \ no\ provienen \ de \ una \ distribución \ normal\)
shapiro.test(datos4$residuales)
Shapiro-Wilk normality test
data: datos4$residuales W = 0.92061, p-value = 0.4741
De acuerdo con el \(p-valor=0.4741\) se puede decir con un 95% de confianza no hay suficiente evidencia para rechazar la hipotesis nula, lo que significa que los residuales Provienen de una población con ditribución normal.
Es decir que cumple la propiedad de la distribucion sea normal.
Todo modelo de regresión lineal tiene supuestos que son: independencia varianza constante normalidad de los residuales relación lineal Para valorar estos supuestos hay un resumen gráfico que nos da un modelo:
plot(modelo4)
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. En la primera gráfica no indica patrones.
La segunda gráfica es un qqplot o bien un gráfico de probabilidad normal y permite valorar la normalidad de los residuales. Solo hay dos valores alejados de la linea de regresión. pero se ve que sigue la distribucion normal
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, lo que indica que podría haber varianza constante
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.
Las gráficas por separado:
Gráfica de Valores Estimados vs Valores observados
plot(datos4$Deducciones_razonables, modelo4$fitted.values, col="yellow", main="Observados vs Estimados", xlab="Valores observados", ylab ="Valores estimados", pch=20); abline(0,1)
Los puntos están cerca de la línea diagonal, lo que indica que el modelo está haciendo buenas predicciones. Tambien que los valores observados son similares a los estimados
En la primera gráfica tiene distribución aleatoria es decir no indica patrones.
Prueba de Homogeneidad de Varianzas de los residuales:
\(H_{0}: \dfrac{\sigma^2_{1}}{\sigma^{2}}=1\)
\(H_{1}: \dfrac{\sigma^2_{1}}{\sigma^{2}}\neq1\)
plot(datos4$Ingreso_bruto_ajustado, datos4$Deducciones_razonables,
main = "Diagrama de Dispersión\n (Homogeneidad de Varianzas)",
xlab = "Ingreso bruto ajustado", ylab = "Monto razonable de las deducciones declaradas ");abline(modelo4)
abline(v = median(datos4$Ingreso_bruto_ajustado), col = "yellow", lty = 3)
var.test(residuals(modelo4)[datos4$Ingreso_bruto_ajustado > median(datos4$Ingreso_bruto_ajustado)],
residuals(modelo4)[datos4$Ingreso_bruto_ajustado< median(datos4$Ingreso_bruto_ajustado)])
##
## F test to compare two variances
##
## data: residuals(modelo4)[datos4$Ingreso_bruto_ajustado > median(datos4$Ingreso_bruto_ajustado)] and residuals(modelo4)[datos4$Ingreso_bruto_ajustado < median(datos4$Ingreso_bruto_ajustado)]
## F = 7.7503, num df = 2, denom df = 2, p-value = 0.2286
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1987254 302.2613319
## sample estimates:
## ratio of variances
## 7.750291
Como \(p-valor > 0.05\), no se tiene evidencia suficiente para rechaza la hipotesis nula, es decir, tenemos varianza constante.
Prueba de bonferroni para datos atipicos
Valorando los puntos influyentes
library(car)
outlierTest(modelo4)
No Studentized residuals with Bonferroni p < 0.05 Largest |rstudent|: rstudent unadjusted p-value Bonferroni p 7 2.886681 0.044712 0.31298
influenceIndexPlot(modelo4) #Valores influyentes en la predicción
Estadisticos para valorar los supuestos
influencePlot(modelo4) #puntos influyentes
StudRes Hat CookD 1 1.272477 0.3030297 0.3132108 5 -1.454574 0.1512253
0.1540961 7 2.886681 0.6618163 3.3056454
Las burbujas 1,5 y 7 son mas grandes lo que indica que son los puntos mas influyentes en el modelo.
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(modelo4) #prueba de homocedosticidad
Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.2179547, Df = 1, p = 0.6406
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 de 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(modelo4)
studentized Breusch-Pagan test
data: modelo4 BP = 0.84206, df = 1, p-value = 0.3588
Como el \(p-valor\) es mayor a 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(Deducciones_razonables~Ingreso_bruto_ajustado, alternative = "two.sided", data=datos4)
Durbin-Watson test
data: Deducciones_razonables ~ Ingreso_bruto_ajustado DW = 0.95149, p-value = 0.01382 alternative hypothesis: true autocorrelation is not 0
Como \(p-valor=0.01382\), entonces se tiene evidencia para rechazar la hipotesis nula, es decir, Existe correlación, entonces no se cumple el supuesto de independencia.
El análisis de la varianza, para las pruebas de hipotesis
Contraste de hipotesis:
\(H_{0}: \beta_{1}=0\)
\(H_{1}: \beta_{1} \neq0\)
anova(modelo4)
Analysis of Variance Table
Response: Deducciones_razonables Df Sum Sq Mean Sq F value
Pr(>F)
Ingreso_bruto_ajustado 1 199.008 199.008 105.79 0.0001493 *** Residuals
5 9.406 1.881
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Vemos que la variabilidad explicada por Ingreso bruto ajustado de los contribuyentes es mayor a la variabilidad de los residuales, es decir, los ingresos brutos ajustados es significativo en el modelo ajustado.
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(modelo4)
Analysis of Variance Table
Response: Deducciones_razonables Df Sum Sq Mean Sq F value
Pr(>F)
Ingreso_bruto_ajustado 1 199.008 199.008 105.79 0.0001493 *** Residuals
5 9.406 1.881
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
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 Monto razonable de las deducciones declaradas y el ingreso bruto ajustado 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 = datos4, mapping = aes(x = Ingreso_bruto_ajustado, y = Deducciones_razonables)) +
geom_point(color = "yellow", size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
labs(title = "Monto razonable de las deducciones declaradas ~ Ingreso bruto ajustado", x = "Ingreso bruto ajustado", y = "Monto razonable de las deducciones declaradas") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
Existe una relación positiva entre el Ingreso bruto ajustado y el Monto razonable de las deducciones declaradas ,en general, los contribuyentes con mas ingreso bruto razonable tienden a dar Montos razonable de las deducciones declaradas más altas. La dispersión indica que el Ingreso bruto ajustado no explica todo; hay otros factores que también afectan Monto razonable de las deducciones declaradas. El modelo puede usarse para hacer predicciones aproximadas.