Grass example

#Factorial version

The grass example as a verification of optimal design potential. The original experiment is a full factorial design with 84 runs, including two levels for the main control factors (type and tratment) and seven blocks for the carbon bioxide concentration (conc). The response was the uptake per area \(\frac{u mol}{m^2/s}\) and six specimens per type.

full <- lm(uptake~Type+Treatment+conc+Type:Treatment+Type:conc,CO2)
summary (full)

Call:
lm(formula = uptake ~ Type + Treatment + conc + Type:Treatment + 
    Type:conc, data = CO2)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6052  -3.0155   0.7975   3.9059  10.6174 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      25.293514   1.834972  13.784  < 2e-16 ***
TypeMississippi                  -4.726924   2.595043  -1.822   0.0724 .  
Treatmentchilled                 -3.580952   1.793560  -1.997   0.0494 *  
conc                              0.023080   0.003049   7.571 6.36e-11 ***
TypeMississippi:Treatmentchilled -6.557143   2.536477  -2.585   0.0116 *  
TypeMississippi:conc             -0.010699   0.004311  -2.482   0.0152 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.812 on 78 degrees of freedom
Multiple R-squared:  0.7286,    Adjusted R-squared:  0.7112 
F-statistic: 41.88 on 5 and 78 DF,  p-value: < 2.2e-16

anova (full)

Analysis of Variance Table

Response: uptake
               Df Sum Sq Mean Sq F value    Pr(>F)    
Type            1 3365.5  3365.5 99.6398 1.366e-15 ***
Treatment       1  988.1   988.1 29.2541 6.741e-07 ***
conc            1 2285.0  2285.0 67.6494 3.448e-12 ***
Type:Treatment  1  225.7   225.7  6.6829   0.01160 *  
Type:conc       1  208.0   208.0  6.1580   0.01523 *  
Residuals      78 2634.6    33.8                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library (sjPlot)
plot_model(full, type = "est", show.values = T) #Forest plot efectos

plot_model(full, type = "int") #Interacciones

[[1]]

[[2]]

plot_model(full, type = "diag") #Diagnóstico del modelo

[[1]]

[[2]]

`geom_smooth()` using formula 'y ~ x'

[[3]]

[[4]]

`geom_smooth()` using formula 'y ~ x'

Optimal version (automatic)

In a first optimal version the algorithm is allowed to select the number of runs automatically. It selects eleven runs only. The results show a similar structure but the values are not similar.

library (AlgDesign)
opt <- optFederov(uptake~Type+Treatment+conc+Type:Treatment+Type:conc,CO2)
opt

$D
[1] 17.79998

$A
[1] 7.045747

$Ge
[1] 0.764

$Dea
[1] 0.734

$design
Grouped Data: uptake ~ conc | Plant
   Plant        Type  Treatment conc uptake
1    Qn1      Quebec nonchilled   95   16.0
7    Qn1      Quebec nonchilled 1000   39.7
8    Qn2      Quebec nonchilled   95   13.6
21   Qn3      Quebec nonchilled 1000   45.5
22   Qc1      Quebec    chilled   95   14.2
28   Qc1      Quebec    chilled 1000   38.7
43   Mn1 Mississippi nonchilled   95   10.6
49   Mn1 Mississippi nonchilled 1000   35.5
50   Mn2 Mississippi nonchilled   95   12.0
64   Mc1 Mississippi    chilled   95   10.5
70   Mc1 Mississippi    chilled 1000   21.9

$rows
 [1]  1  7  8 21 22 28 43 49 50 64 70

min <- opt$design

fit <- lm(uptake~Type+Treatment+conc+Type:Treatment+Type:conc,min)
summary (fit)

Call:
lm(formula = uptake ~ Type + Treatment + conc + Type:Treatment + 
    Type:conc, data = min)

Residuals:
     1      7      8     21     22     28     43     49     50     64     70 
 0.650 -2.350 -1.750  3.450  1.100 -1.100 -2.529  3.657 -1.129  3.657 -3.657 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      12.547238   2.650796   4.733 0.005180 ** 
TypeMississippi                  -1.383149   3.756384  -0.368 0.727789    
Treatmentchilled                 -2.250000   3.266223  -0.689 0.521560    
conc                              0.029503   0.003403   8.670 0.000337 ***
TypeMississippi:Treatmentchilled -4.035714   4.781261  -0.844 0.437140    
TypeMississippi:conc             -0.008824   0.005144  -1.715 0.146955    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.772 on 5 degrees of freedom
Multiple R-squared:  0.9612,    Adjusted R-squared:  0.9224 
F-statistic: 24.76 on 5 and 5 DF,  p-value: 0.001546

anova (fit)

Analysis of Variance Table

Response: uptake
               Df  Sum Sq Mean Sq F value   Pr(>F)    
