Assignment

Assignment 5

Due Thursday October 20 at 11:59pm on Blackboard

As before, the questions without solutions are an assignment: you need to do these questions yourself and hand them in (instructions below).

The assignment is due on the date shown above. An assignment handed in after the deadline is late, and may or may not be accepted (see course outline). My solutions to the assignment questions will be available when everyone has handed in their assignment.

You are reminded that work handed in with your name on it must be entirely your own work.

Assignments are to be handed in on Blackboard. Instructions are at http://www.utsc.utoronto.ca/

~butler/c32/blackboard-assgt-howto.html, in case you forgot since last week. Markers’ comments and grades will be available on Blackboard as well.

1. I said before that obtaining a ggplot normal quantile plot with a line is not automatic, but let’s explore how it might be done, since the ideas are not difficult ones.

(a) This question assesses the normality of a chi-squared distribution. (Don’t worry if you’ve never heard of the chi-squared distribution before. It’s the one that is used to obtain P-values for any of the various chi-squared tests. If you’ve done any of those you might have used a chi-squared table. If not, it’s fine.) R has a function rchisq that generates random values from this distribution. It takes two things as input, the number of random values to generate, and the number of degrees of freedom. A chi-squared distribution with large degrees of freedom is more normal-like, in the same way that a t distribution with large degrees of freedom is indistinguishable from a normal distribution.

Generate 50 random values from a chi-squared distribution with 5 degrees of freedom, and save them into a data frame (containing just that one column of values).

Solution: First, I need to make sure that my answer is reproducible (so I can talk about it), plus I need the “tidyverse” for later:

set.seed(457299)

library(tidyverse)

## Loading tidyverse: ggplot2

## Loading tidyverse: tibble

## Loading tidyverse: tidyr

## Loading tidyverse: readr

## Loading tidyverse: purrr

## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats

## lag(): dplyr, stats

1

The obvious way is to do it in two steps:

z=rchisq(50,5)

df=data.frame(z)

summary(df)

## z

## Min. : 0.9193

## 1st Qu.: 2.6898

## Median : 4.9005

## Mean : 5.9352

## 3rd Qu.: 8.1580

## Max. :20.2781

You can also do it in one shot:

df=data.frame(z=rchisq(50,5))

summary(df)

## z

## Min. : 0.7468

## 1st Qu.: 3.0280

## Median : 4.3933

## Mean : 4.7460

## 3rd Qu.: 5.9517

## Max. :18.5339

The variable names are of no consequence, but it will be less confusing to avoid x and y because of what is coming up.

(A chi-squared distribution has mean equal to its degrees of freedom, so the mean should be 5. That looks OK here.)

(b) On a normal quantile plot such as the one produced by qqnorm-qqline, the line goes through the first and third quartiles of the data (on the y-axis) and the first and third quartiles of a standard normal distribution (on the x-axis). Calculate these, calling them y and x respectively.

Solution: For the data:

y=quantile(df$z,c(0.25,0.75))

y

## 25% 75%

## 3.028047 5.951711

For the standard normal distribution, the function qnorm produces the values of z that have the given probabilities of being less than them (the “inverse CDF”, if you will, or “reading the table backwards”):

x=qnorm(c(0.25,0.75))

x

## [1] -0.6744898 0.6744898

(c) Work out the slope and intercept of the straight line joining these two points. (The slope is rise over run; the intercept is whatever it has to be to make a line with the slope you just calculated pass through one of the points; it doesn’t matter which one.) See if you can find a lazy way of

Page 2

getting the slope and intercept, and justify why it works.

Solution: Let’s suppose the line has the formula y = a + bx. We first get b as rise over run:

b=(y[2]-y[1])/(x[2]-x[1])

b

## 75%

## 2.167315

then we stop and think a bit. Pick one of the points, say the first one (and use that for x and y); we just found b, so the only unknown thing is a. A tiny amount of algebra produces a = y − bx:

a=y[1]-b*x[1]

a

## 25%

## 4.489879

I said there was a lazy way. For that, note that the regression line through two points must actually go through those two points (so that the sum of squares of residuals can be and is zero).1 In R, lm does regression, or you could even do it by hand if you were determined enough:

yy=lm(y~x)

summary(yy)

##

## Call:

## lm(formula = y ~ x)

##

## Residuals:

## ALL 2 residuals are 0: no residual degrees of freedom!

##

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 4.490 NA NA NA

## x 2.167 NA NA NA

##

## Residual standard error: NaN on 0 degrees of freedom

## Multiple R-squared: 1,Adjusted R-squared: NaN

## F-statistic: NaN on 1 and 0 DF, p-value: NA

Same answers. The other things are all missing, because the residual sum of squares is zero (but we don’t worry about that).

(d) Now we make our ggplot normal quantile plot with line. The key is to note that there is also geom_abline that adds to the plot a line with intercept and slope that you specify. Make a ggplot normal quantile plot with the line you calculated above through the points.

Solution: Like this. I’m using my a and b, but you can type in the numbers that you found if you prefer.

ggplot(df,aes(sample=z))+stat_qq()+

geom_abline(intercept=a,slope=b)

Page 3

● ● ●

●

● ●

● ● ●

● ●

● ● ●

● ● ● ●

● ●●

● ●

●●

●●● ●●

● ● ●

● ● ●

● ●

● ● ● ●

● ●

●

●

●

● ●

●

0

5

10

15

−2 −1 0 1 2 theoretical

sa m

p le

(e) Does this chi-squared distribution look normal? If not, describe how it looks relative to the normal.

Solution: With the line, this is easier to judge. If you used a different seed to me (or none at all), you’ll have a different picture, and you can use yours. With mine, I can go a couple of ways: (i) I could say that the distribution is approximately normal apart from the outlier at the high end, or (ii) all the points that are not near the line are above it, at the bottom and top ends, so the distribution is bunched up at the bottom and spread out at the top, which is to say skewed to the right.

