Stat209/Ed260 D Rogosa 2/26/18
Solutions Assignment 8.
Matching review: Randomized blocks
Problem 1.
>
> meds = read.table(file="http://web.stanford.edu/~rag/stat209/depressdata", header = T)
> names(meds)
[1] "Imipramine" "Placebo"
> attach(meds)
> diff = Imipramine - Placebo
> stem(diff)
The decimal point is at the |
-6 | 000
-4 | 00
-2 | 0000000000
-0 | 000
0 | 00000
2 | 000000
4 | 0
> qqnorm(diff) # not enough data to conclude anything but useful to look at
> t.test(diff)
One Sample t-test
data: diff
t = -2.3731, df = 29, p-value = 0.02449
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-2.3583370 -0.1749963
sample estimates:
mean of x
-1.266667
> # reject Ho no differential effect at Type I error .05, (sidenote : appears drug has a positive effect)
>
(b)
> t.test(Imipramine,Placebo) # carry out the two sample test on the coloumns or create a grouping variable
Welch Two Sample t-test
data: Imipramine and Placebo
t = -1.9646, df = 57.628, p-value = 0.05429
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.55746123 0.02412790
sample estimates:
mean of x mean of y
6.300000 7.566667
> # loss of pairing information is quite consequential; reduces precision
> # of comparison so that a difference between drug and placebo can no
> # longer be detected.
> # full credit requires carrying out test and noting the difference from part a.
>
(c)
> # problem asks us to do a sign test.
> sort(diff) #easy way to count for sign test
[1] -6 -6 -6 -5 -5 -3 -3 -3 -3 -3 -3 -3 -2 -2 -2 -1 -1 -1 0 0 0 0 1 2 2
[26] 2 2 3 3 5
> sum(diff > 0); sum(diff < 0) # 18 negs 8 positives (toss 0's), let the computer count
[1] 8
[1] 18
#Refer Treatment of Zeros, page 367 of SW
> binom.test(8,26)
Exact binomial test
data: 8 and 26
number of successes = 8, number of trials = 26, p-value = 0.07552
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1432600 0.5178964
sample estimates:
probability of success
0.3076923
> # can no longer reject at Type I .05 because we had to throw
> # away info about magnitude of paired differences.
>
----------------------------------------------------------
Problem 2.
Matching to increase precision: Randomized blocks designs
> dblock = read.table(file="[path]dental.dat", header = T)
> attach(dblock)
# get 2x2 cell means in either order
> tapply(Pain,list(as.factor(Acup),as.factor(Codeine)), mean)
1 2
1 0.600 1.0625
2 1.175 1.7875
> tapply(Pain,list(as.factor(Codeine),as.factor(Acup)), mean)
1 2
1 0.6000 1.1750
2 1.0625 1.7875
> blkaov = aov(Pain ~ as.factor(Block) + as.factor(Codeine) + as.factor(Acup) +
+ (as.factor(Codeine)*as.factor(Acup)))
> summary(blkaov)
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(Block) 7 5.5987 0.7998 55.2963 4.126e-12 ***
as.factor(Codeine) 1 2.3112 2.3112 159.7901 2.773e-11 ***
as.factor(Acup) 1 3.3800 3.3800 233.6790 7.465e-13 ***
as.factor(Codeine):as.factor(Acup) 1 0.0450 0.0450 3.1111 0.0923 .
Residuals 21 0.3037 0.0145
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> # now ignore blocking variable
> blknoaov = aov(Pain ~ as.factor(Codeine) + as.factor(Acup) +
+ (as.factor(Codeine)*as.factor(Acup)))
> summary(blknoaov)
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(Codeine) 1 2.3112 2.3112 10.9640 0.0025659 **
as.factor(Acup) 1 3.3800 3.3800 16.0339 0.0004155 ***
as.factor(Codeine):as.factor(Acup) 1 0.0450 0.0450 0.2135 0.6476320
Residuals 28 5.9025 0.2108
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> # so the 2x2 with 8 replications per cell (no blocking version) has
> # same SS and MS for the experimental factors but Mean square residual
> # is approx 15 times larger, .2108/.0145. So relative efficiency is
> # around 15 (there are more complex versions of rel eff). Instead of
> # 32 subjects, would need at least 400 to obtain same precision on
> # experimental effects. In this ex blocking is valuable, as an approximate
> # matching (on pain tolerance)
>
------------------------------------------------------------