Tema: AIR QUALITY DATA SET

Curso: Análisis de regresión

Profesor: Jesús Gamboa

Universidad Nacional Agraria La Molina

SOBRE LOS DATOS

Los datos contienen 9358 filas de respuestas promediadas por hora de una matriz de 5 sensores químicos que tenían como finalidad evaluar la calidad del aire. El dispositivo estaba ubicado en el campo en un área significativamente contaminada, al nivel de la carretera, dentro de una ciudad italiana. Los datos se registraron desde marzo de 2004 hasta febrero de 2005 (un año). Se registraron valores de monóxido de carbono (CO), hidrocarburos no metálicos, temperaturas, humedades relativas, etc.

Para más información y adquirir los datos, click en el enlace:https://archive.ics.uci.edu/ml/datasets/Air+Quality .

Fuente:

Saverio De Vito ( ), ENEA - Agencia Nacional de Nuevas Tecnologías, Energía y Desarrollo Económico Sostenible

PAQUETES UTILIZADOS

library(readxl)
library(dplyr)
library(TeachingSampling)
library(broom)
library(ggplot2)
library(kableExtra)
library(DT)

LECTURA DE DATOS

base <- read_xlsx("AirQualityUCI.xlsx")

LIMPIEZA DE DATOS

datos <- base |> 
  filter(base$`CO(GT)`>-200 & base$`NO2(GT)`>0)

SELECCION DE MUESTRA

Realizamos un pequeño muestreo aleatorio simple sin reemplazo, con un tamaño de muestra de 500. Muestra con la que se realizaran las estimaciones.

N<-nrow(datos) 
n<-500
RNGkind(sample.kind="Rounding")
set.seed(20)
m1<-S.SI(N,n)
mues<-datos[m1,]
mues |> datatable(option=list(scrollX=T))

Establecemos como variables estocásticas:

  • Variable X: CO(GT) Concentración real promediada por hora de CO en mg/m^3.
  • Variable Y: NO2(GT) Concentración de NO2 promediada por hora real en microg/m^3.
attach(mues)
x <- `CO(GT)`
y <- `NO2(GT)`

GRAFICA ENTRE LAS VARIABLES ANALIZADAS

p <- ggplot(mues, aes(x, y)) + geom_point() + geom_smooth(method=lm)
p
## `geom_smooth()` using formula 'y ~ x'

DESARROLLO DE PROBLEMAS

1. Para dos variables cuantitativas (pueden ser las mismas que el trabajo corto 1), estimar los coeficientes de la ecuación de regresión utilizando las fórmulas, verificar con el comando coef e interpretar los valores obtenidos. En caso no sea interpretable, explicar.

s2_x <- (sum((x-mean(x))^2)/n) |> print()
## [1] 2.097111
s2_y <- (sum((y-mean(y))^2)/n) |> print()
## [1] 2325.418
s_xy <- (sum((x-mean(x))*(y-mean(y)))/n) |> print()
## [1] 48.72111
B1 <- (s_xy/s2_x) |> print()
## [1] 23.23249
Bo <- (mean(y)-B1*mean(x)) |> print()
## [1] 63.46834
VERIFICACIÓN DE LOS RESULTADOS DE \(\beta_0\) Y \(\beta_1\)
lm(y ~ x) -> modelo
modelo |> tidy()
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     63.5      2.89      22.0 2.51e-75
## 2 x               23.2      1.07      21.7 3.79e-74
modelo$coefficients
## (Intercept)           x 
##    63.46834    23.23249

Los datos coiciden.

Interpretación

Consideramos que es posible la concentracion de monoxido de carbono igual a 0 ya que en los datos no filtrados, encontramos valores muy proximos al valor de 0 e inclusive -200.

Si \(x=0\), entonces \(y=\beta_o\)

Para el valor de \(\beta_0=63.47\),\(\beta_0\) es la media de la concentracion del dioxido de nitrogeno cuando no hay concentracion de monoxido de carbono.

Para el valor de \(\beta_1\), por cada unidad de incremento de monoxido de carbono, la concentracion de dioxido de nitrogeno incrementa en \(\beta_1\) = 23.23.

2. Probar si β1 = 0 vs β1 ≠ 0, usando el comando aov y summary y añadir una gráfica como la realizada en clase:

Planteamos la prueba de hipotesis.

\[y=\overbrace{\beta_0+\beta_1.x}^{regresion}+\overbrace{\epsilon}^{error}\\H_0:\beta_1=0 : No\: hay\: influencia\: lineal\: de\: x\: sobre\: y.\\H_1:\beta_1≠0;\:Hay\: influencia\: lineal\: de\: x\: sobre\: y.\] \[H_0: \beta_1 = 0 (Variabilidad\: explicada\:por\:la\:Regresion \leqslant Variabilidad\:explicada\:por\:el\:error)\\H_1: \beta_1 ≠ 0 (Variabilidad\: explicada\:por\:la\:Regresion > Variabilidad\:explicada\:por\:el\:error) \]

modelo |> aov() |> summary()
##              Df Sum Sq Mean Sq F value Pr(>F)    
## x             1 565956  565956   472.3 <2e-16 ***
## Residuals   498 596753    1198                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qf(0.90,1,498) #cae en la zona de rechazo
## [1] 2.715637
pf(521.8,1,498,lower.tail=FALSE) #El area es mucho menor
## [1] 1.540484e-79

Graficamente