STAT1000 Regression & Nonparametric Inference

http://www.industrytap.com/making-heavens-open-cloud-seeding-drones/22011

(From Ramsey, F.L. and Schafer, D. W. (2013). The Statistical Sleuth, Brooks/Cole: Boston, MA.)

An experiment was conducted in southern Florida between 1968 and 1972 to test a hypothesis that massive injection of silver iodide (AgI) into cumulus clouds could lead to increased rainfall. On each of 52 days that were deemed suitable for cloud seeding, a random mechanism was used to decide whether to seed the target cloud on that day or to leave it unseeded as a control. An airplane flew through the cloud in both cases, since the experimenters and the pilot themselves were aware of whether on any particular day the seeding mechanism in the plane was loaded or not [with AgI]. Precipitation was measured as the total rain volume falling from the cloud base following the airplane seeding run, as measured by radar, in units of acre-feet.

Your task is to determine, using a subset of the data, whether cloud seeding had an effect on rainfall.

The R data file STAT1000Test1.RData contains the data frame RainSeed, which consists of rainfall measurements (Rainfall) for the two treatment conditions (Seeded and Unseeded). I have also extracted the rainfall amounts corresponding to these two conditions in separate vectors of the same name.

rm(list = ls()) # This removes any other objects in your workspace
load("S2_2016_STAT1000Test1.RData")

1. Plot the rainfall data in two different ways that give you some sense of their distribution. Make sure your plots are appropriately labelled.

# Plot 1
boxplot(Rainfall ~ Treatment, data = RainSeed)

# Plot 2
hist(Seeded, ylab= 'rainfall', main = 'seeded rainfall histogram')

hist(Unseeded)

2. Comment on the plots that you produced above. What, if any, noticeable features are present?

Answer: The boxplot optically suggests that seeded days have a higher rainfall, and also has a larger spread.

The histograms also support this, with data recorded of up to 2000 for seeded but only 1400 for unseeded.

There are obvious outliers for both data sets.

3. Your boss asks you to carry out a \(t\)-test to determine whether seeding led to increased rainfall. Explain to your boss why a \(t\)-test might not be appropriate here.

Answer: The data has several obvious outliers which make a t-test innapropriate.

4. Begrudgingly, you carry out the \(t\)-test.

# Carry out t-test here

1. What hypotheses would you be testing?

Answer: The null hypothesis will be that both data sets have the same mean.

The alternative will be that seeded has a greater mean than unseeded.

2. What are the results of the \(t\)-test?

 t.test(Seeded, Unseeded, alternative ='greater')

    Welch Two Sample t-test

data:  Seeded and Unseeded
t = 1.0278, df = 23.637, p-value = 0.1572
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 -108.3462       Inf
sample estimates:
mean of x mean of y 
 347.5615  184.8000 

5. You explain to your boss that a rank-based nonparametric procedure might be more appropriate.

1. Which one would you choose?

Answer:rank-sum test

2. State the hypotheses being tested, carry out the test showing all intermediate steps, and state your conclusion.

Answer:Ho = uS ~= uU

Ha = uS > uU

# Carry out rank-based test here
U <- tapply(rank(RainSeed$Rainfall), RainSeed$Treatment, sum) - 13*7
U

  Seeded Unseeded 
     137       32 

pval <- pwilcox(32,13,13)
pval

[1] 0.003037517

3. Check your result using the built-in R function for the test you’ve chosen.

# Check result here
wilcox.test(Seeded,Unseeded, alternative = 'greater')

    Wilcoxon rank sum test

data:  Seeded and Unseeded
W = 137, p-value = 0.003038
alternative hypothesis: true location shift is greater than 0

Answer:

4. Does it yield a different result from the \(t\)-test? Explain why.

Answer: Yes, as the median is not as susceptible to outliers as the mean is (the rank will remain the same no matter how large or small certain outliers are)

6. You decide to try a permutation test, but because there are many possible permutations, you randomly choose only 10000 of them for the test. The distribution (PermDistMeanSample) is based on the difference ‘Seeded’ - ‘Unseeded’.

1. What is the \(p\)-value of the permutation test?

# Carry out calculations here
meandiff <- mean(Seeded) - mean(Unseeded)
mean(PermDistMeanSample >= meandiff)

[1] 0.1728

2. How might you modify the permutation test so that it yields a \(p\)-value closer to that of the rank-based test you performed above?

Answer: