Javier Santibáñez
1 de febrero de 2018
Se tienen 20 observaciones de un proceso de destilación química de oxígeno a distintos porcentajes de hidrocarburos presentes en el condensador principal de la unidad de destilación. La variable respuesta es la pureza (%) del oxígeno producido.
El objetivo es ajustar un modelo RLS al conjunto de datos.
file_url <- 'http://sigma.iimas.unam.mx/jsantibanez/Cursos/Regresion/2018_2/ejemplos/ejemplo1.csv'
datos <- read.csv(file_url, header = T)
Recordemos que los estimadores de MCO son
\[\begin{align*} \hat{\beta}_1 &= \frac{S_{xy}}{S_{xx}} = \frac{\sum(x_i - \bar{x}_n)(y_i - \bar{y}_n)}{\sum(x_i - \bar{x}_n)^2}\\ \hat{\beta}_0 &= \bar{y}_n - \hat{\beta}_1 \bar{x}_n \end{align*}\]x <- datos$Concentracion
y <- datos$Pureza
Sxx <- sum((x - mean(x))^2)
Sxy <- sum((x - mean(x)) * (y - mean(y)))
b1_h <- Sxy / Sxx
b0_h <- mean(y) - b1_h * mean(x)
b0_h
## [1] 74.28331
b1_h
## [1] 14.94748
El valor de la suma de cuadrados residual es
q2_min <- sum((y - b0_h - b1_h * x)**2)
Aunque el método de mínimos cuadrados ordinarios no proporciona un estimador para \(\sigma^2\), un estimador razonable está dado por
\[\begin{equation*} \hat{\sigma}^2 = \frac{1}{n-2} \sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1x_i)^2 \end{equation*}\]Entonces, la estimación de la varianza para este conjunto de datos es
s_h <- q2_min / (dim(datos)[1] - 2)
s_h
## [1] 1.180545
lm
de R
.La función de R
que podemos utilizar para ajustar modelos de regresión lineal (simple o múltiple) es lm
. Para consultar la documentación relacionada con esta función podemos utilizar el siguiente código:
help('lm')
Una vez revisada la información de la función lm
, ajustamos un modelo RLS al conjunto de datos con el siguiente código
ajuste_lm <- lm(y ~ x)
summary(ajuste_lm)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.83029 -0.73334 0.04497 0.69969 1.96809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 74.283 1.593 46.62 < 2e-16 ***
## x 14.947 1.317 11.35 1.23e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.087 on 18 degrees of freedom
## Multiple R-squared: 0.8774, Adjusted R-squared: 0.8706
## F-statistic: 128.9 on 1 and 18 DF, p-value: 1.227e-09
Podemos oservar que las estimaciones de \(\beta_0\) y \(\beta_1\) coinciden con la solución analítica de la sección anterior. El método que usa la función lm
aproxima numéricamente a la solución de las ecuaciones normales.
Como obtuvimos las mismas estimaciones de \(\hat{\beta}_0\) y \(\hat{\beta}_1\), el valor mínimo de la suma de cuadrados de los residuos y \(\hat{\sigma}^2\) son los mismos que en la sección anterior.
Podemos utilizar el siguiente código para calcular la suma de cuadrados a partir del resultado del ajuste con la función lm
. La función residuals
nos permite obtener los residuos del modelo ajustado.
Finalmente, podemos obtener los estimadores MCO de \(\beta_0\) y \(\beta_1\) optimizando numéricamente la función suma de cuadrados de los residuos. Primero debemos declarar la función que calcula la suma de cuadrados para valores de \(\beta_0\) y \(\beta_1\) dados
q2 <- function(beta)
sum((y - beta[1] - beta[2]*x) ** 2)
El argumento de la función q2
es un vector con dos entradas, la primera corresponde al valor de \(\beta_0\) y la segunda a \(\beta_1\), esto para permitir la optimización numérica utilizando la funicón nlm
de R
. De igual manera, se puede consultar la documentación relacaionada con la función nlm
con el código help('nlm')
.
Por el momento no nos interesan todos los argumentos de la función nlm
nos basta con indicar la función a minimizar y un par de valores iniciales para \(\beta_0\) y \(\beta_1\).
ajuste_nlm <- nlm(q2, c(0.5, 0.5))
ajuste_nlm
## $minimum
## [1] 21.24982
##
## $estimate
## [1] 74.28210 14.94846
##
## $gradient
## [1] 0.00000e+00 -2.85197e-09
##
## $code
## [1] 1
##
## $iterations
## [1] 6
De todos los resultados que proporciona nlm
nos interesan tres:
minimum
: valor mínimo de la función q2
,estimate
: los estimadores de MCO de \(\beta_0\) y \(\beta_1\),code
: el valor 1 nos indica la convergencia del algoritmo iterativo.Observamos nuevamente que los estimadores obtenidos son los mismos que en las secciones anteriores.
¿Cuál es la pureza esperada con una concentración de 1.50%?
x_new <- 1.5
y_new <- b0_h + b1_h * x_new
y_new
## [1] 96.70453