The dataset

The dataset has 6 variables and 28 observations. Next we present the first six observations of the dataset.

Next we present the distinct values of the covariates considered in the experiment.

library(dplyr)
datos %>% select(n_balls, d_balls, freq) %>% distinct()
## # A tibble: 4 x 3
##   n_balls d_balls  freq
##     <dbl>   <dbl> <dbl>
## 1       1      10    25
## 2       8       5    25
## 3      39       3    25
## 4      39       3    30

First Order

En esta sección se ajusta un modelo de primer orden para explicar la variable respuesta yield en función de las variables n_balls, d_balls y freq.

Para crear el modelo se usa el siguiente código.

library(rsm)
m1 <- rsm(yield ~ FO(n_balls, d_balls, freq), data=datos)
summary(m1)
## 
## Call:
## rsm(formula = yield ~ FO(n_balls, d_balls, freq), data = datos)
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) -68.87766   47.15034  -1.4608     0.157    
## n_balls      -3.05260    0.31703  -9.6286 1.024e-09 ***
## d_balls     -17.85697    1.53121 -11.6620 2.253e-11 ***
## freq         11.25000    1.92070   5.8572 4.845e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.8657, Adjusted R-squared:  0.8489 
## F-statistic: 51.58 on 3 and 24 DF,  p-value: 1.295e-10
## 
## Analysis of Variance Table
## 
## Response: yield
##                            Df  Sum Sq Mean Sq F value    Pr(>F)
## FO(n_balls, d_balls, freq)  3 28544.0    9515  51.583 1.295e-10
## Residuals                  24  4426.9     184                  
## Lack of fit                 0     0.0    -Inf                  
## Pure error                 24  4426.9     184                  
## 
## Direction of steepest ascent (at radius 1):
##    n_balls    d_balls       freq 
## -0.1431471 -0.8373757  0.5275517 
## 
## Corresponding increment in original units:
##    n_balls    d_balls       freq 
## -0.1431471 -0.8373757  0.5275517

Usando los resultados de la primera parte de la salida anterior se puede escribir la superficie de respuesta para yield así:

\[ \widehat{yield} = -68.88 - 3.05 n_{balls} - 17.86 d_{balls} + 11.25 freq \]

La expresión anterior depende de 3 variables y es posible dibujar usando curvas de nivel. A continuación se muestra el código para dibujar las curvas de nivel de yield en función de 2 variables. La variable fija en cada panel se fijó en su valor promedio.

par(mfrow=c(1, 3))
contour(m1, ~ n_balls + d_balls + freq, image = TRUE)

Los colores verdes representan valores menores de yield mientras que los colores naranja representan valores mayores de yield.

  • Del panel de la izquierda se observa que valores menores de n_balls y d_balls aumentan a yield.
  • De los páneles central y derecho se observa que freq altas y valores pequeños de n_balls y d_balls aumentan a yield.
  • Estos resultados eran de esperarse por los signos de los coeficientes en la expresión matemática de arriba.

A continuación se muestra como obtener las predicciones \(\hat{y}\) para la variable yield usando el modelo anterior, adicionalmente se calcula el coeficiente de correlación de Pearson entre \(y\) e \(\hat{y}\).

y1 <- predict(m1)
cor(y1, datos$yield)
## [1] 0.9304477

De la salida anterior se observa un valor alto correlación de 93.04%.

A continuación se muestra un gráfico de dispersión entre el valor predicho de yield versus su valor observado en el experimento. La recta a 45 grados es una recta de referencia, un modelo de predicción perfecta mostraría todos los puntos sobre esa línea de referencia.

plot(x=datos$yield, y=y1, pch=19, las=1, xlab='Observed yield', ylab='Predicted yield')
abline(a=0, b=1, lty='dashed')

Pure Quadratic Order

En esta sección se ajusta un modelo de segundo orden para explicar la variable respuesta yield en función de las variables n_balls, d_balls y freq.

Para crear el modelo se usa el siguiente código.

m2 <- rsm(yield ~ PQ(n_balls, d_balls) + freq, data=datos)
## Warning in rsm(yield ~ PQ(n_balls, d_balls) + freq, data = datos): No FO() terms in model; cannot use RSM methods
## An 'lm' object has been returned.
summary(m2)
## 
## Call:
## rsm(formula = yield ~ PQ(n_balls, d_balls) + freq, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.750  -5.979   1.333   2.500  31.250 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -1.556e+02  4.808e+01  -3.237  0.00351 ** 
## PQ(n_balls, d_balls)n_balls^2 -5.085e-02  5.787e-03  -8.787 5.77e-09 ***
## PQ(n_balls, d_balls)d_balls^2 -9.483e-01  8.498e-02 -11.159 5.54e-11 ***
## freq                           1.125e+01  1.921e+00   5.857 4.85e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.58 on 24 degrees of freedom
## Multiple R-squared:  0.8657, Adjusted R-squared:  0.8489 
## F-statistic: 51.58 on 3 and 24 DF,  p-value: 1.295e-10

A continuación se muestra el código para dibujar las curvas de nivel de yield en función de 2 variables. La variable fija en cada panel se fijó en su valor promedio.

par(mfrow=c(1, 3))
contour(m2, ~ n_balls + d_balls + freq, image = TRUE)

A continuación se muestra como obtener las predicciones \(\hat{y}\) para la variable yield usando el modelo anterior, adicionalmente se calcula el coeficiente de correlación de Pearson entre \(y\) e \(\hat{y}\).

y2 <- predict(m2)
cor(y2, datos$yield)
## [1] 0.9304477

A continuación se muestra un gráfico de dispersión entre el valor predicho de yield versus su valor observado en el experimento.

plot(x=datos$yield, y=y2, pch=19, las=1, xlab='Observed yield', ylab='Predicted yield')
abline(a=0, b=1, lty='dashed')