1. Datos.

Los datos son tomados del libro Diseño y análisis de experimentos de Douglas Montgomery. Página 246.

A = c(-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1)
B = c(-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1)
C = c(-1,-1,-1,-1,+1,+1,+1,+1,-1,-1,-1,-1,+1,+1,+1,+1)
D = c(-1,-1,-1,-1,-1,-1,-1,-1,+1,+1,+1,+1,+1,+1,+1,+1)
y = c(45,71,48,65,68,60,80,65,43,100,45,104,75,86,70,96)
datos<-data.frame(A,B,C,D,y);datos
##     A  B  C  D   y
## 1  -1 -1 -1 -1  45
## 2   1 -1 -1 -1  71
## 3  -1  1 -1 -1  48
## 4   1  1 -1 -1  65
## 5  -1 -1  1 -1  68
## 6   1 -1  1 -1  60
## 7  -1  1  1 -1  80
## 8   1  1  1 -1  65
## 9  -1 -1 -1  1  43
## 10  1 -1 -1  1 100
## 11 -1  1 -1  1  45
## 12  1  1 -1  1 104
## 13 -1 -1  1  1  75
## 14  1 -1  1  1  86
## 15 -1  1  1  1  70
## 16  1  1  1  1  96

2. Gráficos de interacciones.

par(mfrow=c(3,2), oma = c(0,0,2,0 ))
interaction.plot(A,B,y, main="Inter. AB", ylab="y")
interaction.plot(A,C,y, main="Inter. AC", ylab="y")
interaction.plot(A,D,y, main="Inter. AD", ylab="y")
interaction.plot(B,C,y, main="Inter. BC", ylab="y")
interaction.plot(B,D,y, main="Inter. BD", ylab="y")
interaction.plot(C,D,y, main="Inter. CD", ylab="y")
title( "Interacciones", outer = TRUE)

3. Análisis de varianza.

modelo=aov(y ~ A+B+C+D+A*B + A*C +A*B*C+ A*D + B*C + B*D +A*B*D+ C*D+B*C*D+A*B*C*D)
summary(modelo)
##             Df Sum Sq Mean Sq
## A            1 1870.6  1870.6
## B            1   39.1    39.1
## C            1  390.1   390.1
## D            1  855.6   855.6
## A:B          1    0.1     0.1
## A:C          1 1314.1  1314.1
## B:C          1   22.6    22.6
## A:D          1 1105.6  1105.6
## B:D          1    0.6     0.6
## C:D          1    5.1     5.1
## A:B:C        1   14.1    14.1
## A:B:D        1   68.1    68.1
## B:C:D        1   27.6    27.6
## A:C:D        1   10.6    10.6
## A:B:C:D      1    7.6     7.6

4. Suma de cuadrados.

sucua<-summary(modelo)[[1]]["Sum Sq"];sucua
##              Sum Sq
## A           1870.56
## B             39.06
## C            390.06
## D            855.56
## A:B            0.06
## A:C         1314.06
## B:C           22.56
## A:D         1105.56
## B:D            0.56
## C:D            5.06
## A:B:C         14.06
## A:B:D         68.06
## B:C:D         27.56
## A:C:D         10.56
## A:B:C:D        7.56
suma<-sum(sucua);suma
## [1] 5730.938
conporc<-sucua/suma*(100);conporc
##                   Sum Sq
## A           32.639729538
## B            0.681607503
## C            6.806259883
## D           14.928840177
## A:B          0.001090572
## A:C         22.929276405
## B:C          0.393696494
## A:D         19.291128197
## B:D          0.009815148
## C:D          0.088336332
## A:B:C        0.245378701
## A:B:D        1.187632913
## B:C:D        0.480942254
## A:C:D        0.184306669
## A:B:C:D      0.131959213

5. Modelo lineal.

