Statistics 2MA3 - Assignment #3 Solutions
2001-04-06
Due 2001-04-06 17:00
Full marks = 100
You can submit answers calculated with a pocket calculator or with
R..
Q1 [15]
Considering the small sample size, the histogram is reasonably
close to a normal distribution.
> mean(hosp$temp1)
[1] 98.308
Confidence interval for mean first temperature, assuming standard
deviation = 1 degree: (97.92, 98.70)
> mean(hosp$temp1)-qnorm(.975)*(1/sqrt(25))
[1] 97.916
> mean(hosp$temp1)+qnorm(.975)*(1/sqrt(25))
[1] 98.7
Test the hypothesis that mean first temperature = 98.6, assuming
standard deviation = 1 degree:
> z0 <- (mean(hosp$temp1)-98.6)/(1/sqrt(25))
> z0
[1] -1.46
> 2*pnorm(z0)
[1] 0.1442901
If standard deviation = 1, there is no evidence (p-value =
0.14) that mean first temperature is not 98.6.
Test the hypothesis that standard deviation = 1 degree:
> var(hosp$temp1)
[1] 0.4641
> chsq0 <- 24*var(hosp$temp1)/1
> chsq0
[1] 11.1384
> 2*pchisq(chsq0,24)
[1] 0.0240196
There is significant evidence (p-value = 0.024) that the
standard deviation is not 1 degree.
Confidence interval for mean first temperature: (98.03, 98.59)
> mean(hosp$temp1)-qt(.975,24)*(sqrt(var(hosp$temp1))/sqrt(25))
[1] 98.0268
> mean(hosp$temp1)+qt(.975,24)*(sqrt(var(hosp$temp1))/sqrt(25))
[1] 98.5892
Test the hypothesis that mean first temperature = 98.6:
> t0 <- (mean(hosp$temp1)-98.6)/(sqrt(var(hosp$temp1))/sqrt(25))
> t0
[1] -2.143123
> 2*pt(t0,24)
[1] 0.04244861
or, using t.test()
> t.test(hosp$temp1,mu=98.6)
One Sample t-test
data: hosp$temp1
t = -2.1431, df = 24, p-value = 0.04245
alternative hypothesis: true mean is not equal to 98.6
95 percent confidence interval:
98.0268 98.5892
sample estimates:
mean of x
98.308
There is significant evidence (p-value = 0.042) that the
true mean is not 98.6.
The calculations assume independent, normal data. We can't test
the assumption of independence. The histogram does not indicate
non-normality. Since the hypothesis that standard deviation = 1 was
rejected, the confidence intervals and tests that assumed standard
deviation = 1 will be invalid. It is safe to assume that the other
confidence intervals and tests are valid.
Q2 [2]
> (0.5*(qnorm(0.95)+qnorm(0.99))/(37-38))^2
[1] 3.94261
Rounding up, we conclude that at least n = 4 observations will be
needed.
Q3 [2]
> (0.5*qnorm(0.975)/0.02)^2
[1] 2400.912
Rounding up, we conclude that at least n = 2401 observations will
be needed. This calculation is based on the "worst case" where p =
0.5 and the variance of the estimate is maximized.
Q4 [2]
> qchisq(0.95,19)/qchisq(0.05,19)
[1] 2.979489
> qchisq(0.95,18)/qchisq(0.05,18)
[1] 3.074324
By trial and error, we see that at least 19 degrees of freedom
are needed, that is, n = 19+1 = 20.
Q5 [8]
You have to "read between the lines"; if the hypothesis is that
the mean body temperature is "unaffected" the alternative would seem
to be 2-sided. However, these patients are fevered and the objective
is obviously to reduce the temperature, so a right-tailed test could
also be justified. But since the p-value is extremely small
(p-value = 0.000023 2-sided, p-value = 0.000012
right-tailed) it doesn't affect the conclusion.
> pharma
before after
1 102.4 99.6
2 103.2 100.1
3 101.9 100.2
4 103.0 101.1
5 101.2 99.8
6 100.7 100.2
7 102.5 101.0
8 103.1 100.1
9 102.8 100.7
10 102.3 101.1
11 101.9 101.3
12 101.4 100.2
> t.test(pharma$before,pharma$after,paired=T)
Paired t-test
data: pharma$before and pharma$after
t = 6.9747, df = 11, p-value = 2.346e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.197756 2.302244
sample estimates:
mean of the differences
1.75
We conclude that aspirin does affect mean body temperature
(p-value = 0.000023) or does reduce mean body temperature
(p-value = 0.000012). The assumptions are that the differences
are independent and normal. We can't test independence. The sample is
too small to test normality but the histogram of differences is
consistent with the assumption of normality.
Q6 [3]
There are 12 positive differences out of 12 non-zero differences,
so the exact p-value for the sign test is easily computed from
the binomial distribution. The exact 2-sided p-value is
2*(0.5^12) = 0.000488. The result is very highly significant, so we
reject the hypothesis.
> 2*(1-pbinom(11,12,0.5))
[1] 0.0004882812
The normal approximation to the binomial does not work very well
in this example, but still gives a highly significant p-value.
> 2*(1-pnorm((12-6-0.5)/sqrt(12/4)))
[1] 0.001496164
Q7 [8]
> obstet
bwt group
1 6.9 treatment
2 7.6 treatment
3 7.3 treatment
4 7.6 treatment
5 6.8 treatment
6 7.2 treatment
7 8.0 treatment
8 5.5 treatment
9 5.8 treatment
10 7.3 treatment
11 8.2 treatment
12 6.9 treatment
13 6.8 treatment
14 5.7 treatment
15 8.6 treatment
16 6.4 control
17 6.7 control
18 5.4 control
19 8.2 control
20 5.3 control
21 6.6 control
22 5.8 control
23 5.7 control
24 6.2 control
25 7.1 control
26 7.0 control
27 6.9 control
28 5.6 control
29 4.2 control
30 6.8 control
> boxplot(split(obstet$bwt,obstet$group),ylab="Birth Weight (lb.)")
> bwtctl <- obstet$bwt[obstet$group=="control"]
> bwttrt <- obstet$bwt[obstet$group=="treatment"]
> t.test(bwtctl,bwttrt,var.equal=T)
Two Sample t-test
data: bwtctl and bwttrt
t = -2.4136, df = 28, p-value = 0.02259
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.5159418 -0.1240582
sample estimates:
mean of x mean of y
6.26 7.08
The 2-sample t-test, assuming equal variances, is significant
(p-value = 0.023), indicating that the means are not equal.
Independence can't be tested. the sample sizes are too small to test
normality but the box plots look reasonably symmetric so there is no
reason to reject normality.
Assuming independence and normality, the assumption of equal
variances (homoscedasticity) can be tested with an F-test.
> sqrt(var(bwtctl))
[1] 0.9605058
> sqrt(var(bwttrt))
[1] 0.8993649
> sp <- sqrt((14*var(bwtctl) + 14*var(bwttrt))/28)
> sp
[1] 0.9304377
> F0 <- var(bwtctl)/var(bwttrt)
> F0
[1] 1.140586
> 2*(1-pf(F0,14,14))
[1] 0.8090466
There is no evidence (p-value = 0.81) from these data that
the variances of the control and treatment groups are unequal.
Q8 [4]
> anova(lm(bwt~group,data=obstet))
Analysis of Variance Table
Response: bwt
Df Sum Sq Mean Sq F value Pr(>F)
group 1 5.0430 5.0430 5.8252 0.02259 *
Residuals 28 24.2400 0.8657
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
The F statistic from the ANOVA is the square of the t statistic
> -2.4136^2
[1] -5.825465
The mean squared residual from the ANOVA is just the pooled
variance estimate
> sp^2
[1] 0.8657143
Q9 [2]
The mean squared residual from the ANOVA, which is the same as
the pooled variance estimate, is on 28 df, so the 95% confidence
interval is (0.545, 1.584).
> c((sp^2)/(qchisq(0.975,28)/28), (sp^2)/(qchisq(0.025,28)/28))
[1] 0.5451995 1.5835002
Q10 [8]
> pulmon
bronch lfg
1 20.8 A
2 4.1 A
3 30.0 A
4 24.7 A
5 13.8 A
6 7.5 B
7 7.5 B
8 11.9 B
9 4.5 B
10 3.1 B
11 8.0 B
12 4.7 B
13 28.1 B
14 10.3 B
15 10.0 B
16 5.1 B
17 2.2 B
18 9.2 C
19 2.0 C
20 2.5 C
21 6.1 C
22 7.5 C
> is.factor(pulmon$lfg)
[1] TRUE
> boxplot(split(pulmon$bronch,pulmon$lfg), ylab="Bronchial Reactivity to Sulphur Dioxide", xlab="Lung Function Group")
Lung function Group A, with FEV1/FVC < 75%, has a
much higher reactivity than group B (75-84%) and group C (> 84%).
Group A appears to have a larger variance than the others and the
distributions appear to be positively skewed.
Taking logs (NOT REQUIRED) seems to over-correct; the three
variances are more nearly equal but the distributions appear to be
negatively skewed. The square-root transformation (NOT REQUIRED and
not shown) seems to be about right.
> boxplot(split(pulmon$bronch,pulmon$lfg), ylab="Bronchial Reactivity to Sulphur Dioxide", xlab="Lung Function Group", log="y")
> anova(lm(bronch~lfg,data=pulmon))
Analysis of Variance Table
Response: bronch
Df Sum Sq Mean Sq F value Pr(>F)
lfg 2 503.55 251.77 4.9893 0.01813 *
Residuals 19 958.80 50.46
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
There is significant evidence (p-value = 0.018) that mean
bronchial reactivity is not the same in each lung-function group.
> msr <- anova(lm(bronch~lfg,data=pulmon))$"Mean Sq"[[2]]
> msr
[1] 50.46329
> c(msr/(qchisq(0.975,19)/19), msr/(qchisq(0.025,19)/19))
[1] 29.18522 107.65180
The 95% confidence interval for the residual variance s2 is (29.2, 107.7).
Q11 [10]
> anesthetics$subject<-as.factor(anesthetics$subject)
> anesthetics$anesthesia<-as.factor(anesthetics$anesthesia)
> anesthetics
subject anesthesia epinephrine
1 1 1 0.28
2 1 1 0.36
3 1 2 0.30
4 1 2 0.88
5 1 3 1.07
6 1 3 1.53
7 2 1 0.51
8 2 1 0.32
9 2 2 0.39
10 2 2 0.39
11 2 3 1.35
12 2 3 0.49
> anova(lm(epinephrine~subject*anesthesia, data=anesthetics))
Analysis of Variance Table
Response: epinephrine
Df Sum Sq Mean Sq F value Pr(>F)
subject 1 0.07841 0.07841 0.7074 0.43255
anesthesia 2 1.26762 0.63381 5.7181 0.04075 *
subject:anesthesia 2 0.11502 0.05751 0.5188 0.61968
Residuals 6 0.66505 0.11084
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Assume that the observations are independent and normal, and that
the variance of plasma epinephrine concentration is the same within
each combination of Subject and Anesthetic.
There is no evidence (p-value = 0.62) of an interaction
between Subject and Anesthesia. That is, both subjects responded to
the three anesthesias in similar ways, so we can test the main
effects. There is significant evidence (p-value = 0.041) of a
difference between anesthesias, but no evidence (p-value =
0.43) of a difference between subjects.
> msr <- anova(lm(epinephrine~subject*anesthesia,data=anesthetics))$"Mean Sq"[[4]]
> msr
[1] 0.1108417
> c(msr/(qchisq(0.975,6)/6), msr/(qchisq(0.025,6)/6))
[1] 0.04602621 0.53748179
The 95% confidence interval for the residual variance s2 is (0.0460, 0.5375).
If the study were to be done again, there should be more
subjects. There should still be replication and a test of the
interaction between Subject and Anesthesia because it is possible
that with more subjects someone will respond differently to the
different anesthesias.
Q12 [6]
> plot(hosp$age[-7],hosp$duration[-7])
> fith1 <- lm(duration~age,data=hosp[-7,])
> summary(fith1)
Call:
lm(formula = duration ~ age, data = hosp[-7, ])
Residuals:
Min 1Q Median 3Q Max
-5.4951 -2.0849 0.1927 1.8022 9.1326
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.17552 1.63336 3.169 0.00445 **
age 0.06260 0.03629 1.725 0.09854 .
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 3.506 on 22 degrees of freedom
Multiple R-Squared: 0.1191, Adjusted R-squared: 0.07911
F-statistic: 2.976 on 1 and 22 degrees of freedom, p-value: 0.09854
> abline(fith1)
> anova(fith1)
Analysis of Variance Table
Response: duration
Df Sum Sq Mean Sq F value Pr(>F)
age 1 36.573 36.573 2.9758 0.09854 .
Residuals 22 270.385 12.290
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Assume that the data are independent, that the true relationship
is linear, and that the conditional distribution of duration, given
age, is normal, with variance independent of age.
There is no evidence (p-value = 0.099) from these data
that duration of stay depends on age. (Compare the corresponding
analysis in Exercise #4 and note that removing the observation makes
very little change to the results; the p-value is increased
from 0.074 to 0.099.)
Q13 [12]
If you don't want to delete an extreme observation, taking logs
can reduce its impact on the analysis. The intention here was that
you would fit log-transformed data and not remove the 7th subject.
> hosp$logdur <- log(hosp$duration)
> dimnames(hosp)[[2]]
[1] "id" "duration" "age" "sex" "temp1"
[6] "wbc1" "antib" "bact" "serv" "logdur"
> cor(hosp[,c(10,3,5,6)])
logdur age temp1 wbc1
logdur 1.00000000 0.3515604 0.1846740 -0.01404079
age 0.35156041 1.0000000 -0.3816794 -0.36932443
temp1 0.18467398 -0.3816794 1.0000000 0.42354776
wbc1 -0.01404079 -0.3693244 0.4235478 1.00000000
> pairs(hosp[,c(10,3,5,6)])
Note that first white blood cell count is correlated with first
temperature but not with log(duration), so if temperature is in the
model then adding white blood cell count will not improve the fit
significantly. Similarly, adding first temperature to the model after
age model does not improve the fit significantly.
> fith2 <- lm(logdur~age+temp1+wbc1, data=hosp)
> anova(fith2)
Analysis of Variance Table
Response: logdur
Df Sum Sq Mean Sq F value Pr(>F)
age 1 0.9599 0.9599 3.4275 0.07824 .
temp1 1 0.9243 0.9243 3.3002 0.08357 .
wbc1 1 0.0011 0.0011 0.0038 0.95116
Residuals 21 5.8815 0.2801
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> plot(fitted(fith2),hosp$logdur)
> abline(0,1)
> plot(fitted(fith2),residuals(fith2))
> abline(h=0)
The ANOVA and the two plots indicate that a multiple linear
regression model is reasonable, but does not give a significant
improvement over the grand mean.
Since the instructions weren't clear, you could also analyze
taking logs AND removing the 7th observation.
> cor(hosp[-7,c(10,3,5,6)])
logdur age temp1 wbc1
logdur 1.00000000 0.2989780 -0.00494539 -0.02268732
age 0.29897801 1.0000000 -0.49544466 -0.37858725
temp1 -0.00494539 -0.4954447 1.00000000 0.45080402
wbc1 -0.02268732 -0.3785873 0.45080402 1.00000000
> pairs(hosp[-7,c(10,3,5,6)])
Note that first white blood cell count is correlated with first
temperature but not with log(duration), so if temperature is in the
model then adding white blood cell count will not improve the fit
significantly. Similarly, adding first temperature to the model after
age is in the model does not improve the fit significantly.
Removing the 7th subject as well as taking log(duration) makes
age and first temperature even less significant.
> fith2a <- lm(logdur~age+temp1+wbc1, data=hosp[-7,])
> anova(fith2a)
Analysis of Variance Table
Response: logdur
Df Sum Sq Mean Sq F value Pr(>F)
age 1 0.5084 0.5084 2.0286 0.1698
temp1 1 0.1545 0.1545 0.6166 0.4415
wbc1 1 0.0122 0.0122 0.0488 0.8274
Residuals 20 5.0123 0.2506
> plot(fitted(fith2a),hosp$logdur[-7])
> abline(0,1)
> plot(fitted(fith2a),residuals(fith2a))
> abline(h=0)
We assume that the data are independent and normal and that the
effects of age, first temperature and first white blood cell count
are linear.
We conclude that the dependence of log(duration of stay) on age,
first temperature and first white blood cell count is too weak to
have any predictive power.
Q14 [4]
If we fit age after fitting first temperature the improvement in
the fit is significant (p-value = 0.026) for predicting
log(duration of stay) if the 7th observation is included, but (NOT
REQUIRED) not if it is excluded.
> anova(lm(logdur~temp1+age+wbc1, data=hosp))
Analysis of Variance Table
Response: logdur
Df Sum Sq Mean Sq F value Pr(>F)
temp1 1 0.2649 0.2649 0.9458 0.3419
age 1 1.6194 1.6194 5.7819 0.0255 *
wbc1 1 0.0011 0.0011 0.0038 0.9512
Residuals 21 5.8815 0.2801
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> anova(lm(logdur~temp1+age+wbc1, data=hosp[-7,]))
Analysis of Variance Table
Response: logdur
Df Sum Sq Mean Sq F value Pr(>F)
temp1 1 0.0001 0.0001 0.0006 0.9814
age 1 0.6628 0.6628 2.6446 0.1196
wbc1 1 0.0122 0.0122 0.0488 0.8274
Residuals 20 5.0123 0.2506
Q15 [10]
> cattle
bodywt mcr bodywtf
1 110 235 110
2 110 198 110
3 110 173 110
4 230 174 230
5 230 149 230
6 230 124 230
7 360 115 360
8 360 130 360
9 360 102 360
10 360 95 360
11 505 122 505
12 505 112 505
13 505 98 505
14 505 96 505
> fitc1 <- lm(mcr~bodywt, data=cattle)
> plot(cattle$bodywt,cattle$mcr)
> abline(fitc1)
> anova(fitc1)
Analysis of Variance Table
Response: mcr
Df Sum Sq Mean Sq F value Pr(>F)
bodywt 1 16634.2 16634.2 27.567 0.0002043 ***
Residuals 12 7241.0 603.4
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> anova(lm(mcr~bodywtf, data=cattle))
Analysis of Variance Table
Response: mcr
Df Sum Sq Mean Sq F value Pr(>F)
bodywtf 3 19514.2 6504.7 14.916 0.0005067 ***
Residuals 10 4361.0 436.1
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> anova(lm(mcr~bodywt+bodywtf, data=cattle))
Analysis of Variance Table
Response: mcr
Df Sum Sq Mean Sq F value Pr(>F)
bodywt 1 16634.2 16634.2 38.143 0.0001047 ***
bodywtf 2 2880.0 1440.0 3.302 0.0792380 .
Residuals 10 4361.0 436.1
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Assume that the data are independent and normal and the variance
in mcr does not depend on body weight.
Although the plot suggests that mcr levels out when body weight
is greater than 350, the non-linearity is not significant
(p-value = 0.079). Hence we can test the slope, which is very
highly significant (p-value = 0.0001).
Interpolating on the fitted line, a body weight of 300 kg would
predict a metabolic clearance rate of 142.1.
> predict(fitc1,data.frame(bodywt=300))
[1] 142.0674
> msrc1 <- anova(fitc1)$"Mean Sq"[[2]]
> msrc1
[1] 603.4178
> c(msrc1/(qchisq(0.975,12)/12), msrc1/(qchisq(0.025,12)/12))
[1] 310.2848 1644.2691
> msrc2 <- anova(lm(mcr~bodywtf,data=cattle))$"Mean Sq"[[2]]
> msrc2
[1] 436.1
> c(msrc2/(qchisq(0.975,10)/10), msrc2/(qchisq(0.025,10)/10))
[1] 212.9064 1343.0972
The 95% confidence interval for the residual variance is (310,
1644) if we assume a linear regression and use the MSR obtained after
fitting the straight line. It is (213, 1343) if we do not assume a
linear regression.
Q16 [4]
> library(ctest)
> std <- matrix(c(40,10,15,30,20,40,130,70,45),nrow=3)
> dimnames(std) <- list(c("Penicillin","Spect.low","Spect.high"),c("Spos","Sneg.Cpos","Sneg.Cneg"))
> std
Spos Sneg.Cpos Sneg.Cneg
Penicillin 40 30 130
Spect.low 10 20 70
Spect.high 15 40 45
> chisq.test(std)$expected
Spos Sneg.Cpos Sneg.Cneg
Penicillin 32.50 45.0 122.50
Spect.low 16.25 22.5 61.25
Spect.high 16.25 22.5 61.25
> chisq.test(std)
Pearson's Chi-squared test
data: std
X-squared = 29.1401, df = 4, p-value = 7.322e-06
There is highly significant evidence (p-value = 0.0000007)
that the response depends on which of the three treatments was used.
In particular, penicillin gave more (- smear, + culture) results than
expected under independence, while spectinomycin (high dose) gave
more (- smear, - culture) than expected.