Introduction

This document contains the design and the analysis of an experiment that compares the mean effective life of four insulating fluids at an accelerated load of 35kV. In the first part, the number of samples of each fluid needed to detect a difference of 2 hours between the mean lives is determined for the cases of minimum, intermediate and maximum variability. In the second part, six observations of each fluid are analyzed with a one-way ANOVA, the adequacy of the model is checked, and Tukey’s test is used to identify which fluids differ. All of the data manipulation, analysis and plotting is done within R.

Question 1: Sample Size

The variance of the fluid life is estimated to be 3.5 hrs, the type 1 error probability is 0.05 and the desired power is 0.80. The smallest mean is 18 hrs and the largest mean is 20 hrs, so the difference between them is 2 hrs. The position of the two remaining means defines the three cases of variability. In the minimum variability case the two remaining means are in the middle (18, 19, 19, 20), in the intermediate case the four means are equally spaced (18, 18.67, 19.33, 20), and in the maximum variability case half of the means are at each extreme (18, 18, 20, 20).

means_min<-c(18,19,19,20)
means_int<-c(18,18+2/3,18+4/3,20)
means_max<-c(18,18,20,20)

The sample size is obtained with the function power.anova.test(), where between.var is the variance of the means and within.var is the variance of the fluid life.

Minimum Variability

p_min<-power.anova.test(groups=4,between.var=var(means_min),within.var=3.5,sig.level=0.05,power=0.80)
p_min
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 4
##               n = 20.08368
##     between.var = 0.6666667
##      within.var = 3.5
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

Intermediate Variability

p_int<-power.anova.test(groups=4,between.var=var(means_int),within.var=3.5,sig.level=0.05,power=0.80)
p_int
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 4
##               n = 18.17867
##     between.var = 0.7407407
##      within.var = 3.5
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

Maximum Variability

p_max<-power.anova.test(groups=4,between.var=var(means_max),within.var=3.5,sig.level=0.05,power=0.80)
p_max
## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 4
##               n = 10.56952
##     between.var = 1.333333
##      within.var = 3.5
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

Summary of Sample Sizes

cat("samples of each fluid (min, intermediate, max variability) \n",ceiling(c(p_min$n,p_int$n,p_max$n)),"\n")
## samples of each fluid (min, intermediate, max variability) 
##  21 19 11

Since the number of samples must be an integer, the values of n are rounded up. We need 21 samples of each fluid in the case of minimum variability, 19 samples in the case of intermediate variability and 11 samples in the case of maximum variability. The minimum variability case requires the most samples because the means are closest together, which makes the difference harder to detect. If we do not know how the means are arranged, the safest choice is the minimum variability case, since it guarantees the desired power for any arrangement.

Question 2: Analysis of the Test Data

Load the Data

The data is entered in a tidy format, with one column for the fluid type and one column for the life, so that each row is one observation.

life<-c(17.6,18.9,16.3,17.4,20.1,21.6,
        16.9,15.3,18.6,17.1,19.5,20.3,
        21.4,23.6,19.4,18.5,20.5,22.3,
        19.3,21.1,16.9,17.5,18.3,19.8)
fluid<-factor(rep(1:4,each=6))
df<-data.frame(fluid,life)
df
##    fluid life
## 1      1 17.6
## 2      1 18.9
## 3      1 16.3
## 4      1 17.4
## 5      1 20.1
## 6      1 21.6
## 7      2 16.9
## 8      2 15.3
## 9      2 18.6
## 10     2 17.1
## 11     2 19.5
## 12     2 20.3
## 13     3 21.4
## 14     3 23.6
## 15     3 19.4
## 16     3 18.5
## 17     3 20.5
## 18     3 22.3
## 19     4 19.3
## 20     4 21.1
## 21     4 16.9
## 22     4 17.5
## 23     4 18.3
## 24     4 19.8

Descriptive Statistics

cat("sample mean of the life of each fluid \n",tapply(df$life,df$fluid,mean),"\n")
## sample mean of the life of each fluid 
##  18.65 17.95 20.95 18.81667
cat("sample standard deviation of the life of each fluid \n",tapply(df$life,df$fluid,sd),"\n")
## sample standard deviation of the life of each fluid 
##  1.952178 1.854454 1.879096 1.554885

The sample means of fluids 1, 2, 3 and 4 are about 18.65, 17.95, 20.95 and 18.82 hrs. Fluid 3 has the largest mean life and fluid 2 has the smallest, and the difference between them is 3 hrs. The sample standard deviations are all between 1.55 and 1.95 hrs, so the four fluids show similar variability.

Side by Side Box Plots

boxplot(life~fluid,data=df,main="Life of the Insulating Fluids at 35kV",xlab="Fluid Type",ylab="Life (hr)",col=c("blue","red","green","orange"))

The box of fluid 3 is clearly higher than the other three boxes, while fluids 1, 2 and 4 are centered at about the same level. The boxes have similar heights, which is consistent with the similar standard deviations, and no observation is plotted as an outlier.

a. Hypothesis Test

The null hypothesis is that the mean life of the four fluids is the same, H0: mu1 = mu2 = mu3 = mu4, and the alternative hypothesis is that at least one mean is different. The test is done at an alpha = 0.10 level of significance.