lineal<-lm(y~ A+B+C+D+A*B + A*C +A*B*C+ A*D + B*C + B*D +A*B*D+ C*D+B*C*D+A*B*C*D)
summary(lineal)
## 
## Call:
## lm.default(formula = y ~ A + B + C + D + A * B + A * C + A * 
##     B * C + A * D + B * C + B * D + A * B * D + C * D + B * C * 
##     D + A * B * C * D)
## 
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  70.0625        NaN     NaN      NaN
## A            10.8125        NaN     NaN      NaN
## B             1.5625        NaN     NaN      NaN
## C             4.9375        NaN     NaN      NaN
## D             7.3125        NaN     NaN      NaN
## A:B           0.0625        NaN     NaN      NaN
## A:C          -9.0625        NaN     NaN      NaN
## B:C           1.1875        NaN     NaN      NaN
## A:D           8.3125        NaN     NaN      NaN
## B:D          -0.1875        NaN     NaN      NaN
## C:D          -0.5625        NaN     NaN      NaN
## A:B:C         0.9375        NaN     NaN      NaN
## A:B:D         2.0625        NaN     NaN      NaN
## B:C:D        -1.3125        NaN     NaN      NaN
## A:C:D        -0.8125        NaN     NaN      NaN
## A:B:C:D       0.6875        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 15 and 0 DF,  p-value: NA

6. Efectos de los factores.

efe<-2*coef(lineal)
efecto<-efe[-1];efecto
##       A       B       C       D     A:B     A:C     B:C     A:D     B:D     C:D 
##  21.625   3.125   9.875  14.625   0.125 -18.125   2.375  16.625  -0.375  -1.125 
##   A:B:C   A:B:D   B:C:D   A:C:D A:B:C:D 
##   1.875   4.125  -2.625  -1.625   1.375
efectos<-abs(efecto);efectos
##       A       B       C       D     A:B     A:C     B:C     A:D     B:D     C:D 
##  21.625   3.125   9.875  14.625   0.125  18.125   2.375  16.625   0.375   1.125 
##   A:B:C   A:B:D   B:C:D   A:C:D A:B:C:D 
##   1.875   4.125   2.625   1.625   1.375
da<-data.frame(efecto,sucua,conporc);da
##          efecto    Sum.Sq     Sum.Sq.1
## A        21.625 1870.5625 32.639729538
## B         3.125   39.0625  0.681607503
## C         9.875  390.0625  6.806259883
## D        14.625  855.5625 14.928840177
## A:B       0.125    0.0625  0.001090572
## A:C     -18.125 1314.0625 22.929276405
## B:C       2.375   22.5625  0.393696494
## A:D      16.625 1105.5625 19.291128197
## B:D      -0.375    0.5625  0.009815148
## C:D      -1.125    5.0625  0.088336332
## A:B:C     1.875   14.0625  0.245378701
## A:B:D     4.125   68.0625  1.187632913
## B:C:D    -2.625   27.5625  0.480942254
## A:C:D    -1.625   10.5625  0.184306669
## A:B:C:D   1.375    7.5625  0.131959213

7. Gráfico de Pareto.

pareto.chart(efectos, ylab = "Efectos", col=rainbow(length(efectos)))

##          
## Pareto chart analysis for efectos
##             Frequency   Cum.Freq.  Percentage Cum.Percent.
##   A        21.6250000  21.6250000  21.7063990   21.7063990
##   A:C      18.1250000  39.7500000  18.1932246   39.8996236
##   A:D      16.6250000  56.3750000  16.6875784   56.5872020
##   D        14.6250000  71.0000000  14.6800502   71.2672522
##   C         9.8750000  80.8750000   9.9121706   81.1794228
##   A:B:D     4.1250000  85.0000000   4.1405270   85.3199498
##   B         3.1250000  88.1250000   3.1367629   88.4567127
##   B:C:D     2.6250000  90.7500000   2.6348808   91.0915935
##   B:C       2.3750000  93.1250000   2.3839398   93.4755332
##   A:B:C     1.8750000  95.0000000   1.8820577   95.3575910
##   A:C:D     1.6250000  96.6250000   1.6311167   96.9887077
##   A:B:C:D   1.3750000  98.0000000   1.3801757   98.3688833
##   C:D       1.1250000  99.1250000   1.1292346   99.4981179
##   B:D       0.3750000  99.5000000   0.3764115   99.8745295
##   A:B       0.1250000  99.6250000   0.1254705  100.0000000

8. Gráfico de Daniel.

