Hypothesis Two-Sample Equal and Unequal Variances (D&B Chp. 10)
Difference Between Two Means Data: trout1.txt, trout2.txt & lead.washington.txt EqualScript: Lec12.TwoSampEqV.ssc UnEqualScript: Lec12.TwoSampUnEqV.ssc
Brown Trout (Salmo trutta)
Trout Length at Stocking
1977 Population : trout1
n = 99 mean = 136.5 mm var = 218.9 mm2
1978 Population : trout2
n = 100 mean = 207.8 mm var = 273.8 mm2
Is there a difference?
Concern that size at stocking will affect survivorship. First we must see if, in fact, they are different in size.
15
0
5
10
Trout 1
12
100
150
200
250
Trout 2
02468 100
150 Length (mm)
200
250
Histogram Script
my.breaks = seq(95,250,by=5) par(mfrow=c(2,1)) par(mar=c(0,4,4,2)) hist(trout1,breaks=my.breaks,cex=1.3) box() par(mar=c(5,4,0,2)) hist(trout2,breaks=my.breaks, xlab="Length (mm)",cex=1.3) box() par(mar=c(5,4,4,2)+.1)
Length (mm) 100 150 200 250
1977 n = 99 1978 n = 100
Boxplot Script
boxplot(trout1,trout2, names=c("1977","1978"), ylab="Length (mm)") text(45,mean(trout1), paste("n =",length(trout1))) text(90,mean(trout2), paste("n =",length(trout2)))
Other Displays of Information
plot(c(0,3),c(0,6), type='n',xlab="Examples",ylab="Length") boxes(x=c(1),width=c(0.25), stats=cbind(c(1,2,3,4,5))) points(2,3,pch=15,col=1) lines(c(2,2),c(2,4),lwd=3)
Length
0 0.0
1
2
3
4
5
6
0.5
1.0
1.5 Examples
2.0
2.5
3.0
Other Statistics
range(trout1) 99.06 175.26 range(trout2) 160.02 246.38 stdev(trout1) 14.79436 stdev(trout2) 16.54675
Length (mm) 100 150 200 250
1977 Year 1978
CI Boxes Script
plot(c(0,3),c(100,250), axes=F,type='n',ylab="Length (mm)",xlab="Year") axis(side=2) axis(side=1,at=c(1,2),labels=c("1977","1978")) box() boxes(c(1,2),width=0.25, cbind( c(range(trout1)[1], mean(trout1)-2*stdev(trout1)/sqrt(length(trout1)), mean(trout1), mean(trout1)+2*stdev(trout1)/sqrt(length(trout1)), range(trout1)[2]), c(range(trout2)[1], mean(trout2)-2*stdev(trout2)/sqrt(length(trout2)), mean(trout2), mean(trout2)+2*stdev(trout2)/sqrt(length(trout2)), range(trout2)[2])))
Length (mm) 100 150 200 250
1977 1978
Boxplot with Notch
boxplot(trout1,trout2,notch=T, names=c("1977","1978"),ylab="Length (mm)")
But is there a difference?
H 0 : 1 = 2 H A : 1 2
Test Statistic
x1 x2 ( 1 2 ) z= SE x1 x2
Difference?
x1 x2
V ( x1 x2 ) = ?
The difference represents a statistic, that under the null hypothesis of no difference should have a mean equal to zero. Variance of sum is sum of the variances.
Assumptions to Check
If the samples came from populations having normal distributions If the two populations have equal variances
Then
x1 x2 t= s x1 x2 where s x1 x2 = s
2 p
n1
+
s
2 p
n2
SS1 + SS 2 s= 1 + 2
2 p
Are they normal populations?
Use order of observations to determine normal quantiles: 50th quantile z = 0.0 2.5th quantile z = -2.0 97.5th quantile z = 2.0
trout1
100
120
140
160
qqnorm(trout1) qqline(trout1)
-2
-1
0 Quantiles of Standard Normal
1
2
240
qqnorm(trout2) qqline(trout2)
trout2
160
180
200
220
-2
-1
0 Quantiles of Standard Normal
1
2
Guidance:
The rejection region will be approximately the correct region even when the population distribution is non-normal provided the sample size is large; in most cases, n 30 is sufficient.
Are the variances the same?
Ratio of variances? We use the F-statistic var(trout2)/var(trout1) 2 * (1 - pf(var(trout2)/var(trout1), 100-1, 99-1)) 0.2685195 Reject if < 0.05 More on F statistics and F distributions later
So pool estimate of pop variance
SS1 + SS 2 s= 1 + 2
2 p
(sum((trout1-mean(trout1))^2)+sum((trout2-mean(trout2))^2))/ (length(trout1)-1 + length(trout2)-1)
246.4733 Compare with var(trout1) 218.873 var(trout2) 273.7948
And calculate standard error
s x1 x2 = s
2 p
n1
+
s
2 p
n2
sp2 = (sum((trout1-mean(trout1))^2)+sum((trout2-mean(trout2))^2))/ (length(trout1)-1 + length(trout2)-1) SEdiff = sqrt((sp2/length(trout1))+(sp2/length(trout2))) SEdiff 2.22584
And compute t-statistic
x1 x2 t= s x1 x2 t0.05( 2 ), = t0.05( 2 ),197 = 1.97
t.stat = (mean(trout1)-mean(trout2))/SEdiff t.stat -32.02334 degrees.freedom = (length(trout1)-1 + length(trout2)-1) degrees.freedom 197 pt(t.stat,degrees.freedom)
0
qt(0.025, 197) -1.972079
Conclusion
The populations are measurably (statistically significantly) different Next: Determine differences in mortality (but we will leave that for another time).
Lead poisoning
It is estimated that lead poisoning resulting from an unnatural craving for substances such as paint may affect as many as a quarter of a million children each year, causing them to suffer from severe, irreversible retardation. Explanations for why children voluntarily consume lead range from Improper parental supervision to childs need to mouth objects.
Lead poisoning (continued)
Some researchers, however, have investigated whether the habit of eating such substances has some nutritional explanation. One such study involved a comparison of a regular diet and a calcium-deficient diet on the ingestion of a lead-acetate solution in rats.
Lead poisoning (continued)
Each rat in a group of 20 rats was randomly assigned to either an experimental or a control group. Those in the control group received a normal diet, while the experimental group received a calcium-deficient diet. Each of the rats occupied a separate cage and was monitored to observed the quantity of a .15% lead-acetate solution consumed during the study period.
Lead poisoning (data)
Rat.Control = c(5.4,6.2,3.1,3.8,6.5,5.8,6.4,4.5,4.9,4.0) Rat.Treatment = c(8.8,9.5,10.6,9.6,7.5,6.9,7.4,6.5,10.5,8.3)
Two-Sample Hypothesis Unequal Variances
One- and Two-sided Tests Testing for Unequal Variances More Lead Poisoning
Comparing Two Populations
t-distributed if
Normal Equal variances
What if unequal variances? What if not normal? An introduction to the F-distribution
Washington State Childhood Blood Lead Levels
Washington State Dept. of Health Should children be screened for lead poisoning? Concern about higher levels in Eastern Washington
Author: Ossiander, E. WSDOH
lead.washington
lead.washington[1:10,] county pcntlead lon 1 Adams 4.2 -118.5447 2 Asotin NA -117.1584 3 Benton 1.0 -119.5364 4 Chelan 9.1 -120.6175 5 Clallam NA -123.9556 6 Clark 1.7 -122.4763 7 Columbia NA -117.9213 8 Cowlitz NA -122.6754 9 Douglas NA -119.6425 10 Ferry NA -118.5298 lat west 46.99545 0 46.18353 0 46.22546 0 47.84803 0 48.02465 1 45.76676 1 46.34998 0 46.20385 1 47.74384 0 48.55577 0
County Adams Benton Chelan
Percent Above 10ug/dL 4.2 1.0 9.1 1.7 0.8 1.6 1.7 4.0 5.0 3.8 3.7 4.7 2.1 1.9 2.8 3.4 0.0 2.6 10.8 2.8 5.3
West 1.6 (1) or East (0) 0 0 0 1 0 0 1 1 1 1 0 1 1 1 1 0 0 1 0 1 0 0
Percent of children with elevated blood lead levels (that is levels above 12ug/dL) by county 19931998.
Clark Franklin Grant Island King Kitsap Lewis Okanogan Pierce Skagit Skamania Snohomish Spokane Stevens Thurston Walla Walla Whatcom Whitman Yakima
Percent with Elevated Blood Lead
0.0 - 1.3 1.4 - 4.1 4.2 - 6.8 6.9 - 9.6 9.7 -11.0
library(maps)
# Data in lead.washington.txt
map("county","washington",fill=T, color=round(2+4*lead.washington$pcntlead/ max(lead.washington$pcntlead,na.rm=T))) id = seq(length(lead.washington$lon) )[lead.washington$west==1] points(lead.washington$lon[id], lead.washington$lat[id]) map("county","washington",lwd=4,add=T) #text(lead.washington$lon, # lead.washington$lat, # as.character(lead.washington$county)) key( corner=c(0,1), text = list( c("0.0 - 1.3","1.4 - 4.1","4.2 - 6.8", "6.9 - 9.6","9.7 -11.0"), col = 1, adj = 0.5), line = list(lty = 1, col = 2:6, lwd = 20) )
Percent Blood Lead Levels > 12 ug/Dl 0 2 4 6 8 10
Percent with Elevated Blood Lead
West East
Script for Boxplot
lead.west = lead.washington$pcntlead[lead.washington$west==1] lead.west = lead.west[!is.na(lead.west)] lead.east = lead.washington$pcntlead[lead.washington$west==0] lead.east = lead.east[!is.na(lead.east)] boxplot(lead.west,lead.east,names=c("West","East"))
Wish to examine if mean percent with elevated blood lead is higher in eastern region of Washington
Suggests a one-sided two-sample hypothesis test mean(lead.west) 3.07 mean(lead.east) 3.79 To use t-test Normal Equal variances
Two-tailed vs. One-tailed Tests
If the choice about the direction of the test is clear a priori then choosing a one-tailed test may allow you greater power.
Are population means equal? Two-tailed Test:
H 0 : 1 = 2 H A : 1 2 Level of significance, x1 x2 Test statistic, t = SE ( x1 x2 )
-3 -2 -1 0 x 1 2 3 0.4 P(x) 0.1 0.2 0.3
Reject H 0 , if t < -t1 / 2,df or t > t1 / 2,df
0.0
Is the mean of population 1 smaller than the mean of 2?
H 0 : 1 = 2 H A : 1 < 2 Level of significance, x1 x2 Test statistic, t = SE ( x1 x2 )
-3 -2 -1 0 x 1 2 3 0.4 P(x) 0.1 0.2 0.3
Reject H 0 , if t < t1 ,df t < t ,df
0.0
Is the mean of population 1 greater than the mean of 2?
H 0 : 1 = 2 H A : 1 > 2 Level of significance, x1 x2 Test statistic, t = SE ( x1 x2 )
-3 -2 -1 0 x 1 2 3 0.4 dnorm(x) 0.0 0.1 0.2 0.3
Reject H 0 , if t > t1 ,df
Now check on assumptions
Independent samples?
(Follows from sampling design.)
Normality? Equal variances?
0.0
1.0
2.0
3.00
Number Observed 1 2 3 4
0 2 4 6 Length (cm) 8 10
West
East
West
5 10
East
4
3
2
-1 0 1 Quantiles of Standard Normal
0 -1 0 1 Quantiles of Standard Normal
2
4
6
8
Script for Quantile Plots
par(mfrow=c(1,2)) qqnorm(lead.west) qqline(lead.west) title("West") qqnorm(lead.east) qqline(lead.east) title("East") par(mfrow=c(1,1))
Try: qqnorm(rpois(100,1)) qqnorm(rpois(100,10)) qqnorm(rpois(100,100))
Quantile Plots
Quantile plots suggest normality might be a problem Why? The sample sizes are not large, so normality of the observations should be considered. We will assume it is ok for the moment, and proceed to checking to see if the variances are equal.
Examine Variances
F-statistic =
2 2 / 2 F (df 2, df 1) = 2 1 / 1
v1 = var(lead.west) v2 = var(lead.east) n1 = length(lead.west) n2 = length(lead.east)
# 1.600182 # 10.96992 # 11 # 12
2*(1-pf(v2/v1,n2-1,n1-1)) # 0.005015549 Accept or reject variances equal?
F Distribution 11 df,10 df : 95% CI
0.6 df(x, 11, 10) 0.4
0.0
0.2
Location of our F-statistic
0
2
4 x
6
8
Computing the Standard Error and Effective Degrees of Freedom
x1 x2 t' = s ' x1 x2 where
2 s12 s2 + n1 n2
s ' x1 x2 =
s s + n n 2 1 '= 22 22 s1 s2 n n 1 + 2 n1 1 n2 1
2 1 2 2
2
T&T Suggest Degrees of Freedom
Df = smaller of n1-1 and n2-1 More conservative
Calculate standard error
s ' x1 x2 = s s + n1 n2
2 1 2 2
SEdiff = sqrt((var(lead.west)/length(lead.west))+ (var(lead.east)/length(lead.east))) SEdiff
1.029384
Calculate degrees of freedom
s s + n n ' = 12 2 2 2 s12 s2 n n 1 + 2 n1 1 n2 1
2 1 2 2 2
S21 = v1/n1 S22 = v2/n2 dfreedom = trunc(((S21+S22)^2)/ ( (S21^2)/(n1-1) + (S22^2)/(n2-1) ))
dff = function(v1,v2,n1,n2) { S21 = v1/n1 S22 = v2/n2 dfreedom = trunc(((S21+S22)^2)/ ( (S21^2)/(n1-1) + (S22^2)/(n2-1) )) }
Calculate degrees of freedom
s s + n n ' = 12 2 2 2 s12 s2 n n 1 + 2 n1 1 n2 2
2 1 2 2 2
dff(100,100,10,10) 18 dff(v1,v2,n1,n2) 14
Test the Hypothesis H0: The two regions are the same
t' = x1 x2 s s + n1 n2
2 1 2 2
mean(lead.west) - mean(lead.east) -0.7189394 SEdiff 1.029384 dfreedom 14 pt((mean(lead.west) - mean(lead.east))/ SEdiff, dfreedom) 0.2481774
Accept or reject that means are equal?
Cannot Reject H0
The two regions are not statistically significantly different at the alpha=0.05 level
t.test Command
t.test(lead.west,lead.east,alternative="less",var.equal=F) Welch Modified Two-Sample t-Test data: lead.west and lead.east t = -0.6984, df = 14.3789039230301, p-value = 0.248 alternative hypothesis: difference in means is less than 0 95 percent confidence interval: NA 1.090755 sample estimates: mean of x mean of y 3.072727 3.791667
t.test
Student's t-Tests DESCRIPTION: Performs a one-sample, two-sample, or paired t-test, or a Welch modified twosample t-test. USAGE: t.test(x, y=NULL, alternative="two.sided", mu=0, paired=F, var.equal=T, conf.level=.95)
Do we feel comfortable about meeting our assumptions?
Variances unequal, the two-sample t-test with unequal variance calculation of the degrees of freedom accounts for that. What about normality? If the assumption of normality is not satisfied what could we do about that?
BOOTSTRAP !!!
Bootstrap Script
boot.stat = rep(0,1000) for(i in seq(1000)) { x1 = sample(lead.west,length(lead.west),replace=T) x2 = sample(lead.east,length(lead.east),replace=T) boot.stat[i] = mean(x1)-mean(x2) }
Histogram of Bootstrap Distribution of Difference in Means
150 0 50 100
-4
-2 boot.stat
0
2
What alpha level would we need to reject zero?
length(boot.stat[boot.stat>0])/ length(boot.stat) 0.240 Compare with result from t-test 0.248
Conclusions
Blood lead levels in Eastern Washington comparable to those in Western Washington Different variances test suggests roughly same statistic with modified degrees of freedom Note variation on the theme of t-tests and ztests, only we are using F-tests Bootstrap Best thing since sliced bread t statistics pretty robust
Consider
Paired sample hypotheses Computing sample sizes