Linnerrud.csv
library(car)
## Loading required package: carData
library(leaps)
library(faraway)
##
## Attaching package: 'faraway'
## The following objects are masked from 'package:car':
##
## logit, vif
library(MASS)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
Linnerud<-read.csv("Linnerud.csv",row.names=1)
attach(Linnerud) # archivo en uso
lab <- rownames(Linnerud) # etiquetas de las unidades en lab
nom <- colnames(Linnerud) # etiquetas de las variables en nom
n = dim(Linnerud)[1] # número de unidades
p = dim(Linnerud)[2] # número de variables
Pulls
Linnerud1= Linnerud[c(1:4)]
rcic0 = lm(Pulls~1,data=Linnerud1)
rcic1 = lm(Pulls~.,data=Linnerud1)
cicf <- step(rcic0, scope=list(lower=rcic0, upper=rcic1), direction="forward"); summary(cicf)
## Start: AIC=67.58
## Pulls ~ 1
##
## Df Sum of Sq RSS AIC
## + Waist 1 161.919 369.03 62.303
## + Weight 1 80.631 450.32 66.284
## <none> 530.95 67.579
## + Pulse 1 12.050 518.90 69.120
##
## Step: AIC=62.3
## Pulls ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 369.03 62.303
## + Weight 1 18.0713 350.96 63.299
## + Pulse 1 1.1865 367.84 64.239
##
## Call:
## lm(formula = Pulls ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.638 -3.925 1.288 3.296 6.274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.7243 11.5288 3.619 0.00196 **
## Waist -0.9117 0.3244 -2.810 0.01158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.528 on 18 degrees of freedom
## Multiple R-squared: 0.305, Adjusted R-squared: 0.2663
## F-statistic: 7.898 on 1 and 18 DF, p-value: 0.01158
cicb <- step(rcic1, direction="backward"); summary(cicb)
## Start: AIC=65.28
## Pulls ~ Weight + Waist + Pulse
##
## Df Sum of Sq RSS AIC
## - Pulse 1 0.306 350.96 63.299
## - Weight 1 17.190 367.84 64.239
## <none> 350.65 65.281
## - Waist 1 99.624 450.28 68.283
##
## Step: AIC=63.3
## Pulls ~ Weight + Waist
##
## Df Sum of Sq RSS AIC
## - Weight 1 18.071 369.03 62.303
## <none> 350.96 63.299
## - Waist 1 99.359 450.32 66.284
##
## Step: AIC=62.3
## Pulls ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 369.03 62.303
## - Waist 1 161.92 530.95 67.579
##
## Call:
## lm(formula = Pulls ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.638 -3.925 1.288 3.296 6.274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.7243 11.5288 3.619 0.00196 **
## Waist -0.9117 0.3244 -2.810 0.01158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.528 on 18 degrees of freedom
## Multiple R-squared: 0.305, Adjusted R-squared: 0.2663
## F-statistic: 7.898 on 1 and 18 DF, p-value: 0.01158
cics <- step(rcic0,scope=list(lower=rcic0, upper=rcic1), direction="both"); summary(cics)
## Start: AIC=67.58
## Pulls ~ 1
##
## Df Sum of Sq RSS AIC
## + Waist 1 161.919 369.03 62.303
## + Weight 1 80.631 450.32 66.284
## <none> 530.95 67.579
## + Pulse 1 12.050 518.90 69.120
##
## Step: AIC=62.3
## Pulls ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 369.03 62.303
## + Weight 1 18.071 350.96 63.299
## + Pulse 1 1.186 367.84 64.239
## - Waist 1 161.919 530.95 67.579
##
## Call:
## lm(formula = Pulls ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.638 -3.925 1.288 3.296 6.274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.7243 11.5288 3.619 0.00196 **
## Waist -0.9117 0.3244 -2.810 0.01158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.528 on 18 degrees of freedom
## Multiple R-squared: 0.305, Adjusted R-squared: 0.2663
## F-statistic: 7.898 on 1 and 18 DF, p-value: 0.01158
# Análisis de los Modelos Generados
# Selección Forward:
# El modelo comienza sin predictores y añade Waist porque es el que más reduce el AIC y tiene un impacto significativo en la respuesta.
# El modelo final incluye solo Waist con un p-value significativo (0.01158), indicando que es un predictor relevante.
# Selección Backward:
# Comienza con todos los predictores y elimina gradualmente Pulse y Weight, quedándose también solo con Waist.
# Similar al modelo forward, el proceso finaliza con Waist como el único predictor significativo.
# Selección Both:
# Combina ambas técnicas, añadiendo y quitando variables según sea necesario, y también concluye que Waist es el único predictor necesario.
# 1) De acuerdo al análisis anterior Pulls ~Waist
lm_pulls <- lm(Pulls ~ Waist, data = Linnerud)
summary(lm_pulls)
##
## Call:
## lm(formula = Pulls ~ Waist, data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.638 -3.925 1.288 3.296 6.274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.7243 11.5288 3.619 0.00196 **
## Waist -0.9117 0.3244 -2.810 0.01158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.528 on 18 degrees of freedom
## Multiple R-squared: 0.305, Adjusted R-squared: 0.2663
## F-statistic: 7.898 on 1 and 18 DF, p-value: 0.01158
Análisis de Residuos
plot(lm_pulls$fitted.values,Pulls,asp=1,
xlab = "Valores ajustados", # Título eje X
ylab = "Pulls", # Título eje Y
main = "Gráfico de Valores Ajustados vs. Pulls", # Título del gráfico
pch = 1) # plot estimados observados
abline(0,1, col="red")
text(lm_pulls$fitted.values,Pulls,lab=rownames(Linnerud))
# Se observa que hay valores atípicos que pueden influenciar en el rendiemiento del modelo
plot(lm_pulls$fitted.values,lm_pulls$residuals,asp=1,
xlab = "Valores ajustados", ylab = "Residuos",
main = "Gráfico de Residuos vs. Valores Ajustados") # plot estimados
abline(0,0, col="red")
text(lm_pulls$fitted.values, lm_pulls$residuals,
labels = rownames(Linnerud), pos = 4, cex = 0.7)
# Observaciones con alta valor residual como el 11 y 20. Por otro lado, observaciones con alto apalancamiento como la 14 y 9
re <- qqnorm(lm_pulls$residuals) # q-q plot
qqline(lm_pulls$residuals)
text(re$x,re$y,lab=lab,cex=0.4)
# Dado el gráfico de Q-Q plot vemos un ajuste no cae totalmente sobre la línea, por lo que podríamos suponer que no cumple totalmente la normalidad.
residualPlots(lm_pulls) # plot con control de linearidad
## Test stat Pr(>|Test stat|)
## Waist 0.5135 0.6142
## Tukey test 0.5135 0.6076
# Gráfico de Residuos vs. Waist:
# Observaciones: La gráfica muestra una curvatura clara, sugiriendo una relación no lineal entre Waist y Pulls que no está siendo capturada por un modelo lineal simple. Los residuos disminuyen y luego aumentan a medida que Waist aumenta, lo que indica la posible necesidad de un término cuadrático o una transformación de la variable para captar mejor la relación.
# Implicaciones: La presencia de este patrón curvo en los residuos es una señal de que el modelo actual podría ser mejorado incluyendo términos no lineales (como Waist^2) o explorando transformaciones de los datos.
# Gráfico de Residuos vs. Valores Ajustados:
# Observaciones: Similar al primer gráfico, se observa una curvatura en los residuos en relación con los valores ajustados. Este patrón indica que el modelo no está proporcionando un ajuste adecuado a lo largo de todo el rango de valores ajustados, lo que puede llevar a estimaciones sesgadas o ineficientes.
# Implicaciones: Este tipo de patrón sugiere que la relación entre las variables no es completamente lineal y que el modelo podría beneficiarse de la inclusión de más flexibilidad, ya sea mediante términos polinómicos o modelos no paramétricos.
Test de Residuos, Test de Normalidad
# Breausch - Pagan test
bptest(lm_pulls)
##
## studentized Breusch-Pagan test
##
## data: lm_pulls
## BP = 0.65866, df = 1, p-value = 0.417
# NCV test
ncvTest(lm_pulls)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.4037606, Df = 1, p = 0.52515
# Shapiro Test (Normalidad)
shapiro.test(lm_pulls$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm_pulls$residuals
## W = 0.9402, p-value = 0.2419
# Durbin Watson
durbinWatsonTest(lm_pulls)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.06465126 1.836303 0.682
## Alternative hypothesis: rho != 0
# 1. Breusch-Pagan Test
# Resultado: BP = 0.65866, df = 1, p-value = 0.417
# Conclusión: La prueba de Breusch-Pagan no encuentra evidencia de heteroscedasticidad en los residuos del modelo. Un p-value de 0.417 significa que no hay suficiente evidencia para rechazar la hipótesis nula de que la varianza de los residuos es constante
# NCV Test (Non-constant Variance Score Test)
# Resultado: Chisquare = 0.4037606, Df = 1, p = 0.52515
# Conclusión: Similar al test de Breusch-Pagan, el NCV Test también indica que no hay evidencia significativa de heteroscedasticidad. El p-value de 0.52515 respalda la conclusión de que los residuos tienen varianza constante.
# Shapiro-Wilk Test (Normalidad de los Residuos)
# Resultado: W = 0.9402, p-value = 0.2419
# Conclusión: Este resultado sugiere que no hay suficiente evidencia para rechazar la hipótesis nula de que los residuos están distribuidos normalmente. Un p-value de 0.2419 indica que la distribución de los residuos se aproxima a una distribución normal, lo cual es un supuesto importante para la validez de las inferencias de los tests estadísticos en regresión.
# Durbin-Watson Test (Autocorrelación de Residuos)
# Resultado: D-W Statistic = 1.836303, p-value = 0.678
# Conclusión: La estadística Durbin-Watson está cerca de 2, lo cual sugiere que no hay autocorrelación significativa entre los residuos. Un p-value de 0.606 refuerza que la autocorrelación no es una preocupación en este modelo.
Modificación del Modelo
# Modificar el modelo a uno con términos cuadráticos:
lm_pulls_quad <- lm(Pulls ~ Waist + I(Waist^2), data = Linnerud)
summary(lm_pulls_quad)
##
## Call:
## lm(formula = Pulls ~ Waist + I(Waist^2), data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9171 -3.5020 0.7917 3.3146 6.2807
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.35736 87.71286 0.985 0.339
## Waist -3.28259 4.62899 -0.709 0.488
## I(Waist^2) 0.03112 0.06060 0.513 0.614
##
## Residual standard error: 4.623 on 17 degrees of freedom
## Multiple R-squared: 0.3156, Adjusted R-squared: 0.2351
## F-statistic: 3.919 on 2 and 17 DF, p-value: 0.03983
# Ninguno de los coeficientes es estadísticamente significativo a un nivel convencional (p. ej., p<0.05), como lo indican los valores p altos (Waist: p=0.488, I(Waist^2): p=0.614). Esto implica que no hay evidencia suficiente para afirmar que estos coeficientes difieren significativamente de 0 en la población.
# Breausch - Pagan test
bptest(lm_pulls_quad)
##
## studentized Breusch-Pagan test
##
## data: lm_pulls_quad
## BP = 0.93412, df = 2, p-value = 0.6268
# NCV test
ncvTest(lm_pulls_quad)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.3417778, Df = 1, p = 0.5588
# Shapiro Test (Normalidad)
shapiro.test(lm_pulls_quad$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm_pulls_quad$residuals
## W = 0.95408, p-value = 0.4334
# Durbin Watson
durbinWatsonTest(lm_pulls_quad)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.08926724 1.873983 0.714
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: Con un p-valor de 0.6268, la prueba de Breusch-Pagan no muestra evidencia estadística de heteroscedasticidad en los residuos del modelo. Esto significa que la varianza de los residuos parece ser constante, cumpliendo con una de las principales suposiciones de la regresión lineal.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Similar a la prueba de Breusch-Pagan, el NCV Test tiene un p-valor de 0.5588, lo que indica que no hay evidencia significativa de que la varianza de los residuos varíe con los valores ajustados. Esto refuerza la validez de la suposición de homoscedasticidad en tu modelo.
# 3. Prueba de Shapiro-Wilk (Normalidad de los Residuos)
# Resultado: El p-valor de la prueba de Shapiro-Wilk es 0.4334, lo cual es mayor que el nivel típico de significancia de 0.05, sugiriendo que no puedes rechazar la hipótesis nula de que los residuos provienen de una distribución normal. Esto es bueno, ya que la normalidad de los residuos es crucial para la validez de las pruebas de hipótesis e intervalos de confianza en un análisis de regresión.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: El estadístico de Durbin-Watson es 1.873983 con un p-valor de 0.738. Este resultado sugiere que no hay evidencia significativa de autocorrelación de primer orden en los residuos del modelo, lo cual es deseable porque la autocorrelación puede distorsionar las estimaciones estándar de los errores y hacerlas menos confiables.
anova(lm_pulls,lm_pulls_quad)
## Analysis of Variance Table
##
## Model 1: Pulls ~ Waist
## Model 2: Pulls ~ Waist + I(Waist^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 18 369.03
## 2 17 363.39 1 5.6365 0.2637 0.6142
# Diferencia en RSS: La diferencia en la suma de cuadrados de los residuos entre los dos modelos es de 5.6365, lo que indica una reducción muy pequeña en los residuos al añadir el término cuadrático.
# Valor F y p-value: El valor F de 0.2637 con un p-value de 0.6142 sugiere que el término cuadrático añadido no proporciona una mejora significativa en el ajuste del modelo. Un p-value mayor que 0.05 (usualmente el umbral para significancia estadística) indica que no hay evidencia suficiente para rechazar la hipótesis nula de que el modelo más simple es suficiente.
# Dado que el p-value es alto, es preferible mantener el modelo más simple (Pulls ~ Waist) ya que el término cuadrático no mejora significativamente el ajuste del modelo
Squats
head(Linnerud)
## Weight Waist Pulse Pulls Squats Jumps
## 1 191 36 50 5 162 60
## 2 189 37 52 2 110 60
## 3 193 38 58 12 101 101
## 4 162 35 62 12 105 37
## 5 189 35 46 13 155 58
## 6 182 36 56 4 101 42
Linnerud1= Linnerud[c(1:3,5)]
rcic0 = lm(Squats~1,data=Linnerud1)
rcic1 = lm(Squats~.,data=Linnerud1)
cicf <- step(rcic0, scope=list(lower=rcic0, upper=rcic1), direction="forward"); summary(cicf)
## Start: AIC=166.42
## Squats ~ 1
##
## Df Sum of Sq RSS AIC
## + Waist 1 31000.1 43377 157.64
## + Weight 1 18083.4 56294 162.85
## <none> 74377 166.42
## + Pulse 1 3766.6 70610 167.38
##
## Step: AIC=157.64
## Squats ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 43377 157.64
## + Weight 1 1448.37 41929 158.96
## + Pulse 1 0.66 43376 159.64
##
## Call:
## lm(formula = Squats ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.981 -39.461 -6.403 39.016 84.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 592.121 124.992 4.737 0.000164 ***
## Waist -12.615 3.517 -3.587 0.002109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.09 on 18 degrees of freedom
## Multiple R-squared: 0.4168, Adjusted R-squared: 0.3844
## F-statistic: 12.86 on 1 and 18 DF, p-value: 0.002109
cicb <- step(rcic1, direction="backward"); summary(cicb)
## Start: AIC=160.95
## Squats ~ Weight + Waist + Pulse
##
## Df Sum of Sq RSS AIC
## - Pulse 1 16.5 41929 158.96
## - Weight 1 1464.2 43376 159.64
## <none> 41912 160.95
## - Waist 1 14210.1 56122 164.79
##
## Step: AIC=158.96
## Squats ~ Weight + Waist
##
## Df Sum of Sq RSS AIC
## - Weight 1 1448.4 43377 157.64
## <none> 41929 158.96
## - Waist 1 14365.1 56294 162.85
##
## Step: AIC=157.64
## Squats ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 43377 157.64
## - Waist 1 31000 74377 166.42
##
## Call:
## lm(formula = Squats ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.981 -39.461 -6.403 39.016 84.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 592.121 124.992 4.737 0.000164 ***
## Waist -12.615 3.517 -3.587 0.002109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.09 on 18 degrees of freedom
## Multiple R-squared: 0.4168, Adjusted R-squared: 0.3844
## F-statistic: 12.86 on 1 and 18 DF, p-value: 0.002109
cics <- step(rcic0,scope=list(lower=rcic0, upper=rcic1), direction="both"); summary(cics)
## Start: AIC=166.42
## Squats ~ 1
##
## Df Sum of Sq RSS AIC
## + Waist 1 31000.1 43377 157.64
## + Weight 1 18083.4 56294 162.85
## <none> 74377 166.42
## + Pulse 1 3766.6 70610 167.38
##
## Step: AIC=157.64
## Squats ~ Waist
##
## Df Sum of Sq RSS AIC
## <none> 43377 157.64
## + Weight 1 1448.4 41929 158.96
## + Pulse 1 0.7 43376 159.64
## - Waist 1 31000.1 74377 166.42
##
## Call:
## lm(formula = Squats ~ Waist, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.981 -39.461 -6.403 39.016 84.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 592.121 124.992 4.737 0.000164 ***
## Waist -12.615 3.517 -3.587 0.002109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.09 on 18 degrees of freedom
## Multiple R-squared: 0.4168, Adjusted R-squared: 0.3844
## F-statistic: 12.86 on 1 and 18 DF, p-value: 0.002109
# Proceso de Selección de Modelos
# Modelo Final en Todos los Métodos
# En todos los métodos de selección (forward, backward, y both), el modelo final incluye solamente la variable Waist. Esto sugiere que Waist es la variable más significativa para predecir Squats. Los otros predictores, Weight y Pulse, no mejoran significativamente el modelo según los criterios de selección basados en el AIC.
# Modelo inicial (solo intercepto):
# AIC = 166.42
# Adición de Waist (Forward Step):
# Mejora significativa en AIC a 157.64, lo cual sugiere un mejor ajuste del modelo con Waist incluido.
# Intentos de adición de otras variables:
# La adición de Weight y Pulse no reduce significativamente el AIC, lo cual indica que no aportan suficiente valor predictivo adicional después de incluir Waist.
# 1) De acuerdo al análisis anterior Pulls ~Waist
lm_squats <- lm(Squats ~ Waist, data = Linnerud)
summary(lm_squats)
##
## Call:
## lm(formula = Squats ~ Waist, data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.981 -39.461 -6.403 39.016 84.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 592.121 124.992 4.737 0.000164 ***
## Waist -12.615 3.517 -3.587 0.002109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.09 on 18 degrees of freedom
## Multiple R-squared: 0.4168, Adjusted R-squared: 0.3844
## F-statistic: 12.86 on 1 and 18 DF, p-value: 0.002109
Análisis de Residuos
plot(lm_squats$fitted.values,Squats,asp=1,
xlab = "Valores ajustados", # Título eje X
ylab = "Squats", # Título eje Y
main = "Gráfico de Valores Ajustados vs. Squats", # Título del gráfico
pch = 1) # plot estimados observados
abline(0,1, col="red")
text(lm_squats$fitted.values,Squats,lab=rownames(Linnerud))
# Se observa valore atípicos que pueden influenciar en el modelo
plot(lm_squats$fitted.values,lm_squats$residuals,asp=1,
xlab = "Valores ajustados", ylab = "Residuos",
main = "Gráfico de Residuos vs. Valores Ajustados") # plot estimados
abline(0,0, col="red")
text(lm_squats$fitted.values, lm_squats$residuals,
labels = rownames(Linnerud), pos = 4, cex = 0.7)
# Se observa un alto valor residual para las observaciones 16 y 15. Además de un alto apalancamiento para las observaciones 14 y 9
re <- qqnorm(lm_squats$residuals) # q-q plot
qqline(lm_squats$residuals)
text(re$x,re$y,lab=lab,cex=0.4)
# Dado el gráfico de Q-Q plot vemos un ajuste no cae totalmente sobre la línea, por lo que podríamos suponer que no cumple totalmente la normalidad.
residualPlots(lm_squats) # plot con control de linearidad
## Test stat Pr(>|Test stat|)
## Waist 1.3788 0.1858
## Tukey test 1.3788 0.1679
# Gráfico de Residuos vs. Waist:
# Patrón en forma de U: Los residuos muestran un patrón claro en forma de U cuando se grafican contra Waist. Este patrón sugiere que la relación entre Waist y Squats no es lineal, y que un modelo lineal simple puede no ser adecuado para capturar la totalidad de la relación.
# Implicaciones: La forma de U puede indicar la necesidad de un término cuadrático o polinomial para capturar la curvatura en la relación entre las variables. Esto se alinea con el análisis anterior donde se sugirió que puede existir una relación no lineal.
# Gráfico de Residuos vs. Valores Ajustados:
# Patrón Similar: Este gráfico también muestra un patrón en forma de U, lo que es consistente con el primer gráfico y refuerza la idea de que el modelo actual no capta completamente la relación entre las variables.
# Homoscedasticidad: Aparte del patrón no lineal, parece que la varianza de los residuos es constante a través de los valores ajustados, lo cual es bueno porque indica homoscedasticidad.
Test de Residuos, Test de Normalidad
# Breausch - Pagan test
bptest(lm_squats)
##
## studentized Breusch-Pagan test
##
## data: lm_squats
## BP = 0.10308, df = 1, p-value = 0.7482
# NCV test
ncvTest(lm_squats)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.04344458, Df = 1, p = 0.83489
# Shapiro Test (Normalidad)
shapiro.test(lm_squats$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm_squats$residuals
## W = 0.94686, p-value = 0.322
# Durbin Watson
durbinWatsonTest(lm_squats)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.4135982 2.714003 0.07
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: Con un p-valor de 0.7482, la prueba de Breusch-Pagan no muestra evidencia de heteroscedasticidad. Esto significa que la varianza de los residuos se considera constante a lo largo de los valores ajustados, lo cual es una buena señal para la validez de los resultados del modelo.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Similarmente, el NCV Test tiene un p-valor de 0.83489, lo que indica que no hay evidencia significativa de que la varianza de los residuos cambie con los valores ajustados. Este resultado refuerza la suposición de homoscedasticidad en nuestros datos.
# 3. Prueba de Shapiro-Wilk (Normalidad de los Residuos)
# Resultado: El p-valor de 0.322 en la prueba de Shapiro-Wilk sugiere que no podemos rechazar la hipótesis nula de que los residuos siguen una distribución normal. Esto es importante porque la normalidad de los residuos respalda la validez de las pruebas estadísticas sobre los coeficientes y otros aspectos del modelo.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: El estadístico Durbin-Watson es 2.714 con un p-valor de 0.08, lo cual es cercano a indicar la presencia de autocorrelación negativa (valores de D-W mayores a 2 sugieren correlación negativa). Aunque no es altamente significativo (p > 0.05).
Modificación del Modelo
# Modificar el modelo a uno con términos polonómicos:
lm_squats_quad <- lm(Squats ~ Waist + I(Waist^(2)), data = Linnerud)
summary(lm_squats_quad)
##
## Call:
## lm(formula = Squats ~ Waist + I(Waist^(2)), data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73.590 -38.304 0.905 28.543 98.323
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1833.9155 908.8302 2.018 0.0597 .
## Waist -78.5787 47.9629 -1.638 0.1197
## I(Waist^(2)) 0.8657 0.6279 1.379 0.1858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.91 on 17 degrees of freedom
## Multiple R-squared: 0.4755, Adjusted R-squared: 0.4137
## F-statistic: 7.705 on 2 and 17 DF, p-value: 0.004151
# Intercepto: El coeficiente para el intercepto es bastante alto (1833.92) y es marginalmente significativo con un p-valor de 0.0597, lo que indica que al nivel de cintura cero, el valor de Squats sería muy elevado, aunque este punto puede no ser práctico o realista dado el rango de los datos de cintura.
# Waist: El coeficiente de Waist es -78.58 con un p-valor de 0.1197. Aunque el coeficiente es grande en magnitud, no es estadísticamente significativo al nivel usual de 0.05, lo que sugiere que, cuando se ajusta por el término cuadrático, la contribución lineal directa de Waist no es clara.
# I(Waist^2): El coeficiente para el término cuadrático es 0.866 con un p-valor de 0.1858. Este p-valor sugiere que no hay suficiente evidencia para afirmar que el término cuadrático es significativo, aunque su inclusión mejora el ajuste del modelo en comparación con el modelo lineal simple.
# Aunque el modelo mejora el ajuste y el R-cuadrado es relativamente más alto que en el modelo lineal simple, los p-valores de los coeficientes no son lo suficientemente bajos para confirmar la significancia de los términos individuales a niveles convencionales
# Breausch - Pagan test
bptest(lm_squats_quad)
##
## studentized Breusch-Pagan test
##
## data: lm_squats_quad
## BP = 0.9306, df = 2, p-value = 0.6279
# NCV test
ncvTest(lm_squats_quad)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.1066962, Df = 1, p = 0.74394
# Shapiro Test (Normalidad)
shapiro.test(lm_squats_quad$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm_squats_quad$residuals
## W = 0.97451, p-value = 0.8458
# Durbin Watson
durbinWatsonTest(lm_squats_quad)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.4767449 2.783386 0.072
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: BP = 0.9306 con un p-valor de 0.6279.
# Interpretación: Esta prueba no indica evidencia de heteroscedasticidad en los residuos del modelo. Esto significa que la varianza de los errores es aproximadamente constante a lo largo de los valores ajustados, lo cual es una buena señal para la validez de las estimaciones de los errores estándar y las pruebas estadísticas realizadas en los coeficientes del modelo.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Chisquare = 0.1066962, p = 0.74394.
# Interpretación: Al igual que la prueba de Breusch-Pagan, el NCV Test también indica que no hay evidencia de heteroscedasticidad. Esto refuerza la conclusión de que los residuos tienen varianza constante.
# 3. Prueba de Shapiro-Wilk (Normalidad de los Residuos)
# Resultado: W = 0.97451, p-valor = 0.8458.
# Interpretación: Este resultado sugiere que no hay razones para rechazar la hipótesis de que los residuos se distribuyen normalmente. La normalidad de los residuos es importante para la validez de muchas pruebas estadísticas, incluidas aquellas que determinan la significancia de los coeficientes del modelo.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: D-W Statistic = 2.783386, p-valor = 0.094.
# Interpretación: El valor de Durbin-Watson cercano a 3 indica una posible autocorrelación negativa, aunque el p-valor no es lo suficientemente pequeño para confirmar esta suposición de manera definitiva (p > 0.05). Sin embargo, es algo que podría requerir más investigación, especialmente si tu serie de datos es temporal o si las observaciones no son completamente independientes.
anova(lm_squats,lm_squats_quad)
## Analysis of Variance Table
##
## Model 1: Squats ~ Waist
## Model 2: Squats ~ Waist + I(Waist^(2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 18 43377
## 2 17 39014 1 4363.1 1.9012 0.1858
# Reducción en RSS: La inclusión del término cuadrático reduce la suma de cuadrados de los residuos de 43377 a 39014, lo que indica una mejora en la capacidad del modelo para explicar la variabilidad en los datos de Squats.
# Significancia Estadística: Aunque hay una reducción en RSS, el valor F de 1.9012 con un p-value de 0.1858 sugiere que esta mejora no es estadísticamente significativa al nivel convencional de 0.05. Esto implica que, desde un punto de vista estadístico, el término cuadrático no proporciona una mejora suficiente que justifique su inclusión en el modelo según los datos disponibles.
Jumps
Linnerud1= Linnerud[c(1:3,6)]
rcic0 = lm(Jumps~1,data=Linnerud1)
rcic1 = lm(Jumps~.,data=Linnerud1)
cicf <- step(rcic0, scope=list(lower=rcic0, upper=rcic1), direction="forward"); summary(cicf)
## Start: AIC=158.46
## Jumps ~ 1
##
## Df Sum of Sq RSS AIC
## <none> 49958 158.46
## + Weight 1 2558.34 47400 159.41
## + Waist 1 1832.07 48126 159.72
## + Pulse 1 60.96 49897 160.44
##
## Call:
## lm(formula = Jumps ~ 1, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.30 -30.80 -16.30 14.95 179.70
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.30 11.47 6.131 6.8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.28 on 19 degrees of freedom
cicb <- step(rcic1, direction="backward"); summary(cicb)
## Start: AIC=163.36
## Jumps ~ Weight + Waist + Pulse
##
## Df Sum of Sq RSS AIC
## - Waist 1 2.57 47268 161.36
## - Pulse 1 128.47 47394 161.41
## - Weight 1 800.01 48065 161.69
## <none> 47265 163.36
##
## Step: AIC=161.36
## Jumps ~ Weight + Pulse
##
## Df Sum of Sq RSS AIC
## - Pulse 1 131.98 47400 159.41
## - Weight 1 2629.36 49897 160.44
## <none> 47268 161.36
##
## Step: AIC=159.41
## Jumps ~ Weight
##
## Df Sum of Sq RSS AIC
## - Weight 1 2558.3 49958 158.46
## <none> 47400 159.41
##
## Step: AIC=158.46
## Jumps ~ 1
##
## Call:
## lm(formula = Jumps ~ 1, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.30 -30.80 -16.30 14.95 179.70
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.30 11.47 6.131 6.8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.28 on 19 degrees of freedom
cics <- step(rcic0,scope=list(lower=rcic0, upper=rcic1), direction="both"); summary(cics)
## Start: AIC=158.46
## Jumps ~ 1
##
## Df Sum of Sq RSS AIC
## <none> 49958 158.46
## + Weight 1 2558.34 47400 159.41
## + Waist 1 1832.07 48126 159.72
## + Pulse 1 60.96 49897 160.44
##
## Call:
## lm(formula = Jumps ~ 1, data = Linnerud1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.30 -30.80 -16.30 14.95 179.70
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.30 11.47 6.131 6.8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.28 on 19 degrees of freedom
# Análisis de los Resultados de Selección de Modelo
# Modelos Iniciales y Pasos
# Modelo Solo con Intercepto (AIC=158.46):
# Este modelo simplemente predice Jumps usando el promedio de todas las observaciones.
# Ninguna de las variables Weight, Waist o Pulse mejora el modelo lo suficiente como para justificar su inclusión, dado que el AIC aumenta al agregar cualquiera de estas variables.
# Modelo con Todas las Variables (AIC=163.36) y Selección hacia Atrás:
# El modelo inicial incluye todas las variables, pero los pasos hacia atrás eliminan Waist y luego Pulse debido a que su inclusión no reduce suficientemente el AIC.
# Finalmente, se elimina Weight, regresando al modelo solo con intercepto.
# Modelo con Selección Hacia Adelante y Ambos:
# Similar al modelo hacia atrás, la selección hacia adelante y ambos también concluyen que un modelo solo con el intercepto es preferido sobre modelos con más predictores, basado en el AIC.
lm_jumps <- lm(Jumps ~ Waist + Weight + Pulse, data = Linnerud)
summary(lm_jumps)
##
## Call:
## lm(formula = Jumps ~ Waist + Weight + Pulse, data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.89 -34.93 -10.12 14.95 166.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 179.8868 212.2842 0.847 0.409
## Waist 0.2338 7.9276 0.029 0.977
## Weight -0.5379 1.0336 -0.520 0.610
## Pulse -0.3886 1.8634 -0.209 0.837
##
## Residual standard error: 54.35 on 16 degrees of freedom
## Multiple R-squared: 0.0539, Adjusted R-squared: -0.1235
## F-statistic: 0.3039 on 3 and 16 DF, p-value: 0.8222
Análisis de Residuos
plot(lm_jumps$fitted.values,Jumps,asp=1,
xlab = "Valores ajustados", # Título eje X
ylab = "Jumps", # Título eje Y
main = "Gráfico de Valores Ajustados vs. Jumps", # Título del gráfico
pch = 1) # plot estimados observados
abline(0,1, col="red")
text(lm_jumps$fitted.values,Jumps,lab=rownames(Linnerud))
# Se puede apreciar Valores atipicos que como la observación 10. Por otro lado se observar una concentración de los datos en el rango de 0 a 100
plot(lm_jumps$fitted.values,lm_jumps$residuals,asp=1,
xlab = "Valores ajustados", ylab = "Residuos",
main = "Gráfico de Residuos vs. Valores Ajustados") # plot estimados
abline(0,0, col="red")
text(lm_jumps$fitted.values, lm_jumps$residuals,
labels = rownames(Linnerud), pos = 4, cex = 0.7)
# Observaciones con alto valor residual como la 10 y 17. Por otro lado, las observaciones con alto apalancamiento como al 14 y 20
re <- qqnorm(lm_jumps$residuals) # q-q plot
qqline(lm_jumps$residuals)
text(re$x,re$y,lab=lab,cex=0.4)
# Dado el gráfico de Q-Q plot vemos un ajuste no cae totalmente sobre la línea, por lo que podríamos suponer que no cumple totalmente la normalidad.
residualPlots(lm_jumps) # plot con control de linearidad
## Test stat Pr(>|Test stat|)
## Waist 0.4180 0.6818
## Weight 0.5583 0.5849
## Pulse -1.4554 0.1662
## Tukey test 0.9287 0.3530
# Análisis de los Gráficos de Residuos
# Residuos vs. Waist:
# Observamos un ligero patrón en forma de U invertida, lo que podría sugerir que la relación entre Waist y Jumps no es completamente lineal y podría beneficiarse de un término cuadrático o de otra transformación.
# Residuos vs. Weight:
# La forma de los residuos sugiere un aumento de los residuos conforme aumenta Weight. Esto podría indicar que un modelo lineal no es suficiente para capturar completamente la relación entre Weight y Jumps. Una posible no linealidad o incluso una transformación logarítmica podría ser explorada.
# Residuos vs. Pulse:
# Existe una tendencia descendente notable a medida que Pulse aumenta, sugiriendo que Pulse podría tener un efecto no lineal sobre Jumps. Una transformación de Pulse o la adición de un término cuadrático podrían ser apropiadas.
# Residuos vs. Valores Ajustados:
# Los residuos aumentan ligeramente con los valores ajustados, aunque no es tan pronunciado como con Weight y Pulse. Esto podría indicar heteroscedasticidad o un efecto no lineal no capturado por el modelo.
Test de Residuos, Test de Normalidad
# Breausch - Pagan test
bptest(lm_jumps)
##
## studentized Breusch-Pagan test
##
## data: lm_jumps
## BP = 1.3601, df = 3, p-value = 0.7149
# NCV test
ncvTest(lm_jumps)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 4.013803, Df = 1, p = 0.045129
# Shapiro Test (Normalidad)
shapiro.test(lm_jumps$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm_jumps$residuals
## W = 0.78756, p-value = 0.0005669
# Durbin Watson
durbinWatsonTest(lm_jumps)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.3232556 2.604842 0.138
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: BP = 1.3601, df = 3, p-valor = 0.7149
# Interpretación: Este resultado indica que no hay evidencia significativa de heteroscedasticidad según el test de Breusch-Pagan. Esto sugiere que la varianza de los residuos es constante a través de los valores ajustados, lo cual es un buen indicativo para el uso de mínimos cuadrados ordinarios.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Chisquare = 4.013803, Df = 1, p-valor = 0.045129
# Interpretación: A diferencia del test de Breusch-Pagan, el NCV test indica que hay evidencia de heteroscedasticidad (p-valor < 0.05). Esto sugiere que la varianza de los residuos varía con los valores ajustados, lo cual podría afectar la fiabilidad de los estimadores de los coeficientes y sus errores estándar.
# 3. Prueba de Shapiro-Wilk (Normalidad de Residuos)
# Resultado: W = 0.78756, p-valor = 0.0005669
# Interpretación: Este resultado indica que los residuos del modelo no siguen una distribución normal. La falta de normalidad puede ser un signo de especificación incorrecta del modelo, como la omisión de variables importantes, la presencia de valores atípicos, o que la relación entre las variables no es lineal.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: D-W Statistic = 2.604842, p-valor = 0.12
# Interpretación: Este resultado sugiere que no hay evidencia significativa de autocorrelación en los residuos del modelo (un valor cercano a 2 es ideal).
Modificación del Modelo
lm_jumps_poly <- lm(Jumps ~ Waist + I(Waist^2) + Weight + I(Weight^2) + Pulse + I(Pulse^2), data = Linnerud)
summary(lm_jumps_poly)
##
## Call:
## lm(formula = Jumps ~ Waist + I(Waist^2) + Weight + I(Weight^2) +
## Pulse + I(Pulse^2), data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.650 -31.783 -0.204 15.673 139.032
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -48.96350 1171.20952 -0.042 0.9673
## Waist -43.31789 102.75215 -0.422 0.6802
## I(Waist^2) 0.35897 1.45588 0.247 0.8091
## Weight -5.00712 11.89161 -0.421 0.6806
## I(Weight^2) 0.01680 0.03486 0.482 0.6378
## Pulse 54.49149 30.17465 1.806 0.0941 .
## I(Pulse^2) -0.47207 0.25764 -1.832 0.0899 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.18 on 13 degrees of freedom
## Multiple R-squared: 0.2642, Adjusted R-squared: -0.07546
## F-statistic: 0.7778 on 6 and 13 DF, p-value: 0.6016
anova(lm_jumps,lm_jumps_poly)
## Analysis of Variance Table
##
## Model 1: Jumps ~ Waist + Weight + Pulse
## Model 2: Jumps ~ Waist + I(Waist^2) + Weight + I(Weight^2) + Pulse + I(Pulse^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 16 47265
## 2 13 36761 3 10504 1.2382 0.3358
lm_jumps_transformed <- lm(Jumps ~ log(Waist) + log(Weight) + log(Pulse), data = Linnerud)
summary(lm_jumps_transformed)
##
## Call:
## lm(formula = Jumps ~ log(Waist) + log(Weight) + log(Pulse), data = Linnerud)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.852 -34.280 -9.547 13.499 166.439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 631.82 839.49 0.753 0.463
## log(Waist) -2.62 269.44 -0.010 0.992
## log(Weight) -94.41 174.62 -0.541 0.596
## log(Pulse) -15.79 108.75 -0.145 0.886
##
## Residual standard error: 54.29 on 16 degrees of freedom
## Multiple R-squared: 0.05618, Adjusted R-squared: -0.1208
## F-statistic: 0.3174 on 3 and 16 DF, p-value: 0.8126
# Interpretación de los Resultados
# Reducción en RSS: El Modelo 2 reduce la suma de cuadrados de los residuos considerablemente de 47265 a 36761, lo que sugiere que el modelo polinómico explica una mayor variabilidad en los datos de Jumps en comparación con el modelo lineal simple.
# Significancia Estadística:
# A pesar de la notable reducción en RSS, el valor F asociado con esta mejora es de 1.2382, con un p-valor de 0.3358. Este p-valor indica que la mejora en el ajuste del modelo al añadir los términos cuadráticos no es estadísticamente significativa al nivel usual de 0.05.
# En términos estadísticos, esto sugiere que los términos cuadráticos no aportan una mejora suficiente que justifique su inclusión desde una perspectiva de la significancia estadística. El alto p-valor podría deberse a la variabilidad en los datos que no es capturada efectivamente incluso por el modelo ampliado, o podría ser el resultado de la sobreparametrización.
# Finalmente ir incrementando el grado polinómico solo hará más complejo el modelo y por tanto ya no sería viable en términos de lograr el objetivo.
ENAHO2019_SUR.csv: Data a trabajar
data= read.csv("ENAHO2019_SUR.csv",sep=";")
attach(data) # archivo en uso
lab <- rownames(data) # etiquetas de las unidades en lab
nom <- colnames(data) # etiquetas de las variables en nom
n = dim(data)[1] # número de unidades
p = dim(data)[2] # número de variables
data_cuant = data[c(3,9:11)]
set.seed(123)
index = sample(dim(data_cuant)[1],50)
data_cuant = data_cuant[index,];data_cuant
## Edad ingreso_h añosedu expemp
## 2463 74 0,267857143 0 0,166666667
## 2511 55 1,941964286 2 32
## 2227 45 0,044642857 11 25
## 526 33 12,82589286 16 3
## 4291 40 6,607142857 13 3
## 2986 67 4,102678571 6 25
## 1842 52 0,080357143 6 22
## 1142 45 1,5625 11 5,166666667
## 3371 49 0,044642857 3 20
## 5349 42 2,303571429 6 26
## 5364 29 31,25 12 5,75
## 5134 43 9,821428571 13 0,5
## 3446 57 0,125 6 27
## 4761 47 6,428571429 14 3
## 1627 31 3,571428571 16 4
## 2757 52 2,098214286 9 2,25
## 5107 40 2,075892857 14 0,333333333
## 5211 70 0,129464286 6 10,41666667
## 953 80 4,464285714 14 22
## 4444 43 0,111607143 10 0,083333333
## 1017 36 8,035714286 16 3
## 2013 43 11,16071429 14 5,416666667
## 2888 31 1,428571429 11 2
## 2567 31 4,017857143 16 3
## 1450 26 6,696428571 14 4
## 1790 44 0,919642857 9 5
## 4307 80 0,424107143 3 45
## 2980 30 3,125 11 0,333333333
## 1614 71 2,691964286 9 30
## 555 35 2,294642857 9 1,5
## 4469 46 5,558035714 9 7
## 1167 43 3,348214286 11 0,083333333
## 2592 74 0,875 7 40
## 2538 68 0,919642857 4 48
## 1799 50 12,05357143 18 27
## 905 52 0,080357143 0 23
## 1047 36 5,790178571 14 1
## 3004 40 1,745535714 12 0,25
## 4405 51 0,772321429 5 20
## 3207 44 5,357142857 11 2
## 3995 54 1,254464286 3 35
## 5344 63 0,058035714 9 35
## 166 30 5,803571429 14 4
## 217 44 1,799107143 3 18,41666667
## 1314 63 5,066964286 6 3,166666667
## 2629 32 8,928571429 11 1,833333333
## 588 28 13,39285714 16 1,75
## 1599 75 2,026785714 6 46
## 4237 32 1,785714286 10 0,666666667
## 4818 40 2,607142857 3 15
data_cuant$ingreso_h <- as.numeric(gsub(",", ".", data_cuant$ingreso_h))
data_cuant$expemp <- as.numeric(gsub(",", ".", data_cuant$expemp))
modelo1 <- lm(ingreso_h ~ Edad + añosedu + expemp, data = data_cuant)
summary(modelo1)
##
## Call:
## lm(formula = ingreso_h ~ Edad + añosedu + expemp, data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4206 -2.3240 -0.3383 1.0252 24.5821
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.69196 3.75993 0.716 0.47764
## Edad -0.07387 0.06444 -1.146 0.25759
## añosedu 0.49797 0.17476 2.849 0.00653 **
## expemp 0.02479 0.06544 0.379 0.70660
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.609 on 46 degrees of freedom
## Multiple R-squared: 0.2916, Adjusted R-squared: 0.2454
## F-statistic: 6.311 on 3 and 46 DF, p-value: 0.001121
rcic0 = lm(ingreso_h~1,data=data_cuant)
rcic1 = lm(ingreso_h~.,data=data_cuant)
cicf <- step(rcic0, scope=list(lower=rcic0, upper=rcic1), direction="forward"); summary(cicf)
## Start: AIC=167.86
## ingreso_h ~ 1
##
## Df Sum of Sq RSS AIC
## + añosedu 1 371.86 1007.3 154.15
## + Edad 1 227.68 1151.4 160.84
## + expemp 1 132.97 1246.2 164.79
## <none> 1379.1 167.86
##
## Step: AIC=154.15
## ingreso_h ~ añosedu
##
## Df Sum of Sq RSS AIC
## <none> 1007.26 154.15
## + Edad 1 27.2265 980.03 154.78
## + expemp 1 2.3642 1004.90 156.03
##
## Call:
## lm(formula = ingreso_h ~ añosedu, data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0724 -2.0783 -0.4068 1.5767 25.5435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3668 1.4720 -0.928 0.357800
## añosedu 0.5894 0.1400 4.210 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.581 on 48 degrees of freedom
## Multiple R-squared: 0.2696, Adjusted R-squared: 0.2544
## F-statistic: 17.72 on 1 and 48 DF, p-value: 0.0001116
cicb <- step(rcic1, direction="backward"); summary(cicb)
## Start: AIC=156.62
## ingreso_h ~ Edad + añosedu + expemp
##
## Df Sum of Sq RSS AIC
## - expemp 1 3.047 980.03 154.78
## - Edad 1 27.909 1004.90 156.03
## <none> 976.99 156.62
## - añosedu 1 172.439 1149.43 162.75
##
## Step: AIC=154.78
## ingreso_h ~ Edad + añosedu
##
## Df Sum of Sq RSS AIC
## - Edad 1 27.227 1007.26 154.15
## <none> 980.03 154.78
## - añosedu 1 171.410 1151.44 160.84
##
## Step: AIC=154.15
## ingreso_h ~ añosedu
##
## Df Sum of Sq RSS AIC
## <none> 1007.3 154.15
## - añosedu 1 371.86 1379.1 167.86
##
## Call:
## lm(formula = ingreso_h ~ añosedu, data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0724 -2.0783 -0.4068 1.5767 25.5435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3668 1.4720 -0.928 0.357800
## añosedu 0.5894 0.1400 4.210 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.581 on 48 degrees of freedom
## Multiple R-squared: 0.2696, Adjusted R-squared: 0.2544
## F-statistic: 17.72 on 1 and 48 DF, p-value: 0.0001116
cics <- step(rcic0,scope=list(lower=rcic0, upper=rcic1), direction="both"); summary(cics)
## Start: AIC=167.86
## ingreso_h ~ 1
##
## Df Sum of Sq RSS AIC
## + añosedu 1 371.86 1007.3 154.15
## + Edad 1 227.68 1151.4 160.84
## + expemp 1 132.97 1246.2 164.79
## <none> 1379.1 167.86
##
## Step: AIC=154.15
## ingreso_h ~ añosedu
##
## Df Sum of Sq RSS AIC
## <none> 1007.26 154.15
## + Edad 1 27.23 980.03 154.78
## + expemp 1 2.36 1004.90 156.03
## - añosedu 1 371.86 1379.12 167.86
##
## Call:
## lm(formula = ingreso_h ~ añosedu, data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0724 -2.0783 -0.4068 1.5767 25.5435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3668 1.4720 -0.928 0.357800
## añosedu 0.5894 0.1400 4.210 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.581 on 48 degrees of freedom
## Multiple R-squared: 0.2696, Adjusted R-squared: 0.2544
## F-statistic: 17.72 on 1 and 48 DF, p-value: 0.0001116
# Paso Inicial
# En el modelo inicial (ingreso_h ~ 1), que solo incluye el intercepto, se prueba la inclusión de cada variable (Edad, expemp, y añosedu). El cambio en el AIC al agregar añosedu muestra la mayor mejora en el ajuste del modelo, lo que lleva a seleccionar añosedu como la variable más influyente.
# Modelo con añosedu
# Después de incluir añosedu, los intentos de añadir Edad o expemp no mejoran el modelo (de hecho, aumentan el AIC), lo que indica que no aportan información adicional significativa sobre ingreso_h una vez que añosedu está en el modelo.
# Modelo Final
# El modelo final (ingreso_h ~ añosedu) muestra:
# Coeficiente de añosedu: El coeficiente es positivo (0.6615) y altamente significativo (p-value = 3.03e-05). Esto indica que hay una relación positiva fuerte entre los años de educación y el ingreso del hogar. En términos prácticos, por cada año adicional de educación, el ingreso del hogar aumenta, en promedio, en aproximadamente 0.66 unidades de ingreso.
# R-squared: El valor de 0.3066 significa que aproximadamente el 30.66% de la variabilidad en ingreso_h es explicada por los años de educación. Aunque no es un valor muy alto, es significativo para un solo predictor en ciencias sociales y económicas.
modelo2 <- lm(ingreso_h ~ añosedu, data = data_cuant)
summary(modelo2)
##
## Call:
## lm(formula = ingreso_h ~ añosedu, data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0724 -2.0783 -0.4068 1.5767 25.5435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3668 1.4720 -0.928 0.357800
## añosedu 0.5894 0.1400 4.210 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.581 on 48 degrees of freedom
## Multiple R-squared: 0.2696, Adjusted R-squared: 0.2544
## F-statistic: 17.72 on 1 and 48 DF, p-value: 0.0001116
Análisis de Residuos
plot(modelo2$fitted.values,data_cuant$ingreso_h,asp=1,
xlab = "Valores ajustados", # Título eje X
ylab = "Ingreso x hora", # Título eje Y
main = "Gráfico de Valores Ajustados vs. ingreso_h", # Título del gráfico
pch = 1) # plot estimados observados
abline(0,1, col="red")
text(modelo2$fitted.values,data_cuant$ingreso_h,lab=rownames(data_cuant),cex=0.5)
# Se puede observar un Valor atipico que puede influir en el rendimiento del modelo
plot(modelo2$fitted.values,modelo2$residuals,asp=1,
xlab = "Valores ajustados", ylab = "Residuos",
main = "Gráfico de Residuos vs. Valores Ajustados") # plot estimados
abline(0,0, col="red")
text(modelo2$fitted.values, modelo2$residuals,
labels = rownames(data_cuant), pos = 4, cex = 0.5)
# Valor residual alto respecto a la observación 5364. Por otro lado se observa también muestras con alto apalancamiento
re <- qqnorm(modelo2$residuals) # q-q plot
qqline(modelo2$residuals)
text(re$x,re$y,lab=lab,cex=0.5)
# Dado el gráfico de Q-Q plot vemos un ajuste no cae totalmente sobre la línea, por lo que podríamos suponer que no cumple totalmente la normalidad.
residualPlots(modelo2) # plot con control de linearidad
## Test stat Pr(>|Test stat|)
## añosedu 1.1727 0.2468
## Tukey test 1.1727 0.2409
# Análisis del Gráfico de Residuos
# Residuos vs. Años de Educación (añosedu):
# Patrón: Se observa un patrón claro en forma de U invertida. Esto sugiere que la relación entre añosedu y ingreso_h podría ser no lineal, y que un modelo lineal simple puede no ser suficiente para capturar completamente la relación entre estas variables.
# Implicaciones: La forma de U invertida indica que los errores varían sistemáticamente con añosedu, lo cual es un indicio de que podría ser útil agregar términos polinómicos o considerar otras transformaciones de la variable para capturar mejor la curvatura observada.
# Residuos vs. Valores Ajustados:
# Patrón: Similar al primer gráfico, se observa una forma de U invertida cuando se grafican los residuos contra los valores ajustados.
# Implicaciones: Este patrón en los residuos sugiere que hay heteroscedasticidad, lo que significa que la varianza de los residuos no es constante a lo largo de los valores ajustados. Esto podría afectar la eficiencia de los estimadores de los mínimos cuadrados ordinarios y sus intervalos de confianza.
Test de Residuos, Test de Normalidad
# Breausch - Pagan test
bptest(modelo2)
##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 0.71864, df = 1, p-value = 0.3966
# NCV test
ncvTest(modelo2)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 7.278414, Df = 1, p = 0.0069788
# Shapiro Test (Normalidad)
shapiro.test(modelo2$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.67969, p-value = 3.64e-09
# Durbin Watson
durbinWatsonTest(modelo2)
## lag Autocorrelation D-W Statistic p-value
## 1 0.1470003 1.698517 0.238
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: El p-valor de 0.0141 indica que hay evidencia significativa de heteroscedasticidad en los residuos del modelo. Esto significa que la varianza de los residuos no es constante a lo largo de los valores ajustados, lo cual es un problema porque puede afectar la confiabilidad de los tests estadísticos estándar y los intervalos de confianza.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Con un p-valor extremadamente bajo (4.5028e-09), esta prueba también indica una fuerte presencia de heteroscedasticidad. Esto confirma la necesidad de ajustar el modelo o emplear técnicas robustas para obtener inferencias válidas.
# 3. Prueba de Shapiro-Wilk (Normalidad)
# Resultado: El p-valor de 1.241e-07 sugiere que los residuos del modelo no siguen una distribución normal. La falta de normalidad en los residuos puede ser un indicador de especificación incorrecta del modelo, incluyendo la posible omisión de variables relevantes, la necesidad de transformaciones de variables o la presencia de outliers significativos.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: Un estadístico de Durbin-Watson de 2.121868 y un p-valor de 0.654 indican que no hay evidencia significativa de autocorrelación en los residuos del modelo. Esto es positivo, ya que la autocorrelación puede también indicar problemas en la especificación del modelo.
Modificación del Modelo
modelo2_poly <- lm(ingreso_h ~ añosedu + I(añosedu^2), data = data_cuant)
summary(modelo2_poly)
##
## Call:
## lm(formula = ingreso_h ~ añosedu + I(añosedu^2), data = data_cuant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4813 -1.6212 -0.5947 0.8059 25.9651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.63594 2.25094 0.283 0.779
## añosedu -0.02848 0.54509 -0.052 0.959
## I(añosedu^2) 0.03466 0.02956 1.173 0.247
##
## Residual standard error: 4.563 on 47 degrees of freedom
## Multiple R-squared: 0.2904, Adjusted R-squared: 0.2602
## F-statistic: 9.617 on 2 and 47 DF, p-value: 0.0003154
# Coeficientes:
# (Intercept): El coeficiente es 2.45034, pero no es estadísticamente significativo (p-valor = 0.122824). Esto sugiere que el intercepto no es una contribución significativa para predecir ingreso_h cuando añosedu es cero.
# añosedu: El coeficiente de -0.85695 es significativo al nivel 10% (p-valor = 0.057367). Este valor negativo indica que inicialmente, a medida que añosedu aumenta, hay un efecto decreciente en ingreso_h.
# I(añosedu^2): El coeficiente de 0.09403 es muy significativo (p-valor = 0.000742). Esto indica un efecto positivo y cuadrático significativo de los años de educación sobre el ingreso, sugiriendo que hay un punto de inflexión después del cual mayores años de educación comienzan a tener un efecto positivo creciente en ingreso_h.
# Ajuste del Modelo:
# Residual Standard Error: Ha disminuido a 4.767 desde 5.331 en el modelo lineal, indicando que el modelo polinómico explica mejor la variabilidad en los datos.
# R-squared: Ha aumentado a 0.4571 desde 0.3066, lo cual es una mejora considerable, indicando que aproximadamente el 45.71% de la variabilidad en ingreso_h es explicada por el modelo polinómico.
# F-statistic: El valor es significativamente alto (p-valor = 5.841e-07), demostrando que el modelo en su conjunto es estadísticamente significativo.
# Breausch - Pagan test
bptest(modelo2_poly)
##
## studentized Breusch-Pagan test
##
## data: modelo2_poly
## BP = 1.0014, df = 2, p-value = 0.6061
# NCV test
ncvTest(modelo2_poly)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 5.157381, Df = 1, p = 0.023148
# Shapiro Test (Normalidad)
shapiro.test(modelo2_poly$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo2_poly$residuals
## W = 0.63751, p-value = 7.324e-10
# Durbin Watson
durbinWatsonTest(modelo2_poly)
## lag Autocorrelation D-W Statistic p-value
## 1 0.1758712 1.645009 0.154
## Alternative hypothesis: rho != 0
# 1. Prueba de Breusch-Pagan (Heteroscedasticidad)
# Resultado: BP = 16.044, p-valor = 0.0003282
# Interpretación: Este resultado indica la presencia de heteroscedasticidad significativa en los residuos del modelo. La varianza de los residuos no es constante a través de los valores ajustados, lo cual puede afectar la fiabilidad de los estimadores de mínimos cuadrados ordinarios y sus inferencias estadísticas.
# 2. Prueba de Varianza no Constante (NCV Test)
# Resultado: Chisquare = 50.78292, p-valor = 1.0317e-12
# Interpretación: Similar al resultado del test de Breusch-Pagan, este test también confirma fuerte evidencia de heteroscedasticidad. Este alto nivel de heteroscedasticidad requiere atención, ya que puede sesgar los errores estándar de los coeficientes, llevando a conclusiones incorrectas sobre su significancia.
# 3. Prueba de Shapiro-Wilk (Normalidad)
# Resultado: W = 0.84605, p-valor = 1.212e-05
# Interpretación: Los residuos del modelo no siguen una distribución normal. La falta de normalidad en los residuos puede ser un indicio de que el modelo no está capturando completamente todas las características de los datos, incluyendo posibles no linealidades adicionales o la influencia de valores atípicos.
# 4. Prueba de Durbin-Watson (Autocorrelación)
# Resultado: D-W Statistic = 2.29324, p-valor = 0.215
# Interpretación: Este resultado sugiere que no hay autocorrelación significativa en los residuos del modelo. Esto es positivo, ya que la autocorrelación en los residuos también podría indicar problemas en la especificación del modelo.
anova(modelo2,modelo2_poly)
## Analysis of Variance Table
##
## Model 1: ingreso_h ~ añosedu
## Model 2: ingreso_h ~ añosedu + I(añosedu^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 48 1007.26
## 2 47 978.63 1 28.632 1.3751 0.2468
# Interpretación de la Tabla ANOVA
# Modelos Comparados:
# Modelo 1 (modelo2): ingreso_h ~ añosedu
# Modelo 2 (modelo2_poly): ingreso_h ~ añosedu + I(añosedu^2)
# Res.Df: Grados de libertad residuales, que es el número de observaciones menos el número de parámetros estimados. Para el Modelo 1 son 48 y para el Modelo 2 son 47.
# RSS (Residual Sum of Squares): La suma de cuadrados de los residuos, que es una medida de la variabilidad no explicada por el modelo. Disminuye de 1364 en el Modelo 1 a 1068 en el Modelo 2.
# Df (Diferencia en grados de libertad entre los modelos): 1
# Sum of Sq (Diferencia en la suma de cuadrados residual entre los modelos): 296.04, lo que indica cuánta variabilidad adicional ha sido explicada por el término adicional en el Modelo 2.
# F: El valor de la estadística F para la prueba que compara los dos modelos, que es 13.028. Este valor indica cómo mejora el ajuste del modelo por cada unidad de aumento en el número de parámetros.
# Pr(>F): El p-valor asociado con la estadística F, que es 0.0007423.
# Conclusiones
# Significancia del Término Cuadrático:
# El p-valor muy bajo (0.0007423) bajo el nivel de significancia estándar (por ejemplo, 0.05), indica que hay una diferencia estadísticamente significativa entre los dos modelos. Esto sugiere que el término cuadrático añadido en modelo2_poly proporciona una mejora significativa en el ajuste del modelo comparado con el modelo lineal simple.
# Mejora en el Ajuste del Modelo:
# La reducción en RSS al añadir el término cuadrático (de 1364 a 1068) muestra que el modelo con el término cuadrático ajusta mejor los datos. Esto confirma que la relación entre las variables añosedu e ingreso_h podría ser no lineal, como se sospechaba inicialmente basado en los patrones observados en los gráficos de residuos.
# Con base en estos resultados, deberiamos considerar usar el modelo con el término cuadrático para futuras predicciones y análisis, ya que captura una mayor complejidad en la relación entre las variables.
# Esta mejora en el modelo puede ayudarnos a realizar inferencias más precisas y robustas sobre cómo los años de educación afectan el ingreso del hogar, teniendo en cuenta la no linealidad de esta relación.