I happen to know that (ii) is actually the right answer, but my random data made it hard to distinguish between (i) and (ii). (A chi-squared distribution never goes below zero, so I expected a bigger straggle of points just above zero than I actually got.)

2. A biologist wished to study the effects of ethanol on sleep time in rats. A sample of 20 rats (all the same age) was selected, and each rat was given an injection having a particular concentration (0, 1, 2 or 4

Page 4

grams per kilogram of body weight) of ethanol. These are labelled e0, e1, e2, e4. The “0” treatment was a control group. The rapid eye movement (REM) sleep time was then recorded for each rat. The data are in http://www.utsc.utoronto.ca/~butler/c32/ratsleep.txt. Use R for this question.

(a) Read the data in from the file. Check that you have four rows of observations and five columns of sleep times.

Solution:

sleep1=read.table("ratsleep.txt",header=T)

sleep1

## treatment obs1 obs2 obs3 obs4 obs5

## 1 e0 88.6 73.2 91.4 68.0 75.2

## 2 e1 63.0 53.9 69.2 50.1 71.5

## 3 e2 44.9 59.5 40.2 56.3 38.7

## 4 e4 31.0 39.6 45.3 25.2 22.7

There are six columns, but one of them labels the groups, and there are correctly five columns of sleep times.

(b) Unfortunately, the data are in the wrong format. All the sleep times for each treatment group are on one row, and we should have one column containing all the sleep times, and the corresponding row should show which treatment group that sleep time came from.

If you prefer to skip this part: read in the data from http://www.utsc.utoronto.ca/~butler/ c32/ratsleep2.txt, and proceed to the boxplots in (c).

First install the package tidyr by typing this in the console window:

install.packages("tidyr")

(you only have to do this once), and then in every R session that you want to use tidyr, you type

library(tidyr)

Or you can load the “tidyverse” as per last week:

library(tidyverse)

which loads tidyr and dplyr and ggplot and other handy things.

The tidyr function gather turns wide format (which we have) into long format (which we want). gather needs four things fed into it: a data frame, what makes the columns different, what makes them the same, and finally which columns are to be gathered together (combined into one column), the first one, a colon, and the last one. See my solution below if you don’t want to try it yourself. Save the result of gather into a data frame, and look at it. Do you have 20 rows of not-very-many variables?

Solution: What makes the columns obs1 through obs5 different is that they are different observation numbers (“replicates”, in the jargon). I’ll call that rep. What makes them the same is that they are all sleep times. Columns obs1 through obs5 are the ones we want to combine, thus. I already did the install.packages thing. In fact, since I did library(tidyverse) above, my library(tidyr) is redundant, since tidyr is part of the tidyverse:

Page 5

library(tidyr)

sleep=gather(sleep1,rep,sleeptime,obs1:obs5)

sleep

## treatment rep sleeptime

## 1 e0 obs1 88.6

## 2 e1 obs1 63.0

## 3 e2 obs1 44.9

## 4 e4 obs1 31.0

## 5 e0 obs2 73.2

## 6 e1 obs2 53.9

## 7 e2 obs2 59.5

## 8 e4 obs2 39.6

## 9 e0 obs3 91.4

## 10 e1 obs3 69.2

## 11 e2 obs3 40.2

## 12 e4 obs3 45.3

## 13 e0 obs4 68.0

## 14 e1 obs4 50.1

## 15 e2 obs4 56.3

## 16 e4 obs4 25.2

## 17 e0 obs5 75.2

## 18 e1 obs5 71.5

## 19 e2 obs5 38.7

## 20 e4 obs5 22.7

We have 20 rows of 3 columns. The column rep is not very interesting: it just says which observation each one was within its group.2 The interesting things are treatment and sleeptime, which are the two variables we’ll need for our analysis of variance.

(c) Using your new data frame, make side-by-side boxplots of sleep time by treatment group.

Solution:

ggplot(sleep,aes(x=treatment,y=sleeptime))+geom_boxplot()

Page 6

20

40

60

80

e0 e1 e2 e4 treatment

sl e

e p

tim e

(d) In your boxplots, how does the median sleep time appear to depend on treatment group?

Solution: It appears to decrease as the dose of ethanol increases, and pretty substantially so (in that the differences ought to be significant, but that’s coming up).

(e) There is an assumption about spread that the analysis of variance needs in order to be reliable. Do your boxplots indicate that this assumption is satisfied for these data, bearing in mind that you have only five observations per group?

Solution: The assumption is that the population SDs of each group are all equal. Now, the boxplots show IQRs, which are kind of a surrogate for SD, and because we only have five observations per group to base the IQRs on, the sample IQRs might vary a bit. So we should look at the heights of the boxes on

Page 7

the boxplot, and see whether they are grossly unequal. They appear to be to be of very similar heights, all things considered, so I am happy.

If you want the SDs themselves, aggregate will get them for you. Don’t forget the name of the data frame before the sd:

aggregate(sleeptime~treatment,sleep,sd)

## treatment sleeptime

## 1 e0 10.180963

## 2 e1 9.344143

## 3 e2 9.461078

## 4 e4 9.558923

Those are very similar, given only 5 observations per group. No problems here.

If you have the dplyr group-by and summarize ideas in your head, that works too (since I loaded dplyr with the “tidyverse” above):

tmp=group_by(sleep,treatment)

summarize(tmp,stddev=sd(sleeptime))

## # A tibble: 4 × 2 ## treatment stddev

## <fctr> <dbl>

## 1 e0 10.180963

## 2 e1 9.344143

## 3 e2 9.461078

## 4 e4 9.558923

(f) Run an analysis of variance to see whether sleep time differs significantly among treatment groups. What do you conclude?