Type            1  264.61  264.61 18.6025 0.007621 ** 
Treatment       1   18.25   18.25  1.2832 0.308684    
conc            1 1422.25 1422.25 99.9873 0.000171 ***
Type:Treatment  1   14.36   14.36  1.0097 0.361093    
Type:conc       1   41.85   41.85  2.9422 0.146955    
Residuals       5   71.12   14.22                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library (sjPlot)
plot_model(fit, type = "est", show.values = T) #Forest plot efectos

plot_model(fit, type = "int") #Interacciones

[[1]]

[[2]]

plot_model(fit, type = "diag") #Diagnóstico del modelo

[[1]]

[[2]]

`geom_smooth()` using formula 'y ~ x'

[[3]]

[[4]]

`geom_smooth()` using formula 'y ~ x'

A second optimal design with fixed runs

This time the algorithm was set twenty three runs.

library (AlgDesign)
opt <- optFederov(~Type+Treatment+conc+Type:Treatment+Type:conc,CO2, nTrials = 23)
opt

$D
[1] 18.46528

$A
[1] 7.540847

$Ge
[1] 0.783

$Dea
[1] 0.757

$design
Grouped Data: uptake ~ conc | Plant
   Plant        Type  Treatment conc uptake
1    Qn1      Quebec nonchilled   95   16.0
7    Qn1      Quebec nonchilled 1000   39.7
8    Qn2      Quebec nonchilled   95   13.6
14   Qn2      Quebec nonchilled 1000   44.3
15   Qn3      Quebec nonchilled   95   16.2
21   Qn3      Quebec nonchilled 1000   45.5
22   Qc1      Quebec    chilled   95   14.2
28   Qc1      Quebec    chilled 1000   38.7
29   Qc2      Quebec    chilled   95    9.3
35   Qc2      Quebec    chilled 1000   42.4
36   Qc3      Quebec    chilled   95   15.1
43   Mn1 Mississippi nonchilled   95   10.6
49   Mn1 Mississippi nonchilled 1000   35.5
50   Mn2 Mississippi nonchilled   95   12.0
56   Mn2 Mississippi nonchilled 1000   31.5
57   Mn3 Mississippi nonchilled   95   11.3
63   Mn3 Mississippi nonchilled 1000   27.8
64   Mc1 Mississippi    chilled   95   10.5
70   Mc1 Mississippi    chilled 1000   21.9
71   Mc2 Mississippi    chilled   95    7.7
77   Mc2 Mississippi    chilled 1000   14.4
78   Mc3 Mississippi    chilled   95   10.6
84   Mc3 Mississippi    chilled 1000   19.9

$rows
 [1]  1  7  8 14 15 21 22 28 29 35 36 43 49 50 56 57 63 64 70 71 77 78 84

min <- opt$design

f23 <- lm(uptake~Type+Treatment+conc+Type:Treatment+Type:conc,min)
summary (f23)

Call:
lm(formula = uptake ~ Type + Treatment + conc + Type:Treatment + 
    Type:conc, data = min)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.1250 -2.0069  0.6852  2.0296  6.6917 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      12.396194   1.927770   6.430 6.22e-06 ***
TypeMississippi                   0.150629   2.695535   0.056 0.956088    
Treatmentchilled                 -2.496296   2.134088  -1.170 0.258251    
conc                              0.030722   0.002358  13.028 2.83e-10 ***
TypeMississippi:Treatmentchilled -4.787037   2.941637  -1.627 0.122056    
TypeMississippi:conc             -0.014461   0.003250  -4.449 0.000352 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.507 on 17 degrees of freedom
Multiple R-squared:  0.9427,    Adjusted R-squared:  0.9258 
F-statistic: 55.92 on 5 and 17 DF,  p-value: 5.899e-10

anova (f23)

Analysis of Variance Table

Response: uptake
               Df  Sum Sq Mean Sq  F value    Pr(>F)    
Type            1  465.89  465.89  37.8872 1.058e-05 ***
Treatment       1  229.32  229.32  18.6493 0.0004661 ***
conc            1 2475.40 2475.40 201.3061 7.461e-11 ***
Type:Treatment  1   23.93   23.93   1.9464 0.1809407    
Type:conc       1  243.38  243.38  19.7927 0.0003522 ***
Residuals      17  209.04   12.30                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library (sjPlot)
plot_model(f23, type = "est", show.values = T) #Forest plot efectos

plot_model(f23, type = "int") #Interacciones

[[1]]

[[2]]

plot_model(f23, type = "diag") #Diagnóstico del modelo

[[1]]

[[2]]

`geom_smooth()` using formula 'y ~ x'

[[3]]

[[4]]

`geom_smooth()` using formula 'y ~ x'