library('latex2exp')
library('lsmeans')
## Loading required package: emmeans
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
library('multcomp')
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## 
## Attaching package: 'multcomp'
## The following object is masked from 'package:emmeans':
## 
##     cld
library('nCDunnett')

assumptions:

misc topics to be covered later:

1 Example

An engineer is interested in the relationship between Radio-Frequency power setting and the etch rate. The engineer will test 5 wafer at each of the different RF power setting 160,180,200 and 220 Watt.

Ensure that your randomize the test to prevent the effect of unknow nuisance variables from contaminating the data. We are randomizing the experiment unit under which the treatment is applied. This ensure that the design is completely randomized designed

Input = ("
Order Power Response
13 '160' 575
14 '160' 542
8 '160' 530
5 '160' 539
4 '160' 570

18 '180' 565
9 '180' 593
6 '180' 590
16 '180' 579
17 '180' 610


7 '200' 600
19 '200' 651
10 '200' 610
20 '200' 637
1 '200' 629


2 '220' 725 
3 '220' 700 
15 '220' 715 
11 '220' 685 
12 '220' 710 
")

Data = read.table(textConnection(Input),header=TRUE)
Data = Data[order(Data$Order),]
print(Data)
##    Order Power Response
## 15     1   200      629
## 16     2   220      725
## 17     3   220      700
## 5      4   160      570
## 4      5   160      539
## 8      6   180      590
## 11     7   200      600
## 3      8   160      530
## 7      9   180      593
## 13    10   200      610
## 19    11   220      685
## 20    12   220      710
## 1     13   160      575
## 2     14   160      542
## 18    15   220      715
## 9     16   180      579
## 10    17   180      610
## 6     18   180      565
## 12    19   200      651
## 14    20   200      637
par(mfrow=c(1,2))
plot(Data$Power, Data$Response, main="Power(w) scatter plot", xlab = "Power(w)", ylab = "Etch Rate A/min")
boxplot(Response~Power,data=Data,main="Power(w) Boxplot", xlab = "Power(w)", ylab = "Etch Rate A/min")

2 The Analysis of Variance

Means Model: \[y_{ij} = \mu_i + \epsilon_{ij} ~~ i=1,2,...,a~ and~~j=1,2,...,n\]

effects Model: \[y_{ij} = \mu + \tau_{i} + \epsilon_{ij}~~i=1,2,...,a~and~~j=1,2,...,n\]

\(\epsilon_{ij}\) is a random error that encposses all forms of variations inherited from the experiments such as variabilities that arises from both controllered and uncontrolled factors, the regular noise of the process so on and so forth. The errors are assumed to be normally and independently distributed with mean 0 and variance \(\sigma^2\)

So observations \(y_{ij} \sim N(\mu + \tau_i, \sigma^2)\) where the observations are independent.

3 The Analysis of Fixed Effect Model

\[y_i. = \sum_{j=1}^ny_{ij} ~~ i=1,2,...,a~ and~~j=1,2,...,n ~~~~~~~ \bar{y}_{i.} = y_{i.}/n ~~ i=1,2,...,n \]

\(y_{..} = \sum_{i=1}^a\sum_{j=1}^by_{ij} ~~~~~ \bar{y}_{..}= y_{..}/N\)

\(H_0: \mu_1 = \mu_2 = \dots = \mu_a\)

\(H_1: \mu_i \neq \mu_j\) for at least one pair (i,j)

overall mean: \(\frac{\sum_{i=1}^a\mu_i}{a} = \mu\) which implies that \(\tau\) are deviation from the overall mean that are equal to zero under the \(H_0\)

so we can re-write the hypothesis as:

\(H_0: \tau_1 = \tau_2 = \dots = \tau_a\)

\(H_1: \tau_i \neq \tau_j\) for at least one pair (i,j)

Now we need to decomponse the total variability into its sum of squares components. In general \(SS_Total = SS_{Treatements}+SS_E\) and can be written as:

\(\sum_{i=1}^a\sum_{j=1}^n(y_{ij} - \bar{y}_{..})^2 = n\sum_{i=1}^a(\bar{y}_{i.} - \bar{y}_{..})^2 + \sum_{i=1}^a\sum_{j=1}^n(y_{ij}-\bar{y}_{i.})^2\)

Another appraoch to calculate Sum of Squares:

\(SS_T = \sum_{i=1}^a\sum_{j=1}^ny_{ij}^2 - y_{..}^2/N\)

\(SS_{Treatments} = \frac{1}{n}\sum_{i=1}^a\sum_{j=1}^ny_{i.}^2 - y_{..}^2/N\)

Anova for Plasma Etching Experiment

Data$Power <- as.factor(Data$Power)
power_response_model = aov(Response~Power,data=Data)
power_response_model_summary = summary(power_response_model)
power_response_model_summary
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Power        3  66871   22290    66.8 2.88e-09 ***
## Residuals   16   5339     334                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
x <- seq(0,30,by=0.001)
plot(1,type='n', xlim = c(0,30), ylim = c(0,0.9),xlab = TeX("$F_0$"), ylab = 'Probability Density',axes=F)
axis(side=1,at=seq(0,30,2))
axis(side=2,at=seq(0,1,0.1))
lines(x,df(x,3,16),col=1,lwd=2)
abline(v=power_response_model_summary[[1]]$`F value`,col='green')
abline(v=qf(0.95,3,16),col='red')
text(qf(0.95,3,16)+0.5, 0.8, TeX("$F_{0.05,3,16}"), col = "red",srt=90) 
abline(v=qf(0.99,3,16),col='red')
text(qf(0.99,3,16)+0.5, 0.8, TeX("$F_{0.01,3,16}"), col = "red",srt=90)

Our test statistic \(F_0\) is 66.8. So assuming \(H_0\) is true, the probability of obtaining a random sample whose F-value is at least as far as 66.8 is 2.88e-9. This is a very very rare event. this p-value of 2.88e-9 is also the risk of wrongly rejecting the null hypothesis which is also very small compared to the risk that we are willing to take on of p-value 0.05. So we reject the \(H_0\)

3.1 Esimation of Model Parameters and Confidence Interval

\(\hat{\mu} = \bar{y}_{..}\) \(\hat{\tau}_i = \bar{y}_{i.} - \bar{y}_{..}\) for i = 1,2,…,a

Assuming that the errors are normally distributed, then \(\bar{y}_{i.} \sim N(\mu,\sigma^2/n)\). If \(\sigma^2\) is known then use the normal districution otherwise use \(MS_E\) as an estimator of \(\sigma^2\) with t-distribution in the confidence interval below.\(1-\alpha\) applies only to one particular estimate.

\(P(-t_{\alpha/2} \leq \frac{\bar{y}_{i.}-\mu_i}{s/\sqrt{n}})=1- \alpha\)

\(\bar{y}_{i.} - t_{\alpha/2,N-a}*\sqrt{MS_E/n} \leq \mu_{i.} \leq \bar{y}_{i.} + t_{\alpha/2,N-a}*\sqrt{MS_E/n}\)

\(\bar{y}_{i.} - \bar{y}_{j.} - t_{\alpha/2,N-a}*\sqrt{2MS_E/n} \leq \mu_{i.} - \mu_{j.} \leq \bar{y}_{i.} - \bar{y}_{j.} + t_{\alpha/2,N-a}*\sqrt{2MS_E/n}\)

power_response_lm_model = lm(Response~Power,data=Data)
lsmeans(power_response_lm_model,"Power")
##  Power lsmean       SE df lower.CL upper.CL
##  160    551.2 8.169455 16 533.8815 568.5185
##  180    587.4 8.169455 16 570.0815 604.7185
##  200    625.4 8.169455 16 608.0815 642.7185
##  220    707.0 8.169455 16 689.6815 724.3185
## 
## Confidence level used: 0.95

Simultaneous Confidence Intervals

If you wish to calculate several confidence intervals r at 100(1-\(\alpha\)) percent confidendece, then the probability that the r intervals will simultaneously be correct is at least 1 - \(r\alpha\). \(r\alpha\) is called the experimentwise errr rate. Bonferroni method suggests to replace \(\alpha/2\) with \(\alpha/(2r)\)

Contrasts = list( "1v2" = c(1,-1,0,0), 
                  "1v3"=c(1,0,-1,0), 
                  "1v4"=c(1,0,0,-1),
                  "2v3"=c(0,1,-1,0),
                  "2v4"=c(0,1,0,-1),
                  "3v4"=c(0,0,1,-1))
power_response_lsmeans = lsmeans(power_response_lm_model,~Power)
contrast(power_response_lsmeans,Contrasts,adjust = "bonferroni")
##  contrast estimate       SE df t.ratio p.value
##  1v2         -36.2 11.55335 16  -3.133  0.0385
##  1v3         -74.2 11.55335 16  -6.422  0.0001
##  1v4        -155.8 11.55335 16 -13.485  <.0001
##  2v3         -38.0 11.55335 16  -3.289  0.0277
##  2v4        -119.6 11.55335 16 -10.352  <.0001
##  3v4         -81.6 11.55335 16  -7.063  <.0001
## 
## P value adjustment: bonferroni method for 6 tests

Another way of calling the above bonferroni’s method is using pairwise function

lsmeans(power_response_lm_model, pairwise ~ Power, adjust="bonferroni")
## $lsmeans
##  Power lsmean       SE df lower.CL upper.CL
##  160    551.2 8.169455 16 533.8815 568.5185
##  180    587.4 8.169455 16 570.0815 604.7185
##  200    625.4 8.169455 16 608.0815 642.7185
##  220    707.0 8.169455 16 689.6815 724.3185
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate       SE df t.ratio p.value
##  160 - 180    -36.2 11.55335 16  -3.133  0.0385
##  160 - 200    -74.2 11.55335 16  -6.422  0.0001
##  160 - 220   -155.8 11.55335 16 -13.485  <.0001
##  180 - 200    -38.0 11.55335 16  -3.289  0.0277
##  180 - 220   -119.6 11.55335 16 -10.352  <.0001
##  200 - 220    -81.6 11.55335 16  -7.063  <.0001
## 
## P value adjustment: bonferroni method for 6 tests
contrast(power_response_lsmeans,Contrasts,adjust = "tukey")
##  contrast estimate       SE df t.ratio p.value
##  1v2         -36.2 11.55335 16  -3.133  0.0379
##  1v3         -74.2 11.55335 16  -6.422  0.0001
##  1v4        -155.8 11.55335 16 -13.485  <.0001
##  2v3         -38.0 11.55335 16  -3.289  0.0274
##  2v4        -119.6 11.55335 16 -10.352  <.0001
##  3v4         -81.6 11.55335 16  -7.063  <.0001
## 
## P value adjustment: sidak method for 6 tests

Another way of calling the above tukey’s method is using pairwise function

lsmeans(power_response_lm_model, pairwise ~ Power, adjust="tukey")
## $lsmeans
##  Power lsmean       SE df lower.CL upper.CL
##  160    551.2 8.169455 16 533.8815 568.5185
##  180    587.4 8.169455 16 570.0815 604.7185
##  200    625.4 8.169455 16 608.0815 642.7185
##  220    707.0 8.169455 16 689.6815 724.3185
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate       SE df t.ratio p.value
##  160 - 180    -36.2 11.55335 16  -3.133  0.0294
##  160 - 200    -74.2 11.55335 16  -6.422  <.0001
##  160 - 220   -155.8 11.55335 16 -13.485  <.0001
##  180 - 200    -38.0 11.55335 16  -3.289  0.0216
##  180 - 220   -119.6 11.55335 16 -10.352  <.0001
##  200 - 220    -81.6 11.55335 16  -7.063  <.0001
## 
## P value adjustment: tukey method for comparing a family of 4 estimates

3.2 Model Adequacy Checking

#qqnorm implementation
N_residuals = length(power_response_lm_model$residuals)
residuals_ = power_response_lm_model$residuals
residuals_ = residuals_[order(residuals_)]


quantile_list = seq(1, N_residuals )
theoretical_quantile = (quantile_list-0.5)/N_residuals
theoretical_quantile = qnorm(theoretical_quantile)
theoretical_quantile = theoretical_quantile[order(theoretical_quantile)]
#or ppoints(N_residuals)




#plot(x=theoretical_quantile, y=residuals_)
#qqline(residuals_)

#or
plot(x=theoretical_quantile, y=residuals_)
#qqline implementation
slope = diff(quantile(residuals_,probs = c(0.25,0.75)))/diff(qnorm(c(0.25,0.75)))
x = qnorm(c(0.25,0.75))
y = quantile(residuals_,probs = c(0.25,0.75))
intercept = y[1] - slope*x[1]
abline(intercept,slope)

#qqnorm(power_response_lm_model$residuals)
#qqline(residuals(aov(power_response_lm_model)))

So the circle data seems to follow a straight line to indicate that it follows a normal distribution. But the tails of the distribution seems to be thinner. Ploting normal data to other distribution may help understand the patterns in the data. for example normal vs t-distribution

x = seq(-4,4,length.out = 2000)
y = dt(x,df=3)
plot(x,dnorm(x),type='l',lwd=3)
lines(x,dt(x,1), col='green', type = 'l' )
lines(x,dt(x,2), col='blue', type = 'l')
lines(x,dt(x,3), col='red', type = 'l')
lines(x,dt(x,4), col='orange', type = 'l')

x = qnorm(ppoints(30)) 
qqplot(x,x)
lines(x, qt(ppoints(30),df=1), col='green', type = 'l')
lines(x, qt(ppoints(30),df=2), col='blue', type = 'l')
lines(x, qt(ppoints(30),df=3), col='red', type = 'l')
lines(x, qt(ppoints(30),df=4), col='orange', type = 'l')

we need to check for outliers and one way to check that is to standardize the residuals. Shown below the maximum standardized residual is 1.52

\[d_{ij} = \frac{\epsilon_{ij}}{\sqrt{MS_E}}\]

standardized_residuals = residuals_/sd(residuals_)
max_error = max(abs(standardized_residuals))
#max_error
#pnorm(1.52)
sprintf("the max standardized residual is %s where %s percentage of the data fall",round(max_error,3),round(pnorm((max_error))*100,3))
## [1] "the max standardized residual is 1.527 where 93.664 percentage of the data fall"

Plot of Residuals in Time Sequence

In first plot below watch for positive and negative trends which violate the independence assumption on the errors. Proper randomization is key to help avoid dependencies between residuals.The residuals should be structureless nad unreleated to any other variables/treatments including predicted response. A simple check can be made by the second plot below.This plot was if the variance grows with the maganitude of the observations. This could indicate a nonconstant variance (funnel effect).

Note: if homogeneity of variance is violated F test is only slightly affected when data is balanced in a fixed effect model.always attempt to chose equal sample sizes.

To deal with nonconstant variance you may need to deal with transformation to bring the error distribution closer to normal.

if you know the distribution of the observation here are some suggestions: 1- if observations follow Poisson, then use square root transformation 2- if observations follow lognormal distribution, then use logarithmic transformation 3- if observations follow a binomial expressed as fractions, then use arcsin sqrt(Yij)

par(mfrow=c(2,2))
xlab1 = 'Run Order'
ylab = 'Residuals'
xlab2 = 'Fitted/Predicted'
plot(seq(1, length(Data$Order)), power_response_lm_model$residuals,xlab = xlab1,ylab = ylab)
plot(power_response_lm_model$fitted.values, power_response_lm_model$residuals,xlab = xlab2,ylab = ylab)

There are other statistical tools for testing for equality of variance such as Bartlett’s test if K-squared is bigger than K-squared > \(\chi_{\alpha,a-1}^2\) then we reject the null hypothesis. This test is sensitive to normality assumption

bartlett.test(Data$Response, Data$Power)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Data$Response and Data$Power
## Bartlett's K-squared = 0.43349, df = 3, p-value = 0.9332

Use modified Levene test in case of normality assumption doesn’t hold, use Levene test. To test the hypothesis of equal variances in all treatments, Levene test using the absolute deviation of the observation \(y_{ij}\)

\[ d_{ij} = | y_{ij} - median(y_{i.})|\] for i=1,2,…,a j=1,2…,\(n_i\)

let’s take an example where we already calculate the \(d_{ij}\) as shown below

Input = ("
Procedure Flow
1 0.18
1 0.40
1 0.71
1 0.18
1 1.23
1 0.4

2 1.7
2 0.33
2 0.47
2 0.25
2 0.25
2 1.94


3 1.495
3 0.565
3 1.945
3 1.715
3 2.015
3 0.565

4 1.56
4 3.77
4 4.64
4 1.61
4 1.24
4 1.23


")

Flow_Data = read.table(textConnection(Input),header=TRUE)
Flow_Data$Procedure = as.factor(Flow_Data$Procedure)
Flow_Data
##    Procedure  Flow
## 1          1 0.180
## 2          1 0.400
## 3          1 0.710
## 4          1 0.180
## 5          1 1.230
## 6          1 0.400
## 7          2 1.700
## 8          2 0.330
## 9          2 0.470
## 10         2 0.250
## 11         2 0.250
## 12         2 1.940
## 13         3 1.495
## 14         3 0.565
## 15         3 1.945
## 16         3 1.715
## 17         3 2.015
## 18         3 0.565
## 19         4 1.560
## 20         4 3.770
## 21         4 4.640
## 22         4 1.610
## 23         4 1.240
## 24         4 1.230

The anova method shows \(F_0\) = 4.55 for which p-value is p=0.0135, meaning that we reject the null hypothesis of equal variance. That indicates that a transformation is needed

summary(aov(Flow~Procedure, data = Flow_Data))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Procedure    3  11.57   3.857   4.554 0.0137 *
## Residuals   20  16.94   0.847                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

watch for the funnel effect

flow_model = aov(Flow~Procedure, data = Flow_Data)
plot(flow_model$fitted.values, flow_model$residuals,pch=16)

Input = ("
Procedure Flow
1 0.34
1 0.12
1 1.23
1 0.70
1 1.75
1 0.12

2 0.91
2 2.94
2 2.14
2 2.36
2 2.86
2 4.55


3 6.31
3 8.37
3 9.75
3 6.09
3 9.82
3 7.24

4 17.15
4 11.82
4 10.95
4 17.20
4 14.35
4 16.82


")


Flow_Data_Raw = read.table(textConnection(Input),header=TRUE)
Flow_Data_Raw$Procedure = as.factor(Flow_Data_Raw$Procedure)
Flow_Data_Raw
##    Procedure  Flow
## 1          1  0.34
## 2          1  0.12
## 3          1  1.23
## 4          1  0.70
## 5          1  1.75
## 6          1  0.12
## 7          2  0.91
## 8          2  2.94
## 9          2  2.14
## 10         2  2.36
## 11         2  2.86
## 12         2  4.55
## 13         3  6.31
## 14         3  8.37
## 15         3  9.75
## 16         3  6.09
## 17         3  9.82
## 18         3  7.24
## 19         4 17.15
## 20         4 11.82
## 21         4 10.95
## 22         4 17.20
## 23         4 14.35
## 24         4 16.82

empirical method for finding alpha value

we can empirical estimate apha if we think that $ _{y_i} ^$, where theta is constant of proportionality. so \(log \sigma_{y_i} = \log \theta + \alpha log\mu_i\), we can estimate \(\alpha\) from the data as shown below

#empirical method for finding alpha value
t1 = aggregate(x=Flow_Data_Raw$Flow, by=list(Flow_Data_Raw$Procedure), function(x) { c(M=mean(x),S=sd(x))})
plot( log(t1$x[1:4,2]), log( t1$x[1:4,1]))

Fitting the line using lm, we find that the slope is close to 0.5 (0.4465)

lm(log(t1$x[1:4,2])~log(t1$x[1:4,1]))
## 
## Call:
## lm(formula = log(t1$x[1:4, 2]) ~ log(t1$x[1:4, 1]))
## 
## Coefficients:
##       (Intercept)  log(t1$x[1:4, 1])  
##           -0.2781             0.4465
mean(diff(log(t1$x[1:4,2]))/diff(log(t1$x[1:4,1])))
## [1] 0.5340247

we can now power flow_rate by alpha = 0.5 and look for residuation vs fitted values .

Flow_Data_Raw$Flow_powerby0_5 = Flow_Data_Raw$Flow^(1/2)

m1 = aov(Flow_powerby0_5~Procedure, data=Flow_Data_Raw)
plot(m1$fitted.values, m1$residuals)

You can use Levene’s test instead of Bartlett’s test if normality is at risk by running aov on the \(d_ij\)

Flow_Data_Raw$Flow_powered_Median = median(Flow_Data_Raw$Flow_powerby0_5)
Flow_Data_Raw$Flow_powered_distance = abs( Flow_Data_Raw$Flow_powerby0_5 - median(Flow_Data_Raw$Flow_powerby0_5))
summary(aov(Flow_powered_distance~Procedure, data=Flow_Data_Raw))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Procedure    3  5.231  1.7437   12.97 6.25e-05 ***
## Residuals   20  2.688  0.1344                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.3 Practical Interpretation Results

3.3.1 A Regression Model

As one can see from the plot below there is a linear trend, so if factors are quantitative, perhaps someone is interested in predicting values between the two levels. If data for temperatures 10 and 30 were collectect, you could predict the response at 15 using a regressing model. \(regression~model: y = \beta_0 + \beta_1x + \epsilon\) If you look closely at the data, you might see some quadtratic effect, so quadratic fit maybe important here. check code below

Data$Power1 = as.numeric(as.character(droplevels(Data$Power)))
pairs(Data$Response~Data$Power1,data=Data)

lm(Response~ Power1,data=Data)
## 
## Call:
## lm(formula = Response ~ Power1, data = Data)
## 
## Coefficients:
## (Intercept)       Power1  
##     137.620        2.527
power_lm_model = lm(Response ~ Power1 + I(Power1^2),data=Data)
power_lm_model
## 
## Call:
## lm(formula = Response ~ Power1 + I(Power1^2), data = Data)
## 
## Coefficients:
## (Intercept)       Power1  I(Power1^2)  
##  1147.77000     -8.25550      0.02838
sprintf("let's predict the response value")
## [1] "let's predict the response value"
predict(power_lm_model,data.frame(Power1=c(200)))
##      1 
## 631.67
plot(Response~Power1, data=Data)
Power1_range = range(Data$Power1)
prd = data.frame(Power1 = seq(Power1_range[1],Power1_range[2],length.out = 100))
err = predict(power_lm_model,newdata=prd)
lines(seq(Power1_range[1],Power1_range[2],length.out = 100),err)

** sliding t-distribution for means comparison**

3.3.2 Contrasts

A contrast is a linear combination of parameters \[\Gamma = \sum_{i=1}^a c_i\mu_i ~~~where~ \sum_{i=1}^a c_i=0\]

You can also standardize the contrast by dividing each \(c_i\) by \(\sqrt{\frac{1}{n}\sum_{i=1}^ac_i^2}\)

in case the data is of unequal sample size, the contrast now requires that \(\sum_{i=1}^an_ic_i = 0\)

if I suspect that the average of the lowest two levels is equal to average of highest levels of power in the etching plasma example, I can construct a hypothesis as follow: \[H_0: \mu_1 + \mu_2 = \mu3 + \mu4\\H_1:\mu_1 + \mu_2 \ne \mu3 + \mu4 \] it can also be expressed as: \[ H_0: \sum_{i=1}^a c_i\mu_i = 0 \\ H_1: \sum_{i=1}^a c_i\mu_i \ne 0\]

so let C=\(\sum_{i=1}^a c_i\bar{y}_{i.}\), V(C) = \(\sigma^2/n\sum_{i=1}^a c_i^2\), if null hypothesis is true then \(\frac{\sum_{i=1}^a c_i\bar{y}_{i.}}{\sigma^2/n\sum_{i=1}^a c_i^2}\) has \(\sim N(0,1)\). If we don’t know \(\sigma^2\) we can replce it by its \(MS_E\) estimate and use the statistics \[ \frac{\sum_{i=1}^a c_i\bar{y}_{i.}}{\sqrt{\frac{MS_E}{n}\sum_{i=1}^a c_i^2}}\] The null hypothesis would be reject if \(|t_0|\) exceeds \(t_{\alpha/2,N-a}\)

for unequal sample size $ $

and the 100(1-alpha)% confidence interval \[\sum_{i=1}^a c_i\bar{y}_{i.}-t_{\alpha/2,N-a}*\sqrt{\frac{MS_E}{n}\sum_{i=1}^a c_i^2} \leq \sum_{i=1}^a c_i\mu_i \leq \sum_{i=1}^a c_i\bar{y}_{i.} + t_{\alpha/2,N-a}*\sqrt{\frac{MS_E}{n}\sum_{i=1}^a c_i^2} \]

A second approach is to use \(F_0 = t_0^2 = \frac{(\sum_{i=1}^a c_i\bar{y}_{i.})^2}{\frac{MS_E}{n}\sum_{i=1}^a c_i^2}\) and reject the null hypothese if \(F_0 > F_{\alpha,1,N-a}\)

Input = ("
Order Power Response C1v2 C12v34 C3v4
13 '160' 575 1 1 0
14 '160' 542 1 1 0
8 '160' 530 1 1 0
5 '160' 539 1 1 0
4 '160' 570 1 1 0

18 '180' 565 -1 1 0
9 '180' 593 -1 1 0
6 '180' 590 -1 1 0
16 '180' 579 -1 1 0
17 '180' 610 -1 1 0


7 '200' 600 0 -1 1
19 '200' 651 0 -1 1
10 '200' 610 0 -1 1
20 '200' 637 0 -1 1
1 '200' 629 0 -1 1


2 '220' 725  0 -1 -1
3 '220' 700  0 -1 -1
15 '220' 715  0 -1 -1
11 '220' 685  0 -1 -1
12 '220' 710  0 -1 -1
")

#added the orthogonal contrasts into the data frame.orthogonality completely partitions the treatment sum of squares with 1 df each.

Data2 = read.table(textConnection(Input),header=TRUE)
summary(aov(Response~C1v2+C12v34+C3v4,data=Data2))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## C1v2         1   3276    3276   9.818  0.00642 ** 
## C12v34       1  46948   46948 140.689 2.43e-09 ***
## C3v4         1  16646   16646  49.884 2.68e-06 ***
## Residuals   16   5339     334                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#another way of getting the estimates
#-------------------------------------

#Contrasts = list( "1v2" = c(1,-1,0,0), "12v34"=c(1,1,-1,-1), "3v4"=c(0,0,1,-1))
#power_response_lsmeans = lsmeans(power_response_lm_model,~Power)
#contrast(power_response_lsmeans,Contrasts,adjust = "bonferroni")
Input = ("
Order Power Response C1v2 C12v34 C3v4 C1v4
13 '160' 575 1 1 0 1
14 '160' 542 1 1 0 1
8 '160' 530 1 1 0 1
5 '160' 539 1 1 0 1
4 '160' 570 1 1 0 1

18 '180' 565 -1 1 0 0
9 '180' 593 -1 1 0 0
6 '180' 590 -1 1 0 0
16 '180' 579 -1 1 0 0
17 '180' 610 -1 1 0 0


7 '200' 600 0 -1 1 0
19 '200' 651 0 -1 1 0
10 '200' 610 0 -1 1 0
20 '200' 637 0 -1 1 0
1 '200' 629 0 -1 1 0


2 '220' 725  0 -1 -1 -1
3 '220' 700  0 -1 -1 -1
15 '220' 715  0 -1 -1 -1
11 '220' 685  0 -1 -1 -1
12 '220' 710  0 -1 -1 -1
")
Data3 = read.table(textConnection(Input),header=TRUE)
Data3$Power = as.factor(Data3$Power)

fit = aov(Response~Power, data=Data3)
#contrast
contrast_list = list(c(1,1,-1,-1))
fit.scheffe = glht(fit, linfct = mcp(Power=c(1,1,-1,-1)))
summary(fit.scheffe)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = Response ~ Power, data = Data3)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  -193.80      16.34  -11.86 2.43e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
pf(70.32,3,16,lower.tail = F)
## [1] 1.971775e-09

for multipel pairwise comparison use Tukey

fit = aov(Response~Power, data=Data3)
fit.tukey = glht(fit,linfct=mcp(Power=c(1,-1,1,1)))
summary(fit.tukey)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = Response ~ Power, data = Data3)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0   193.80      16.34   11.86 2.43e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Tukey??

k = matrix(c(1,1,-1,-1),nrow = 1)
k = rbind(k,c(1,-1,0,0))
fit.tukey = glht(fit,linfct=mcp(Power= k))
summary(fit.tukey)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = Response ~ Power, data = Data3)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  -193.80      16.34 -11.861 4.86e-09 ***
## 2 == 0   -36.20      11.55  -3.133   0.0126 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Dunnett’s test

One might be interested in comparing a-1 treatement levels to a control level. Dunnett procedure address this hypothesis. it is recommended to have more observations for control treatments \(n_a\) than for other treatments. if other treatments have equal numbers of observations, then it is recommended to chose sample size such as $ = $

\[H_0: \mu_i = \mu_a ~~~ for~i = 1,2,\dots,a-1\\ H_1: \mu_i \ne \mu_a\]

The null hypothesis is rejected if

\[|\bar{y}_{i.} - \bar{y}_{a.}| \ge d_\alpha(a-1,f)\sqrt{MS_E(1/n_i + 1/n_a)}\] the \(d_\alpha\) can be caculated in using nCDunnett library as shown below

library(nCDunnett)
qNCDun(1-0.05, nu=16, rho= (rep(0.5,times=3)), delta=rep(0,times=3),two.sided = T )
## [1] 2.592548

The d statistics is plotted

dun = dNCDun(seq(0,10,length.out = 100),16, rho=(rep(0.5,times=3)),delta=rep(0,times=3))
plot(seq(0,10,length.out = 100),dun,type='l')

so manually we can caculate the critical difference as shown below. you can see all the \(|\bar{y}_{i.} - \bar{y}_{a.}|\) are greater than the critical difference value.

qNCDun(1-0.05, nu=16, rho= (rep(0.5,times=3)), delta=rep(0,times=3),two.sided = T )*sqrt( (2*333.7)/5)
## [1] 29.95263

you can alos using glht method, but you need to order the levels accordingly since this method assume that the control group is the first group in the data frame.

Data3$Power2 = Data3$Power
Data3$Power2 = factor(Data3$Power2,levels = c("220","160","180","200"))

fit = aov(Response~Power2, data=Data3)
glht(fit,linfct=mcp(Power2="Dunnet"))
## 
##   General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Linear Hypotheses:
##                Estimate
## 160 - 220 == 0   -155.8
## 180 - 220 == 0   -119.6
## 200 - 220 == 0    -81.6

3.4 Sample Size

3.4.1 Operating Characteristic Curves

There are 4 ways to determine sample size:

  1. Use exiting OC chart by caculating \(\Phi\) by using a sample size ‘n’ and compute the corresponding df1=a and df2=a(n-1), then use this information to read the probability of type II error (\(\beta\)). or you can use the power.anova.test in R. from the method’s output, clearly a sample size of 4 or 5 will be adequate.
groupmeans <- c(575,600,650,675)
power.anova.test(groups=4, between.var=var(groupmeans),within.var=25^2,n=seq(2,20),sig.level=0.01)
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 4
##               n = 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
##     between.var = 2083.333
##      within.var = 625
##       sig.level = 0.01
##           power = 0.2361515, 0.7534019, 0.9621239, 0.9963525, 0.9997422, 0.9999854, 0.9999993, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000
## 
## NOTE: n is number in each group
  1. Using the difference between any two treatments is as large as D. The same as above but using delta instead

\[\Phi^2 = \frac{nD^2}{2a\sigma^2}\]

we can also using power.t.test in R. to achieve a power of 0.999, we would need to collect at least 5 to 6 samples

power.t.test(delta=75,sd=25,sig.level = 0.05,type='one.sample',alternative = "two.sided", power=0.998)
## 
##      One-sample t test power calculation 
## 
##               n = 5.051791
##           delta = 75
##              sd = 25
##       sig.level = 0.05
##           power = 0.998
##     alternative = two.sided
  1. Standard deviation Increase

If the treatments means are different, then standard deviation = \(\sqrt{\sigma^2 + (\sum_{i=1}^a\tau_i^2/a)}\) so \(\frac{\sqrt{\sigma^2 + (\sum_{i=1}^a\tau_i^2/a)}}{\sigma} = 1 + 0.01P\) (P = percent)

So the way we comoputer Phi is \[\Phi = \frac{\sqrt{\sum_{i=1}^a\tau_i^2/a}}{\sigma/\sqrt{n}} = \sqrt{(1+0.01P)^2 -1} (\sqrt{n})\]

  1. Using confidence Interval Estimation Method

Here we are using confidences interval. You specify how wide the C.I to be at 95% on the difference in treatment means. for example if 95% CI on difference in mean etch rate for any two settings to be +/- 30 with an estimate of sigma=25, then using

\[\pm t_{\alpha/2, N-a} \sqrt{ \frac{2MS_E}{n}} \]

n = 5
a = 4
sd = 25
df = a*(n-1) ##a(n-1) = 4(6-1) = 20
qt(0.975,df)*sqrt(2*sd^2/n)
## [1] 33.51865
# this doesn't meet the of 95% CI on delta of 

n = 6
a = 4
sd = 25
df = a*(n-1) ##a(n-1) = 4(6-1) = 20
qt(0.975,df)*sqrt(2*sd^2/n)
## [1] 30.10829
# this doesn't meet the of 95% CI on delta of 30


n = 7
a = 4
sd = 25
df = a*(n-1) ##a(n-1) = 4(6-1) = 20
qt(0.975,df)*sqrt(2*sd^2/n)
## [1] 27.58
# this doesn't meet the of 95% CI on delta of