Solution: I use aov here, because I might be following up with Tukey in a minute:

sleep.1=aov(sleeptime~treatment,data=sleep)

summary(sleep.1)

## Df Sum Sq Mean Sq F value Pr(>F)

## treatment 3 5882 1961 21.09 8.32e-06 ***

## Residuals 16 1487 93

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is a very small P-value, so my conclusion is that the mean sleep times are not all the same for the treatment groups. Further than that I am not entitled to say (yet).

The technique here is to save the output from aov in something, look at that (via summary), and then that same something gets fed into TukeyHSD later.

(g) Would it be a good idea to run Tukey’s method here? Explain briefly why or why not, and if you think it would be a good idea, run it.

Solution: Tukey’s method is useful when (a) we have run an analysis of variance and got a significant result and (b) when we want to know which groups differ significantly from which. Both (a) and (b) are

Page 8

true here. So, noting the uppercase-lowercase:

TukeyHSD(sleep.1)

## Tukey multiple comparisons of means

## 95% family-wise confidence level

##

## Fit: aov(formula = sleeptime ~ treatment, data = sleep)

##

## $treatment

## diff lwr upr p adj

## e1-e0 -17.74 -35.18636 -0.2936428 0.0455781

## e2-e0 -31.36 -48.80636 -13.9136428 0.0005142

## e4-e0 -46.52 -63.96636 -29.0736428 0.0000056

## e2-e1 -13.62 -31.06636 3.8263572 0.1563545

## e4-e1 -28.78 -46.22636 -11.3336428 0.0011925

## e4-e2 -15.16 -32.60636 2.2863572 0.1005398

(h) What do you conclude from Tukey’s method? (This is liable to be a bit complicated.) Is there a treatment that is clearly best, in terms of the sleep time being largest?

Solution: All the differences are significant except treatment e2 vs. e1 and e4. All the differences involving the control group e0 are significant, and if you look back at the boxplots in (c), you’ll see that the control group e0 had the highest mean sleep time. So the control group is best (from this point of view), or another way of saying it is that any dose of ethanol is significantly reducing mean sleep time.

The other comparisons are a bit confusing, because the 1-4 difference is significant, but neither of the differences involving 2 are. That is, 1 is better than 4, but 2 is not significantly worse than 1 nor better than 4. This seems like it should be a logical impossibility, but the story is that we don’t have enough data to decide where 2 fits relative to 1 or 4. If we had 10 or 20 observations per group, we might be able to conclude that 2 is in between 1 and 4 as the boxplots suggest.

3. Does the age of a sweet white wine affect its acidity? A winemaker took random samples of Riesling (a sweet white wine) from casks that are 15, 20 and 25 years old and measured the acidity of each one. The results are in http://www.utsc.utoronto.ca/~butler/c32/wine-acidity.txt. The two columns are the age of the wine and the acidity.

Note that the data are in the correct format, so we don’t need to do the tidyr stuff (or SAS equivalent thereof).

(a) Read the data into SAS and confirm, without listing all the data values, that all of the acidity values are between 0.7 and 0.85.

Solution: The trick here is to remember that proc means also produces the maximum and minimum, which is exactly what we need here:

data wine;

infile '/home/ken/wine-acidity.txt';

input age acidity;

with this output:

proc means;

Page 9

var acidity;

The MEANS Procedure

Analysis Variable : acidity

N Mean Std Dev Minimum Maximum

------------------------------------------------------------------

39 0.8058615 0.0151667 0.7735000 0.8298000

------------------------------------------------------------------

Page 10

The smallest value is 0.77 and the largest 0.83, so all the acidity values must be between 0.7 and 0.85.

An alternative here is to use proc univariate and look at the “extreme values” down near the bottom of the output. It works, but you need to make sure you only grab the table of extreme values (and not the rest of it), because the output from proc univariate is rather long:

proc univariate;

var acidity;

The UNIVARIATE Procedure

Variable: acidity

Extreme Observations

-----Lowest----- -----Highest----

Value Obs Value Obs

0.7735 27 0.8251 13

0.7789 34 0.8262 21

0.7801 38 0.8265 20

0.7813 28 0.8291 3

0.7843 39 0.8298 5

Again, the smallest value is bigger than 0.7 and the largest is smaller than 0.85, so all the values must be between 0.7 and 0.85.

Listing all the data values (with proc print) is smart until you have convinced yourself that the values are correct, but it is not smart to make the marker read it every time (unless the data set is small).

I have 39 data values altogether (which is getting a bit long to list all of), so I imagine there are 13 per group. You don’t need to, but I want to check:

proc freq;

table age;

The FREQ Procedure

Cumulative Cumulative

age Frequency Percent Frequency Percent

--------------------------------------------------------

15 13 33.33 13 33.33

20 13 33.33 26 66.67

25 13 33.33 39 100.00

I do. Or, of course, proc means with class age would have told me the same thing.

(b) Obtain side by side boxplots of acidity levels for each age of wine.

Solution: This is just like the one you did in the tutorial:

proc sgplot;

Page 11

vbox acidity / category=age;

(c) What do you learn from the boxplots? (Two things, one about means/medians and one about spreads. For the second, why is what you learn important?).

Solution: The mean and median acidity is much lower for the oldest wines (age 25). The age-15 and age-20 wines have about the same median acidity.

The groups each have similar spreads, and there are no outliers and more or less symmetric shapes in all cases. (The age 20 wines are the most questionable in terms of symmetry.) This is important because an assumption of ANOVA is that you have (approximately) normal data with (approximately) equal SDs, and if we have that, as I think we do, we can trust the results of the ANOVA. As ever, you can ask for perfection, but you will generally be disappointed. If it’s not bad enough to cause problems, that’s all we need.

