Construcción de modelo lineal aplicado a los datos LandRent
Estudiantes Maestría en estadística: Univalle
8/12/2020
Introducción
El estudio a continuación es basado bajo el reporte económico realizado por los autores Hasbargen, Thomas, and Rolfes (1977), los cuales buscaban proporcionar una actualización en los patrones de alquiler de terrenos dedicados a la agricultura en el estado de Minnesota. Este conjunto de datos ha sido ampliamente trabajado en cursos de modelos lineales, utilizando el conjunto de variables que proporciona para realizar estimaciones sobre los mismos. En específico, partiremos del ejercicio planteado por (Weisberg 2005) , donde se interesa estudiar la relación entre el precio de terrenos destinados para agricultura y la densidad de vacas en la zona. El conjunto de datos landrent de la librería alr4 contiene información sobre la renta pagada en 1977 por las tierras agrícolas plantadas con alfalfa en 67 condados con tierras agrícolas alquiladas en Minnesota (Estados Unidos). Para esto los autores asumen la siguiente propuesta.
- Hipótesis de trabajo: Se cree que la renta de la tierra planteada con alfalfa en relación con la renta para otros fines agrícolas es mayor en áreas con una alta densidad de vacas lecheras (Weisberg 2005)
Aspectos importantes del problema
Es importante señalar que, dado el alto contenido proteínico, el alfalfa es un alimento adecuado para las vacas lechera. Esta alcanza a cumplir con cerca del \(90\%\) de los requerimientos nutricionales de la vaca lechera para la producción de leche.
Muestra del conjunto de datos
Método de recolección de datos
“Gran parte de los datos presentados en las siguientes paginas están basadas en una encuesta de 1979, aunque las comparaciones están hechas con los tarifas de alquiler determinadas en la encuesta de 1978 y también con las encuestas tempranas realizadas por la Universidad de Minnesota” (Hasbargen, Thomas, and Rolfes 1977)
Variables del conjunto de datos
\(x_{1}\): Renta promedio paga por toda la tierra cultivable (dolares)
\(x_{2}\): Densidad de vacas lecheras (número por milla cuadrada).
\(x_{3}\): proporción de tierra agrícola utilizada como pasto.
\(x_{4}\): Si el condado está limitado solo al cultivo de alfalfa: \(x4 = 1\) o no \(x4 = 0\)
\(y\): Renta promedio por acre sembrado de alfalfa (dolares).
Estadísticas descriptivas
| Rpromedio | Densidad | Proporcion | Respuesta | Limitacion | |
|---|---|---|---|---|---|
| Min. | 6,1700 | 1,5300 | 0,0200 | 5,0000 | |
| Q.1st | 24,4200 | 7,1100 | 0,0650 | 23,5000 | |
| Median | 44,5600 | 16,1200 | 0,1200 | 39,1700 | SI: 34 |
| Mean | 43,6399 | 20,5633 | 0,1697 | 42,1661 | |
| Q.3rd | 59,6800 | 31,2150 | 0,2350 | 57,0850 | NO: 33 |
| Max. | 83,9000 | 58,6000 | 0,7200 | 99,1700 | |
| CV% | 48% | 75% | 85% | 54% |
par(mfrow=c(2,2))
for (v in names(data)[c(1,2,3,5)]) {
boxplot(data[v][,1]~data$`condado limitado`,ylab = v,xlab = "condado limitado")
}
par(mfrow=c(2,3))
for (v in names(data)[c(1,2,3,5)]) {
hist(data[v][,1],xlab=v,ylab ="Frecuencia",main = "")
}
barplot(table(data$`condado limitado`),xlab ="condado limitado")
a<-data %>% ggplot(aes(x=`Renta promedio tierra`,y=`Renta promedio acre`,color=factor(`condado limitado`)))+
geom_point(size=2)+
#scale_color_brewer()+
labs(color="Limitado")+
theme_light()
b<-data %>% ggplot(aes(x=`densidad vacas lecheras`,y=`Renta promedio acre`,color=factor(`condado limitado`)))+
geom_point(size=2)+
#scale_color_brewer()+
labs(color="Limitado")+
theme_light()
c<-data %>% ggplot(aes(x=`proporcion de tierra pasto`,y=`Renta promedio acre`,color=factor(`condado limitado`)))+
geom_point(size=2)+
#scale_color_brewer()+
labs(color="Limitado")+
theme_light()
ggarrange(a,b,c)
corr<-cor(data[-4])
corrplot(corr, type="upper", order="hclust", tl.col="black", tl.srt=45,
method = "number")
Selección de variables en el modelo
El ANOVA del modelo completo con todas las variables y todas las interacciónes es la siguiente
| Df | Sum Sq M | ean Sq | F value | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Rpromedio | 1 | 25824,2 | 25824,2 | 308,3679 | < 2.2e-16 | * |
| Densidad | 1 | 2386,3 | 2386,3 | 28,4951 | 1,57E-06 | * |
| Proporcion | 1 | 73,9 | 73,9 | 0,882 | 0,35149 | |
| I(Limitacion) | 1 | 10,9 | 10,9 | 0,1305 | 0,71916 | |
| Rpromedio:Limitacion | 1 | 119,2 | 119,2 | 1,4228 | 0,23772 | |
| Densidad:Limitacion | 1 | 307,5 | 307,5 | 3,672 | 0,06018 | . |
| Proporcion:Limitacion | 1 | 7,2 | 7,2 | 0,0859 | 0,77047 | |
| Residuals | 59 | 4941,0 | 83, 7 |
Con el objetivo de efectuar el proceso de selección de variables utilizamos la función regsubsets la cual nos permite especificar un modelo óptimo de regresión lineal dado el conjunto de datos utilizado. Así, se definió como método todas las regresiones posibles
| Paso | Rpromedio | Densidad | Densidad * Limitación | I (Limitación)1 | Proporción | Rpromedio*Limitación | Proporción |
|---|---|---|---|---|---|---|---|
| Step 1 | * | ||||||
| Step 2 | * | * | |||||
| Step 3 | * | * | * | ||||
| Step 4 | * | * | * | * | |||
| Step 5 | * | * | * | * | * | ||
| Step 6 | * | * | * | * | * | * | |
| Step 7 | * | * | * | * | * | * | * |
| Step 8 | * | * | * | * | * | * | * |
Posteriormente, realizamos una evaluación de los resultados obtenidos graficando los criterios estadísticos definidos para la comparación de modelos. Si bien se observa que un modelo adecuado podría ser trabajar con 4 covariables, decidimos enfocarnos en los resultados del BIC y AIC. Aún así, para el modelo de 4 covariables, 2 de ellas no son significativas.
data_orig<-landrent
modelo2 = regsubsets(
Y~X1+X2+X3+I(X4)+X1:X4+X2:X4+X3:X4,
data=data_orig
)
Propuesta del modelo
Con base en los resultados el modelo propuesto es
\[ Y= \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2} + \varepsilon \] * \(\beta_0:\) Es el intercepto del modelo pero no tiene una interpretación práctica. * \(\beta_1:\) Por cada dolar adicional en la renta promedio de la tierra cultivable, se espera que la renta promedio por acre de alfalfa cultivado aumente 0.92 dolares * \(\beta_2:\) Por cada vaca adicional en una milla cuadrada se espera que la renta promedio por acre de alfalfa cultivado aumente 0.39 dolares
mod<-lm(data$`Renta promedio acre`~data$`Renta promedio tierra`+
data$`densidad vacas lecheras`,data = data)
anova(mod)## Analysis of Variance Table
##
## Response: data$`Renta promedio acre`
## Df Sum Sq Mean Sq F value Pr(>F)
## data$`Renta promedio tierra` 1 25824.2 25824.2 302.724 < 2.2e-16 ***
## data$`densidad vacas lecheras` 1 2386.3 2386.3 27.974 1.593e-06 ***
## Residuals 64 5459.6 85.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod)##
## Call:
## lm(formula = data$`Renta promedio acre` ~ data$`Renta promedio tierra` +
## data$`densidad vacas lecheras`, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.4827 -5.8720 0.3321 4.3855 28.6007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.11433 2.96123 -2.065 0.043 *
## data$`Renta promedio tierra` 0.92137 0.05382 17.121 < 2e-16 ***
## data$`densidad vacas lecheras` 0.39255 0.07422 5.289 1.59e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.236 on 64 degrees of freedom
## Multiple R-squared: 0.8379, Adjusted R-squared: 0.8328
## F-statistic: 165.3 on 2 and 64 DF, p-value: < 2.2e-16Evaluacion de los supuestos del modelo
Media de los errores igual a 0
\(H_0:E(\varepsilon)=0\)
\(H_1:E(\varepsilon)\neq0\)
res<-residuals(mod)
t.test(res,mu = 0)##
## One Sample t-test
##
## data: res
## t = 2.8521e-17, df = 66, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -2.218474 2.218474
## sample estimates:
## mean of x
## 3.169107e-17Supuesto de normalidad
\(H_0:\) Los errores se distribuyen normal
\(H_1:\) Los errores no tieneen distribucion normal
par(mfrow=c(1,2))
res_estud<- rstudent(mod)
qqPlot(res, col="green",pch=20)
shapiro.test(res)
text(-1, 20,"Spiro Wilk: 0.3855")
qqPlot(res_estud, col="green",pch=20)
Con base en el análisis gráfico y el valor-p de la prueba shapiro- Wilk, no se rechaza \(H_0\) por tanto se cumple el supuesto.
Supuesto de homogeneidad de varianzas
\(H_0:\) Los errores son homocedasticos
\(H_1:\) Los errores no son homocedasticos
pred=predict(mod)
plot(pred,res,pch=20,xlab="Valores predichos",ylab="Residuales")
abline(h=0,col=2,lty=2)
text(20, 25,"Breush Pagan: 0.00006")
lines(lowess(res~pred),col=4,lty=2)
Basado en el análisis grafico y la prueba de Breush-Pagan, se rechaza \(H_0\) y no se cumple el supuesto ## Supuesto de independencia
Para validar este supuesto se podrían tener en cuenta correlaciones espaciales y estudiar algunas pruebas que tengan en cuenta este problema.
Como no se cumple el supuesto de homocedasticidad entonces se aplica una transformacion de box-cox
Transformación de box-cox
library(car)
bm<-boxCox(mod)
valor_de_lambda<-bm$ lambdahat
lamvalue=bm$x[which(bm$y==max(bm$y))]
lamvalue## [1] 0.4646465Con base en el valor de \(\lambda\) encontrado, se decide realizar una transformacion de raiz cuadrada a la variable de respuesta \(Y\)
Ajustando el modelo transformado
\[ \sqrt{Y}= \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2} \]
mod2=lm(sqrt(Y)~X1+X2,data = data_orig)
par(mfrow=c(2,2))
plot(mod2)
bptest(mod2)##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 3.742, df = 2, p-value = 0.154Dagnóstico de valores atípicos
Regresión robusta
Aplicando una regresion robusta nos damos cuenta si las estimaciones de los parámetros cambian y con esto se tiene una idea de la existencia de posibles observaciones atípicas
mod4=rlm(sqrt(Y)~X1+X2,data = data_orig,method = "MM")
summary(mod4)##
## Call: rlm(formula = sqrt(Y) ~ X1 + X2, data = data_orig, method = "MM")
## Residuals:
## Min 1Q Median 3Q Max
## -1.64187 -0.43413 0.01955 0.37771 1.46387
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 2.3508 0.2164 10.8652
## X1 0.0752 0.0039 19.1135
## X2 0.0304 0.0054 5.6126
##
## Residual standard error: 0.6248 on 64 degrees of freedomsummary(mod2)##
## Call:
## lm(formula = sqrt(Y) ~ X1 + X2, data = data_orig)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5963 -0.4298 0.0418 0.3468 1.5083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.370239 0.210117 11.281 < 2e-16 ***
## X1 0.073798 0.003819 19.326 < 2e-16 ***
## X2 0.031991 0.005266 6.075 7.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6554 on 64 degrees of freedom
## Multiple R-squared: 0.8685, Adjusted R-squared: 0.8644
## F-statistic: 211.4 on 2 and 64 DF, p-value: < 2.2e-16Como no cambian las estimaciones, se decide de profundizar sobre posibles datos que ejercen una influencia muy fuerte o desproporcionada sobre los coeficientes estimados y las propiedades del modelo decidimos utilizar otros estadísticos. Así, La figura XX (bolas) muestran los valores hii de la matriz hat, la cual nos indica las observaciones que son potencialmente influyentes, pero que están lejos del espacio x a 2p/n, de las cuales se logran evidenciar x. Además se utilizan la D de Cook como alternativa más fuerte cuando queremos identificar influyentes teniendo en cuenta el espacio de x y también la variable de respuesta, utilizando el criterio de Montgomery, los atípicos son aquellos que tienen valores superiores de 1.
#cooks.distance(mod2)
corte=4/(67-3-2)
par(mfrow=c(1,2))
plot(mod2,which = 4,cook.levels = corte,las=1)
abline(h=corte,lty="dashed",col="dodgerblue2")
influencePlot(mod2,id.method="identify",main="influence Plot")
## StudRes Hat CookD
## 5 -2.4738833 0.08750174 0.18113293
## 32 -1.7487399 0.09805668 0.10736954
## 56 -0.4817041 0.12105854 0.01078246
## 67 -2.5859896 0.03392631 0.07189268Detectamos como atipicos o influyentes los puntos 5, 32, 36 y 67 los cuales serán eliminados del conjunto de datos
mod3=lm(sqrt(Y)~X1+X2,data = data_orig[-c(5, 32, 36, 67),])
summary(mod)##
## Call:
## lm(formula = data$`Renta promedio acre` ~ data$`Renta promedio tierra` +
## data$`densidad vacas lecheras`, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.4827 -5.8720 0.3321 4.3855 28.6007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.11433 2.96123 -2.065 0.043 *
## data$`Renta promedio tierra` 0.92137 0.05382 17.121 < 2e-16 ***
## data$`densidad vacas lecheras` 0.39255 0.07422 5.289 1.59e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.236 on 64 degrees of freedom
## Multiple R-squared: 0.8379, Adjusted R-squared: 0.8328
## F-statistic: 165.3 on 2 and 64 DF, p-value: < 2.2e-16par(mfrow=c(1,2))
corte=4/(63-3-2)
plot(mod3,which = 4,cook.levels = corte,las=1)
abline(h=corte,lty="dashed",col="dodgerblue2")
influencePlot(mod3,id.method="identify",main="influence Plot")
## StudRes Hat CookD
## 19 2.5745805 0.03604030 0.07552317
## 20 -1.6382433 0.09936006 0.09600105
## 33 2.2542219 0.06922191 0.11794729
## 56 -0.5474743 0.13880310 0.01629302
## 66 -1.8897135 0.09348777 0.11771446Detectando multicolinealidad
Basado en la matriz de correlación y en el factor de inflacíon de varianza, se puede ver que no existe una dependencia lineal entre las covariables, puesto que el vif es menor que 3, indicando que no existe este problema en los datos y el supuesto 2 del teorema de gauss markov no se esta violando
vif(mod2)## X1 X2
## 1.002379 1.002379Conclusiones
Basado en la pregunta de investigación y los resultados presentados, se llega a un modelo que cumple a cabalidad con los supuestos de interés, en esa medida, podemos afirmar que efectivamente existe una relacion entre la densidad de vacas lecheras por milla cuadrada y la renta promedio por acre sembrado de alfalfa, pues esta variable es significativa en dicho modelo. Pero además la construccion del modelo permite evidenciar que por cada vaca adicional en una milla cuadrada, el impacto sobre la renta es de 0.39 dolares en promedio.
Referencias
Hasbargen, Paul R, Kenneth H Thomas, and Nick Rolfes. 1977. “Land Rental Arrangements in Minnesota.”
Weisberg, Sanford. 2005. Applied Linear Regression. Vol. 528. John Wiley & Sons.