Regresion lineal simple
Kenneth Roy Cabrera Torres
3 de octubre de 2018
Fundamentación de la regresión lineal simple.
Se tiene una variable \(x\) denominada explicativa, regresora o independiente y una variable \(y\) denominada respuesta, regresada o dependiente.
Al tener un conjunto de pares de elementos de \(x\) y \(y\) se pretende construir un modelo que relacione de manera lineal y mediante una línea recta el comportamiento medio de \(y\) dado los valores de \(x\).
Es decir se trata de hallar los valores de \(\beta_0\) y \(\beta_1\) tal que puedan representar el siguiente sistema de ecuaciones:
\[ \begin{eqnarray} y_1 & = & \beta_0 + \beta_1 x_1 + \epsilon_1 \\ y_2 & = & \beta_0 + \beta_1 x_2 + \epsilon_2 \\ \vdots & = & \vdots \\ y_i & = & \beta_0 + \beta_1 x_i + \epsilon_i \\ \vdots & = & \vdots \\ y_n & = & \beta_0 + \beta_1 x_n + \epsilon_n \\ \end{eqnarray} \]
Donde el supesto más utilizado sobre los errores es:
\[ e_i \sim \mathcal{N}(0, \sigma^2) \] Donde \(\sigma^2\) el la varianza de los errores aleatorios.
En forma matricial se puede escribir:
\[ Y = X\beta + \epsilon \]
Donde:
\[ Y = \left[ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_i \\ \vdots \\ y_n \end{matrix} \right], \qquad X = \left[ \begin{matrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_i \\ \vdots \\ 1 & x_n \end{matrix} \right], \qquad \beta = \left[ \begin{matrix} \beta_0 \\ \beta_1 \\ \end{matrix} \right], \qquad \epsilon = \left[ \begin{matrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_i \\ \vdots \\ \epsilon_n \end{matrix} \right] \]
Recordemos que la solución mínimo cuadrática para obtener los valores de \(\beta\) es:
\[ X'Y = X'X\beta \]
Al despejar el valor de \(\beta\) se tiene que:
\[ \begin{eqnarray} X'Y & = & X'X\beta \\ (X'X)^{-1}X'Y & = & (X'X)^{-1}(X'X)\beta \\ (X'X)^{-1}X'Y & = & I\beta \\ (X'X)^{-1}X'Y & = & \beta \\ \end{eqnarray} \]
Es decir que el estimador mínimo cuadrático para \(\beta\) es:
\[ \hat{\beta} = (X'X)^{-1}X'Y \]
Luego el valor estimado de \(Y\), es decir \(\hat{Y}\) es:
\[ \hat{Y} = X\hat{\beta} = X(X'X)^{-1}X'Y \]
A la matrix que toma los valores de \(Y\) y obtiene los valores de \(\hat{Y}\) se le denomina la matriz “hat” o sombrero y se denota por \(H\).
\[ H = X(X'X)^{-1}X' \]
Esta matriz tiene propiedades para identificar valores influenciales de los datos de \(X\).
Ejemplo con el pez Longnose dace (Rhinichthys cataractae).
El pecesito Langnose dace (Rhinichthys cataractae) habita en la corrientes de agua de Norteamérica.
En la página del John H. McDonald se encuentra una base de datos que el investigador menciona que fue tomada del sito del Maryland Biological Stream Suervey.
El MBSS realiza unos muestreos siguiendo un protocolo establecida para determinar la abundancia de peces en un tramo de la quebrada de 75 metros.
La base de datos tiene las siguientes columanas:
- quebrada: Nombre de la quebrada muestreada.
- num_indiv: Número de individuos contados en un tramo de 75 m.
- area: El área (en acres) drenada por la quebrada hasta el punto de muestreo.
- oxi_dis: El oxigeno disuelto (en mg/litro).
- prof_max: Profundidad máxima (en cm) del tramo de 75 metros.
- nitrato: Concentración de nitrato en (mg/litro).
- sulfato: Concentración de sulfato en (mg/litro).
- temperatura: La temperatura (en grados centígrados) de la quebrada en el fecha de muestreo.
El interés es conocer e identificar las variables que influencian en la abundancia del pecesito.
Y constuir un modelo que permita estimar la abundancia esperada dada las características ambientales de la quebrada.
library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.0.0 ✔ purrr 0.2.5
## ✔ tibble 1.4.2 ✔ dplyr 0.7.6
## ✔ tidyr 0.8.1 ✔ stringr 1.3.1
## ✔ readr 1.1.1 ✔ forcats 0.3.0## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()Lectura de la base de datos.
longnose <- read.table("datos_Longnose.txt", header = TRUE)Construcción del área drenada en Km\(^2\).
Hay 0.00404685642 Km\(^2\) en un acre.
longnose <- longnose %>%
mutate(areaKm2 = area * 0.00404685642)longnose## quebrada num_indiv area oxi_dis prof_max nitrato sulfato temperatura areaKm2
## 1 BASIN_RUN 13 2528 9.6 80 2.28 16.75 15.3 10.230453
## 2 BEAR_BR 12 3333 8.5 83 5.34 7.74 19.4 13.488172
## 3 BEAR_CR 54 19611 8.3 96 0.99 10.92 19.5 79.362901
## 4 BEAVER_DAM_CR 19 3570 9.2 56 5.44 16.53 17.0 14.447277
## 5 BEAVER_RUN 37 1722 8.1 43 5.66 5.91 19.3 6.968687
## 6 BENNETT_CR 2 583 9.2 51 2.26 8.81 12.9 2.359317
## 7 BIG_BR 72 4790 9.4 91 4.10 5.65 16.7 19.384442
## 8 BIG_ELK_CR 164 35971 10.2 81 3.20 17.53 13.8 145.569472
## 9 BIG_PIPE_CR 18 25440 7.5 120 3.53 8.20 13.7 102.952027
## 10 BLUE_LICK_RUN 1 2217 8.5 46 1.20 10.85 14.3 8.971881
## 11 BROAD_RUN 53 1971 11.9 56 3.25 11.12 22.2 7.976354
## 12 BUFFALO_RUN 16 12620 8.3 37 0.61 18.87 16.8 51.071328
## 13 BUSH_CR 32 19046 8.3 120 2.93 11.31 18.0 77.076427
## 14 CABIN_JOHN_CR 21 8612 8.2 103 1.57 16.09 15.0 34.851527
## 15 CARROLL_BR 23 3896 10.4 105 2.77 12.79 18.4 15.766553
## 16 COLLIER_RUN 18 6298 8.6 42 0.26 17.63 18.2 25.487102
## 17 CONOWINGO_CR 112 27350 8.5 65 6.95 14.94 24.1 110.681523
## 18 DEAD_RUN 25 4145 8.7 51 0.34 44.93 23.0 16.774220
## 19 DEEP_RUN 5 1175 7.7 57 1.30 21.68 21.8 4.755056
## 20 DEER_CR 26 8297 9.9 60 5.26 6.36 19.1 33.576768
## 21 DORSEY_RUN 8 7814 6.8 160 0.44 20.24 22.6 31.622136
## 22 FALLS_RUN 15 1745 9.4 48 2.19 10.27 14.3 7.061764
## 23 FISHING_CR 11 5046 7.6 109 0.73 7.10 19.0 20.420437
## 24 FLINTSTONE_CR 11 18943 9.2 50 0.25 14.21 18.5 76.659601
## 25 GREAT_SENECA_CR 87 8624 8.6 78 3.37 7.51 21.3 34.900090
## 26 GREENE_BR 33 2225 9.1 41 2.30 9.72 20.5 9.004256
## 27 GUNPOWDER_FALLS 22 12659 9.7 65 3.30 5.98 18.0 51.229155
## 28 HAINES_BR 98 1967 8.6 50 7.71 26.44 16.8 7.960167
## 29 HAWLINGS_R 1 1172 8.3 73 2.62 4.64 20.5 4.742916
## 30 HAY_MEADOW_BR 5 639 9.5 26 3.53 4.46 20.1 2.585941
## 31 HERRINGTON_RUN 1 7056 6.4 60 0.25 9.82 24.5 28.554619
## 32 HOLLANDS_BR 38 1934 10.5 85 2.34 11.44 12.0 7.826620
## 33 ISRAEL_CR 30 6260 9.5 133 2.41 13.77 21.0 25.333321
## 34 LIBERTY_RES 12 424 8.3 62 3.49 5.82 20.2 1.715867
## 35 LITTLE_ANTIETAM_CR 24 3488 9.3 44 2.11 13.37 24.0 14.115435
## 36 LITTLE_BEAR_CR 6 3330 9.1 67 0.81 8.16 14.9 13.476032
## 37 LITTLE_CONOCOCHEAGUE_CR 15 2227 6.8 54 0.33 7.60 24.0 9.012349
## 38 LITTLE_DEER_CR 38 8115 9.6 110 3.40 9.22 20.5 32.840240
## 39 LITTLE_FALLS 84 1600 10.2 56 3.54 5.69 19.5 6.474970
## 40 LITTLE_GUNPOWDER_R 3 15305 9.7 85 2.60 6.96 17.5 61.937138
## 41 LITTLE_HUNTING_CR 18 7121 9.5 58 0.51 7.41 16.0 28.817665
## 42 LITTLE_PAINT_BR 63 5794 9.4 34 1.19 12.27 17.5 23.447486
## 43 MAINSTEM_PATUXENT_R 239 8636 8.4 150 3.31 5.95 18.1 34.948652
## 44 MEADOW_BR 234 4803 8.5 93 5.01 10.98 24.3 19.437051
## 45 MILL_CR 6 1097 8.3 53 1.71 15.77 13.1 4.439401
## 46 MORGAN_RUN 76 9765 9.3 130 4.38 5.74 16.9 39.517553
## 47 MUDDY_BR 25 4266 8.9 68 2.05 12.77 17.0 17.263889
## 48 MUDLICK_RUN 8 1507 7.4 51 0.84 16.30 21.0 6.098613
## 49 NORTH_BR 23 3836 8.3 121 1.32 7.36 18.5 15.523741
## 50 NORTH_BR_CASSELMAN_R 16 17419 7.4 48 0.29 2.50 18.0 70.492192
## 51 NORTHWEST_BR 6 8735 8.2 63 1.56 13.22 20.8 35.349291
## 52 NORTHWEST_BR_ANACOSTIA_R 100 22550 8.4 107 1.41 14.45 23.0 91.256612
## 53 OWENS_CR 80 9961 8.6 79 1.02 9.07 21.8 40.310737
## 54 PATAPSCO_R 28 4706 8.9 61 4.06 9.90 19.7 19.044506
## 55 PINEY_BR 48 4011 8.3 52 4.70 5.38 18.9 16.231941
## 56 PINEY_CR 18 6949 9.3 100 4.57 17.84 18.6 28.121605
## 57 PINEY_RUN 36 11405 9.2 70 2.17 10.17 23.6 46.154397
## 58 PRETTYBOY_BR 19 904 9.8 39 6.81 9.20 19.2 3.658358
## 59 RED_RUN 32 3332 8.4 73 2.09 5.50 17.7 13.484126
## 60 ROCK_CR 3 575 6.8 33 2.47 7.61 18.0 2.326942
## 61 SAVAGE_R 106 29708 7.7 73 0.63 12.28 21.4 120.224011
## 62 SECOND_MINE_BR 62 2511 10.2 60 4.17 10.75 17.7 10.161656
## 63 SENECA_CR 23 18422 9.9 45 1.58 8.37 20.1 74.551189
## 64 SOUTH_BR_CASSELMAN_R 2 6311 7.6 46 0.64 21.16 18.5 25.539711
## 65 SOUTH_BR_PATAPSCO 26 1450 7.9 60 2.96 8.84 18.6 5.867942
## 66 SOUTH_FORK_LINGANORE_CR 20 4106 10.0 96 2.62 5.45 15.4 16.616392
## 67 TUSCARORA_CR 38 10274 9.3 90 5.45 24.76 15.0 41.577403
## 68 WATTS_BR 19 510 6.7 82 5.25 14.19 26.5 2.063897Matriz de diagramas de dispersión.
## Función auxiliar para constuir un histograma en la diagonal de
## la matriz de dispersión tomado de la ayuda de comando `pairs()`.
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}pairs(longnose[, -1], panel = panel.smooth,
cex = 1.5, pch = 24, bg = "light blue",
diag.panel = panel.hist, cex.labels = 2, font.labels = 2)
En particula se examinará el comportamiento de la abundancia del pecesito tomando en cuenta el área drenada hasta el sitio de muestreo.
Diagrama de dispersión de la abundancia vs el área drenada.
ggplot(longnose, aes(areaKm2, num_indiv)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_smooth(method = lm, se = FALSE, col = "red")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Estimación de los parámetros de regresión lineal simple.
Y <- as.matrix(longnose[, 2])
X <- as.matrix(cbind(1, longnose[, "areaKm2"]))
beta <- solve(t(X) %*% X) %*% t(X) %*% Y
beta ## [,1]
## [1,] 22.818155
## [2,] 0.522332Estimación del valor esperado de abundancia por cada quebrada.
Y_est <- X %*% beta
Y_est## [,1]
## [1,] 28.16185
## [2,] 29.86346
## [3,] 64.27194
## [4,] 30.36443
## [5,] 26.45812
## [6,] 24.05050
## [7,] 32.94327
## [8,] 98.85375
## [9,] 76.59330
## [10,] 27.50446
## [11,] 26.98446
## [12,] 49.49435
## [13,] 63.07764
## [14,] 41.02222
## [15,] 31.05353
## [16,] 36.13088
## [17,] 80.63066
## [18,] 31.57987
## [19,] 25.30187
## [20,] 40.35638
## [21,] 39.33541
## [22,] 26.50674
## [23,] 33.48440
## [24,] 62.85992
## [25,] 41.04759
## [26,] 27.52137
## [27,] 49.57678
## [28,] 26.97600
## [29,] 25.29553
## [30,] 24.16887
## [31,] 37.73315
## [32,] 26.90625
## [33,] 36.05056
## [34,] 23.71441
## [35,] 30.19110
## [36,] 29.85712
## [37,] 27.52559
## [38,] 39.97166
## [39,] 26.20024
## [40,] 55.16991
## [41,] 37.87054
## [42,] 35.06553
## [43,] 41.07295
## [44,] 32.97075
## [45,] 25.13700
## [46,] 43.45944
## [47,] 31.83564
## [48,] 26.00366
## [49,] 30.92670
## [50,] 59.63848
## [51,] 41.28222
## [52,] 70.48441
## [53,] 43.87374
## [54,] 32.76571
## [55,] 31.29662
## [56,] 37.50697
## [57,] 46.92607
## [58,] 24.72903
## [59,] 29.86135
## [60,] 24.03359
## [61,] 85.61501
## [62,] 28.12591
## [63,] 61.75863
## [64,] 36.15836
## [65,] 25.88317
## [66,] 31.49743
## [67,] 44.53536
## [68,] 23.89619Simulación de regresión lineal simple.
Para tener una estimación del temaño de muestra adecuado se requiere suponder loa valores “verdaderos” de \(\beta_0\) y \(\beta_1\) como también un valor de la varianza del error aleatorio \(\sigma^2\).
También se necesita determinar si el valor de \(X\) es aleatorio o son valores fijos.
En este contexto en particular los valores de \(x\) son aleatorios.
Distribución de los valores aleatorios de X.
hist(X[, 2], nclass = 18,
main = "Histograma de frecuencia del área drenada",
freq = FALSE)
lines(density(X[, 2]), col = "blue", lwd = 2)
Para la simulación se tomarán valores aleatorios de los datos originales con reemplazamiento de las áreas drenadas.
Determinación de los valores de la desviación estándar de los errores.
segmentosX <- cut(X[, 2], seq(0, 160, 10))
desv_estandard <- tapply(Y, segmentosX, sd)
desv_estandard## (0,10] (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,80] (80,90] (90,100]
## 27.152590 54.263401 19.737112 77.230522 24.846194 4.242641 NA 16.932218 NA NA
## (100,110] (110,120] (120,130] (130,140] (140,150] (150,160]
## NA NA NA NA NA NAdesv_est_estim <- mean(desv_estandard, na.rm = TRUE)
desv_est_estim## [1] 32.05781Parámetros de simulación.
Se tiene entonces que:
- Los valores de \(x\) son aleatorios.
- \(\beta_0\) = 23
- \(\beta_1\) = 0.5
- \(\sigma\) = 32, luego \(\sigma^2\) = 1024.
beta0 <- 23
beta1 <- 0.5
sigma <- 32
sigma2 <- sigma ^ 2# Para hacer la simulación reproducible.
# set.seed(0)tam_muestra <- 50x <- sort(sample(X[, 2], tam_muestra, replace = TRUE))
y <- beta0 + beta1 * x
e <- rnorm(tam_muestra, 0, sigma)
y <- y + e
y <- as.integer(round(y, 0))
# Los valores de y < 0 se remplazarán por 0
y[y < 0] <- 0
datos <- data.frame(x = x, y = y)plot(datos, las = 1)
Y <- as.matrix(datos$y)
X <- as.matrix(cbind(1, datos$x))betas_est <- solve(t(X) %*% X) %*% t(X) %*% Y
betas_est## [,1]
## [1,] 23.6837518
## [2,] 0.4606606est_param <- function(datos) {
Y <- as.matrix(datos[, "y"])
X <- as.matrix(cbind(1 , datos[, "x"]))
beta_est <- solve(t(X) %*% X) %*% t(X) %*% Y
return(beta_est)
}est_param(datos)## [,1]
## [1,] 23.6837518
## [2,] 0.4606606gen_datos <- function(tam_muestra, beta0, beta1, sigma, x_orig){
x <- sort(sample(x_orig, tam_muestra, replace = TRUE))
y <- beta0 + beta1 * x
e <- rnorm(tam_muestra, 0, sigma)
y <- y + e
y <- as.integer(round(y, 0))
# Los valores de y < 0 se remplazarán por 0
y[y < 0] <- 0
datos <- cbind(x = x, y = y)
return(datos)
}numsim <- 5000
tam_muestra <- 10
vector_tam_m <- rep(tam_muestra, numsim)
datosT <- lapply(vector_tam_m, gen_datos,
beta0 = beta0,
beta1 = beta1,
sigma = sigma,
x_orig = longnose$areaKm2)datosT[1:3]## [[1]]
## x y
## [1,] 1.715867 0
## [2,] 13.476032 0
## [3,] 16.231941 34
## [4,] 19.044506 13
## [5,] 19.384442 17
## [6,] 28.121605 68
## [7,] 31.622136 60
## [8,] 32.840240 46
## [9,] 32.840240 14
## [10,] 145.569472 76
##
## [[2]]
## x y
## [1,] 4.439401 1
## [2,] 8.971881 0
## [3,] 10.230453 32
## [4,] 19.437051 6
## [5,] 35.349291 99
## [6,] 46.154397 66
## [7,] 61.937138 45
## [8,] 76.659601 63
## [9,] 102.952027 46
## [10,] 102.952027 52
##
## [[3]]
## x y
## [1,] 6.968687 6
## [2,] 7.061764 34
## [3,] 7.826620 51
## [4,] 16.774220 3
## [5,] 23.447486 93
## [6,] 25.487102 0
## [7,] 34.900090 0
## [8,] 46.154397 13
## [9,] 76.659601 77
## [10,] 120.224011 106parametros <- sapply(datosT, est_param)
parametros[, 1:5]## [,1] [,2] [,3] [,4] [,5]
## [1,] 17.0578387 21.0705372 13.3631139 25.0368123 7.1233477
## [2,] 0.4618549 0.4248598 0.6822603 0.4955086 0.7652985hist(parametros[1, ])
hist(parametros[2, ])
Determinación del grado de precisión.
nivel_significancia <- 0.05
IC <- quantile(parametros[2, ],
c(nivel_significancia/2, (1-nivel_significancia/2)))
IC## 2.5% 97.5%
## -0.4599897 1.3014837IC - beta1## 2.5% 97.5%
## -0.9599897 0.8014837