model<-aov(life~fluid,data=df)
summary(model)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## fluid        3  30.17   10.05   3.047 0.0525 .
## Residuals   20  65.99    3.30                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The F statistic is 3.047 with 3 and 20 degrees of freedom, and the p-value is 0.0525. Since the p-value is smaller than alpha = 0.10, we reject the null hypothesis. There is evidence that the mean life is not the same for all four fluids.

b. Model Adequacy

The residuals of the model are used to check the assumptions of normality, constant variance and independence.

res<-residuals(model)
qqnorm(res,main="Normal Probability Plot of the Residuals")
qqline(res)

The residuals fall close to the straight line, with only small departures in the tails and no systematic curvature, so the normality assumption is reasonable.

plot(fitted(model),res,main="Residuals vs Fitted Values",xlab="Fitted Values (hr)",ylab="Residuals")
abline(h=0)

The residuals are spread about zero in a similar way for all four fitted values, and there is no funnel shape, so the constant variance assumption is reasonable.

plot(res,type="b",main="Residuals vs Observation Order",xlab="Observation Order",ylab="Residuals")
abline(h=0)

The residuals do not show any trend or pattern along the observation order, so there is no evidence against independence.

shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.95671, p-value = 0.376
bartlett.test(life~fluid,data=df)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  life by fluid
## Bartlett's K-squared = 0.26691, df = 3, p-value = 0.9661

The Shapiro-Wilk test has a p-value of about 0.38 and the Bartlett test has a p-value of about 0.97. Both are larger than 0.10, so we do not reject normality of the residuals or equality of the variances. The plots and the tests agree, so the model is adequate.

c. Tukey’s Test

Tukey’s test is done with a familywise error rate of alpha = 0.10, which corresponds to 90% confidence intervals.

tukey<-TukeyHSD(model,conf.level=0.90)
tukey
##   Tukey multiple comparisons of means
##     90% family-wise confidence level
## 
## Fit: aov(formula = life ~ fluid, data = df)
## 
## $fluid
##           diff        lwr       upr     p adj
## 2-1 -0.7000000 -3.2670196 1.8670196 0.9080815
## 3-1  2.3000000 -0.2670196 4.8670196 0.1593262
## 4-1  0.1666667 -2.4003529 2.7336862 0.9985213
## 3-2  3.0000000  0.4329804 5.5670196 0.0440578
## 4-2  0.8666667 -1.7003529 3.4336862 0.8413288
## 4-3 -2.1333333 -4.7003529 0.4336862 0.2090635
plot(tukey,las=1)

A pair of fluids differs significantly when its confidence interval does not contain zero. The only interval that does not contain zero is the one of fluids 3 and 2, with a difference of 3 hrs, an interval from about 0.43 to 5.57 hrs and an adjusted p-value of 0.044. So fluid 3 has a significantly longer mean life than fluid 2. The pairs 3-1 and 4-3 have intervals that almost exclude zero, but they are not significant at the 0.10 level, and all other pairs are clearly not different.

Complete Code

means_min<-c(18,19,19,20)
means_int<-c(18,18+2/3,18+4/3,20)
means_max<-c(18,18,20,20)
p_min<-power.anova.test(groups=4,between.var=var(means_min),within.var=3.5,sig.level=0.05,power=0.80)
p_min
p_int<-power.anova.test(groups=4,between.var=var(means_int),within.var=3.5,sig.level=0.05,power=0.80)
p_int
p_max<-power.anova.test(groups=4,between.var=var(means_max),within.var=3.5,sig.level=0.05,power=0.80)
p_max
cat("samples of each fluid (min, intermediate, max variability) \n",ceiling(c(p_min$n,p_int$n,p_max$n)),"\n")
life<-c(17.6,18.9,16.3,17.4,20.1,21.6,
        16.9,15.3,18.6,17.1,19.5,20.3,
        21.4,23.6,19.4,18.5,20.5,22.3,
        19.3,21.1,16.9,17.5,18.3,19.8)
fluid<-factor(rep(1:4,each=6))
df<-data.frame(fluid,life)
df
cat("sample mean of the life of each fluid \n",tapply(df$life,df$fluid,mean),"\n")
cat("sample standard deviation of the life of each fluid \n",tapply(df$life,df$fluid,sd),"\n")
boxplot(life~fluid,data=df,main="Life of the Insulating Fluids at 35kV",xlab="Fluid Type",ylab="Life (hr)",col=c("blue","red","green","orange"))
model<-aov(life~fluid,data=df)
summary(model)
res<-residuals(model)
qqnorm(res,main="Normal Probability Plot of the Residuals")
qqline(res)
plot(fitted(model),res,main="Residuals vs Fitted Values",xlab="Fitted Values (hr)",ylab="Residuals")
abline(h=0)
plot(res,type="b",main="Residuals vs Observation Order",xlab="Observation Order",ylab="Residuals")
abline(h=0)
shapiro.test(res)
bartlett.test(life~fluid,data=df)
tukey<-TukeyHSD(model,conf.level=0.90)
tukey
plot(tukey,las=1)