Planteamiento

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)

Solución analítica

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

Ajuste con la función `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.

Aproximación numérica

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.

Ejemplo: estimación por Mínimos Cuadrados Ordinarios

Planteamiento

Solución analítica

Ajuste con la función `lm` de `R`.

Aproximación numérica

Aplicaciones del modelo

Ejemplo: estimación por Mínimos Cuadrados Ordinarios

Planteamiento

Solución analítica

Ajuste con la función lm de R.

Aproximación numérica

Aplicaciones del modelo

Ajuste con la función `lm` de `R`.