AI2UC1_3
Santana Celaya Alec Demian
9/6/2021
Regresion lineal simple
La regresión lineal es una técnica de modelado estadístico que se emplea para describir una variable de respuesta continua como una función de una o varias variables predictoras. Puede ayudar a comprender y predecir el comportamiento de sistemas complejos o a analizar datos experimentales, financieros y biológicos.
Fuente: https://la.mathworks.com/discovery/linear-regression.html
Importar datos
#grasas <- read.table('http://verso.mat.uam.es/~joser.berrendero/datos/EdadPesoGrasas.txt', header = TRUE)
#names(grasas)
library("readr")## Warning: package 'readr' was built under R version 3.6.3Prodmiel2018muni <- read_csv("Prodmiel2018muni.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## Nomestado = col_character(),
## Nommunicipio = col_character(),
## Volumen = col_double(),
## Valor = col_double()
## )Prodmiel2018muni <- Prodmiel2018muni[,3:4]
#View(Prodmiel2018muni)
#names(Prodmiel2018muni)
#typeof(Prodmiel2018muni)Matriz de diagramas de dispersión
#pairs(grasas)
pairs(Prodmiel2018muni)
Aqui podemos observar que hay una estrecha relación entre el volumen y el valor de la miel en los diferentes estados y municipios de México, se ve claramente que van en una diagonal en aumento lo que nos deja en claro algo obvio: cuando uno aumenta tambien el otro, esto lo podemos comprobar con lo siguiente.
Matriz de coeficientes de correlacion
En estadística, el coeficiente de correlación de Pearson es una medida de dependencia lineal entre dos variables aleatorias cuantitativas. A diferencia de la covarianza, la correlación de Pearson es independiente de la escala de medida de las variables. De manera menos formal, podemos definir el coeficiente de correlación de Pearson como un índice que puede utilizarse para medir el grado de relación de dos variables siempre y cuando ambas sean cuantitativas y continuas.
#cor(grasas)
cor(Prodmiel2018muni)## Volumen Valor
## Volumen 1.0000000 0.9920269
## Valor 0.9920269 1.0000000Se puede observar que tienen una corrrelación de caso el 100% con un 99% de correlación.
Estimacion y representacion de la recta de minimos cuadrados
El comando básico es lm (linear models). El primer argumento de este comando es una fórmula y ~ x en la que se especifica cuál es la variable respuesta o dependiente (y ) y cuál es la variable regresora o independiente (x). El segundo argumento, llamado data especifica cuál es el fichero en el que se encuentran las variables. El resultado lo guardamos en un objeto llamado regresion. Este objeto es una lista que contiene toda la información relevante sobre el análisis. Mediante el comando summary obtenemos un resumen de los principales resultados:
#regresion <- lm(grasas ~ edad, data = grasas )
#summary(regresion)
regresionMiel <- lm(Valor ~ Volumen, data = Prodmiel2018muni)
summary(regresionMiel)##
## Call:
## lm(formula = Valor ~ Volumen, data = Prodmiel2018muni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6163.6 -168.9 -121.8 20.0 7667.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 186.117 15.677 11.87 <2e-16 ***
## Volumen 38.384 0.124 309.60 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 582.9 on 1547 degrees of freedom
## Multiple R-squared: 0.9841, Adjusted R-squared: 0.9841
## F-statistic: 9.585e+04 on 1 and 1547 DF, p-value: < 2.2e-16Los parámetros de la ecuación de la recta de mínimos cuadrados que relaciona la cantidad de volumen en función del valor que vienen dados por la columna ´Estimate´ de la tabla ´Coefficients´ de la salida anterior. Por lo tanto, en este ejemplo la ecuación de la recta de mínimos cuadrados es:
\[ y = 186.117 + 38.384 x \]
Para tomar el primer valor tomamos el ´Intercept´ y X lo tomamos de la columna ´Estimate´ en la fila de la variable en este caso Volumen
Representacion de la recta de minimos cuadrados
#plot(grasas$edad, grasas$grasas, xlab="Edad", ylab = "Grasas")
#abline(regresion)
plot(Prodmiel2018muni$Volumen, Prodmiel2018muni$Valor, xlab="Volumen", ylab = "Valor")
abline(regresionMiel)
El coeficiente de determinación (es decir, el coeficiente de correlación al cuadrado) mide la bondad del ajuste de la recta a los datos. A partir de la salida anterior, vemos que su valor en este caso es Multiple R-squared: 0.9841.
Estimacion de predicciones
#nuevas.edades <- data.frame(edad = seq(30, 50))
#predict(regresion, nuevas.edades)
nuevos.volumenes <- data.frame(Volumen = seq(100, 500))
predict(regresionMiel, nuevos.volumenes)## 1 2 3 4 5 6 7 8
## 4024.491 4062.875 4101.258 4139.642 4178.026 4216.410 4254.793 4293.177
## 9 10 11 12 13 14 15 16
## 4331.561 4369.945 4408.328 4446.712 4485.096 4523.480 4561.863 4600.247
## 17 18 19 20 21 22 23 24
## 4638.631 4677.014 4715.398 4753.782 4792.166 4830.549 4868.933 4907.317
## 25 26 27 28 29 30 31 32
## 4945.701 4984.084 5022.468 5060.852 5099.236 5137.619 5176.003 5214.387
## 33 34 35 36 37 38 39 40
## 5252.771 5291.154 5329.538 5367.922 5406.306 5444.689 5483.073 5521.457
## 41 42 43 44 45 46 47 48
## 5559.841 5598.224 5636.608 5674.992 5713.376 5751.759 5790.143 5828.527
## 49 50 51 52 53 54 55 56
## 5866.911 5905.294 5943.678 5982.062 6020.445 6058.829 6097.213 6135.597
## 57 58 59 60 61 62 63 64
## 6173.980 6212.364 6250.748 6289.132 6327.515 6365.899 6404.283 6442.667
## 65 66 67 68 69 70 71 72
## 6481.050 6519.434 6557.818 6596.202 6634.585 6672.969 6711.353 6749.737
## 73 74 75 76 77 78 79 80
## 6788.120 6826.504 6864.888 6903.272 6941.655 6980.039 7018.423 7056.807
## 81 82 83 84 85 86 87 88
## 7095.190 7133.574 7171.958 7210.341 7248.725 7287.109 7325.493 7363.876
## 89 90 91 92 93 94 95 96
## 7402.260 7440.644 7479.028 7517.411 7555.795 7594.179 7632.563 7670.946
## 97 98 99 100 101 102 103 104
## 7709.330 7747.714 7786.098 7824.481 7862.865 7901.249 7939.633 7978.016
## 105 106 107 108 109 110 111 112
## 8016.400 8054.784 8093.168 8131.551 8169.935 8208.319 8246.703 8285.086
## 113 114 115 116 117 118 119 120
## 8323.470 8361.854 8400.238 8438.621 8477.005 8515.389 8553.772 8592.156
## 121 122 123 124 125 126 127 128
## 8630.540 8668.924 8707.307 8745.691 8784.075 8822.459 8860.842 8899.226
## 129 130 131 132 133 134 135 136
## 8937.610 8975.994 9014.377 9052.761 9091.145 9129.529 9167.912 9206.296
## 137 138 139 140 141 142 143 144
## 9244.680 9283.064 9321.447 9359.831 9398.215 9436.599 9474.982 9513.366
## 145 146 147 148 149 150 151 152
## 9551.750 9590.134 9628.517 9666.901 9705.285 9743.668 9782.052 9820.436
## 153 154 155 156 157 158 159 160
## 9858.820 9897.203 9935.587 9973.971 10012.355 10050.738 10089.122 10127.506
## 161 162 163 164 165 166 167 168
## 10165.890 10204.273 10242.657 10281.041 10319.425 10357.808 10396.192 10434.576
## 169 170 171 172 173 174 175 176
## 10472.960 10511.343 10549.727 10588.111 10626.495 10664.878 10703.262 10741.646
## 177 178 179 180 181 182 183 184
## 10780.030 10818.413 10856.797 10895.181 10933.565 10971.948 11010.332 11048.716
## 185 186 187 188 189 190 191 192
## 11087.099 11125.483 11163.867 11202.251 11240.634 11279.018 11317.402 11355.786
## 193 194 195 196 197 198 199 200
## 11394.169 11432.553 11470.937 11509.321 11547.704 11586.088 11624.472 11662.856
## 201 202 203 204 205 206 207 208
## 11701.239 11739.623 11778.007 11816.391 11854.774 11893.158 11931.542 11969.926
## 209 210 211 212 213 214 215 216
## 12008.309 12046.693 12085.077 12123.461 12161.844 12200.228 12238.612 12276.996
## 217 218 219 220 221 222 223 224
## 12315.379 12353.763 12392.147 12430.530 12468.914 12507.298 12545.682 12584.065
## 225 226 227 228 229 230 231 232
## 12622.449 12660.833 12699.217 12737.600 12775.984 12814.368 12852.752 12891.135
## 233 234 235 236 237 238 239 240
## 12929.519 12967.903 13006.287 13044.670 13083.054 13121.438 13159.822 13198.205
## 241 242 243 244 245 246 247 248
## 13236.589 13274.973 13313.357 13351.740 13390.124 13428.508 13466.892 13505.275
## 249 250 251 252 253 254 255 256
## 13543.659 13582.043 13620.426 13658.810 13697.194 13735.578 13773.961 13812.345
## 257 258 259 260 261 262 263 264
## 13850.729 13889.113 13927.496 13965.880 14004.264 14042.648 14081.031 14119.415
## 265 266 267 268 269 270 271 272
## 14157.799 14196.183 14234.566 14272.950 14311.334 14349.718 14388.101 14426.485
## 273 274 275 276 277 278 279 280
## 14464.869 14503.253 14541.636 14580.020 14618.404 14656.788 14695.171 14733.555
## 281 282 283 284 285 286 287 288
## 14771.939 14810.323 14848.706 14887.090 14925.474 14963.857 15002.241 15040.625
## 289 290 291 292 293 294 295 296
## 15079.009 15117.392 15155.776 15194.160 15232.544 15270.927 15309.311 15347.695
## 297 298 299 300 301 302 303 304
## 15386.079 15424.462 15462.846 15501.230 15539.614 15577.997 15616.381 15654.765
## 305 306 307 308 309 310 311 312
## 15693.149 15731.532 15769.916 15808.300 15846.684 15885.067 15923.451 15961.835
## 313 314 315 316 317 318 319 320
## 16000.219 16038.602 16076.986 16115.370 16153.753 16192.137 16230.521 16268.905
## 321 322 323 324 325 326 327 328
## 16307.288 16345.672 16384.056 16422.440 16460.823 16499.207 16537.591 16575.975
## 329 330 331 332 333 334 335 336
## 16614.358 16652.742 16691.126 16729.510 16767.893 16806.277 16844.661 16883.045
## 337 338 339 340 341 342 343 344
## 16921.428 16959.812 16998.196 17036.580 17074.963 17113.347 17151.731 17190.115
## 345 346 347 348 349 350 351 352
## 17228.498 17266.882 17305.266 17343.650 17382.033 17420.417 17458.801 17497.184
## 353 354 355 356 357 358 359 360
## 17535.568 17573.952 17612.336 17650.719 17689.103 17727.487 17765.871 17804.254
## 361 362 363 364 365 366 367 368
## 17842.638 17881.022 17919.406 17957.789 17996.173 18034.557 18072.941 18111.324
## 369 370 371 372 373 374 375 376
## 18149.708 18188.092 18226.476 18264.859 18303.243 18341.627 18380.011 18418.394
## 377 378 379 380 381 382 383 384
## 18456.778 18495.162 18533.546 18571.929 18610.313 18648.697 18687.081 18725.464
## 385 386 387 388 389 390 391 392
## 18763.848 18802.232 18840.615 18878.999 18917.383 18955.767 18994.150 19032.534
## 393 394 395 396 397 398 399 400
## 19070.918 19109.302 19147.685 19186.069 19224.453 19262.837 19301.220 19339.604
## 401
## 19377.988Inferencia en modelo de regresion lineal simple
Los datos vienen de esta ecuacion:
\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \ \ \ \ i=1,\ldots,n, \] donde los errores aleatorios ϵi son independientes con distribución normal de media 0 y varianza σ2
Intervalos de confianza
#confint(regresion)
confint(regresionMiel)## 2.5 % 97.5 %
## (Intercept) 155.36564 216.86758
## Volumen 38.14056 38.62692#confint(regresion, level=0.90)
confint(regresionMiel, level=0.90)## 5 % 95 %
## (Intercept) 160.3143 211.91889
## Volumen 38.1797 38.58779Representacion grafica de los intervalos de confianza
Bajo este modelo
Los errores típicos de los estimadores de los parámetros β0 y β1 se encuentran en la columna Std Error de la salida anterior. En el ejemplo, sus valores son 15.677 y 0.124 respectivamente.
La columna t value contiene el estadístico t, es decir, cociente entre cada estimador y su error típico. Estos cocientes son la base para llevar a cabo los contrastes H0:β0=0 y H0:β1=0
. Los correspondientes p-valores aparecen en la columna Pr(>|t|). En este caso son muy pequeños por lo que se rechazan ambas hipótesis para los niveles de significación habituales.
El estimador de la desviación típica de los errores σ aparece como Residual standard error y su valor en el ejemplo es 582.9
Los intervalos de confianza para los parámetros se obtienen con el comando confint. El parámetro level permite elegir el nivel de confianza (por defecto es 0.95):
Los intervalos de confianza para la respuesta media y los intervalos de predicción para la respuesta se pueden obtener usando el comando predict. Por ejemplo, el siguiente código calcula y representa los dos tipos de intervalos para el rango de volumen que van de 0 a 2000 años (los de predicción en rojo):
nuevos.volumenes <- data.frame(Volumen = seq(0, 2000))
# Grafico de dispersion y recta
plot(Prodmiel2018muni$Volumen, Prodmiel2018muni$Valor, xlab="Volumen", ylab = "Valor")
abline(regresionMiel)
# Intervalos de confianza de la respuesta media
# ic es una matriz con tres columnas: prediccion, las otras 2 son los extremos del intervalo
ic <- predict(regresionMiel, nuevos.volumenes, interval = "confidence")
lines(nuevos.volumenes$Volumen, ic[, 2], lty = 2) #lwr
lines(nuevos.volumenes$Volumen, ic[, 3], lty = 3) #upr
# Intervalos de prediccion
ic <- predict(regresionMiel, nuevos.volumenes, interval = "prediction")
lines(nuevos.volumenes$Volumen, ic[, 2], lty = 2, col="red") #lwr
lines(nuevos.volumenes$Volumen, ic[, 3], lty = 3, col="red") #upr
- La tabla de analisis de Varianza
#anova(regresion)
anova(regresionMiel)## Analysis of Variance Table
##
## Response: Valor
## Df Sum Sq Mean Sq F value Pr(>F)
## Volumen 1 3.2566e+10 3.2566e+10 95855 < 2.2e-16 ***
## Residuals 1547 5.2558e+08 3.3974e+05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnostico del modelo
residuos <- rstandard(regresionMiel) #y
valores.ajustados <- fitted(regresionMiel) #x
plot(valores.ajustados, residuos)
qqnorm(residuos)
qqline(residuos)