Statistical Sleuth Chapter 3 - Case 1 & 2 : A Closer Look at Assumptions
Although modern statistical softwares are able to supply confidence intervals and p-values whenever an analysis is conducted, the mathematical assumptions behind the model for two independent samples come from normal populations with the same standard deviation are almost never strictly met in in real world.
At this stage an understanding of actual conditions to evaluate appropriateness of the tools for a particular problem becomes a important issue. Cases in Chapter 3, mainly focuses on deflections from these assumptions, severity of these deviations and the solution…
Case 1 : Does Cloud Seeding Increases Rainfall or No?
The data used in this case were collected between 1968 and 1972 to test injection of silver iodide into cumulus clouds can lead to increased rainfall
case31<- Sleuth2::case0301
tapply(case31$Rainfall,case31$Treatment,summary) #Basic Desciptive
$Unseeded
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 24.82 44.20 164.60 159.20 1203.00
$Seeded
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.10 98.13 221.60 442.00 406.00 2746.00 require(lattice)
bwplot(case0301$Rainfall ~ case0301$Treatment,ylab='Volume of Rainfall per acres')
As it could be seen from the basic descriptives and box plots seeded clouds created substantially higher volume of rain fall. Mean volume of rainfall for seeded days is almost 3 times as large as the unseeded days.
It’s also a good alternative to inspect histogram to see whether two separate treatment groups have similar distributions or not ? In addition a histogram can be a effective tool to observe level of skewness in the data
require(lattice) #We will use histogram function from lattice package to produce a spilt histogram to view both factors.
histogram(~ case0301$Rainfall | case0301$Treatment, xlab='Rainfall volume per acres')
The box plots and histograms of the data in both displays show that both distributions are excessively skewed. But there is more variability occurred in seeded group. Moreover Shapiro test whose results given bellow also confirms that there is no enough evidence to say both of these two distribution is normal.
tapply(case0301$Rainfall,case0301$Treatment,shapiro.test)$Unseeded
Shapiro-Wilk normality test
data: X[[i]]
W = 0.60216, p-value = 3.131e-07
$Seeded
Shapiro-Wilk normality test
data: X[[i]]
W = 0.65626, p-value = 1.411e-06var.test(case0301$Rainfall ~ case0301$Treatment)
F test to compare two variances
data: case0301$Rainfall by case0301$Treatment
F = 0.18304, num df = 25, denom df = 25, p-value = 6.695e-05
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.0820690 0.4082315
sample estimates:
ratio of variances
0.1830387 Moreover, equality of variances test also points that variances are substantially differs from each other. At first look this phenomena means a strong deflection from the assumptions and statistical distribution of our statistics will no longer converges to t distribution. But while stating assumptions, the statistical theory lets some approximations such as the one valid here.
If the skewness in the two populations have the same standart deviations and approximately the same shapes and if the sample sizes are about equal, then validity of t tools is affected moderately by long tailedness (kurtosis) and very little by skewness.
But in order to state this approximation and apply t test, one should initially remove adverse outcomes printed by validity of assumptions tests (variance equality and normality test). At this stage logarithmic calculation offers a very practical solution. By nature any logarithmic transformation series will keep same ordering but smaller observations get spread out more, while large numbers squeezed more closely. Could these distribution be more less skewed in logarithms? Let’s see…
##Transfromation of the data into logarithcmic scale
Log_Rainfall <- log(case0301$Rainfall)
# In R log() operator calculates natural logarithms in default but if you wish to use base you are able to assign it via base = is.numeric argument. tapply(Log_Rainfall,case31$Treatment,summary) #Basic stats for transformed values.$Unseeded
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 3.211 3.786 3.990 5.069 7.092
$Seeded
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.411 4.581 5.396 5.134 6.001 7.918 bwplot(log(case0301$Rainfall) ~ case0301$Treatment, ylab='Log. Transformation of Volume of Rainfall per Acres')
It is clearly seen that logarithmic transformation data suffer less from skewness problem. Maybe a histogram would be a better sketch to compare log transformed volume of rainfalls originated from seeded and unseeded clouds with skewness of a perfectly normal distribution.
library('lattice')
histogram(~ log(case0301$Rainfall) | case0301$Treatment, xlab='Logarithmic Transformed Volume of Rainfall per Acres')
This is graphical presentation, but we can also inspect improvement after logarithmic transformation with t-tools. If we calculate skewness ratios for our raw and logarithmic transformed values it could easily seen that transformation causes substantial improvement. As you can see from the plots, a skewness ratio near “0” is the closest figure to represent a normally distributed data.
# The e1071 package
install.packages('e1071')Error in install.packages : Updating loaded packagesrequire(e1071)
#Skewness of raw data
Skewness<- tapply(case0301$Rainfall,case0301$Treatment,skewness)
#Skewness of log transformed data
log_Skewness<- tapply(log(case0301$Rainfall),case0301$Treatment,skewness)
Table_Skewness <- t(rbind(Skewness,log_Skewness))
print("Improvement in Skewness (after log-transformation)")[1] "Improvement in Skewness (after log-transformation)"Table_Skewness Skewness log_Skewness
Unseeded 2.475611 -0.2036855
Seeded 2.161456 -0.4589744Now we are ready to test our assumptions once more…
shapiro.test(log(case0301$Rainfall))
Shapiro-Wilk normality test
data: log(case0301$Rainfall)
W = 0.98464, p-value = 0.7351var.test(log(case0301$Rainfall) ~ case0301$Treatment)
F test to compare two variances
data: log(case0301$Rainfall) by case0301$Treatment
F = 1.0536, num df = 25, denom df = 25, p-value = 0.8971
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.4724176 2.3499219
sample estimates:
ratio of variances
1.053634 t.test(log(case0301$Rainfall,exp(1)) ~ case0301$Treatment,var.equal= TRUE,alternative='less' ) #order is Unseeded - Seeded so we state 'less'
Two Sample t-test
data: log(case0301$Rainfall, exp(1)) by case0301$Treatment
t = -2.5444, df = 50, p-value = 0.007041
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf -0.3904045
sample estimates:
mean in group Unseeded mean in group Seeded
3.990406 5.134187 diff <- log(case0301$Rainfall[case0301$Treatment=='Seeded']) - log(case0301$Rainfall[case0301$Treatment=='Unseeded'])
mean(diff)[1] 1.143781Regarding on Shapiro and var.test now we can confidently say that our data is suitable for any application of t.tools and we can confidently say that average volume of rainfall differs significantly between seeded and unseeded days.
On our first observation on basic statistics for seeded and unseeded days, we said that mean rainfall volume of seed days is about 3.1 times of unseeded days averages. Taking ratios in consideration we could require to estimate confidence interval for this ratio…
Ratio<- (case0301$Rainfall[case0301$Treatment=='Seeded']/case0301$Rainfall[case0301$Treatment=='Unseeded']) #Computes the ratio of seeded days vs unseeded days
ttest<-t.test(Ratio) #calculates interval in ttest
ttest$conf.int #prints confidence interval estimates[1] 2.874151 3.791230
attr(,"conf.level")
[1] 0.95Case 2 Effects of Agent Orange on Troops in Vietnam - An observational Study
Many Vietnam veterans blames Agent Orange, a herbicide sprayed in South Vietnam between 1962 to 1970 on their declining health status. Recent medical studies have shown that high levels of this dioxin can be detected 20 or more years after heavy exposure. Consequently, CDC has compared 1987 dioxin levels in living Vietnam and non- Vietnam veterans who entered the Army between 1965 and 1987.
Blood sample data from 97 non-Veterans and 646 Veterans are collected. Dioxin concentrations for both factors could be seen from following Box-plot and summary
tapply(case0302$Dioxin,case0302$Veteran,summary)$Vietnam
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 3.00 4.00 4.26 5.00 45.00
$Other
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 3.000 4.000 4.186 5.000 15.000 require(lattice) #Lattice package enables R to draw split box plots.
bwplot(case0302$Dioxin ~ case0302$Veteran,ylab='Dioxin Consentration in Parts per Tril.')
As it can be seen from summary statistics median dioxin levels of non-Vietnam and Vietnam veterans is exactly the same, while means differs minimally (this is caused by presence of outliers, but as number of samples for Vietnam Veterans is 646 we see a much more stabilized effect). So it’s not possible to say data provides enough evidence that Vietnam Veterans has increased level of dioxin in their blood compared Veterans served in USA and West Germany at the same period.
It’s a good alternative to employ a t.test whether these assessments are also valid in scope of a linear model. Here it can be seen that normality assumption doesn’t holds while equality of variances assumption is satisfied.
shapiro.test(case0302$Dioxin)
Shapiro-Wilk normality test
data: case0302$Dioxin
W = 0.65326, p-value < 2.2e-16var.test(case0302$Dioxin ~ case0302$Veteran)
F test to compare two variances
data: case0302$Dioxin by case0302$Veteran
F = 1.318, num df = 645, denom df = 96, p-value = 0.09203
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.9551141 1.7566600
sample estimates:
ratio of variances
1.317991 t.test(case0302$Dioxin[case0302$Veteran=='Vietnam'],case0302$Dioxin[case0302$Veteran=='Other'],alternative='greater',var.equal=TRUE,conf.level=0.95)
Two Sample t-test
data: case0302$Dioxin[case0302$Veteran == "Vietnam"] and case0302$Dioxin[case0302$Veteran == "Other"]
t = 0.26302, df = 741, p-value = 0.3963
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-0.391951 Inf
sample estimates:
mean of x mean of y
4.260062 4.185567 T.test results also confirms our inference. There is not enough evidence to state group means differs from each other. On the other hand, scope of inference for this analysis can also be criticized. Since samples are not random, inference to the population is speculative. For example, non-participating Vietnam veterans may have been failed to participate because of dioxin-related illnesses. If so, statistical inferences could be seriously biased.
Putting aside inference discussions, observing distortions caused by outliers can be a interesting effort for now. As it could be seen from Shapiro’s test our data doesn’t fits to a normal distribution and box plots reflects significant presence of outliers in data. Here we can list outlier values among Vietnam Veterans…
boxplot.stats(case0302$Dioxin[case0302$Veteran=='Vietnam'])$out [1] 9 9 9 9 9 9 9 9 9 9 10 10 11 12 13 15 16 25 4525 and 45 are the obs.645 and obs.646 respectively and they seem as the most dominant outliers. Let us exclude them and calculate two sample t.test statistics once more..
NO_Outlier_Dioxin <- subset(case0302,case0302$Dioxin[case0302$Veteran=='Vietnam']<25) #Data excluding both two outliers
Error in subset(case0302, case0302$Dioxin[case0302$Veteran == "Vietnam"] < :
object 'case0302' not foundBox plot of the data after removing outliers,
And consequently t test statistics after removing both outliers and the largest outlier…
t.test(NO_Outlier_Dioxin$Dioxin~ NO_Outlier_Dioxin$Veteran, alternative='greater' )
Welch Two Sample t-test
data: NO_Outlier_Dioxin$Dioxin by NO_Outlier_Dioxin$Veteran
t = -0.0853, df = 117.33, p-value = 0.5339
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-0.4285707 Inf
sample estimates:
mean in group Vietnam mean in group Other
4.164596 4.185567 t.test(Dioxin_Excl_646$Dioxin~ Dioxin_Excl_646$Veteran, alternative='greater' )
Welch Two Sample t-test
data: Dioxin_Excl_646$Dioxin by Dioxin_Excl_646$Veteran
t = 0.045709, df = 121.27, p-value = 0.4818
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-0.3996048 Inf
sample estimates:
mean in group Vietnam mean in group Other
4.196899 4.185567 It could clearly seen that after removing outliers tendency to reject initial hypothesis lowers substantially. Here we didn’t any change in our decision but it’s a good exhibit to see effects of outlier on t.tools.