TALLER 1 MODELACIÓN LINEAL 02-2020
Mariana, Gisela Agudelo, Yonatan Moran , Elkin Danilo
21/8/2020
Partiendo de los datos registrados en la Tabla 1 correspondiente a una parcela permanente de 1/40 ha se realizaron los análisis propuestos para la solución al taller 1 de la asignatura Modelación Lineal, utilizando un archivo de valores cuyas columnas están nombradas como “dap”, “alt” y “vol” para el diámetro, la altura y el volumen medidos, respectivamente. Para los cálculos se utilizó el software R en su versión 3.6.3
## X dap alt vol
## 1 1 12.01 10.4 0.054
## 2 2 17.40 20.0 0.253
## 3 3 25.01 18.3 0.367
## 4 4 22.62 22.2 0.247
## 5 5 32.12 25.9 0.857
## 6 6 17.36 19.8 0.179
## 7 7 18.68 20.1 0.281
## 8 8 30.44 22.3 0.697
## 9 9 33.86 27.7 1.424
## 10 10 15.70 18.1 0.149
## 11 11 18.88 21.2 0.325
## 12 12 33.23 27.0 1.676
## 13 13 25.59 23.9 0.501
## 14 14 13.81 41.2 0.276
## 15 15 116.01 19.2 0.156
## 16 16 34.62 24.0 0.894
## 17 17 12.86 13.3 0.075
## 18 18 14.17 97.6 0.119
## 19 19 22.91 22.7 0.398
## 20 20 11.91 10.4 0.054
## 21 21 16.65 19.1 0.176
## 22 22 26.05 24.6 0.549
## 23 23 26.40 21.8 0.523
## 24 24 30.39 24.1 0.731
## 25 25 11.63 9.6 0.048
## 26 26 28.53 25.1 0.671Para un análisis general de los datos se hace lectura del documento en R, observamos la estructura de cada variable y se obtiene un resumen general de los datos.
summary(datos)## X dap alt vol
## Min. : 1.00 Min. : 11.63 Min. : 9.60 Min. :0.0480
## 1st Qu.: 7.25 1st Qu.: 15.94 1st Qu.:19.12 1st Qu.:0.1610
## Median :13.50 Median : 22.77 Median :22.00 Median :0.3030
## Mean :13.50 Mean : 25.72 Mean :24.22 Mean :0.4492
## 3rd Qu.:19.75 3rd Qu.: 29.93 3rd Qu.:24.48 3rd Qu.:0.6405
## Max. :26.00 Max. :116.01 Max. :97.60 Max. :1.6760str(datos)## 'data.frame': 26 obs. of 4 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ dap: num 12 17.4 25 22.6 32.1 ...
## $ alt: num 10.4 20 18.3 22.2 25.9 19.8 20.1 22.3 27.7 18.1 ...
## $ vol: num 0.054 0.253 0.367 0.247 0.857 ...1. Proponga, lógica y justificadamente, con estas variables un modelo de regresión lineal simple, sin ignorar ninguna de las variables, luego de una revisión exhaustiva de los datos.
Inicialmente, se plantea el comportamiento de los datos que componen cada una de las variables en estudio, y se relaciona su comportamiento con respecto a las demás. Para ello utilizamos diagramas de cajas para el DAP y para la altura obteniendo la media y desviación estándar de cada uno de ellos de una manera gráfica, lo que implica que cualquier dato por fuera de dicho diagrama debería ser posteriormente analizado. Luego, se utiliza el test de Dotchart el cual nos plantea una nube de puntos con el comportamiento de los datos correspondientes a la variable estudiada. Cuando un dato que queda por fuera del diagrama de cajas y bigotes coincide con un comportamiento aislado en el test de Dotchart, queda más cerca de la posibilidad de ser descartado, pero antes de ello se analiza la ontología de los individuos, y se logra establecer de manera lógica la consecuencia de los datos de DAP con respecto a los de altura.
windows()
par(mfrow = c(1,2))
## Datos DAP
dotchart(datos$dap, main = "Test de Dorchart para DAP", pch = 1,col=2)
boxplot(datos$dap, main = "Boxplot para DAP", col= 2)
## Datos Altura
dotchart(datos$alt, main = "Test de Dotchart altura", col= 4)
dotchart(datos$alt, main = "Boxplot para altura", col= 4)
De esta manera, se obtiene que las observaciones 15 y 18 presentes en la Tabla 1, deben ser descartadas para los análisis posteriores, debido a que no sólo fueron datos extraordinarios gráficamente, sino que también la observación 15 presentaba un valor muy extremo de diámetro a la altura del pecho para una altura relativamente pequeña, análogo a la observación 18 en la que obtenemos un valor demasiado alto para la altura con respecto al diámetro del individuo. Después de analizar los datos y tomar las decisiones planteadas con anterioridad, se requiere plantear un modelo que incluya una variable dependiente y una variable independiente. Partiendo de la premisa de que el fuste en los individuos arbóreos podría ser representado como un conjunto de conos integrados para los cuales existe un volumen conocido dado por la Ecuación 1, pero se sabe que la conicidad sería un caso ideal, por lo que se multiplica la anterior ecuación por un factor de forma según el individuo (FF) obteniendo la Ecuación 2, de esta manera logramos plantear el Volumen como una función de DAP^2 por la Altura, el cual dependerá de un valor desconocido para FF * (π/4), y de esta manera logramos establecer un modelo lineal simple cuya variable dependiente será el volumen y su variable independiente será una variable combinada dada por DAP^2 * Altura, obteniendo la Ecuación 3 para el modelo de regresión donde V es Volumen en metros cúbicos, d es el diámetro a la altura del pecho (DAP) en centimetros, H es la altura en metros y los β_0 y B_1 son parámetros del modelo intercepto y pendiente, respectivamente.
\[ V = \frac{pi}{4}*d^2*h \]
\[ V = \{beta__0} + {beta_1}*(DAP^2*H) \] \[ V = \frac{Pi}{4} * d^2 * h * FF\] Ecuacion 3
Para proseguir con los procesos eliminamos de la base de datos las observaciones que planteamos con anterioridad y añadimos a la misma la variable combinada de esta manera:
Posteriormente, graficamos un posible modelo lineal para la base de datos original (recta verde) y uno al retirar los datos seleccionados (recta azul). Gráfica 3
with(datos, plot(dap2alt, vol))
abline(lm(vol ~ dap2alt, datos), col = "green")
abline(lm(vol ~ dap2alt, datos[-c(15, 18),]), col = "blue")
###1.1 Encuentre de manera empírica el modelo para estos datos y explique sus procedimientos. Para hallar el modelo de esta manera, inicialmente se hizo un estimado burdo con la observación de los datos con los que contamos en la tabla 1 excluyendo las dos observaciones necesarias. Como intercepto, se buscó el dato Y para el cual X \[(DAP^2 * Altrua)\] tuviera un valor de 0, pero ningún valor lo cumplió, por lo tanto, se asumió que el modelo no cuenta con un intercepto. La pendiente se obtuvo con la división del cambio en Y entre el cambio en X, del valor máximo y mínimo de cada variable; de la siguiente manera: \[ b = \frac{Δy}{Δx}\].
datos1 <- datos[-c(15, 18,9),]
(b1 <- (max(datos1$vol) - min(datos1$vol))/(max(datos1$dap2alt) - min(datos1$dap2alt)))## [1] 5.709111e-05range(datos1$dap2alt)## [1] 1298.466 29814.288(b0 <- 0)## [1] 0Posteriormente, se graficó el modelo empírico resultante, pero, como se puede observar en la Ilustración 6, los datos no se ajustan de la mejor manera a dicho modelo, quedando muy alejados de la recta de regresión hallada, por lo que se obtendrían residuales muy altos.
with(datos1, plot(dap2alt, vol, xlab= "DAP^2 * Altrua",
ylab = "Volumen", main = "Grafica Modelo Empirico"))
abline(b0, b1, col = 4)
1.2 Encuentre sus parámetros de acuerdo con las definiciones salidas para ellos del método de los mínimos cuadrados
Teniendo en cuenta que los mínimos cuadrados hacen alusión a los residuales obtenidos por un modelo determinado y que a partir de ellos podemos obtener los mejores parámetros para el modelo; utilizamos las definiciones establecidas por el método de los mínimos cuadrados para adquirir un valor de b0 y b1 aplicando las fórmulas de la Ecuación 4 y 5. \[ b_1 =\sum_{i=1}^{n} x_i \frac {-b \pm \sqrt {b^2 - 4ac}}{2a}\] \[ b_1= \frac{\sum X_iY_i - \frac{\ (\sum X_i \sum Y_i)}{n}}{\sum X_i^2 - \frac{\ (\sum X)^2}{n}} = \frac{\sum (X_i - \overline{X})(Y_i - \overline{Y})}{\sum (X_i - \overline{X})^2} = \frac{SPCXY}{SCCX}\] Ecuacion 4
\[ b_o= \frac{(\sum Y_i - b_i \sum X_i)}{n} = \overline{Y}-b_i \overline{X} \] Ecuacion 5
datos1## X dap alt vol dap2alt
## 1 1 12.01 10.4 0.054 1500.097
## 2 2 17.40 20.0 0.253 6055.200
## 3 3 25.01 18.3 0.367 11446.652
## 4 4 22.62 22.2 0.247 11358.950
## 5 5 32.12 25.9 0.857 26720.885
## 6 6 17.36 19.8 0.179 5967.118
## 7 7 18.68 20.1 0.281 7013.742
## 8 8 30.44 22.3 0.697 20663.037
## 10 10 15.70 18.1 0.149 4461.469
## 11 11 18.88 21.2 0.325 7556.833
## 12 12 33.23 27.0 1.676 29814.288
## 13 13 25.59 23.9 0.501 15650.870
## 14 14 13.81 41.2 0.276 7857.503
## 16 16 34.62 24.0 0.894 28765.066
## 17 17 12.86 13.3 0.075 2199.549
## 19 19 22.91 22.7 0.398 11914.506
## 20 20 11.91 10.4 0.054 1475.220
## 21 21 16.65 19.1 0.176 5294.950
## 22 22 26.05 24.6 0.549 16693.622
## 23 23 26.40 21.8 0.523 15193.728
## 24 24 30.39 24.1 0.731 22257.606
## 25 25 11.63 9.6 0.048 1298.466
## 26 26 28.53 25.1 0.671 20430.419datos1$xiyi <- with(datos1, dap2alt*vol)
datos1$xi2 <- with(datos1, dap2alt^2)
n1 <- length(datos1$dap)
b1_mod <- with(datos1,
((sum(xiyi)-((sum(dap2alt)*sum(vol))/n1))/
(sum(xi2)-((sum(dap2alt))^2)/n1)))
b0_mod <- with(datos1, mean(vol) - (b1_mod * mean(dap2alt)))Finalmente, se obtuvo el valor para b0 y b1, correspondientes a -0.0558 y 4.067e-05, respectivamente. Logrando plantear el siguiente modelo para el Volumen y se elabora la Gráfica 5 en la cual se evidencia un mayor ajuste de los datos al modelo obtenido y por ende valores menores para los residuales con respecto a los del modelo obtenido de manera empírica.
\[ vol= -0.0558 + 4.067_e-0.5 * (dap^2 * alt) \]
with(datos1,plot(dap2alt,vol, xlab = "DAP^2 * Altura",
ylab = "Volumen", main = " Grafica Modelo por metodo de minimos cuadrados"))
abline(b0_mod,b1_mod, col=4)
2. Ajuste por medio de R el modelo elegido en 1 y analícelo profusamente conceptual y gráficamente
Una vez realizado el análisis preliminar de los datos y eliminando el dato que se consideró que afectaba notoriamente el ajuste del modelo se procedió a establecer el modelo con la ayuda de R, el cual se describe mediante los siguientes parámetros: \[ y = b0 * b_1 (X) \] donde: Y= Volumen estimado por el modelo X = Variable combinada \[ (DAP^2 * Alt)\] \[b_0 y b_1 =\] parámetros del modelo
Se realizó el procedimiento en el Software para obtener el ajuste del modelo planteado
modelo <- lm(vol~dap2alt, datos1)
summary(modelo)##
## Call:
## lm(formula = vol ~ dap2alt, data = datos1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.17778 -0.06326 -0.01268 0.03525 0.56371
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.868e-02 5.147e-02 -0.752 0.461
## dap2alt 3.861e-05 3.414e-06 11.306 2.17e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.144 on 21 degrees of freedom
## Multiple R-squared: 0.8589, Adjusted R-squared: 0.8522
## F-statistic: 127.8 on 1 and 21 DF, p-value: 2.166e-10Una vez realizado el procedimiento en R para ajustar el modelo se encontraron los siguientes resultados
\[ vol (m^3)= -5.6_e - 02 + 4.07_e -05 *(dap^2 * alt) \] 
A partir del modelo que se ajustó con la ayuda de R se pueden hacer los siguientes análisis:
-Si bien los parámetros del modelo no fueron significativos los dos como fue el caso del b0 el cual arrojó un resultado de no significancia, se dice que el modelo es válido porque b1 es significativo, esto quiere decir que al b1 ser significativo es diferente de 0, lo cual establece una pendiente para la regresión independientemente de si sea una relación fuerte o débil entre las variables. -Al analizar el valor de R del modelo se aprecia un valor de R correspondiente a 0,89, dicho valor al realizar el análisis con la tabla de Snedecor para 22 grados de libertad y una variable independiente la cual es combinada, con un valor p de 0.05 se obtiene que para valores superiores a un R de 0.41 será significativo dicho valor. Para el caso que estamos analizando se obtiene entonces que el valor de R es significativo; dado este resultado se dice que la variable independiente tiene un efecto significativo al lograr la reducción en la variabilidad de la variable dependiente. -Con respecto a los residuales del modelo (Gráfica 7), si bien se espera que se comporten de una forma normal se observa cómo con la prueba de Shapiro - Wilk se presenta una no normalidad, esto puede ser el resultado de pocos datos utilizados en el ajuste del modelo si se tiene en cuenta la sensibilidad de dicha prueba. Por otro lado, cuando realizamos la prueba gráfica se obtiene que se distribuyen a lo largo de la recta de la normalidad (Gráfica 8) lo que nos lleva a aceptar dicho supuesto.
plot(resid(modelo),main = "Residuales del Modelo", pch=1, col= 10, ylab = "Residuales")
qqnorm(resid(modelo))
qqline(resid(modelo),col="red")
3. Verifique las propiedades de los residuales usando las fórmulas dadas para ellos.
1. La media de los residuales es cero para verificar esta propiedad se trabajo con la ecuacion 6: \[ \overline{e} = \frac{\sum e_i }{n} \] \[ E ( e_i)=0 \]
pred1 <- predict(modelo)
resid1 <- datos1$vol - pred1
mean(resid1)## [1] 3.288697e-17Una vez obtenidos los residuales y posteriormente calculada la media se obtuvo un resultado que fue de 2.573145e-17 el cual asumimos como 0, verificando así el cumplimiento de la primera propiedad. 2. La viaranza de los residuales es igual a MSE
var(resid1)## [1] 0.01978959summary(aov(modelo))## Df Sum Sq Mean Sq F value Pr(>F)
## dap2alt 1 2.6502 2.6502 127.8 2.17e-10 ***
## Residuals 21 0.4354 0.0207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MSE <- anova(modelo)$"Mean Sq"[2]
MSE## [1] 0.02073195En esta segunda propiedad se calculó la varianza de los residuales del modelo trabajado obteniendo un valor de 0.0208657. Posteriormente, se extrajo el valor del Cuadrado Medio de los Residuales (MSE) del modelo cuyo cálculo parte de la Ecuación 7, obteniendo un valor de 0.02181414, concluyendo que los valores no son exactamente iguales, pero que al aproximar a tres cifras significativas se cumple la propiedad. \[ \frac{\displaystyle\sum_{i=1}^{n} (e_i - \overline e)^2}{n-2}= \frac{SSE}{n-2}= MSE \]
Con un resultado de MSE= 0.02073195, donde se concluye que los valores si bien no son exactamente iguales si se aproximan cumpliendo esta segunda propiedad. 3. normalización de residuales Como lo planteamos en el punto 2, al realizar la prueba de Shapiro – Wilk obtenemos no normalidad, pues el valor p es menor a 0.05 lo que conlleva a rechazar la hipótesis nula del test (Los datos se distribuyen de manera normal), pero al observar la Gráfica 8 podemos deducir que se cumple esta propiedad \[ e_{iest} = e_{is} = \frac{e_i}{\sqrt{(MSE)}}= \frac{(e_i - \overline{e})}{S} \]
shapiro.test(resid1 / sqrt(MSE))##
## Shapiro-Wilk normality test
##
## data: resid1/sqrt(MSE)
## W = 0.6878, p-value = 9.764e-064. Los residuales no son aleatorios
4 -Aconstruya (y explique) con R una tabla de residuales con la función: tabderes <- cbind(rstandard(modelo),(y-ypredicho))
En este numeral obtuvimos la Tabla 2 en la que se tabulan los residuales estandarizados hallados con una función propia del Software y los residuales obtenidos mediante la expresión (Y_i-Y ̂_i). Es evidente que los residuales de la segunda columna son relativamente más pequeños que los de la primera y oscilan en torno a 0, mientras que los estandarizados llegan a valores más alto, incluso encontrando un valor de 3.865 en contraste a un valor de 0.519, por lo que dicha observación se puede considerar como un valor atípico en el ajuste del modelo, puesto que los residuales estandarizados no sólo plantean la distancia entre un valor predicho y uno observado de la variable dependiente con respecto a la variable explicatoria, sino que además se relacionan con el número de desviaciones estándar que puede presentar una observación con respecto a la media de los residuales que por propiedades es igual a 0. En este caso, no eliminamos el dato para los análisis posteriores.
tabderes<-cbind(rstandard(modelo),(datos1$vol - pred1))
write.csv2(tabderes, "residuales.cesv")B)Use la función de la library(MASS): stdres(modelo) y compárela con studres(modelo)
Los residuales estudentizados y estandarizados, son normalizados a la variación de la unidad, pero en la función de esta segunda se tiene en cuenta el punto de datos actual y en la primera no. Se encuentran los resultados propuestos a continuación. Observamos que se obtienen valores muy cercanos entre los residuales estandarizados y los estudentizados para cada observación. Además, los residuales estudentizados pretenden realizar un ajuste a los estandarizados, y al igual que en el punto anterior obtenemos valores muy descontextualizados para la observación 12, esto debido a que los residuales estudentizados son más sensibles a la hora de encontrar valores atípicos.
library(MASS)
stdres(modelo)## 1 2 3 4 5 6
## 0.25573453 0.41602512 -0.25720106 -1.08550536 -1.03039706 -0.09108778
## 7 8 10 11 12 13
## 0.35020254 -0.44979765 0.11170005 0.51427761 4.42471820 -0.45970282
## 14 16 17 19 20 21
## 0.08102159 -1.37786300 0.21062554 -0.16527783 0.26284243 0.07399454
## 22 23 24 25 26
## -0.40549957 -0.17703749 -0.65563875 0.26922658 -0.57258156Al graficar los residuales estandarizados vs los estudentizados Gráfica 9, logramos observar que tienen una relación casi 1:1 y que sólo a medida que aumenta la desviación estándar de los mismos aumenta la diferencia entre ellos.
plot(stdres(modelo), studres(modelo), main = "Residuales estandarizados vs Residuales Estudentizados",
xlim = c(-2,7), ylim = c(-2,7), xlab = "Estandarizados", ylab = "Estudentizados", col = 4)
C)Encuentre los intervalos de confianza para los parámetros del modelo ajustado
Los valores más probables de los intervalos de confianza para los parámetros B_0 y B_1 del modelo ajustado se hallaron con ayuda de la función confint, obteniendo así, que para el ajuste de modelos lineales simples satisfactorios, lo más adecuado sería encontrar valores para el intercepto entre -1.624e-01 y 0.051, y para la pendiente entre 3.405e-05 y 0.00005..
confint(modelo)## 2.5 % 97.5 %
## (Intercept) -1.457128e-01 0.0683525075
## dap2alt 3.150378e-05 0.0000457053###5a)Calcule los intervalos de confianza de su variable dependiente, para los valores de las observaciones 5, 10, 15, 20 y, 25 y haga sus gráficas
Cuando se obtiene un modelo de regresión es posible predecir el valor de la variable dependiente para un determinado valor de la variable independiente, pero cada una de estas predicciones tendrá una probabilidad de error permitida la cual determina el investigador, además cada predicción tendrá intervalos de confianza que nos definen valores extremos para un intervalo que contendrá el valor real de la observación; de esta manera el ejercicio propone calcular dichos intervalos para 5 valores. Para el cálculo de dichos intervalos utilizamos las Ecuaciones 9 y 10.
\(\hat{Y_h} +- t_{1-\frac{\alpha}{2};(n-2)gl} s(\hat{Y_h})\) $s^2()=MSE[+] $
b) Encuentre gráficamente los intervalos de confianza para las posibles rectas de regresión y explique lo obtenido por dos procesos diferentes
###Calcule con sus intervalos de confianza para y escoja algunos valores para los cuales desee hacer algunas predicciones, explicando bien cada proceso
Al calcular los intervalos de confianza para las observaciones propuestas, obtenemos la Tabla 3 y luego ilustramos dichos intervalos para cada observación en la Gráfica 10.
Para encontrar los intervalos de confianza de las rectas de regresión gráficamente, utilizamos dos métodos, el primero consiste en instruir al software un algoritmo donde se plantean cada uno de los pasos necesarios para completar la ecuación que nos permite encontrar los límites superior e inferior del intervalo y otro en el que mediante el uso de un paquete el algoritmo queda de manera interna en la función “ggplot” y se hace automáticamente. Para ambos métodos obtenemos gráficas semejantes, que sólo difieren en cuestiones estéticas (Ver Gráfica 11 y 12), pero que nos dan la misma información. Cualquier recta de regresión modelada correctamente quedaría entre el intervalo que proponen dichas ilustraciones.
plot(datos1$dap2alt,datos1$vol,pch=20,main="Intervalos de confianza para las rectas de regresión",
cex.lab=1.5, cex.main=1.5, cex.axis =1.5, xlab = "DAP^2 * H", ylab = "Volumen")
abline(modelo,col="red")
ci.lines(modelo)
library(ggplot2)
ggplot(datos1, aes(x=dap2alt, y=vol)) + geom_point() + ggtitle ("Diagrama de Dispersión y limites de
confianza recta") + xlab("Dap^2 * H") + ylab("Volumen") + geom_smooth(method=lm)## `geom_smooth()` using formula 'y ~ x'
6 Calcule con sus intervalos de confianza para escoja algunos valores para los cuales desee hacer algunas predicciones, explicando bien cada proceso ```