Hypotheses Multisample Analysis of Variance
Sources of Variation Multiple Comparisons D&B Chapter 11
Motivation
Thus far weve examined measuring differences between one sample and another Examples:
Experimental
Treatment and Control
Observational
Alternate Conditions Alternate Populations
Motivation
Under those circumstances we conducted a t-test to examin the difference between two means. But what if we wish to examine differences across a variety of treatments or conditions?
Examples
Continuous or Discrete Data
Agricultural production vs. amount of rain Intellectual impairment vs. level of lead in the blood
Nominal Data
Agricultural production vs. variety Physiological response vs. drug type
Response relative to predictor (factor: level)
Continuous-Discrete Data
Stanford-Binet Intelligence Score 100 60 0 70 80 90 110
5
10
15
20
Blood Lead Concentration
Nominal Data
30 Foot size (cm) 22 24 26 28
Male Gender
Female
Nominal Data Script
plot(c(0,3),c(min(females),max(males)),type='n', xlab="Gender",ylab="Foot size (cm)",axes=F) axis(side=2) axis(side=1,at=c(1,2),labels=c("Male","Female")) box() points(jitter(rep(1,length(males))), jitter(males),pch=1) points(jitter(rep(2,length(females))), jitter(females),pch=1)
Types of Analysis
Linear Regression
Gaussian Response with ContinuousDiscrete Predictors
Analysis of Variance
Gaussian Response with Nominal Predictors
Example: Guayule (wy-oo-lee)
Rubber Plant
Parthenium argentatum
Produces a kind of latex that is hypoallergenic Federer (1955) discusses four crop treatments on eight varieties to see the effect on seedling germination Help available on guayule data in Splus, but not in R Data: guayule.txt http://en.wikipedia.org/wiki/Guayule
Guayule Data
guayule variety treatment reps plants flats 1 V1 T1 1 66 1.V1 2 V2 T1 1 77 1.V2 3 V3 T1 1 51 1.V3 4 V4 T1 1 52 1.V4 5 V5 T1 1 45 1.V5 6 V6 T1 1 59 1.V6 7 V7 T1 1 56 1.V7 8 V8 T1 1 49 1.V8 9 V1 T2 1 12 1.V1 10 V2 T2 1 26 1.V2
Germination vs. Variety
80 Percent germination 20 40 60
V1
V2
V3
V4 Variety
V5
V6
V7
V8
Variety Script
boxplot(split(guayule$plants, guayule$variety), style.bxp="att") mtext(side=2,"Percent germination",line=3) mtext(side=1,"Variety",line=3)
mean(guayule$plants[guayule$variety=="V1"]) 24.66667 median(guayule$plants[guayule$variety=="V1"]) 12.5
Germination vs. Treatment
80 Percent germination 20 40 60
T1
T2 Treatment
T3
T4
Is there a treatment effect?
Previously
H 0 : 1 = 2
Now
H 0 : 1 = 2 = 3 = 4
How would we approach such a problem?
Compute the grand mean and then the means for each treatment Examine the variation around each set of means and see which does the better job relative to the number of parameters (means) used Analysis of Variance
Male and Female Foot Size
30 my.foot 22 0 24 26 28 32
20
40
60
80
100
Grand Mean : Sums Squares Total
32 30
my.foot
X ij X
22 0
24
26
28
20
40
60
80
100
Group Means : Sums Squares Error
30 32
my.foot
X ij X i
22 0
24
26
28
20
40
60
80
100
Group - Grand: Sums Squares Treatment
32 30
my.foot
Xi X
22 0
24
26
28
20
40
60
80
100
Sums of Squares Total
1 X= N
X
i =1 j =1 ni k i =1 j =1
k
ni
ij
SST = (X ij X ) dfTotal = N 1
2
Sums of Squares Error
1 Xi = ni
X
j =1 k ni i =1 j =1
ni
ij
SSE = (X ij X i ) df Error = N k MSE = SSE / df Error
2
Sums of Squares Treatment
SSTr = ni (X i X )
k i =1 2
dfTreatment = k 1 MSTr = SSTr / dfTreatment Note : May see SSTr expressed as SSR
Analysis of Variance
SST SSTr SSE dfTotal df Treatment df Error MSTr MSE
MSTr = SSTr / dfTreatment
SST = SSTr + SSE dfTotal = dfTreatment + df Error
MSE = SSE / df Error
Sums of Squares
SSTr / = ni (X i X ) / ~ (k 1)
k 2 2 2 2 i =1
SSE / = (X ij X i ) / 2 ~ 2 ( N k )
k nj 2 2 i =1 j =1
Mean Square
2 (k 1) /(k 1) MSTr =2 = F (k 1, N k ) MSE ( N k ) /( N k )
When null hypothesis is true MSTr=MSE, so F=1. When null hypothesis is false MSTr >MSE, so F>1 This is a one-tailed test!
One-tailed F distribution, =.05
1.0 df(x, 30, 30) 0.0 0.0 0.2 0.4 0.6 0.8 1.2
0.5
1.0
1.5 x
2.0
2.5
3.0
F Df Distribution, = 10,20,30,40
1.2 df(x, 40, 40) 0.0 0.0 0.2 0.4 0.6 0.8 1.0
0.5
1.0
1.5 x
2.0
2.5
3.0
Analysis of Variance Table
SST SSTr SSE dfTotal df Treatment df Error MSTr MSE Fstat P value
Assumptions
All populations have normal distribution. The variances are all equal. Measurements are independent samples.
Guayule Analysis
Lec17.ANOVA1.ssc
Using the attach command
names(guayule) [1] "variety" [4] "plants" attach(guayule) mean(plants) 25.30208 "treatment" "reps" "flats"
Using the attach command
search() [1] "C:\\NTRES 3130" [2] "guayule" [3] "splus" [4] "stat" [5] "data" [6] "trellis" [7] "nlme3" [8] "menu" [9] "sgui" [10] "main" detach(2) # when you are done
Grand Summary Statistics
mean(plants) 25.3 var(plants) 386.6 length(plants) 96 ss.func = function(x) sum((x - mean(x))^2) ss.func(plants) 36736.2
Error Summary Statistics
tapply(plants, treatment, mean) T1 T2 T3 T4 55.8 13.9 20.0 11.4 tapply(plants, treatment, var) T1 T2 T3 T4 129.9 28.6 85.3 15.2 tapply(plants, treatment, length) T1 T2 T3 T4 24 24 24 24 tapply(plants, treatment, ss.func) T1 T2 T3 T4 2989.3 657.8 1962.9 351.8
Set up for ANOVA
SST = ss.func(plants) MST = SST/(length(plants)-1) SSE = sum(tapply(plants,treatment,ss.func)) MSE = SSE/(length(plants)-4) SSTr = SST-SSE MSTr = SSTr/(4-1)
Set up for ANOVA
SST 36736.2 SSE 5961.9 SSTr 30774.2 MST 386.6 MSE 64.8 MSTr 10258.0
ANOVA
MSTr/MSE 158.2 # Our ANOVA F statistic
1-pf(MSTr/MSE, 4-1,length(plants)-4) 0 # P-value
ANOVA
guayule.aov = aov(plants ~ treatment) summary(guayule.aov) Df Sum of Sq Mean Sq F Value Pr(F) treatment 3 30774.28 10258.09 158.29 0 Residuals 92 5961.96 64.80
Conclusion: Reject H0 that treatments are the same
Multiple Comparisons
Multiple Comparisons
Ok, weve determined that treatments account for a significant part of the variation in guayule germination But, to what degree do they differ from one another?
Percent germination 20 40 60 80
T1 T2 Treatment T3 T4
Tukey: All Pairwise Comparisons
H 0 : B = A XB XA q= SE s 2 SE = s = MSE , n = samples per group n critical value q ( , , k )
2
Multicomparison Script
SEtukey = sqrt(ErrorMS/24) my.means = tapply(plants,treatment,mean)
(my.means["T1"]-my.means["T2"])/SEtukey (my.means["T1"]-my.means["T3"])/SEtukey (my.means["T1"]-my.means["T4"])/SEtukey (my.means["T2"]-my.means["T3"])/SEtukey (my.means["T2"]-my.means["T4"])/SEtukey (my.means["T3"]-my.means["T4"])/SEtukey
Multicomparison Results
T1-T2 25.5089 T1-T3 21.78146 T1-T4 27.03031 T2-T3 -3.727444 T2-T4 1.521406 T3-T4 -3 -2 5.24885 qtukey(0.95, 4, 92)/sqrt(2) 2.616614
P(x) 0.0 0.1 0.2 0.3 0.4
-1
0 x
1
2
3
Splus: Tukey Pairwise Differences
multicomp(guayule.aov) 95 % simultaneous confidence intervals for specified linear combinations, by the Tukey method critical point: 2.6166 response variable: plants intervals excluding 0 are flagged by '****' Estimate Std.Error Lower Bound Upper Bound 41.90 2.32 35.80 48.0000 35.80 2.32 29.70 41.9000 44.40 2.32 38.30 50.5000 -6.13 2.32 -12.20 -0.0444 2.50 2.32 -3.58 8.5800 8.63 2.32 2.54 14.7000
T1-T2 T1-T3 T1-T4 T2-T3 T2-T4 T3-T4
**** **** **** **** ****
plot(multicomp(guayule.aov))
Multiple Comparison Figure Treatment 1 Definitely Different
T1-T2 T1-T3 T1-T4 T2-T3 T2-T4 T3-T4 -15
( ( ( ( ( ( ) ) ) )
) )
-5 0 5 10 15 20 25 30 35 40 45 simultaneous 95 % confidence limits, Tukey method response variable: plants
50
55
R: Tukey Pairwise Differences
TukeyHSD(guayule.aov) Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = plants ~ treatment) $treatment diff T2-T1 -41.9 T3-T1 -35.8 T4-T1 -44.4 T3-T2 6.1 T4-T2 -2.5 T4-T3 -8.6
lwr upr p adj -47.997 -35.8 0.00 -41.872 -29.7 0.00 -50.497 -38.3 0.00 0.044 12.2 0.05 -8.581 3.6 0.71 -14.706 -2.5 0.00
plot(TukeyHSD(guayule.aov))
95% family-wise confidence level
T2-T1 T4-T3 -50 T4-T2 T3-T2 T4-T1 T3-T1
-40
-30
-20
-10
0
10
Differences in mean levels of treatment