Una compañía observa que la demanda de uno de sus productos cambió debido a una variación rápida de su precio por unidad. Se observa la demanda del producto (unidades por cada mil habitantes) en una región en particular sobre un intervalo bastante amplio de precios de venta (dolares). Los datos que se encuentran en el archivo company.dat.
Se trata de realizar un análisis de dichos datos, indicando una ecuación que marque la relación entre el precio del producto por unidad y la demanda del mismo.
Vamos a importar y chequear los datos.
company<-read.table("./datos/company.dat",head=T)
colnames(company)
## [1] "unidades" "dolares"
Ajustamos el modelo de regresión lineal de las ventas del producto, \(unidades\), frente al precio, \(dolares\):
\(unidades=\beta_0+\beta_1 dolares +\varepsilon\)
ml1<-lm(unidades~dolares,data=company)
summary(ml1)
##
## Call:
## lm(formula = unidades ~ dolares, data = company)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53.827 -36.979 -3.641 31.720 77.730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 497.156 60.852 8.170 1.87e-05 ***
## dolares -24.419 4.594 -5.315 0.000484 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.53 on 9 degrees of freedom
## Multiple R-squared: 0.7584, Adjusted R-squared: 0.7315
## F-statistic: 28.25 on 1 and 9 DF, p-value: 0.0004841
En este caso, podemos dibujar los datos y la recta ajustada:
plot(unidades~dolares, data=company)
abline(ml1)
Vamos a realizar un diagnóstico del modelo, particularmente la linealidad mediante los gráficos de diagnóstico:
layout(matrix(1:4,2))
plot(ml1)
Comentarios:
Introducimos en primer lugar el término de grado 2, ajustando el modelo:
\(unidades=\beta_0+\beta_1 dolares +\beta_2 dolares^2 +\varepsilon\)
ml2<-lm(unidades~dolares+I(dolares^2),data=company)
summary(ml2)
##
## Call:
## lm(formula = unidades ~ dolares + I(dolares^2), data = company)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.213 -10.067 -5.116 22.734 33.192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1330.407 179.565 7.409 7.55e-05 ***
## dolares -155.467 27.868 -5.579 0.000523 ***
## I(dolares^2) 4.866 1.031 4.722 0.001499 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.91 on 8 degrees of freedom
## Multiple R-squared: 0.9362, Adjusted R-squared: 0.9202
## F-statistic: 58.69 on 2 and 8 DF, p-value: 1.657e-05
Comentarios:
Comprobamos la significación del término de grado 3, ajustando el modelo cúbico:
\(unidades=\beta_0+\beta_1 dolares +\beta_2 dolares^2 + \beta_3 dolares^3 +\varepsilon\)
ml3<-lm(unidades~dolares+I(dolares^2)+I(dolares^3),data=company)
summary(ml3)
##
## Call:
## lm(formula = unidades ~ dolares + I(dolares^2) + I(dolares^3),
## data = company)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.019 -14.183 3.279 12.058 31.413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2555.7171 913.5503 2.798 0.0266 *
## dolares -444.3646 213.2536 -2.084 0.0757 .
## I(dolares^2) 26.8745 16.1498 1.664 0.1400
## I(dolares^3) -0.5426 0.3974 -1.365 0.2144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.61 on 7 degrees of freedom
## Multiple R-squared: 0.9496, Adjusted R-squared: 0.928
## F-statistic: 43.97 on 3 and 7 DF, p-value: 6.553e-05
Comentarios:
Antes de dar por bueno el modelo cuadrático, realicemos el análisis de residuos:
#gráficos de diagnóstico
layout(matrix(1:4,2))
plot(ml2)
#normalidad
shapiro.test(rstandard(ml2))
##
## Shapiro-Wilk normality test
##
## data: rstandard(ml2)
## W = 0.91339, p-value = 0.2674
#homocedasticidad
library(lmtest)
bptest(ml2)
##
## studentized Breusch-Pagan test
##
## data: ml2
## BP = 1.5018, df = 2, p-value = 0.472
#errores incorrelados
library(lmtest)
dwtest(ml2)
##
## Durbin-Watson test
##
## data: ml2
## DW = 2.498, p-value = 0.5372
## alternative hypothesis: true autocorrelation is greater than 0
Comentarios:
summary(ml2)
##
## Call:
## lm(formula = unidades ~ dolares + I(dolares^2), data = company)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.213 -10.067 -5.116 22.734 33.192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1330.407 179.565 7.409 7.55e-05 ***
## dolares -155.467 27.868 -5.579 0.000523 ***
## I(dolares^2) 4.866 1.031 4.722 0.001499 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.91 on 8 degrees of freedom
## Multiple R-squared: 0.9362, Adjusted R-squared: 0.9202
## F-statistic: 58.69 on 2 and 8 DF, p-value: 1.657e-05
Comentarios:
Vamos a estimar la demanda para un precio del producto de 14 dólares.
predict(ml2,data.frame(dolares=14), interval="confidence")
## fit lwr upr
## 1 107.6297 77.48593 137.7735
predict(ml2,data.frame(dolares=14), interval="prediction")
## fit lwr upr
## 1 107.6297 40.71657 174.5429
Las puntuaciones de cierto examen pueden depender de las horas de estudio y las horas de ejercicio físico previas al examen. Para establecer una ecuación que relacione estas variables, se realiza un experimento con 100 estudiantes, que registran cuántas horas de estudio y ejercicio tuvieron en los días previos al examen. Posteriormente se registró también la puntuación del examen. Los datos se encuentran en el archivo examen.dat.
Vamos a importar y chequear los datos.
examen<-read.table("./datos/examen.dat",head=T)
colnames(examen)
## [1] "Horas_Estudio" "Horas_Ejercicio" "Puntuacion"
Ajustamos el modelo
\[Puntuacion=\beta_0 + \beta_1 \mbox{Horas_Estudio} + \beta_2 \mbox{Horas_Ejercicio} + \varepsilon\]
ml4 <- lm(Puntuacion ~ Horas_Estudio+Horas_Ejercicio, data = examen)
summary(ml4)
##
## Call:
## lm(formula = Puntuacion ~ Horas_Estudio + Horas_Ejercicio, data = examen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.4923 -4.0776 -0.2336 3.8730 12.4726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.2730 1.5833 12.804 <2e-16 ***
## Horas_Estudio 5.2482 0.1727 30.385 <2e-16 ***
## Horas_Ejercicio 0.5062 0.2328 2.175 0.0321 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.978 on 97 degrees of freedom
## Multiple R-squared: 0.905, Adjusted R-squared: 0.9031
## F-statistic: 462.1 on 2 and 97 DF, p-value: < 2.2e-16
Dibujemos los gráficos de residuos:
layout(matrix(1:4,2))
plot(ml4)
Comentarios:
Ajustamos el modelo con los términos de grado 2:
\[ Puntuacion=\beta_0+\beta_1 \mbox{Horas_Estudio}+\beta_2 \mbox{Horas_Ejercicio}+\beta_3 \mbox{Horas_Estudio}^2+\beta_4 \mbox{Horas_Ejercicio}^2+\beta_5 \mbox{Horas_Estudio}*\mbox{Horas_Ejercicio}+\varepsilon \]
ml5 <- lm(Puntuacion ~ Horas_Estudio+Horas_Ejercicio+I(Horas_Estudio^2)+I(Horas_Ejercicio^2)+I(Horas_Estudio*Horas_Ejercicio), data = examen)
Veamos si los nuevos términos introducidos son globalmente significativos:
#significación de los términos introducidos
anova(ml4,ml5)
## Analysis of Variance Table
##
## Model 1: Puntuacion ~ Horas_Estudio + Horas_Ejercicio
## Model 2: Puntuacion ~ Horas_Estudio + Horas_Ejercicio + I(Horas_Estudio^2) +
## I(Horas_Ejercicio^2) + I(Horas_Estudio * Horas_Ejercicio)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 97 3466.9
## 2 94 1651.9 3 1815 34.427 4.203e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Comentarios:
Realicemos un análisis del modelo:
summary(ml5)
##
## Call:
## lm(formula = Puntuacion ~ Horas_Estudio + Horas_Ejercicio + I(Horas_Estudio^2) +
## I(Horas_Ejercicio^2) + I(Horas_Estudio * Horas_Ejercicio),
## data = examen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9064 -2.6521 -0.1569 3.1261 8.3559
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.00290 2.17453 15.637 < 2e-16 ***
## Horas_Estudio 0.52992 0.51552 1.028 0.3066
## Horas_Ejercicio -1.32109 0.73158 -1.806 0.0741 .
## I(Horas_Estudio^2) 0.35013 0.03601 9.723 7.14e-16 ***
## I(Horas_Ejercicio^2) 0.20388 0.07296 2.795 0.0063 **
## I(Horas_Estudio * Horas_Ejercicio) 0.02464 0.04771 0.517 0.6067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.192 on 94 degrees of freedom
## Multiple R-squared: 0.9547, Adjusted R-squared: 0.9523
## F-statistic: 396.6 on 5 and 94 DF, p-value: < 2.2e-16
Comentarios:
Uno de los términos de grado 2 no es significativo, ajustemos el modelo sin dicho término y hagamos su análisis.
#ajuste del modelo de grado 2 sin el término producto
ml6 <- lm(Puntuacion ~ Horas_Estudio+Horas_Ejercicio+I(Horas_Estudio^2)+I(Horas_Ejercicio^2), data = examen)
summary(ml6)
##
## Call:
## lm(formula = Puntuacion ~ Horas_Estudio + Horas_Ejercicio + I(Horas_Estudio^2) +
## I(Horas_Ejercicio^2), data = examen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7368 -2.7177 -0.1745 3.1674 8.4390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.33447 1.74082 19.149 < 2e-16 ***
## Horas_Estudio 0.61431 0.48706 1.261 0.21031
## Horas_Ejercicio -1.11254 0.60772 -1.831 0.07029 .
## I(Horas_Estudio^2) 0.35112 0.03582 9.802 4.37e-16 ***
## I(Horas_Ejercicio^2) 0.19784 0.07174 2.758 0.00698 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.176 on 95 degrees of freedom
## Multiple R-squared: 0.9546, Adjusted R-squared: 0.9527
## F-statistic: 499.5 on 4 and 95 DF, p-value: < 2.2e-16
Comentarios:
Realicemos un diagnóstico de este último modelo
#Diagnóstico del modelo
layout(matrix(1:4,2))
plot(ml6)
Comentarios:
#predicciones
#Estimar la puntuación de un alumno con 7 horas de estudio y 5 horas de ejercicio
predict(ml6,data.frame(Horas_Estudio=7, Horas_Ejercicio=5), interval="confidence")
## fit lwr upr
## 1 54.22295 52.73281 55.71309
predict(ml6,data.frame(Horas_Estudio=7, Horas_Ejercicio=5), interval="prediction")
## fit lwr upr
## 1 54.22295 45.79994 62.64596
En ciertas empresas se realizan pruebas de selección para conocer la aptitud para el trabajo de los aspirantes a obtener un empleo. Resulta por tanto clave saber qué pruebas pueden predecir la aptitud para el trabajo de una persona. En el archivo TRABAJO.DAT se encuentran las puntuaciones de un test psicotécnico (psico) y de un test de adecuación (adec) que obtuvieron 20 individuos que fueron contratados por la compañía. También aparece el resultado de la evaluación a que fueron sometidos estos 20 empleados tras un periodo de dos años (eval). Analiza estos datos mediante el modelo teórico adecuado determinando si existe relación entre la evaluación y el resultado de los dos test.
Se pretende estudiar la relación existente entre la producción de sulfato y de nitrato en determinadas reacciones químicas. Para ello se realizan 20 de estas reacciones midiendo las cantidades de dichos compuestos producidas. Los datos registrados se encuentran en el archivo nitratosulfato.dat. Analiza estos datos utilizando para ello el modelo de regresión adecuado.