2
$\begingroup$

I am trying to implement and analyse a full factorial experiment in R but I don't understand why the results presented in the book are different. Here are the problem details:

enter image description here

I tried to use the least square model to estimate the effects of different factors such as gap, power and flow rate but the effect sizes mentioned in the book are completely different :

enter image description here

My implementation of the problem in R and the results are as follows:

et_rate = c(550, 669, 633, 642, 1037, 749, 1075, 729,
            604, 650, 601, 635, 1052, 868, 1063, 860)

gap = factor(rep(1:2, times = 8))

flw_rate = factor(rep(1:2, each = 2, times = 4))

pwr = factor(rep(1:2, each = 4, times= 2))

df <- data.frame(gap, flw_rate, pwr, et_rate)


md3 <- lm(et_rate ~ .^3, data = df)
summary(md3)

And my results are:

Call:
lm(formula = et_rate ~ .^3, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-65.50 -11.12   0.00  11.12  65.50 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           577.00      33.56  17.193 1.33e-07 ***
gap2                   82.50      47.46   1.738  0.12036    
flw_rate2              40.00      47.46   0.843  0.42382    
pwr2                  467.50      47.46   9.850 9.50e-06 ***
gap2:flw_rate2        -61.00      67.12  -0.909  0.39000    
gap2:pwr2            -318.50      67.12  -4.745  0.00145 ** 
flw_rate2:pwr2        -15.50      67.12  -0.231  0.82317    
gap2:flw_rate2:pwr2    22.50      94.92   0.237  0.81859    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 47.46 on 8 degrees of freedom
Multiple R-squared:  0.9661,    Adjusted R-squared:  0.9364 
F-statistic: 32.56 on 7 and 8 DF,  p-value: 2.896e-05

Show in New WindowClear OutputExpand/Collapse Output

Call:
lm(formula = et_rate ~ gap * pwr, data = df)

Coefficients:
(Intercept)         gap2         pwr2    gap2:pwr2  
      597.0         52.0        459.7       -307.2  

I was expecting the coefficients of my model to be equal to the effect size estimate but they are completely different than mentioned in the solution in the book. Am I mistaken in my approach to get the effect sizes?

$\endgroup$
5
  • $\begingroup$ How does your book define "effect estimate"? $\endgroup$ – gung - Reinstate Monica Nov 18 '20 at 21:13
  • $\begingroup$ Note that your book is using effect coding for the factors, not level means coding, which R uses by default. You can easily reproduce the sums of squares, etc, w/ lm(..., contrasts=list(gap="contr.sum", flw_rate="contr.sum", pwr="contr.sum")); anova(md3). $\endgroup$ – gung - Reinstate Monica Nov 18 '20 at 21:16
  • $\begingroup$ The book doesn't mention the coding type or gives the details on that. That was so confusing. Regression modeling is new to me so I wasn't sure how to search explanation for that. Thanks for pointing me in the right direction :) $\endgroup$ – Zaffresky Nov 18 '20 at 22:58
  • $\begingroup$ There's a limit to what I can tell you w/o more information. The contrast type, etc, I can just recognize from what's there. $\endgroup$ – gung - Reinstate Monica Nov 18 '20 at 23:09
  • $\begingroup$ If i am not mistaken then the book uses orthogonal contrasts using treatment averages. In fact the book focuses designing and analyzing experiments using jmp but I don't have access to jmp. so I am trying to implement it in R $\endgroup$ – Zaffresky Nov 19 '20 at 10:01
2
$\begingroup$

The effect estimate is the difference in the result from when your factor changes from the high value to the low value.

In your problem statement above the effect estimate for factor A is the average etch rate when factor A=1 minus the average etch rate when A=-1.

When A= 1: Average( 669, 650, 642, 635, 749, 868, 729, 860)= 725.25

When A= -1: Average( 550, 604, 633, 601, 1037, 1052, 1075, 1063)= 826.875

Thus the effect estimate for factor A = 727.25 - 826.875 = -101.625.

Now repeat for the remaining factors and for the interaction factors.

In R, here is the script:

et_rate = c(550, 669, 633, 642, 1037, 749, 1075, 729,
            604, 650, 601, 635, 1052, 868, 1063, 860)

gap = (rep(c(-1, 1), times = 8))  
flw_rate = (rep(c(-1, 1), each = 2, times = 4)) 
pwr = (rep(c(-1, 1), each = 4, times= 2))

df <- data.frame(gap, flw_rate, pwr, et_rate)

df$AB <- df$gap*df$flw_rate
df$AC <- df$gap*df$pwr
df$BC <- df$flw_rate*df$pwr
df$ABC <- df$gap*df$flw_rate*df$pwr

variables <- names(df)[names(df)!="et_rate"]
sapply(variables, function(var) {
   mean(df$et_rate[df[[var]]==1] - df$et_rate[df[[var]]==-1])
})

#Output:
     gap flw_rate      pwr       AB       AC       BC      ABC 
-101.625    7.375  306.125  -24.875 -153.625   -2.125    5.625 

#perform the ANOVA
summary(aov(et_rate ~ gap*flw_rate*pwr, data=df))

                Df Sum Sq Mean Sq F value   Pr(>F)    
gap               1  41311   41311  18.339 0.002679 ** 
flw_rate          1    218     218   0.097 0.763911    
pwr               1 374850  374850 166.411 1.23e-06 ***
gap:flw_rate      1   2475    2475   1.099 0.325168    
gap:pwr           1  94403   94403  41.909 0.000193 ***
flw_rate:pwr      1     18      18   0.008 0.930849    
gap:flw_rate:pwr  1    127     127   0.056 0.818586    
Residuals         8  18020    2253                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This the reason why the levels are generally defined as -1 and 1!

By the way: And being an etch engineer the results do make sense, gas flow has little effect after a minimum threshold, smaller gap - more etch rate and power has the greatest affect on etch rate but worst selectivity.

$\endgroup$
2
  • $\begingroup$ So I need to calculate the effect sizes manually as shown in your code. So how would you interpret the coefficients in my original model. Are they just comparisons between different treatments? $\endgroup$ – Zaffresky Nov 19 '20 at 10:21
  • 1
    $\begingroup$ Since you defined your experimental factors and R factors, the coefficients are the expected increase from the base case when the factors are change from 1 to 2. In the above example, the model says if you change the gap from the low value to the high value, one could expect an increase of 82.5. (this value is incorrect since it is not considering the interactions). It get confusing when the interaction parameters are included. DOE has been around long before computers thus some of concepts can be hand calculated. $\endgroup$ – Dave2e Nov 19 '20 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.