En el conjunto de datos CLOTHING tenemos información sobre ventas, tamaño y otras características de 400 tiendas de moda masculina en Holanda. El objetivo es explicar las ventas por metro cuadrado (sales) a partir de las características de la tienda (número de propietarios, trabajadores a tiempo completo y tiempo parcial, número de horas trabajadas, tamaño de la tienda, etc.).
1. Estima un modelo lineal (modelo A) que explique sales a partir del número total de horas trabajadas (hoursw), el tamaño de la tienda en metros cuadrados (ssize) y una constante. Interpreta los resultados.
modA <- lm(sales ~ hoursw + ssize)
modA.s <- summary(modA)
modA.s
Call:
lm(formula = sales ~ hoursw + ssize)
Residuals:
Min 1Q Median 3Q Max
-5658.5 -1874.1 -466.8 937.6 18076.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5133.590 321.693 15.96 <2e-16 ***
hoursw 37.528 2.837 13.23 <2e-16 ***
ssize -22.145 1.625 -13.63 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2985 on 397 degrees of freedom
Multiple R-squared: 0.3658, Adjusted R-squared: 0.3626
F-statistic: 114.5 on 2 and 397 DF, p-value: < 2.2e-16Se cumple la significacion global del modelo con un p-valor de < 2.2e-16 y también la significación individual para cada una de las variables del modelo, no hay ninguna variable redundante en el modelo.
El R2 tiene un valor de 0.3658 por lo que las variables hoursw y ssize explican solamente el 36.58% de la variablidad total de sales.
El coeficiente de hoursw de 37.528 nos indica que por cada hora extra total trabajada aumenta en 37.528 el valor de sales.
El coeficiente de ssize de -22.145 nos indica que el valor de sales disminuye 22.145 cada que aumenta una unidad ssize, aunque a priori puede parecer antintuitivo esto tiene sentido ya que sales representa las ventas por metro cuadrado y no las ventas totales.
Los valores minimos de hoursw y ssize son 32 y 16 respectivamente por lo que la estimación de la constante no tiene una interpretación interesante en este contexto.
2. Realiza un contraste RESET con Q=2.
FUNCIÓN RESET
reset(modA,2)
RESET test
data: modA
RESET = 2.0913, df1 = 1, df2 = 396, p-value = 0.1489PROGRAMANDO
yhat1 <- modA$fit
yhat2 <- yhat1^2
modaux <- lm(sales ~ hoursw + ssize + yhat2)
summary(modaux)
Call:
lm(formula = sales ~ hoursw + ssize + yhat2)
Residuals:
Min 1Q Median 3Q Max
-8672.3 -1791.4 -533.8 915.1 18141.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.025e+03 3.299e+02 15.235 < 2e-16 ***
hoursw 2.938e+01 6.307e+00 4.659 4.35e-06 ***
ssize -1.864e+01 2.915e+00 -6.395 4.53e-10 ***
yhat2 1.252e-05 8.658e-06 1.446 0.149
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2981 on 396 degrees of freedom
Multiple R-squared: 0.3691, Adjusted R-squared: 0.3644
F-statistic: 77.24 on 3 and 396 DF, p-value: < 2.2e-16Al haber unicamente un término que agregar con Q=2, basta con fijarse en el summary y quedarnos con el p-valor de yhat2 que es de 0.149, este coincidira con el p-valor del test RESET.
Haciendo el contraste formal:
ssr <- crossprod(modaux$resid)
ssrr <- crossprod(modA$resid)
num <- (ssrr-ssr)/1
den <- ssr/modaux$df
fstat <- num/den
pvalor <- pf(fstat,1,modaux$df, lower.tail=F)
pvalor [,1]
[1,] 0.1489321De cualquier manera el p-valor obtenido en los tres casos es de 0.1489 y por tanto con un nivel de significación del 0.05 no podemos rechazar la H0 de que las combinaciones no lineales de los valores ajustados del modelo hasta el grado 2 no ayudan a explicar la variable sales.
No tenemos indicios de problemas de especificación con el contraste RESET con Q=2.
3. Comprueba si el número de propietarios (nown) afecta a las ventas de la tienda, condicional a los niveles de hoursw y ssize.
modA2 <- lm(sales ~ hoursw + ssize + nown + nown:hoursw + nown:ssize)
modA2.s <- summary(modA2)
modA2.s
Call:
lm(formula = sales ~ hoursw + ssize + nown + nown:hoursw + nown:ssize)
Residuals:
Min 1Q Median 3Q Max
-5700 -1864 -479 1016 17767
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5231.140 829.274 6.308 7.60e-10 ***
hoursw 43.225 3.786 11.416 < 2e-16 ***
ssize -26.725 4.036 -6.622 1.16e-10 ***
nown -221.114 614.537 -0.360 0.719
hoursw:nown -1.544 1.106 -1.396 0.164
ssize:nown 2.265 2.332 0.972 0.332
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2975 on 394 degrees of freedom
Multiple R-squared: 0.3749, Adjusted R-squared: 0.3669
F-statistic: 47.26 on 5 and 394 DF, p-value: < 2.2e-16Vemos como el único término que se puede considerar significativo es la interacción de nown con hoursw, esto quiere decir que el número de propietarios o managers influye en la variable sales en función del número de horas totales trabajadas.
Si quitamos la interacción de nown con ssize podemos ver como el p-valor de hoursw:nown mejora aun más.
modA2 <- lm(sales ~ hoursw + ssize + nown + nown:hoursw)
modA2.s <- summary(modA2)
modA2.s
Call:
lm(formula = sales ~ hoursw + ssize + nown + nown:hoursw)
Residuals:
Min 1Q Median 3Q Max
-5884.8 -1875.7 -467.1 936.3 17849.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4685.941 610.569 7.675 1.31e-13 ***
hoursw 42.914 3.773 11.375 < 2e-16 ***
ssize -23.162 1.685 -13.746 < 2e-16 ***
nown 195.728 439.976 0.445 0.657
hoursw:nown -1.757 1.085 -1.619 0.106
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2975 on 395 degrees of freedom
Multiple R-squared: 0.3734, Adjusted R-squared: 0.367
F-statistic: 58.84 on 4 and 395 DF, p-value: < 2.2e-164. Comprueba también si la inclusión del número de trabajadores a tiempo parcial (npart) mejora el modelo.
modA3 <- lm(sales ~ hoursw + ssize + npart)
modA3.s <- summary(modA3)
modA3.s
Call:
lm(formula = sales ~ hoursw + ssize + npart)
Residuals:
Min 1Q Median 3Q Max
-6194.1 -1929.8 -476.4 1080.8 18239.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4175.017 411.993 10.13 < 2e-16 ***
hoursw 37.019 2.798 13.23 < 2e-16 ***
ssize -23.855 1.668 -14.30 < 2e-16 ***
npart 816.719 224.372 3.64 0.000309 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2940 on 396 degrees of freedom
Multiple R-squared: 0.3863, Adjusted R-squared: 0.3817
F-statistic: 83.1 on 3 and 396 DF, p-value: < 2.2e-16Al incluir npart en el modelo el resto de variables siguen siendo significativas. npart tiene un p-valor de 0.000309 con lo que también es significativa en el modelo y lo mejora. El R2 ajustado vemos que aumenta de 0.3626 a 0.3817 por lo que al incluir npart se explica mejor la variabilidad de sales.
5. Estima un modelo lineal (modelo B), que explique sales en función del número de trabajadores a tiempo completo (nfull), el número de trabajadores a tiempo parcial y la superficie de la tienda.
modB <- lm(sales ~ nfull + npart + ssize)
modB.s <- summary(modB)
modB.s
Call:
lm(formula = sales ~ nfull + npart + ssize)
Residuals:
Min 1Q Median 3Q Max
-6988.6 -1989.2 -318.8 1207.5 18793.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4924.568 468.928 10.502 < 2e-16 ***
nfull 1348.898 176.447 7.645 1.6e-13 ***
npart 605.046 255.581 2.367 0.0184 *
ssize -15.415 1.639 -9.405 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3296 on 396 degrees of freedom
Multiple R-squared: 0.2289, Adjusted R-squared: 0.223
F-statistic: 39.18 on 3 and 396 DF, p-value: < 2.2e-16Al igual que en el modelo A se cumple la significatividad conjunta e individual.
6. Compara los modelos A y B a partir del R2 ajustado.
Modelo A: 0.362609
Modelo B: 0.2230226El R2 ajustado del modelo A es bastante mayor que el del modelo B, si hubiese que elegir entre estos dos modelos sería preferible el modelo A ya que explica mayor parte de la variabilidad de sales.
7. Realiza un contraste F no anidado del modelo A frente al modelo B. ¿Llegas a alguna conclusión?
encomptest(modA, modB)Encompassing test
Model 1: sales ~ hoursw + ssize
Model 2: sales ~ nfull + npart + ssize
Model E: sales ~ hoursw + ssize + nfull + npart
Res.Df Df F Pr(>F)
M1 vs. ME 395 -2 12.018 8.58e-06 ***
M2 vs. ME 395 -1 114.516 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Lo que hace encomptest es comparar modelos no anidados, para ello primero genera un modelo que contiene todas las variables (Model E). Luego ese Model E lo compara con el modelo A (Model 1) y con el Modelo B (Model 2).
Vemos que el p-valor de la comparación de M1 con ME es de 8.58e-06 y por tanto existen diferencias significativas entre los dos modelos, en este caso quiere decir que nfull y npart son significativas en el Model E.
Lo mismo ocurre cuando comparamos M2 con ME, el p-valor es de < 2.2e-16 y por tanto hoursw es significativo en el Model E.
Como en ambos casos los p-valores son significativos, con este contraste no podemos sacar ninguna conclusión sobre que modelo es mejor entre el A o B.
Otra manera de obtener los p-valores:
modC <- lm(sales ~ hoursw + nfull + npart + ssize)
modC.s <- summary(modC)
modC.s
Call:
lm(formula = sales ~ hoursw + nfull + npart + ssize)
Residuals:
Min 1Q Median 3Q Max
-6372.2 -1847.0 -514.8 1235.6 17846.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3751.311 427.696 8.771 < 2e-16 ***
hoursw 32.765 3.062 10.701 < 2e-16 ***
nfull 557.308 172.248 3.236 0.00132 **
npart 684.998 225.443 3.038 0.00254 **
ssize -23.908 1.649 -14.502 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2906 on 395 degrees of freedom
Multiple R-squared: 0.4022, Adjusted R-squared: 0.3961
F-statistic: 66.43 on 4 and 395 DF, p-value: < 2.2e-16waldtest(modA, modC)Wald test
Model 1: sales ~ hoursw + ssize
Model 2: sales ~ hoursw + nfull + npart + ssize
Res.Df Df F Pr(>F)
1 397
2 395 2 12.018 8.58e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1waldtest(modB, modC)Wald test
Model 1: sales ~ nfull + npart + ssize
Model 2: sales ~ hoursw + nfull + npart + ssize
Res.Df Df F Pr(>F)
1 396
2 395 1 114.52 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 18. Repite el apartado anterior mediate el contraste J. ¿Mantienes la conclusión anterior?
jtest(modA, modB)J test
Model 1: sales ~ hoursw + ssize
Model 2: sales ~ nfull + npart + ssize
Estimate Std. Error t value Pr(>|t|)
M1 + fitted(M2) 0.51704 0.113017 4.5749 6.384e-06 ***
M2 + fitted(M1) 0.87307 0.081586 10.7012 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1En este caso el contraste J lo que hace es incluir a cada modelo los valores ajustados del otro modelo.
Para el Modelo A vemos que si le agregamos los valores ajustados del modelo B estos son significativos con un p-valor de 6.384e-06.
Lo mismo ocurre para el modelo B cuando agregamos los valores ajustados del modelo A, son significativos con un p-valor de < 2.2e-16.
Esto quiere decir que ninguno de los dos modelos incluye las variables regresoras adecuadas y por tanto puede haber algun problema de especificación.
Las conclusiones respecto a que modelo es mejor son las mismas que en el apartado anterior.
9. Al modelo A le añadimos las variables explicativas del número de trabajadores a tiempo parcial y a tiempo completo. El modelo resultante la llamamos modelo C. Estima el modelo, interpreta los resultados y realiza un contraste RESET. ¿ Estás satisfecho con esta especificación?
Call:
lm(formula = sales ~ hoursw + nfull + npart + ssize)
Residuals:
Min 1Q Median 3Q Max
-6372.2 -1847.0 -514.8 1235.6 17846.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3751.311 427.696 8.771 < 2e-16 ***
hoursw 32.765 3.062 10.701 < 2e-16 ***
nfull 557.308 172.248 3.236 0.00132 **
npart 684.998 225.443 3.038 0.00254 **
ssize -23.908 1.649 -14.502 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2906 on 395 degrees of freedom
Multiple R-squared: 0.4022, Adjusted R-squared: 0.3961
F-statistic: 66.43 on 4 and 395 DF, p-value: < 2.2e-16
RESET test
data: modC
RESET = 17.464, df1 = 2, df2 = 393, p-value = 5.418e-08Se puede ver que en este modelo todas las variables son significativas, se cumple la significatividad individual y conjunta.
El R2 ajustado tiene un valor de 0.3961 por lo que las variables independientes del modelo explican el 39.61% de la variablidad total de sales.
reset(modC)
RESET test
data: modC
RESET = 17.464, df1 = 2, df2 = 393, p-value = 5.418e-08Con un p-valor de 5.418e-08 rechazamos la H0 de que las combinaciones no lineales de los valores ajustados del modelo (de grado 2 y 3) no ayudan a explicar la variable sales.
Esto significa que existe algun problema de especificación en nuestro modelo, puede deberse tanto por la falta de inclusión de variables en el modelo o por la no linealidad de alguna de las variables.