STA6246 - DOE

RCBD - ANOVA Examples

Author

Dr. Cohen

Vascular Graft Experiment

A medical manufacturer wants to investigate a=4 different levels of extrusion pressure on flicks (defects) by using randomized complete block design.

The experimenter considers the batches of resin (raw material) as blocks b=6.

N= 24 runs or observations.

Data description

library(tidyverse)
library(gtsummary)
# enter ER
dataV=read.csv(file = "../datasets/VascularData.csv",row.names = 1)
dataV$Press=factor(dataV$Press, labels=c("8500","8700","8900","9100"))
dataV$Batch=as.factor(dataV$Batch)


dataV %>%
  ggplot(aes(x=Press,y=Yield)) +
  geom_point()+
  theme_bw()+
  labs(x="Extrusion pressure (PSI)",
       y="Yiled (%)")

dataV %>%
  tbl_summary(by=Press)
Characteristic 8500, N = 61 8700, N = 61 8900, N = 61 9100, N = 61
Batch
    1 1 (17%) 1 (17%) 1 (17%) 1 (17%)
    2 1 (17%) 1 (17%) 1 (17%) 1 (17%)
    3 1 (17%) 1 (17%) 1 (17%) 1 (17%)
    4 1 (17%) 1 (17%) 1 (17%) 1 (17%)
    5 1 (17%) 1 (17%) 1 (17%) 1 (17%)
    6 1 (17%) 1 (17%) 1 (17%) 1 (17%)
Yield 92.1 (89.5, 96.9) 91.6 (89.8, 94.2) 88.8 (86.7, 90.5) 86.5 (83.3, 89.0)
1 n (%); Median (IQR)

  • The yield seems to decrease with the extrusion pressure.

  • The variability seems to be similar among the groups (extrusion pressure)

ANOVA in R

Fitting
library(sjPlot)

# ANOVA one factor and one blocking factor
results = aov(Yield ~ Press + Batch ,data = dataV)
summary(results)
            Df Sum Sq Mean Sq F value  Pr(>F)   
Press        3  178.2   59.39   8.107 0.00192 **
Batch        5  192.2   38.45   5.249 0.00553 **
Residuals   15  109.9    7.33                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Individual CIs - need to add the CI for the reference level
confint(results)

The ANOVA was rejected (p-value \(= 0.0019\)). Therefore, the data is compatible with the means difference in the yields population among the 4 extrusion pressures. Next, we will diagnose the model.

Model Adequacy
par(mfrow=c(2,2))
plot(results)

We could consider a transformation because the Scale-Location plots is showing a trend.

Multiple Comparison

Considering 5% level of significance.

#Tukey Test 
THSD=TukeyHSD(results,which = "Press")
plot(THSD,las=1)

# Fisher LSD Test
library(agricolae)
L=LSD.test(results,"Press",console=TRUE,group = F)

Study: results ~ "Press"

LSD t Test for Yield 

Mean Square Error:  7.32575 

Press,  means and individual ( 95 %) CI

        Yield      std r      LCL      UCL  Min  Max
8500 92.81667 4.577081 6 90.46148 95.17185 87.4 98.2
8700 91.68333 3.304189 6 89.32815 94.03852 87.0 95.8
8900 88.91667 2.966760 6 86.56148 91.27185 85.5 93.4
9100 85.76667 4.445072 6 83.41148 88.12185 78.9 90.7

Alpha: 0.05 ; DF Error: 15
Critical Value of t: 2.13145 

Comparison between treatments means

            difference pvalue signif.        LCL       UCL
8500 - 8700   1.133333 0.4795         -2.1974047  4.464071
8500 - 8900   3.900000 0.0247       *  0.5692620  7.230738
8500 - 9100   7.050000 0.0004     ***  3.7192620 10.380738
8700 - 8900   2.766667 0.0970       . -0.5640714  6.097405
8700 - 9100   5.916667 0.0018      **  2.5859286  9.247405
8900 - 9100   3.150000 0.0621       . -0.1807380  6.480738

Conclusion:

The Extrusion Pressure setting does affect the Yield (defects) in the Vascular Graft Experiment. The Lower the Pressure the higher the Yield.

The 9100 pressure was significantly different from 8500 and 8700 pressures. However, the 9100 pressure statistically resulted in no difference in the yield compared to 8900 pressure.

Random Effects RCBD

library(lme4)
anovaRE <- lmer(Yield ~ 1 + Press + (1 | Batch), data = dataV)
summary(anovaRE)
Linear mixed model fit by REML ['lmerMod']
Formula: Yield ~ 1 + Press + (1 | Batch)
   Data: dataV

REML criterion at convergence: 112

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.32253 -0.64269  0.01068  0.54202  1.62882 

Random effects:
 Groups   Name        Variance Std.Dev.
 Batch    (Intercept) 7.781    2.789   
 Residual             7.326    2.707   
Number of obs: 24, groups:  Batch, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)   92.817      1.587  58.494
Press8700     -1.133      1.563  -0.725
Press8900     -3.900      1.563  -2.496
Press9100     -7.050      1.563  -4.512

Correlation of Fixed Effects:
          (Intr) Pr8700 Pr8900
Press8700 -0.492              
Press8900 -0.492  0.500       
Press9100 -0.492  0.500  0.500
# Confidence Intervals Variance Components and Overall Mean
confint(anovaRE,oldNames=FALSE)
                          2.5 %     97.5 %
sd_(Intercept)|Batch   1.136322  5.5401925
sigma                  1.840242  3.5643700
(Intercept)           89.726101 95.9072331
Press8700             -4.085262  1.8185951
Press8900             -6.851928 -0.9480715
Press9100            -10.001928 -4.0980715
# Estimates the random effect
ranef(anovaRE)
$Batch
  (Intercept)
1 -1.69652554
2 -0.03710096
3  0.97474330
4  0.61047937
5 -3.61902963
6  3.76743346

with conditional variances for "Batch" 
anova(anovaRE)
Analysis of Variance Table
      npar Sum Sq Mean Sq F value
Press    3 178.17   59.39  8.1071
icc = 7.781/(7.781+7.326)
icc
[1] 0.5150592

There is a significant effect of the variability coming from the batch to batch.