An alternative (and thus acceptable) answer is to say that the age-15 wines have “clearly smaller” IQR, and therefore we should be cautious about the results of the ANOVA. I don’t like this as much, but if that’s your reasoning, that’s a valid conclusion.

(d) Obtain an analysis of variance. What, as precisely as possible, are you able to conclude?

Solution: It looks as if there ought to be a var line in there, but you don’t need one (what would be on the var line is on the left side of the equals in the model line):

proc anova;

class age;

Page 12

model acidity=age;

Page 13

The ANOVA Procedure

Class Level Information

Class Levels Values

age 3 15 20 25

Number of Observations Read 39

Number of Observations Used 39

The ANOVA Procedure

Dependent Variable: acidity

Sum of

Source DF Squares Mean Square F Value Pr > F

Model 2 0.00492012 0.00246006 23.18 <.0001

Error 36 0.00382102 0.00010614

Corrected Total 38 0.00874113

R-Square Coeff Var Root MSE acidity Mean

0.562870 1.278432 0.010302 0.805862

Source DF Anova SS Mean Square F Value Pr > F

age 2 0.00492012 0.00246006 23.18 <.0001

Page 14

The P-value is less than 0.0001, so we reject our null that all three ages of wine have the same mean acidity. All we can conclude at this point is that the three means are not all the same, but beyond that we cannot go yet.

If you thought the three groups had unequal spreads above, you ought to be a bit cautious about this, but : the P-value is so small that even if our assumptions are slightly off, the conclusions are unlikely to change. (That is, our somewhat unequal spreads would have to change the P-value from less than 0.0001 to greater than 0.05 to affect our conclusions, and I don’t see how our assumptions are “off” by enough to make this happen.)

If you really insist on caution, you should run something like Mood’s Median Test, as per last assignment. We did this before in R. SAS can also do Mood’s median test. It lives in proc npar1way, but not under the name you think: there is a mood there, which is not what we need. We want median, which, according to the documentation, is the “Brown-Mood median test”, the right thing:

proc npar1way median;

var acidity;

class age;

with this output:

The NPAR1WAY Procedure

Median Scores (Number of Points Above Median) for Variable acidity

Classified by Variable age

Sum of Expected Std Dev Mean

age N Scores Under H0 Under H0 Score

15 13 10.0 6.333333 1.490712 0.769231

20 13 8.0 6.333333 1.490712 0.615385

25 13 1.0 6.333333 1.490712 0.076923

Median One-Way Analysis

Chi-Square 13.4000

DF 2

Pr > Chi-Square 0.0012

The P-value of 0.0012 is not quite as small as from the ANOVA, but is easily small enough to reject a null hypothesis that the population medians are all equal.

In the grand scheme of things, our concerns were not enough to invalidate the ANOVA F-test, since the conclusions were the same both ways.

(e) Is it appropriate to run Tukey’s method? Explain briefly why or why not. If appropriate, run Tukey’s method. What do you conclude?

Solution: The F-test was significant, so there are some differences to find. In SAS, to get Tukey, we have to run the whole ANOVA again (so if I had been smart, I would have included the Tukey in the first run of proc anova and ignored it had I not needed it):

proc anova;

Page 15

class age;

model acidity=age;

means age / tukey;

Here’s the extra bit of output:

Page 16

The ANOVA Procedure

Tukey's Studentized Range (HSD) Test for acidity

NOTE: This test controls the Type I experimentwise error rate, but it generally

has a higher Type II error rate than REGWQ.

Alpha 0.05

Error Degrees of Freedom 36

Error Mean Square 0.000106

Critical Value of Studentized Range 3.45675

Minimum Significant Difference 0.0099

Means with the same letter are not significantly different.

Tukey Grouping Mean N age

A 0.814738 13 15

A

A 0.812831 13 20

B 0.790015 13 25

Page 17

The last line of code says “get the mean value of acidity (implied) for each group defined by age, and compare using Tukey.”

The results say that the 15 and 20 means are not significantly different (they are in fact very close), but the age-25 mean is significantly different (less, in fact). (The minimum significant difference, almost exactly 0.01, says that any means differing by this much or more are significantly different. The age 25 group mean differs by more than 0.01 from the others, but the age 15 and 20 groups differ in mean by only 0.002.)

This was exactly the same conclusion that we got (informally) from the boxplots, but now we have P-values to attach to it.

To do something like Tukey with Mood’s median test, what I’ve seen suggested is that you do Mood’s median test on all the pairs of groups (ie. 15 vs. 20, 15 vs. 25, 20 vs. 25). Then you adjust for the fact that you’ve done three tests at once. A simple but “conservative” adjustment is to multiply all those P-values by 3.3

I think I could even do that, with a bit of cleverness:

proc npar1way median;

where age=15 | age=20;

var acidity;

class age;

proc npar1way median;

where age=15 | age=25;

var acidity;

class age;

proc npar1way median;

where age=20 | age=25;

var acidity;

class age;

That has three Mood’s median tests as output:

Page 18

The NPAR1WAY Procedure

Median Scores (Number of Points Above Median) for Variable acidity

Classified by Variable age

Sum of Expected Std Dev Mean

age N Scores Under H0 Under H0 Score

15 13 7.0 6.50 1.30 0.538462

20 13 6.0 6.50 1.30 0.461538

Median Two-Sample Test

Statistic 7.0000

Z 0.3846

One-Sided Pr > Z 0.3503

Two-Sided Pr > |Z| 0.7005

Median One-Way Analysis

Chi-Square 0.1479

DF 1

Pr > Chi-Square 0.7005

The NPAR1WAY Procedure

Median Scores (Number of Points Above Median) for Variable acidity

Classified by Variable age

Sum of Expected Std Dev Mean

age N Scores Under H0 Under H0 Score