names(efectos)<-c("A","B","C","D","AB","AC","BC","AD","BD","CD","ABC","ABD","BCD","ACD","ABCD" )
DanielPlot(lineal,code=TRUE,half=FALSE)

DanielPlot(lineal, half = TRUE, main = "Half-Normal Plot")

9. Gráfico de Length.

LenthPlot(lineal,main="Length Plot")

##     alpha       PSE        ME       SME 
##  0.050000  2.625000  6.747777 13.698960

10. Modelo reducido.

Solo se consideran los efectos significativos: A, C, D, AC y AD.

modelo3<-aov(y ~ A+C+D+A*C+A*D)
summary(modelo3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1 1870.6  1870.6   95.86 1.93e-06 ***
## C            1  390.1   390.1   19.99   0.0012 ** 
## D            1  855.6   855.6   43.85 5.92e-05 ***
## A:C          1 1314.1  1314.1   67.34 9.41e-06 ***
## A:D          1 1105.6  1105.6   56.66 2.00e-05 ***
## Residuals   10  195.1    19.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11. Modelo lineal reducido.

lineal2<-lm(y ~ A+C+D+A*C+A*D)
summary(lineal2)
## 
## Call:
## lm.default(formula = y ~ A + C + D + A * C + A * D)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3750 -1.5000  0.0625  2.9062  5.7500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   70.062      1.104  63.444 2.30e-14 ***
## A             10.812      1.104   9.791 1.93e-06 ***
## C              4.938      1.104   4.471   0.0012 ** 
## D              7.313      1.104   6.622 5.92e-05 ***
## A:C           -9.063      1.104  -8.206 9.41e-06 ***
## A:D            8.312      1.104   7.527 2.00e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.417 on 10 degrees of freedom
## Multiple R-squared:  0.966,  Adjusted R-squared:  0.9489 
## F-statistic: 56.74 on 5 and 10 DF,  p-value: 5.14e-07

12. Ecuación del modelo lineal reducido.

coe<-coef(lineal2);coe
## (Intercept)           A           C           D         A:C         A:D 
##     70.0625     10.8125      4.9375      7.3125     -9.0625      8.3125

\[\begin{equation*} \hat{y}=70.0625+10.8125 \cdot x_{1}+4.9375 \cdot x_{3}+7.3125\cdot x_{4}-9.0625\cdot x_{1}x_{3}+8.3125\cdot x_{1}x_{4} \end{equation*}\]

13. Valores estimados y errores.

estima<-predict(modelo3)
error<-resid(modelo3)
da2<-data.frame(y,estima,error);da2
##      y  estima  error
## 1   45  46.250 -1.250
## 2   71  69.375  1.625
## 3   48  46.250  1.750
## 4   65  69.375 -4.375
## 5   68  74.250 -6.250
## 6   60  61.125 -1.125
## 7   80  74.250  5.750
## 8   65  61.125  3.875
## 9   43  44.250 -1.250
## 10 100 100.625 -0.625
## 11  45  44.250  0.750
## 12 104 100.625  3.375
## 13  75  72.250  2.750
## 14  86  92.375 -6.375
## 15  70  72.250 -2.250
## 16  96  92.375  3.625
qqnorm(error)
qqline(error)

14. Superficie de respuesta.

Haciendo \(x_{4}=1\).

\[\begin{equation*} \hat{y}=77.375 +19.125\cdot x_{1}+4.9375\cdot x_{3}-9.0625x_{4} \end{equation*}\]

x<-seq(-1,1,length=50)
y<-x
f<-function(x,y){
-77.375+19.125*x+4.9375*y-9.0625*x*y
}
z<-outer(x,y,f)
persp(x,y,z,theta=45,phi=45)

contour(x,y,z)

Haciendo \(x_{1}=1\).

\[\begin{equation*} \hat{y}= 80.125-4.125 \cdot x+ 15.625 \cdot y \end{equation*}\]

x<-seq(-1,1,length=50)
y<-x
f<-function(x,y){
80.125-4.125*x+ 15.625*y
}
z<-outer(x,y,f)
persp(x,y,z,theta=45,phi=45)

contour(x,y,z)

|—|

O.M.F.