REGRESIÓN LINEAL MÚLTIPLE
Jorge Méndez González (UAAAN)
2020-05-09
Su estructura
La regresión lineal múltiple permite generar un modelo lineal en el que el valor de la variable dependiente o respuesta (\(Y\)) se determina a partir de un conjunto de variables independientes llamadas predictores (\(X1\), \(X2\), \(X3\),…, \(Xn\)). Es una extensión de la regresión lineal simple.
Los modelos lineales múltiples tiene la siguiente forma: \(Y_{i}=(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i}+\cdots+\beta_{n}X_{ni})+e_{i}\)
Donde \(Y\) es la variablde dependiente.
\(\beta_{0}\) es la ordenada en el origen (o intercepta), el valor de la variable dependiente \(Y\) cuando todos los predictores son cero.
\(\beta_{i}\) es el efecto promedio que tiene el incremento en una unidad de la variable predictora \(Yi\) sobre la variable dependiente \(Xi\), manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.
\(e_{i}\): es el residuo o error. Calculado como la diferencia entre el valor observado y el estimado por el modelo.
El modelo de Regresion Lineal Multiple (MRM) tiene la siguiente estructura matricial. \[Y_{nx1}=X_{nx(k+1)}\beta _{k+1}+\varepsilon _{nx1}\]
Modelo \[\begin{bmatrix}Y_{1} \\Y_{2} \\\vdots \\ Y_{n} \end{bmatrix}=\begin{bmatrix} 1 & X_{11} & \cdots & X_{1k}\\ 1 & X_{12} & \cdots & \cdots\\ \vdots & \vdots & \ddots & \ddots \\ 1 & X_{n1} & \cdots & X_{nk} \end{bmatrix}\begin{bmatrix} \beta _{0}\\ \beta _{1}\ \\\vdots \ \\\beta _{k} \ \end{bmatrix}+\begin{bmatrix} \varepsilon _{1} \\ \varepsilon _{2} \\ \vdots \\ \varepsilon _{n} \end{bmatrix}\]
Un ejemplo
Vamos a construir un modelo lineal múltiple para predecir la biomasa (kg) a partir del volumen maderable \((m^{3})\), utilizando los datos disponibles en las librerías de r.
Para la estimación biomasa de árboles existen los métodos directos e indirectos: i) el primero consiste en derribar los árboles y separarlos por componentes para obtener el peso verde y seco (Picard et al., 2012) y ii) el segundo, por medio de ecuaciones de regresión (Avendaño et al., 2009; Návar, 2009). Cuando no es posible derribar los árboles, es conveniente utilizar el volumen del árbol, el cual es obtenido a través de métodos directos o por troceo simulado, posteriormente el volumen puede ser transformado a biomasa usando su densidad básica (Silva y Návar, 2010).
En 1855 Maximilian Robert Pressler, postuló un método para obtener el volumen del fuste, en el que solo se necesita medir el diámetro normal (1.3 m) y la altura a la que se encuentra la mitad del diámetro normal, denominado de esta manera como ´ecuación de Pressler´ y puede aplicarse a cualquier forma geométrica que presente el fuste siempre que sea no cilíndrico (Romahn, 1999).
Las librerías
El primer paso para hacer este ejercicio será cargar las siguientes librerías en r.
rm(list=ls(all=TRUE))
setwd("G:/UAAAN/MATERIAS/2020/EPIDOMETRIA/Modelos de volumen en r")
library(PerformanceAnalytics)
library(corrplot)
library(ecospat)
library(ggplot2)
library(dplyr)
library(tidyr)
library(sampler)
library(reshape)
library(tidyverse)Enseguida cargamos los datos del archivo llamado “Muestra.csv”. Al escribir attach(Muestra[,2:5], le estamos indicando que agregue solo de las columnas 2 a la 5. La columna uno no porque es el ID de los datos y no lo necesitamos. Asi mismo hacemos una gráfica de Diámetro (D) y Biomasa (B).
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendMuestra <- read.csv("Muestra.csv") # Cargar los datos
attach(Muestra[,2:5]) # Adjuntar columnas de la 2 a la 5
head(Muestra[,2:5]) # Mostrar encabezados## D H V B
## 1 5.570423 24.9936 1.5772458 914.8025
## 2 4.615493 22.5552 1.0278998 596.1819
## 3 3.437747 25.2984 0.5578410 323.5478
## 4 3.628733 23.1648 0.6059795 351.4681
## 5 3.342254 21.9456 0.4643955 269.3494
## 6 2.737465 19.8120 0.2916630 169.1646summary(Muestra[,2:5]) # Resumir estadsiticas de los datos## D H V B
## Min. :2.642 Min. :19.20 Min. :0.2888 Min. : 167.5
## 1st Qu.:3.533 1st Qu.:21.95 1st Qu.:0.5578 1st Qu.: 323.5
## Median :4.106 Median :23.16 Median :0.6853 Median : 397.5
## Mean :4.254 Mean :23.06 Mean :0.8677 Mean : 503.3
## 3rd Qu.:5.093 3rd Qu.:24.38 3rd Qu.:1.0845 3rd Qu.: 629.0
## Max. :6.557 Max. :26.52 Max. :2.1804 Max. :1264.6plot(D, B) # Graficar x=D, y=B
Correlacion múltiple
La correlación múltiple cuantifica el grado de asociación entre una variable dependiente y dos o más independientes, tomadas en conjunto, Al igual que el coeficiente de correlación de Pearson, del que se puede considerar una extensión, encuentra su fundamento en el análisis de la regresión, en este caso múltiple. Varía entre 0 y 1, y se presenta separando la variable dependiente de las independientes con un punto.
El coeficiente de correlación de Pearson es la covarianza estandarizada, y su ecuación difiere dependiendo de si se aplica a una muestra, Coeficiente de Pearson muestral (r), o si se aplica la población Coeficiente de Pearson poblacional (ρ).
\[r_{xy}=\dfrac {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) \left( y_{i}-\overline {y}\right) } {\sqrt {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) ^{2}\sum _{i=1}^{n}}\overline {\left( y_{i}-\overline {y}\right) ^{2}}}\] \[\rho=\dfrac {Cov(X,Y)}{\sigma_{x}\sigma_{y}}\] Hacemos un análisis de correlación entre las variables. Se muestra de forma gráfica así como el valor de la correlación. Los “*” indican la significancia estadística que hay entre las variables.
chart.Correlation(Muestra[,2:5]) # Correlacionar y ver figuras
cor(Muestra[,2:5]) # Correlacionar y ver tablas## D H V B
## D 1.0000000 0.5813875 0.9676677 0.9676677
## H 0.5813875 1.0000000 0.6422566 0.6422566
## V 0.9676677 0.6422566 1.0000000 1.0000000
## B 0.9676677 0.6422566 1.0000000 1.0000000En la Figura se indica que la correlación entre Diámetro (D) y Altura (H) es de 0.52. Entre Altura (H) y volumen (V) la correlación es de 0.62.. La correlación entre volumen (V) y Biomasa (B) es de 1.00, pero es debido a que la biomasa se derivó directamente de volumen.
Correlacion parcial
La correlación estudia la relación (lineal o monotónica) existente entre dos variables. Puede ocurrir que la relación que muestran dos variables se deba a una tercera variable que influye sobre las otras dos, a este fenómeno se le conoce como confounding. La correlación parcial permite estudiar la relación lineal entre dos variables bloqueando el efecto de una tercera (o más) variables. Si el valor de correlación de dos variables es distinto al valor de correlación parcial de esas mismas dos variables cuando se controla una tercera, significa que la tercera variable influye en las otras dos.
Primero hagamos la correlación entre diámetro (D) y biomasa (B). Nota, aquí no importa el orden de la variables.
library(MASS)
cor.test(x = D, y = B, method = "pearson")##
## Pearson's product-moment correlation
##
## data: D and B
## t = 18.399, df = 23, p-value = 2.956e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9269681 0.9858526
## sample estimates:
## cor
## 0.9676677Se observa que la correlación entre diámetro (D) y biomasa (B) es de 0.9632961, y acorde a su intervalo de confianza (95 percent confidence interval:) varía de 0.9173248 a 0.9839196 Por su p-value < 1.244e-14, se concluye que es altamente significativo.
Ahora hagamos la correlación parcial entre estas mismas variables (entre diámetro (D) y biomasa (B)), pero dejando como tercera variable la Altura (H).
library(ppcor)
pcor.test(x = D, y = B, z =H,
method = "pearson")## estimate p.value statistic n gp Method
## 1 0.9529072 6.9973e-13 14.73817 25 1 pearsonCuando bloqueamos la variable Altura (H), la correlación entre diámetro (D) y biomasa (B) es de r = 0.9556755 y su p-value < 3.640216e-13. Como las correlaciones no son iguales, se concluye que la relación entre Diámetro (D) y Biomasa (B), está influida por la Altura (H). Lo mismo puede hacer con otras variables para determinar la influencia de cada una.
Ajuste del modelo
Ajustemos el modelo lineal múltiple siguiente: \(B = \beta_{0} + \beta_{1}D+ \beta_{2}H\). Se trata de ajustar el modelo anterior para predecir Biomasa (B) en función de las variables Diámetro y Altura (H), a este “objeto”" le llamaremos Modelo_1. Con “summary(Modelo_1)”, podemos ver los coeficientes de regresión, y sus estadísticos.
Modelo_1<-lm(B ~ D + H, data = Muestra) # Ajustar el moddelo
summary(Modelo_1) # Resumir el modleo##
## Call:
## lm(formula = B ~ D + H, data = Muestra)
##
## Residuals:
## Min 1Q Median 3Q Max
## -103.152 -60.234 -4.571 55.793 137.767
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -937.796 172.151 -5.448 1.8e-05 ***
## D 244.616 16.597 14.738 7.0e-13 ***
## H 17.372 8.793 1.976 0.0609 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 70.78 on 22 degrees of freedom
## Multiple R-squared: 0.946, Adjusted R-squared: 0.9411
## F-statistic: 192.6 on 2 and 22 DF, p-value: 1.146e-14Segun los resultados, el valor de \(\beta_{0}\) = -1051.833, \(\beta_{1}\) = 242.473 y \(\beta_{2 }\) = 22.531. De la misma forma, la R cuadrada es = 0.9461 y la R-cuadrada ajustada es = 0.9412. Si quisiéramos ver el análisis de varianza basta con escribir “Anova(Modelo_1)” y a ese objeto lo llamamos por ejemplo….Anova
library(car)## Loading required package: carDataAnova<-Anova(Modelo_1)
Anova## Anova Table (Type II tests)
##
## Response: B
## Sum Sq Df F value Pr(>F)
## D 1088280 1 217.2138 6.997e-13 ***
## H 19558 1 3.9037 0.06085 .
## Residuals 110224 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Volviendo… Lo primero es ver si los coeficientes de regresión son diferentes de “cero”, esto es; si son estadísticamente significativos. Según el Pr(>|t|), se tiene “***”, “***” y “*”, para cada coeficiente, es decir, son significativos.
MUY IMPORTANTE. En caso de que \(\beta_{0}\) NO SEA SIGNIFICATIVO, se puede ajustar el modelo SIN ese coeficiente, quedando su estructura así: \(B = \beta_{1}D+ \beta_{2}H\). Si los coeficientes \(\beta_{1}\) y \(\beta_{2}\) no son significativos, NO SE DEBEN ELIMINAR, esto indicaría que el modelo no es bueno.
NO se trata de solo eliminar el coeficiente \(\beta_{1}\), sino de ajustar nuevamente el modelo. Observe el modelo original se escribe así: ** Modelo_1<-lm(B ~ D + H, data = Muestra)** y sin el coeficiente \(\beta_{0}\) se escribiría así: Modelo_2<-lm(B ~ D + H -1, data = Muestra), lo que se agregó al modelo fue -1.
Modelo_2<-lm(B ~ D + H -1, data = Muestra) # Ajustar el moddelo
summary(Modelo_2) ##
## Call:
## lm(formula = B ~ D + H - 1, data = Muestra)
##
## Residuals:
## Min 1Q Median 3Q Max
## -236.02 -57.05 -21.43 92.40 230.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## D 269.432 23.923 11.263 7.75e-11 ***
## H -27.608 4.529 -6.095 3.23e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 106.1 on 23 degrees of freedom
## Multiple R-squared: 0.9691, Adjusted R-squared: 0.9664
## F-statistic: 360.4 on 2 and 23 DF, p-value: < 2.2e-16Observe que ya no aparece el coeficiente de la intercepta ((Intercept), solo los coeficientes de Diámetro (D) y Altura (H). que ahora son: 260.543 Y -25.699, que son diferentes cuando el modelo aún tenía la intercepta.
Regresando… veamos la estructura del objeto “Modelo_1”, al escribir “str(Modelo_1”) le dará todas las características que contiene el objeto. Así, por ejemplo, si deseamos ver los residuales solo escribiríamos Modelo_1$ residuals, si deseamos ver los valores estimados, basta con escribir Modelo_1$ fitted.values. Veamos también los intervalos de confianza de los coeficientes de regresión.
str(Modelo_1)## List of 12
## $ coefficients : Named num [1:3] -937.8 244.6 17.4
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "D" "H"
## $ residuals : Named num [1:25] 55.793 13.124 -19.071 -0.803 8.335 ...
## ..- attr(*, "names")= chr [1:25] "1" "2" "3" "4" ...
## $ effects : Named num [1:25] -2516.45 -1382.1 -139.85 -8.08 -1.16 ...
## ..- attr(*, "names")= chr [1:25] "(Intercept)" "D" "H" "" ...
## $ rank : int 3
## $ fitted.values: Named num [1:25] 859 583 343 352 261 ...
## ..- attr(*, "names")= chr [1:25] "1" "2" "3" "4" ...
## $ assign : int [1:3] 0 1 2
## $ qr :List of 5
## ..$ qr : num [1:25, 1:3] -5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:25] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "D" "H"
## .. ..- attr(*, "assign")= int [1:3] 0 1 2
## ..$ qraux: num [1:3] 1.2 1.03 1.36
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-07
## ..$ rank : int 3
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 22
## $ xlevels : Named list()
## $ call : language lm(formula = B ~ D + H, data = Muestra)
## $ terms :Classes 'terms', 'formula' language B ~ D + H
## .. ..- attr(*, "variables")= language list(B, D, H)
## .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "B" "D" "H"
## .. .. .. ..$ : chr [1:2] "D" "H"
## .. ..- attr(*, "term.labels")= chr [1:2] "D" "H"
## .. ..- attr(*, "order")= int [1:2] 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(B, D, H)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
## .. .. ..- attr(*, "names")= chr [1:3] "B" "D" "H"
## $ model :'data.frame': 25 obs. of 3 variables:
## ..$ B: num [1:25] 915 596 324 351 269 ...
## ..$ D: num [1:25] 5.57 4.62 3.44 3.63 3.34 ...
## ..$ H: num [1:25] 25 22.6 25.3 23.2 21.9 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language B ~ D + H
## .. .. ..- attr(*, "variables")= language list(B, D, H)
## .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:3] "B" "D" "H"
## .. .. .. .. ..$ : chr [1:2] "D" "H"
## .. .. ..- attr(*, "term.labels")= chr [1:2] "D" "H"
## .. .. ..- attr(*, "order")= int [1:2] 1 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(B, D, H)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
## .. .. .. ..- attr(*, "names")= chr [1:3] "B" "D" "H"
## - attr(*, "class")= chr "lm"Modelo_1$ residuals # Ver los residuales## 1 2 3 4 5 6
## 55.7931869 13.1235397 -19.0707227 -0.8032535 8.3353785 93.1565660
## 7 8 9 10 11 12
## 77.9395641 -4.5711207 -7.3727511 -41.5282520 -103.1519206 -93.8695410
## 13 14 15 16 17 18
## 137.7666962 -11.9651905 37.0844813 86.5315480 -49.7401240 90.0404612
## 19 20 21 22 23 24
## -70.8368226 -82.7921866 21.0841696 71.7338125 -80.0009708 -60.2342550
## 25
## -66.6522931Modelo_1$ fitted.values # Ver los valores estimados## 1 2 3 4 5 6 7
## 859.00935 583.05837 342.61848 352.27138 261.01402 76.00800 879.56471
## 8 9 10 11 12 13 14
## 331.40363 352.27138 887.35107 553.16251 458.47666 1126.86159 268.17560
## 15 16 17 18 19 20 21
## 360.37012 80.99064 887.35107 79.12410 591.46951 504.88241 350.09244
## 22 23 24 25
## 838.14161 393.69448 689.26365 475.60352confint(Modelo_1) # ver los intervalos de confianza de los coeficienets## 2.5 % 97.5 %
## (Intercept) -1294.8153156 -580.77663
## D 210.1946802 279.03655
## H -0.8625057 35.60683Puede ver que el residual de la primera observación es 99.671294, el de la segunda es 110.853830 y así sucesivamente. Estos son obtenidos de la resta de los valores observados de Biomasa (B) menos los valores de Biomasa (B)estimados con el modelo.
De la misma forma, el valor de \(\beta_{0}\) = -1051.833, va desde -1399.956928 hasta -703.70859. Lo mismo para el resto de los coeficientes de regresión.
Los observados y estimados
Para ver como se distribuyen los datos observados y estimados de Biomasa (B) en función de la variable independiente, es necesario hacer una figura. Es importante que notar que aquí se están usando dos variables independientes. Diámetro (D) y Altura (H), por lo que esta figura puede hacerla solo con solo variable, o puede hacerlas de todas separadamente. Aquí lo hacemos para ambas.
# Graficar los valores predichos y observados con la variable D
plot(D, B, xlab = "Diámetro normal (cm)", ylab = "Biomasa (kg)", cex=2)
lines(D, fitted.values(Modelo_1), lwd=3, col="red")
Al ser una regresión de más de una variable independiente, NO OBSERVARÁ UNA LÍNEA PASANDO POR EL CENTRO DE TODOS LOS DATOS COMO EN REGRESIÓN LINEAL SIMPLE, pero si puede ver que los datos observados de Biomasa (B) son muy similares a los datos estimados por el modelo.
plot(H, B, xlab = "Altura (m)", ylab = "Biomasa (kg)", cex=2)
points(H, fitted.values(Modelo_1), col="red", pch=19, cex=1.2)
Los círculos blancos representan la Biomasa observada (B), los círculos rojos la Biomasa estimada. En la primera figura se hizo un plot de “líneas”, y en ésta, un plot de “scatter” usando “points”. Puede ver que la biomasa observada y estimada son similares.
Ojo se hiciera esta última figura usando líneas verá un desorden, porque? Porque la columna de Altura (H) no está ordenada. Vea lo que pasa.
plot(H, B, xlab = "Altura (m)", ylab = "Biomasa (kg)", cex=2)
lines(H, fitted.values(Modelo_1), col="red", pch=19, cex=1.2)
Condiciones para la regresión
Todo modelo de regresión ya sea lineal simple, múltiple o no lineal, debe cumplir con una serie de “supuestos” “condicionantes” o “requisitos”. Sin este cumplimiento, las predicciones poder ser sesgadas. En muchas de las veces, los modelos no cumplen, o solo se cumplen algunos de ellos, sin embargo, pueden resolverse de diferentes formas. En este ejemplo, solo se verán si estos supuestos cumplen o no, en otros ejemplos les mostrare como se corrigen en caso de que no sean cumplidos.
1.- Varianza de residuos constante (homocedasticidad):
La homocedasticidad implica que la varianza se mantenga constante. Puede analizarse de forma gráfica representando las observaciones en un diagrama de dispersión y viendo si mantiene una homogeneidad en su dispersión a lo largo del eje X. Una forma cónica es un claro indicativo de falta de homocedasticidad. En algunos libros se menciona el test de Goldfeld-Quandt o el de Breusch-Pagan como test de hipótesis para la homocedasticidad en correlación y regresión.
library(ggplot2)
library(gridExtra)
plot1 <- ggplot(data = Muestra, aes(D, Modelo_1$ residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = Muestra, aes(H, Modelo_1$ residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2) #Graficar los dos plot, en caso de ser mas, agrgarlas## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
####
library(lmtest) # Cargar la liberria
bptest(Modelo_1) # Calcular la prueba de Breusch-Pagan##
## studentized Breusch-Pagan test
##
## data: Modelo_1
## BP = 1.8941, df = 2, p-value = 0.3879Como en este modelo son dos variables independientes, se muestran dos figuras (con D y H) para verificar la homogeneidad de varianza de los errores. En el eje de las “x” se representa la variable independiente y en el eje de las “y” los errores.
Las figuras muestran que los errores a lo largo del eje “x”, se mantiene constante.
La hipótesis nula Ho = Homogeneidad de varianza (homocedasticidad)
La hipótesis alterna Ha = No homogeneidad de varianza (heterosedasticidad)
Para probar estadísticamente Ho, aplicamos la prueba de Breusch-Pagan. El p-value = 0.6919 es muy alto, es decir; la probabilidad es muy alta que los errores se distribuyan homogéneamente.
2.- Linealidad:
La relación entre ambas variables debe ser lineal. Para comprobarlo se puede recurrir a: Graficar ambas variables a la vez (scatterplot o diagrama de dispersión), superponiendo la recta del modelo generado por regresión lineal. Calcular los residuos para cada observación acorde al modelo generado y graficarlos (scatterplot). Deben distribuirse de forma aleatoria en torno al valor 0.
plot(D, B, xlab = "Diámetro normal (cm)", ylab = "Biomasa (kg)", cex=2)
abline(lm(B~D), col="red", col.axis=1) # Colocar la linea de regresion
plot(H, B, xlab = "Altura total (m)", ylab = "Biomasa (kg)", cex=2)
abline(lm(B~H), col="red", col.axis=1) # Colocar la linea de regresion
En las figuras anteriores se observa que existe una relación lineal entre las variables predictoras (independientes) y la variable de respuesta (dependiente). La relación (lineal) más fuerte es Diámetro (D) con Biomasa (B). Se espera que el modelo tenga mejor respuesta con Diámetro (D). Una relación no lineal seria por ejemplo una tendencia de los puntos en forma de curva.
3.- Distribución Normal de los residuos:
Los residuos se tiene que distribuir de forma normal, con media igual a 0. Esto se puede comprobar con un histograma, con la distribución de cuantiles (qqnorm() + qqline()) o con un test de hipótesis de normalidad. Los valores extremos suelen ser una causa frecuente por la que se viola la condición de normalidad.
par(mfrow = c(1, 3))
hist(Modelo_1$residuals, col="#00B2EE", cex=1.8, main="Histograma de frecuencias", ylab="Frecuencia")
qqnorm(Modelo_1$residuals, col="#00B2EE")
qqline(Modelo_1$residuals)
boxplot(Modelo_1$residuals, col="#00B2EE", main="Box plot")
El histograma (Figura izquierda) muestra una distribución normal. En la figura del centro (Q-Q plot) se observa que los puntos se distribuyen alrededor de la línea, entre más cercanos a ella, los errores poseen una distribución normal y viceversa. La tercera figura (Box plot), es también utilizada para verificar la normalidad de los datos. La línea gruesa del centro representa la mediana y también el segundo cuantil (Q3), no se observan datos atípicos o “outliers”, es decir, los errores se distribuyen normal.
Para probar estadísticamente la normalidad de los errores se aplicaron tres pruebas: Shapiro-Wilk, Kolmogorov-Smirnov y la de Jarque Bera.
La hipótesis nula Ho = Normalidad de los errores
La hipótesis alterna Ha = No Normalidad de los errores
Test de Shapiro-Wilk Este test se emplea para contrastar normalidad cuando el tamaño de la muestra es menor de 50. Para muestras grandes es equivalente al test de kolmogorov-Smirnov.
library("nortest")
shapiro.test(Modelo_1$residuals) # Probar estadsiticamente la normalidad##
## Shapiro-Wilk normality test
##
## data: Modelo_1$residuals
## W = 0.95807, p-value = 0.3773Test de Kolmogorov-Smirnov (no solo para normalidad) y modificación de Lillefors El test de Kolmogorov-Smirnov permite estudiar si una muestra procede de una población con una determinada distribución (media y desviación típica), no está limitado únicamente a la distribución normal. Se ejecuta con la función ks.test().
library("nortest")
lillie.test(Modelo_1$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Modelo_1$residuals
## D = 0.095087, p-value = 0.8104Test de normalidad de Jarque-Bera Las medidas de asimetría, sobre todo el coeficiente de asimetría de Fisher, junto con las medidas de apuntamiento o curtosis se utilizan para contrastar si se puede aceptar que una distribución estadística sigue la distribución normal. El test de Jarque-Bera no requiere estimaciones de los parámetros que caracterizan la normal. El estadístico de Jarque-Bera cuantifica que tanto se desvían los coeficientes de asimetría y curtosis de los esperados en una distribución normal.
library("tseries")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoojarque.bera.test(Modelo_1$residuals)##
## Jarque Bera Test
##
## data: Modelo_1$residuals
## X-squared = 1.2658, df = 2, p-value = 0.5311Tanto el análisis gráfico como es test de hipótesis confirman la normalidad. Los p-value fueron de 0.7003 (Shapiro-Wilk), 0.6619 (Kolmogorov-Smirnov), y 0.6233 (Jarque Bera). El p-value es muy alto, es decir; la probabilidad es muy alta que los errores se distribuyan normal.
4.- Observaciones influyentes y atípicas
Hay que estudiar con detenimiento los valores atípicos o extremos ya que pueden generar una falsa correlación que realmente no existe, u ocultar una existente. (Ver descripción detallada en la sección de apuntes varios).
Veamos primero el Box plot de los residuales en sus unidades originales.
boxplot(Modelo_1$residuals)
Como puede observarse en la figura de arriba, las unidades son las mismas que la variable dependiente (kg, en este caso). Es decir; existen diferencias del observado y del estimado de hasta más menos 100 kg. Este tipo de residuales no es muy útil.
ei=Modelo_1$residuals/summary(Modelo_1)$sigma
boxplot(ei)
Aquí se muestra el “Box plot” de los residuales estandarizados No son muy importantes tampoco, pero debe hacerse en análisis.
Aqui, vamos a calcular los residuales estudentizados.
ri=rstudent(Modelo_1) # Calcular los residuales estudentizados
boxplot(ri) # Graficar los residuales estudentizados
data.frame(ei,ri) # Hacerlos en dataframe## ei ri
## 1 0.78823268 0.82806999
## 2 0.18540620 0.18670902
## 3 -0.26942657 -0.29797270
## 4 -0.01134817 -0.01146007
## 5 0.11776022 0.11937239
## 6 1.31609348 1.47582289
## 7 1.10111135 1.18226315
## 8 -0.06457969 -0.06515767
## 9 -0.10416045 -0.10521488
## 10 -0.58670112 -0.61675325
## 11 -1.45730543 -1.69283663
## 12 -1.32616621 -1.38315792
## 13 1.94633463 2.49426277
## 14 -0.16904132 -0.17714251
## 15 0.52392060 0.54417277
## 16 1.22249682 1.39169368
## 17 -0.70271646 -0.74163596
## 18 1.27206990 1.39666309
## 19 -1.00076553 -1.03334240
## 20 -1.16966803 -1.23142963
## 21 0.29787206 0.31241866
## 22 1.01343799 1.07071652
## 23 -1.13023440 -1.16797650
## 24 -0.85097501 -0.90830483
## 25 -0.94164751 -1.09358493par(new=T)
abline(h=3, col="red");abline(h=-3, col="red") # Colocar linea para identificar aquellos que ña sobrepasen, son datos atípicos.
Los residuales estudentizados son muy importantes. Este análisis debe hacerse al ajustar un modelo de regresión (vea la literatura para su cálculo). Su interpretación es la siguiente: un valor (absoluto) estudentizados por arriba de tres, indica que el dato fue capturado mal, pudo ser error de medición del mismo instrumento o de la persona, o simplemente, está fuera de lo normal.
Cuando se detecta este tipo de observaciones (en la columna llamada “ri”), es conveniente eliminarlas porque ejercen influencia en la estimación de los parámetros, en su varianza y provocan sesgos en las estimaciones. Si existirán, r indica que observación sería, ejemplo, si la línea indica que 20, entonces de la base de datos, eliminar la observación 20 y despues VOLVER A AJUSTAR EL MODELO DESDE EL PRINCIPIO HACER TODO EL PROCESO
En ambas figuras de arriba, se observa que no existen datos fuera de los “bigotes” del Box plot. Por lo tanto, no existen observaciones atípicas.
Ahora veamos, si existen datos influyentes.
library(lmtest)
library(car)
summary(influence.measures(Modelo_1))## Potentially influential observations of
## lm(formula = B ~ D + H, data = Muestra) :
##
## dfb.1_ dfb.D dfb.H dffit cov.r cook.d hat
## 3 0.10 0.12 -0.13 -0.16 1.45_* 0.01 0.22
## 13 -0.69 0.83 0.33 1.43_* 0.70 0.55 0.25influencePlot(Modelo_1)
## StudRes Hat CookD
## 11 -1.692837 0.1960621 0.2147478
## 13 2.494263 0.2465793 0.5485236
## 25 -1.093585 0.2519643 0.1330914El análisis de datos influyentes, demuestra que la observación 15 muestra influencia en la covarianza, mientras que la observación 25 influye el en ajuste del modelo (ver los “*”). Sin embargo, ninguna de estas dos observaciones rebasa el valor de 3 (“ri”) para que pudieran ser eliminados de la base de datos.
De la misma forma, en la figura de arriba, se señalan esas mismas dos observaciones. Conforme más grande es el “circulo”, la observación es más influyente. Puede notar que ninguna de estas dos observaciones (15 y 25) rebasan el valor estudentizados de 3 (eje “y”).
5.- Independencia, Autocorrelación:
Las observaciones deben ser independientes unas de otras. Esto es importante tenerlo en cuenta cuando se trata de mediciones temporales (de tiempo, no es el casos de este tipo de datos). Puede detectarse estudiando si los residuos siguen un patrón o tendencia. Otro caso frecuente es el de tener varias mediciones para un mismo sujeto. En estos casos, primero se obtiene la media de cada uno y después se ajusta el modelo empleando las medias.
Vamos a calcular el estadístico de Durbin Whatson. De forma muy general, la prueba establece que un valor cercano a 2, demuestra ausencia de autocorrelación o independencia de los residuales (errores).
La hipótesis nula Ho = Dependencia de los residuales
La hipótesis alterna Ha = Independencia de los residuales
library(car)
dwt(Modelo_1, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.0383582 1.854738 0.768
## Alternative hypothesis: rho != 0acf(Modelo_1$ residuals, col="red")
Los resultados demuestran, que el valor de Durbin Watson (DW) D-W Statistic es 1.143613, muy cercano a 2, esto demuestra la ausencia de independencia de los residuales. Es decir, se rechaza la hipotesis nula.
En forma gráfica, vemos el comportamiento de la autocorrelación de los residuales hasta después de 13 retardos (lags). Interprételo de la siguiente manera: La línea vertical en el valor de cero (eje “x”) indica la correlación de los residuales con ellos mismos, como era de esperarse, la correlación es de 1. La segunda línea vertical (en 1) es la correlación de los residuales después de un lag (retardo), esto es; bajando los valores una línea hacia abajo y haciendo correlación. Si la línea 2, rebasa la línea horizontal (puntedada), se demuestra la presencia de autocorrelación entre los residuales.
6.- Análisis de Inflación de Varianza (VIF):
Cuando en un modelo de regresion lineal existe una fuerte relación lineal entre sus variables independientes se dice que existe multicolinealidad. En esta situacion, el estimador por Mínimos Cuadrados Ordinarios (MCO) puede ofrecer resultados inestables, por lo que no se recomienda su uso y se hace necesario disponer de herramientas que permitan detectar este problema de forma adecuada.
Como medida de la misma hay varios estadísticos propuestos, los más sencillos son los coeficientes de determinación de cada variable independiente con todas las demás y, relacionados con ellos, el factor de inflación de la varianza (FIV) dada por \[VIF_{\hat{\beta }_{j}}=\frac{1}{1-R^{2}}\] y la tolerancia (T). \[Tolerancia_{\hat{\beta}_j} = \frac{1}{VIF_{\hat{\beta}_j}}\] Donde \(R^{2}\) se obtiene de la regresión del predictor \(X_j\) sobre los otros predictores. Esta es la opción más recomendada, los límites de referencia que se suelen emplear son:
VIF = 1: Ausencia total de colinealidad.
1 < VIF < 5: La regresión puede verse afectada por cierta colinialidad.
5 < VIF < 10: Causa de preocupación.
El termino tolerancia es 1/VIF por lo que los límites recomendables están entre 1 y 0.1.
Una regla empírica, citada por Kleinbaum, consiste en considerar que existen problemas de colinealidad si algún FIV es superior a 10, que corresponde a algún \(R^{2}i\) 0,9 y Ti < 0,1.
library(corrplot)## corrplot 0.84 loadedlibrary(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:car':
##
## recode## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:xts':
##
## first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncorrplot(cor(dplyr::select(Muestra, D, H)),
method = "number", tl.col = "black")
library(car)
vif(Modelo_1)## D H
## 1.5106 1.5106sqrt(vif(Modelo_1)) #Cuando se prefiere usar la desviación estándar## D H
## 1.229065 1.229065La figura anterior, demuestra que existe una baja correlación entre variables independientes, es decir; existe colinealidad, pero baja, puesto que la correlación es de solo 0.52. Al calcular el FIV se obtienen valores de 1.36894, es decir, según la regla empírica, al ser valores menores a 10, demuestran la ausencia de colinealidad.
No hay predictores que muestren una correlación lineal muy alta ni inflación de varianza.
Modelos a ajustar:
Como primer paso, vaya a https://rpubs.com/jorge_mendez/601482 para que corra el script y obtenga una muestra aleatoria de un conjunto de datos, posteriormente, ajuste los siguientes modelos.
1.- \(B = \beta_{0} + \beta_{1}D\)
2.- \(B = \beta_{0} + \beta_{1}DH\)
3.- \(B = \beta_{0} + \beta_{1}D^{2}H\)
4.- \(B = \beta_{0} + \beta_{1}lnD^{2}H\)
5.- \(B = \beta_{0} + \beta_{1}lnD^{2}H^{2}\)
6.- \(B = \beta_{0} + \beta_{1}D+ \beta_{2}D^{2}+\beta_{3}D^{3}\)
7.- \(B = \beta_{0} + \beta_{1}lnD+ \beta_{2}lnH+\beta_{3}lnD^{2}H\)
8.- \(B = \beta_{0} + \beta_{1}DH+ \beta_{2}lnH\)
9.- \(B = \beta_{0} + \beta_{1}lnD^{}+ \beta_{2}lnH\)
10.-\(lnB = \beta_{0} + \beta_{1}D^{2}H\)
11.- \(lnB = \beta_{0} + \beta_{1}lnD^{2}H\)
Bibliografia
https://www.bioestadistica.uma.es/baron/apuntes/ficheros/cap06.pdf
https://www.codecogs.com/latex/eqneditor.php?lang=es-es
Picard N., L. Saint-André y M. Henry (2012). FAO. 223 p.
Avendaño H. D. M., M. Acosta M., F. Carrillo A. y J. D. Etchevers B. (2009). Revista Fitotecnia Mexicana 32:233-238.
Návar J. (2009). Forest Ecology and Management 257:427-434.
Silva-Arredondo F. y J. Návar-Cháidez (2010). Revista Mexicana de Ciencias Forestales 1:55-62.
Romahn De la V. C. F. (1999) Relascopía. Una técnica de medición forestal. UACh. Texcoco, Estado de México. https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/trees.html https://www.delta-intkey.com/wood/es/www/rospravi.htm
https://rpubs.com/Joaquin_AR/223351
https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/trees.html datos
https://www.codecogs.com/latex/eqneditor.php?lang=es-esLatex
https://www.cienciadedatos.net/documentos/25_regresion_lineal_multiple