15 13 11.0 6.50 1.30 0.846154

25 13 2.0 6.50 1.30 0.153846

Median Two-Sample Test

Statistic 11.0000

Z 3.4615

One-Sided Pr > Z 0.0003

Two-Sided Pr > |Z| 0.0005

Median One-Way Analysis

Chi-Square 11.9822

DF 1

Pr > Chi-Square 0.0005

The NPAR1WAY Procedure

Median Scores (Number of Points Above Median) for Variable acidity

Classified by Variable age

Sum of Expected Std Dev Mean

age N Scores Under H0 Under H0 Score

20 13 11.0 6.50 1.30 0.846154

25 13 2.0 6.50 1.30 0.153846

Page 19

Median Two-Sample Test

Statistic 11.0000

Z 3.4615

One-Sided Pr > Z 0.0003

Two-Sided Pr > |Z| 0.0005

Median One-Way Analysis

Chi-Square 11.9822

DF 1

Pr > Chi-Square 0.0005

Page 20

For the first one, I only take the acidity measurements that have ages 15 or 20 and do a Mood’s median test on them: that is to say, comparing just age 15 and 20. This has a very large P-value of 0.7005 (two-sided, since we are looking for any differences). Ages 15 and 20 are not significantly different.

The second one compares ages 15 and 25, and produces a very small P-value of 0.0005. Hang on to that one a moment.

The third one compares ages 20 and 25, and likewise produces a P-value of 0.0005.

Now we have just done three tests at once, which we have to account for. The Bonferroni procedure says that, with three tests, we multiply each P-value by 3:

Test P-value Adjusted P-value 15 vs. 20 0.7005 “large” 15 vs. 25 0.0005 0.0015 20 vs. 25 0.0005 0.0015

In this case, the adjustment hasn’t changed anything (since the decisions were very clear-cut: age 25 is significantly different in terms of acidity from each of the other ages, which are not significantly different from each other.

Since we are comparing medians, we should look at the medians for each group. proc means can be cajoled into providing these:4

proc means median;

var acidity;

class age;

The MEANS Procedure

Analysis Variable : acidity

N

age Obs Median

-----------------------------------

15 13 0.8133000

20 13 0.8109000

25 13 0.7882000

-----------------------------------

We see, exactly as we did from ANOVA plus Tukey, that the acidity at age 25 is significantly lower than at the other two ages. In other words, the conclusions are identical from the two strategies. They won’t always be, but in this case they are.

4. We are planning a study to estimate a population mean. The population standard deviation is believed to be 20, and the population distribution is believed to be approximately normal. We will be testing the null hypothesis that the population mean is 100. Suppose the population mean is actually 110, and we want to determine how likely we are to (correctly) reject the null hypothesis in this case, using a two-sided (but one-sample) test with α = 0.05.

(a) We will take a sample of size n = 30. What is the power of this test? Use R.

Page 21

Solution: power.t.test. Fill in: sample size n, difference in means delta (10 = 110−100), population SD sd, type of test type (one.sample) and kind of alternative hypothesis alternative (two.sided). Leave out power since that’s what we want:

power.t.test(n=30,delta=10,sd=20,type="one.sample",alternative="two.sided")

##

## One-sample t test power calculation

##

## n = 30

## delta = 10

## sd = 20

## sig.level = 0.05

## power = 0.7539627

## alternative = two.sided

(b) Use SAS to find this power, and verify that you get an answer consistent with the one you just got from R. You will need to specify the situation (onesamplemeans), the test test=t, the true mean (mean), the null mean nullmean, the population SD (stddev), the “total” sample size ntotal and the power that you are trying to find.

Solution: proc power, suitably modified. Guess or search for these:

proc power;

onesamplemeans

test=t

mean=110

nullmean=100

stddev=20

ntotal=30

power=.;

The POWER Procedure

One-Sample t Test for Mean

Fixed Scenario Elements

Distribution Normal

Method Exact

Null Mean 100

Mean 110

Standard Deviation 20

Total Sample Size 30

Number of Sides 2

Alpha 0.05

Computed Power

Power

0.754

Page 22

The power is the same 0.754 as R (to 3 decimals rather than R’s 7).

(c) Use R to find the sample size necessary to obtain a power of at least 0.80 under these conditions. What sample size do you need? Explain briefly how your answer is consistent with (a).

Solution: This time, in power.t.test, put in 0.80 for power and leave out n. The order of things doesn’t matter (since I have named everything that’s going into power.t.test):

power.t.test(delta=10,power=0.80,sd=20,type="one.sample",alternative="two.sided")

##

## One-sample t test power calculation

##

## n = 33.3672

## delta = 10

## sd = 20

## sig.level = 0.05

## power = 0.8

## alternative = two.sided

To get sample size for power at least 0.80, we have to round 33.36 up to the next whole number, ie. n = 34 is needed. (A sample of size 33 wouldn’t quite have enough power.)

This answer is consistent with (a) because a sample size of 30 gave a power a bit less than 0.80, and so to increase the power by a little (0.75 to 0.80), we had to increase the sample size by a little (30 to 34).

(d) Use SAS to find the sample size necessary to obtain a power of (i) 0.80 and (ii) 0.90. (This is doable in one step.) What sample sizes do you get?

Solution: Use the same SAS code as (b), and modify it to specify the power values (with a space between) and leave the sample size blank:

proc power;

onesamplemeans

test=t

mean=110

nullmean=100

stddev=20

ntotal=.

power=0.80 0.90;

The POWER Procedure

One-Sample t Test for Mean

Fixed Scenario Elements

Distribution Normal

Method Exact

Null Mean 100

Mean 110

Standard Deviation 20

Number of Sides 2

Alpha 0.05

Page 23

Computed N Total

Nominal Actual N

Index Power Power Total

1 0.8 0.808 34

2 0.9 0.900 44

Sample sizes of 34 (for power 0.80) and 44 (for power 0.90).

The first of these is what we got from R.5

5. I’ve mentioned several times that the sign test has less power than the t-test. Let’s investigate this with a specific example.

Let’s suppose we are testing H0 : µ = 40 against Ha : µ 6= 40, where µ is the population mean (and median, as we shall see). Our population actually has a normal distribution with mean 50 and SD 15, so that the null hypothesis is wrong and we want to reject it most of the time. On the other hand, the population actually is normally-distributed and so the t-test is the right one to use.

(a) Use proc power in SAS to find the probability that a t-test correctly rejects the null hypothesis using a sample size of n = 10. You’ll have to specify the right things to set up a one-sample t-test, and give values (possibly missing) for nullmean, mean, stddev, ntotal and power.

Solution: This is a one-sample t-test, so the appropriate code for proc power is this. We need to set power equal to missing, since that’s what we want to find, and supply values for all the others:

proc power;

onesamplemeans

test=t

nullmean=40

mean=50

stddev=15

ntotal=10

power=.;

with output

The POWER Procedure

One-Sample t Test for Mean

Fixed Scenario Elements

Distribution Normal

Method Exact

Null Mean 40

Mean 50

Standard Deviation 15

Total Sample Size 10

Number of Sides 2

Alpha 0.05

Page 24

Computed Power

Power

0.469

The power is 0.469. Not great, but we’ll see how this stacks up against the sign test.

(b) What code in R would draw a random sample of size 10 from the true population distribution and save the sample in a variable?

Solution: The data actually have a normal distribution with mean 50 and SD 15, so we use rnorm with this mean and SD, obtaining 10 values:

x=rnorm(10,50,15)

x

## [1] 44.35918 63.08701 52.72680 47.05784 55.42996 32.01897 48.99054

## [8] 67.66784 46.84664 66.15004

(c) What code would count how many of the sampled values are less than 40 and how many are greater (or equal)? Save the result in a variable again.

Solution: Just this:

tab=table(x<40)

tab

##

## FALSE TRUE

## 9 1

(d) It turns out the sign test would reject H0 : M = 40 against Ha : M 6= 40 at α = 0.05 if the smaller of the numbers in the last part is 1 or less. (M is the population median.) Write a line of code (or possibly two lines) that gives you TRUE if you should reject the null for your data and FALSE otherwise.

Solution: This is actually easier than you might think, because testing whether the smaller table is 1 or less gives you a true or a false:

is.reject=(min(tab)<=1)

is.reject

## [1] TRUE

Or you can obtain the smallest value first, and then test whether it is less or equal to 1. That would be two lines.

You might be wondering where the “1 or less” came from. Getting a P-value for the sign test involves the binomial distribution: if the null is correct, each data value is independently either above or below 40, with probability 0.5 of each, so the number of values below 40 (say) is binomial with n = 10 and p = 0.5. The P-value for 1 observed value below 40 and the rest above is

Page 25

2*pbinom(1,10,0.5)

## [1] 0.02148438

which is less than 0.05; the P-value for 2 values below 40 and the rest above is

2*pbinom(2,10,0.5)

## [1] 0.109375

which is bigger than 0.05.

You might have encountered the term “critical region” for a test. This is the values of your test statistic that you would reject the null hypothesis for. In this case, the critical region is 1 and 0 observations below 40, along with 1 and 0 observations above 40.

When you’re thinking about power, I think it’s easiest to think in terms of the critical region (rather than directly in terms of P-values) since you have a certain α in mind all the way through, 0.05 in the power examples that I’ve done. The steps are then:

• Work out the critical region for your test, the values of the test statistic (or sample mean or sample count) that would lead to rejecting the null hypothesis.

• Under your particular alternative hypothesis, find the probability of falling into your critical region.

When I say “work out”, I mean either calculating (along the lines of STAB57), or simulating, as we have done here.

(e) Put your code from parts (b), (c) and (d) into a function that accepts no input and returns a True or a False according to whether the simulated data led to a rejection of H0 : M = 40 by the sign test or not. (A function with no input has a () after its name.)

Solution: Just copy and paste the code above to form the body of the function. The thing I obtained as is.reject is the thing that needs to be returned.

random.sign.test=function() { x=rnorm(10,50,15)

tab=table(x<40)

is.reject=(min(tab)<=1)

return(is.reject)

}

I am expecting this to reject sometimes, so I’d better try it a few times:

replicate(12,random.sign.test())

## [1] FALSE FALSE TRUE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE

## [12] FALSE

Not exactly exhaustive testing, but it seems to work. (My confidence mainly comes from having built the function one line at a time in (b)-(d) above.)

(f) Run your function 1000 times using replicate, and count how many times it rejects the null median of 40. What is the power of the sign test under these circumstances?

Page 26

Solution: I’m using set.seed to make my results here reproducible (so that I can talk about them), but you don’t need to. The function has to have the () brackets on it, to show that you are calling it with no input:

set.seed(457299)

r=replicate(1000,random.sign.test())

table(r)

## r

## FALSE TRUE

## 799 201

The (simulated) power of the sign test is the number of times it correctly rejected the null median of 40, which is the number of TRUE values in the table. In my case, that is 201 out of 1000 or 0.201. In your case, it will be different (but it should be somewhere close to that).

(g) Which is more powerful in this case, the sign test or the t-test? How do you know?

Solution: The power of the sign test is estimated as 0.201, which is quite a bit less than the power of the t-test, which we found back in (a) to be 0.469. So the t-test, in this situation where it is valid, is the right test to use: it is (a) valid and (b) more powerful.

So the t-test is more powerful. One way to think about how much more powerful is to ask “how much smaller of a sample size would be needed for the t-test to have the same power as this sign test?” The power of my sign test was 0.201, so in SAS’s proc ttest we set power equal to that and leave the sample size ntotal blank:

proc power;

onesamplemeans

test=t

nullmean=40

mean=50

stddev=15

ntotal=.

power=0.201;

which has output

The POWER Procedure

One-Sample t Test for Mean

Fixed Scenario Elements

Distribution Normal

Method Exact

Null Mean 40

Mean 50

Standard Deviation 15

Nominal Power 0.201

Number of Sides 2

Alpha 0.05

Page 27

Computed N Total

Actual N

Power Total

0.211 5

Page 28

A sample of size 5 gives the same power for the t-test that a sample of size 10 does for the sign test. The ratio of these two sample sizes is called the relative efficiency of the two tests: in this case, the t-test is 10/5 = 2 times more efficient. The data that you have are being used “more efficiently” by the t-test.

It is possible to derive6 the limiting relative efficiency of the t test relative to the sign test when the data are actually normal, as the sample size gets larger. This turns out not to depend on how far wrong the null is (as long as it is the same for both the t-test and the sign test). This “asymptotic relative efficiency” is π/2 = 1.57. Our relative efficiency for power 0.201, namely 2, wasn’t especially close to this, but then our sample sizes 10 and 5 are not especially close to infinity either.

This says that, if your data are actually from a normal distribution, you do a lot better to use the t-test than the sign test, because the sign test is wasteful of data (it only uses above/below rather than the actual values).

If your data are not from a normal distribution, then the story can be very different.

Of course you knew I would investigate this. There is a distribution called the “Laplace” or “double exponential” distribution, that has very long tails.7 The distribution is not in base R, but there is a package called smoothmest that contains a function rdoublex to generate random values from this distribution. So we’re going to do a simulation investigation of the power of the sign test for Laplace data, by the same simulation technique that we did above. Like the normal, the Laplace distribution is symmetric, so its mean and median are the same (which makes our life easier).8

Let’s test the hypothesis that the median is zero. We’ll suppose that the true median is 0.5 (this is called mu in rdoublex). The first problem we run into is that we can’t use proc power or power.t.test because they assume normal data, which we are far from having. So we have to do two simulations: one to simulate the power of the t test, and one to simulate the power of the sign test.

To simulate the t test, we first have to generate some Laplace data with the true mean of 0.5. We’ll use a sample size of 50 throughout these simulations.

library(smoothmest)

## Loading required package: MASS

##

## Attaching package: ’MASS’

## The following object is masked from ’package:dplyr’:

##

## select

rl=rdoublex(50,mu=0.5)

rl

## [1] -0.33323285 0.70569291 -1.22513053 0.68517708 0.12778518

## [6] 0.50749949 0.26700527 1.90236874 0.53288312 -0.37374732

## [11] 0.27256566 0.53365929 0.43581431 -0.01545866 0.18594908

## [16] -0.40403202 1.13540289 0.16137306 -0.23360644 -0.74050354

## [21] 2.92089551 -2.72173880 0.48428815 1.23636045 0.17078618

## [26] 1.72456334 0.07903058 0.25210411 0.09512810 2.52310082

## [31] -2.13629814 0.81851434 0.74615575 -0.26068744 2.70683355

## [36] 1.46981530 1.45646489 -0.20232517 6.65249860 1.51575026

## [41] -0.07606399 -1.11338640 -1.20427995 -0.70986104 -1.66466321

## [46] 0.55346854 0.66091469 0.72100677 0.92025176 0.98922656

This seems to have some unusual values, far away from zero:

Page 29

qqnorm(rl)

qqline(rl)

●

●

●

●

●

● ●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●●

●

●

●●

●

●

●

●

● ●

●

●

● ● ●

●●

−2 −1 0 1 2

− 2

0 2

4 6

Normal Q−Q Plot

Theoretical Quantiles

S a

m p

le Q

u a

n til

e s

You see the long tails compared to the normal.

Now, we feed these values into t.test and see whether we reject a null median of zero (at α = 0.05):

Page 30

tt=t.test(rl)

tt

##

## One Sample t-test

##

## data: rl

## t = 2.2556, df = 49, p-value = 0.0286

## alternative hypothesis: true mean is not equal to 0

## 95 percent confidence interval:

## 0.04959344 0.85981931

## sample estimates:

## mean of x

## 0.4547064

Or we can just pull out the P-value and even compare it to 0.05:

pval=tt$p.value

pval

## [1] 0.02859596

is.reject=(pval<=0.05)

is.reject

## [1] TRUE

This one has a small P-value and so the null median of 0 should be (correctly) rejected.

Let’s bear in mind what we’re going to do: we’re going to run this many times, and collect up the results. So we assemble this into a function. I’m going to let the mu parameter (the true median) be an input to the function, since I might want to change that. Same for the sample size, which I’ll call n:

laplace.t.test=function(mu.true,n) { rl=rdoublex(n,mu=mu.true)

tt=t.test(rl)

pval=tt$p.value

is.reject=(pval<0.05)

return(is.reject)

}

The first time I did this, I wrote the first line as function laplace.t.test(mu.true), which might have worked in any number of other languages, but not in R!9

I might also have let the α value be input to the function, since I might want to run simulations with that different. But we’ll leave it like this for now.

Run it a few times for testing:

replicate(10,laplace.t.test(0.5,50))

## [1] FALSE TRUE FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE

That looks reasonable. And then 1000 times for the actual simulation:

Page 31

t.sim=replicate(1000,laplace.t.test(0.5,50))

table(t.sim)

## t.sim

## FALSE TRUE

## 303 697

The t-test rejects 697 times out of 1000, so the power is 0.697.

Now we need to think about the sign test. I’ll write a slightly different version of the sign test function, one that generates the random sample and returns TRUE if we reject and FALSE otherwise. The lines to tabulate the number of values above and below the null mean (0) and to work out the binomial P-value should be familiar by now:

laplace.sign.test=function(mu.true,n) { rl=rdoublex(n,mu=mu.true)

tab=table(rl<0)

stat=min(tab)

pval=2*pbinom(stat,n,0.5)

is.reject=(pval<0.05)

return(is.reject)

}

Does that work? Try it a few times:

replicate(20,laplace.sign.test(0.5,50))

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

## [12] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE

and then do the simulation:

rr=replicate(1000,laplace.sign.test(0.5,50))

table(rr)

## rr

## FALSE TRUE

## 238 762

The power for the sign test is 0.762. And here’s the kicker: for Laplace-distributed data, the sign test is more powerful than the t-test.

This is not to say that you will ever run into data that comes from the Laplace distribution. But the moral of the story is that the sign test can be more powerful than the t-test, under the right circumstances (and the above simulation is the “proof” of that statement). So a blanket statement like “the sign test is not very powerful” needs to be qualified a bit: when your data come from a sufficiently long-tailed distribution, the sign test can be more powerful relative to the t-test than you would think.

6. (Hand this one in.) A biology graduate student exposed each of 32 tomato plants to one of four different colours of light (8 plants to each colour). The growth rate of each plant, in millimetres per week, was recorded. The data are in http://www.utsc.utoronto.ca/~butler/c32/tomatoes.txt.

(a) (2 marks) Read the data into R and confirm that you have 8 rows and 5 columns of data.

(b) (3 marks) Use something from tidyr to re-arrange the data so that you have one column containing all the growth rates, and another column saying which colour light each plant was exposed to. (The aim here is to produce something suitable for feeding into aov later.) One of the earlier problems

Page 32

in this assignment contains an introduction to tidyr.

(c) (2 marks) Save the data in the new format to a text file. This is most easily done using write.table, which is the opposite of read.table. It requires two things: a data frame, and the name of a file to save in. A third thing, quote=F, stops R from adding extra quotation marks to text (we don’t need that for what we’re doing).

(d) (4 marks) Get the data file into a SAS dataset. This will mean something like (i) opening the file in R Studio, (ii) selecting all the lines in it, (iii) pasting it into SAS Studio, (iv) making sure it has the right filename, (v) writing code to read in the file. (You do not need to show that the file was read properly, but if it wasn’t, you’ll get the rest of the question wrong, so it’s up to you to check that it was read properly.)

(e) (3 marks) Find the mean growth rate for each colour. Would you expect an analysis of variance to come out significant? Explain briefly.

(f) (3 marks) Run an analysis of variance to see whether there are any differences in mean growth rate among the different colours. What do you conclude?

(g) (3 marks) Explain briefly whether or not you need to obtain the output for Tukey’s method. If you do need it, obtain it and explain briefly what it tells you.

7. Hand this one in. On Assignment 4, you carried out a sign test to determine whether students could throw a baseball farther than a softball. This time, we will calculate a confidence interval for the median difference baseball minus softball, using the results of sign tests.

(a) (1 mark) Read the data into R from http://www.utsc.utoronto.ca/~butler/c32/throw.txt, giving appropriate names to the columns, and calculating a column of differences. (This is exactly the same as you did in Assignment 4, so only one mark for it this time.)

(b) (3 marks) Make a function to carry out a two-sided sign test, based on a null median called med and data in a variable called z. (There are lots of models scattered around the course that you can borrow from, but you will learn the most by adapting your own code from Assignment 4.) Make the null median be the first thing input to your function (this will make your life easier below).

(c) (2 marks) Test your function on the data of this question, and check that you get two times the P-value that you got from R on Assignment 4.

(d) (2 marks) Based on your P-value in assignment 4 (look at my solutions if you are unsure), do you think 0 is inside the confidence interval or not? Explain briefly.

(e) (4 marks) Obtain a 95% confidence interval for the population median difference, baseball minus softball, using the procedure shown in class, with sapply and refinement as necessary.

Notes

1With three points, the regression line won’t necessarily go through all of them. In fact, it only will if they happened to be on a line in the first place, which will be the regression line. With two points, though, the regression line is the line through the two points.

2Sometimes the column playing the role of “rep” is interesting to us, but not here.

3“Conservative” in statistics means that it may fail to declare some actually significant results significant: “playing it safe”, you might say, thinking of rejecting the null as something you only want to do if you’re sure it’s the right thing. That is to say, the correct P-value is smaller than the one given, but we don’t know how much smaller. This particular adjustment is called a Bonferroni correction, and applies any time you are doing several tests at once. There are better but more complicated adjustments; indeed, one of the reasons Tukey developed his procedure is that he wanted to avoid using something like Bonferroni in an ANOVA context where it is possible to do better.

Page 33

4“proc means” knows about a number of statistics: see http://support.sas.com/documentation/cdl/en/proc/61895/HTML/ default/viewer.htm#a000146729.htm. Look for “statistic-keywords” near the bottom.

5Rounded up: SAS has done the proper rounding for you.

6Meaning, I forget how to do it. But it has something to do with looking at alternatives that are very close to the null.

7If you’ve ever run into the exponential distribution, you’ll recall that this is right skewed with a very long tail. The Laplace distribution looks like two of these glued back to back.

8This is about the only way in which the normal and Laplace distributions are alike.

9I seem to do this one a lot. And then I gaze at it and wonder why it’s not working.

10There are also boxplots in with the proc anova output.

11If you go too far below zero, you’ll get some P-values of 2 (!). For example, if the null median is −10, all the data values are bigger than the null median, so you get a table that says just TRUE with a frequency of 24, and no FALSE. Then small gets to be 24 and p gets to be 1, because I didn’t account for the table only having one thing in it.

Page 34