Exploring assumptions

FIGURE 5.1
I came first in
the competition
for who has the
smallest brain

5.1. What will this chapter tell me? 1

When we were learning to read at primary school, we used to read versions of stories by
the famous storyteller Hans Christian Andersen. One of my favourites was the story of
the ugly duckling. This duckling was a big ugly grey bird, so ugly that even a dog would
not bite him. The poor duckling was ridiculed, ostracized and pecked by the other ducks.
Eventually, it became too much for him and he flew to the swans, the royal birds, hoping
that they would end his misery by killing him because he was so ugly. As he stared into the

166
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 167

water, though, he saw not an ugly grey bird but a beautiful swan. Data are much the same.
Sometimes they’re just big, grey and ugly and don’t do any of the things that they’re sup-
posed to do. When we get data like these, we swear at them, curse them, peck them and
hope that they’ll fly away and be killed by the swans. Alternatively, we can try to force our
data into becoming beautiful swans. That’s what this chapter is all about: assessing how
much of an ugly duckling of a data set you have, and discovering how to turn it into a swan.
Remember, though, a swan can break your arm.1

5.2. What are assumptions? 1

Some academics tend to regard assumptions as rather tedious things about which
no one really need worry. When I mention statistical assumptions to my fellow Why bother with
assumptions?
psychologists they tend to give me that raised eyebrow, ‘good grief, get a life’
look and then ignore me. However, there are good reasons for taking assump-
tions seriously. Imagine that I go over to a friend’s house, the lights are on and
it’s obvious that someone is at home. I ring the doorbell and no one answers.
From that experience, I conclude that my friend hates me and that I am a ter-
rible, unlovable person. How tenable is this conclusion? Well, there is a reality
that I am trying to tap (i.e., whether my friend likes or hates me), and I have
collected data about that reality (I’ve gone to his house, seen that he’s at home,
rung the doorbell and got no response). Imagine that in reality my friend likes me (he’s a
lousy judge of character); in this scenario, my conclusion is false. Why have my data led me
to the wrong conclusion? The answer is simple: I had assumed that my friend’s doorbell
was working and under this assumption the conclusion that I made from my data was accu-
rate (my friend heard the bell but chose to ignore it because he hates me). However, this
assumption was not true – his doorbell was not working, which is why he didn’t answer
the door – and as a consequence the conclusion I drew about reality was completely false.
It pays to check assumptions and your doorbell batteries.
Enough about doorbells, friends and my social life: the point to remember is that when
assumptions are broken we stop being able to draw accurate conclusions about reality.
Different statistical models assume different things, and if these models are going to reflect
reality accurately then these assumptions need to be true. This chapter is going to deal with
some particularly ubiquitous assumptions so that you know how to slay these particular
beasts as we battle our way through the rest of the book. However, be warned: some tests
have their own unique two-headed, fire-breathing, green-scaled assumptions and these will
jump out from behind a mound of blood-soaked moss and try to eat us alive when we least
expect them to. Onward into battle …

5.3. Assumptions of parametric data 1

What are the

Many of the statistical procedures described in this book are paramet- assumptions of parametric
ric tests based on the normal distribution (which is described in section data?
1.7.4). A parametric test is one that requires data from one of the large
catalogue of distributions that statisticians have described, and for data to
be parametric certain assumptions must be true. If you use a parametric
test when your data are not parametric then the results are likely to be
inaccurate. Therefore, it is very important that you check the assump-
tions before deciding which statistical test is appropriate. Throughout
1
Although it is theoretically possible, apparently you’d have to be weak boned, and swans are nice and wouldn’t
do that sort of thing.
168 D I S C O VE R I N G STAT I ST I C S US I N G R

this book you will become aware of my obsession with assumptions and checking them.
Most parametric tests based on the normal distribution have four basic assumptions that
must be met for the test to be accurate. Many students find checking assumptions a pretty
tedious affair, and often get confused about how to tell whether or not an assumption has
been met. Therefore, this chapter is designed to take you on a step-by-step tour of the
world of parametric assumptions. Now, you may think that assumptions are not very excit-
ing, but they can have great benefits: for one thing, you can impress your supervisor/
lecturer by spotting all of the test assumptions that they have violated throughout their
careers. You can then rubbish, on statistical grounds, the theories they have spent their
lifetime developing – and they can’t argue with you,2 but they can poke your eyes out. The
assumptions of parametric tests are:

1 Normally distributed data: This is a tricky and misunderstood assumption because it

means different things in different contexts. For this reason I will spend most of the
chapter discussing this assumption. In short, the rationale behind hypothesis test-
ing relies on having something that is normally distributed (in some cases it’s the
sampling distribution, in others the errors in the model), and so if this assumption
is not met then the logic behind hypothesis testing is flawed (we came across these
principles in Chapters 1 and 2).
2 Homogeneity of variance: This assumption means that the variances should be the
same throughout the data. In designs in which you test several groups of participants
this assumption means that each of these samples comes from populations with the
same variance. In correlational designs, this assumption means that the variance of
one variable should be stable at all levels of the other variable (see section 5.7).
3 Interval data: Data should be measured at least at the interval level. This assumption
is tested by common sense and so won’t be discussed further (but do read section
[Link] again to remind yourself of what we mean by interval data).
4 Independence: This assumption, like that of normality, is different depending on the
test you’re using. In some cases it means that data from different participants are inde-
pendent, which means that the behaviour of one participant does not influence the
behaviour of another. In repeated-measures designs (in which participants are mea-
sured in more than one experimental condition), we expect scores in the experimental
conditions to be non-independent for a given participant, but behaviour between dif-
ferent participants should be independent. As an example, imagine two people, Paul
and Julie, were participants in an experiment where they had to indicate whether they
remembered having seen particular photos earlier on in the experiment. If Paul and
Julie were to confer about whether they’d seen certain pictures then their answers
would not be independent: Julie’s response to a given question would depend on Paul’s
answer, and this would violate the assumption of independence. If Paul and Julie were
unable to confer (if they were locked in different rooms) then their responses should be
independent (unless they’re telepathic): Julie’s should not influence Paul’s responses.
In regression, however, this assumption also relates to the errors in the regression
model being uncorrelated, but we’ll discuss that more in Chapter 7.

We will, therefore, focus in this chapter on the assumptions of normality and homogeneity
of variance.

2
When I was doing my Ph.D., we were set a task by our statistics lecturer in which we had to find some published
papers and criticize the statistical methods in them. I chose one of my supervisor’s papers and proceeded to slag
off every aspect of the data analysis (and I was being very pedantic about it all). Imagine my horror when my
supervisor came bounding down the corridor with a big grin on his face and declared that, unbeknownst to me,
he was the second marker of my essay. Luckily, he had a sense of humour and I got a good mark.-
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 169

5.4. Packages used in this chapter 1

Some useful packages for exploring data are car, ggplot2 (for graphs), pastecs (for descrip-
tive statistics) and psych. Of course, if you plan to use R Commander then you need the
Rcmdr package installed too (see section 3.6). If you do not have these packages installed,
you can install them by executing the following commands:
[Link]("car"); [Link]("ggplot2");

[Link]("pastecs"); [Link]("psych")

You then need to load these packages by executing the commands:

library(car); library(ggplot2); library(pastecs); library(psych);
library(Rcmdr)

5.5. The assumption of normality 1

We encountered the normal distribution back in Chapter 1, we know what it looks like
and we (hopefully) understand it. You’d think then that this assumption would be easy to
understand – it just means that our data are normally distributed, right? Actually, no. In
many statistical tests (e.g., the t-test) we assume that the sampling distribution is normally
distributed. This is a problem because we don’t have access to this distribution – we can’t
simply look at its shape and see whether it is normally distributed. However, we know
from the central limit theorem (section 2.5.1) that if the sample data are approximately
normal then the sampling distribution will be also. Therefore, people tend to look at their
sample data to see if they are normally distributed. If so, then they have a little party to
celebrate and assume that the sampling distribution (which is what actually matters) is also.
We also know from the central limit theorem that in big samples the sampling distribu-
tion tends to be normal anyway – regardless of the shape of the data we actually collected
(and remember that the sampling distribution will tend to be normal regardless of the
population distribution in samples of 30 or more). As our sample gets bigger, then, we
can be more confident that the sampling distribution is normally distributed (but see Jane
Superbrain Box 5.1).
The assumption of normality is also important in research using regression (or general
linear models). General linear models, as we will see in Chapter 7, assume that errors in the
model (basically, the deviations we encountered in section 2.4.2) are normally distributed.
In both cases it might be useful to test for normality, and that’s what this section is
dedicated to explaining. Essentially, we can look for normality visually, look at values that
quantify aspects of a distribution (i.e., skew and kurtosis) and compare the distribution we
have to a normal distribution to see if it is different.

5.5.1. Oh no, it’s that pesky frequency distribution again:

checking normality visually 1
We discovered in section 1.7.1 that frequency distributions are a useful way to look at
the shape of a distribution. In addition, we discovered how to plot these graphs in sec-
tion 4.4.8. Therefore, we are already equipped to look for normality in our sample using
a graph. Let’s return to the Download Festival data from Chapter 4. Remember that a
170 D I S C O VE R I N G STAT I ST I C S US I N G R

biologist had visited the Download Festival (a rock and heavy metal festival in the UK) and
assessed people’s hygiene over the three days of the festival using a standardized technique
that results in a score ranging between 0 (you smell like a rotting corpse that’s hiding up a
skunk’s anus) and 4 (you smell of sweet roses on a fresh spring day). The data file can be
downloaded from the companion website ([Link]) – remember to use the
version of the data for which the outlier has been corrected (if you haven’t a clue what I
mean, then read section 4.4.8 or your graphs will look very different from mine!).

SELF-TEST

9 Using what you learnt in Chapter 4, plot histograms

for the hygiene scores for the three days of the
Download Festival. (For reasons that will become
apparent, use geom_histogram(aes(y = ..density..)
rather than geom_histogram().)

When you drew the histograms, this gave you the distributions. It might be nice to also
have a plot of what a normal distribution looks like, for comparison purposes. Even better
would be if that we could put a normal distribution onto the same plot. Well, we can using
the power of ggplot2. First, load in the data:
dlf <- [Link]("[Link]", header=TRUE)

To draw the histogram, you should have used code something like:
hist.day1 <- ggplot(dlf, aes(day1)) + opts([Link] = "none") +
geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
labs(x = "Hygiene score on day 1", y = "Density")

hist.day1

To see what this function is doing we can break down the command:

G ggplot(dlf, aes(day1)): This tells R to plot the day1 variable from the dlf dataframe.
G opts([Link] = “none”): This command gets rid of the legend of the graph.
G geom_histogram(aes(y=..density..), colour = “black”, fill=”white”): This command
plots the histogram, sets the line colour to be black and the fill colour to be white.
Notice that we have asked for a density plot rather than frequency because we want
to plot the normal curve.
G labs(x = “Hygiene score on day 1”, y = “Density”): this command sets the labels for
the x- and y-axes.

We can add another layer to the chart, which is a normal curve. We need to tell ggplot2
what mean and standard deviation we’d like on that curve though. And what we’d like is
the same mean and standard deviation that we have in our data. To add the normal curve,
we take the existing histogram object (hist.day1) and add a new layer that uses stat_func-
tion() to produce a normal curve and lay it on top of the histogram:
hist.day1 + stat_function(fun = dnorm, args = list(mean = mean(dlf$day1,
[Link] = TRUE), sd = sd(dlf$day1, [Link] = TRUE)), colour = "black", size = 1)

The stat_function() command draws the normal curve using the function dnorm(). This
function basically returns the probability (i.e., the density) for a given value from a normal
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 171

distribution of known mean and standard deviation. The rest of the command specifies the
mean as being the mean of the day1 variable after removing any missing values (mean =
mean (dlf$day1, [Link] = TRUE)), and the standard deviation as being that of day1 (), sd =
sd(dlf$day1, [Link] = TRUE)). We also set the line colour as black and the line width as 1.3

SELF-TEST

9 Add normal curves to the histograms that you drew

for day2 and day3.

There is another useful graph that we can inspect to see if a distribution is normal called a
Q-Q plot (quantile–quantile plot; a quantile is the proportion of cases we find below a certain
value). This graph plots the cumulative values we have in our data against the cumulative
probability of a particular distribution (in this case we would specify a normal distribu-
tion). What this means is that the data are ranked and sorted. Each value is compared to
the expected value that the score should have in a normal distribution and they are plotted
against one another. If the data are normally distributed then the actual scores will have the
same distribution as the score we expect from a normal distribution, and you’ll get a lovely
straight diagonal line. If values fall on the diagonal of the plot then the variable is normally
distributed, but deviations from the diagonal show deviations from normality.
To draw a Q-Q plot using the ggplot2 package, we can use the qplot() function in con-
junction with the qq statistic. Execute the following code:
qqplot.day1 <- qplot(sample = dlf$day1, stat="qq")

qqplot.day1

(Note that by default ggplot2 assumes you want to compare your distribution with a nor-
mal distribution – you can change that if you want to, but it’s so rare that we’re not going
to worry about it here.)

SELF-TEST

9 Create Q-Q plots for the variables day2 and day3.

Figure 5.2 shows the histograms (from the self-test task) and the corresponding Q-Q
plots. The first thing to note is that the data from day 1 look a lot more healthy since we’ve
removed the data point that was mistyped back in section 4.7. In fact the distribution is
amazingly normal looking: it is nicely symmetrical and doesn’t seem too pointy or flat –
these are good things! This is echoed by the Q-Q plot: note that the data points all fall very
close to the ‘ideal’ diagonal line.
3
I have built up the histogram and normal plot in two stages because I think it makes it easier to understand what
you’re doing, but you could build the plot in a single command:

hist.day1 <- ggplot(dlf, aes(day1)) + opts([Link] = "none") + geom_

histogram(aes(y = ..density..), colour = "black", fill = "white") + labs(x =
"Hygiene score on day 1", y = "Density") + stat_function(fun = dnorm, args =
list(mean = mean(dlf$day1, [Link] = TRUE), sd = sd(dlf$day1, [Link] = TRUE)), colour
= "black", size = 1)

hist.day1
172 D I S C O VE R I N G STAT I ST I C S US I N G R

FIGURE 5.2
Histograms (left) 0.7 3.5
and Q-Q plots 0.6 3.0
(right) of the
hygiene scores 0.5 2.5
Density

sample
over the three 0.4 2.0
days of the
Download Festival 0.3 1.5

0.2 1.0

0.1 0.5

0.0
0 1 2 80 −3 −2 −1 0 1 2 3
Hygiene score on day 1 theoretical

0.8
3.0

2.5
0.6
2.0
Density

sample
0.4 1.5

1.0
0.2
0.5

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 −3 −2 −1 0 1 2 3
Hygiene score on day 2 theoretical

0.8 3.5

3.0
0.6 2.5
Density

sample

2.0
0.4
1.5

1.0
0.2
0.5

0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 −3 −2 −1 0 1 2 3
Hygiene score on day 3 theoretical

However, the distributions for days 2 and 3 are not nearly as symmetrical. In fact, they
both look positively skewed. Again, this can be seen in the Q-Q plots by the data val-
ues deviating away from the diagonal. In general, what this seems to suggest is that by
days 2 and 3, hygiene scores were much more clustered around the low end of the scale.
Remember that the lower the score, the less hygienic the person is, so this suggests that
generally people became smellier as the festival progressed. The skew occurs because a
substantial minority insisted on upholding their levels of hygiene (against all odds!) over
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 173

the course of the festival (I find baby wet-wipes are indispensable). However, these skewed
distributions might cause us a problem if we want to use parametric tests. In the next sec-
tion we’ll look at ways to try to quantify the skew and kurtosis of these distributions.

5.5.2. Quantifying normality with numbers 1

It is all very well to look at histograms, but they are subjective and open to abuse (I can
imagine researchers sitting looking at a completely distorted distribution and saying ‘yep,
well Bob, that looks normal to me’, and Bob replying ‘yep, sure does’). Therefore, having
inspected the distribution of hygiene scores visually, we can move on to look at ways to quan-
tify the shape of the distributions and to look for outliers. To further explore the distribution
of the variables, we can use the describe() function, in the psych package.
describe(dlf$day1)

We can also use the [Link]() function of the pastecs package,4 which takes the general
form:
[Link](variable name, basic = TRUE, norm = FALSE)

In this function, we simply name our variable and by default (i.e., if we simply name a vari-
able and don’t include the other commands) we’ll get a whole host of statistics including
some basic ones such as the number of cases (because basic = TRUE by default) but not
including statistics relating to the normal distribution (because norm = FALSE by default).
To my mind the basic statistics are not very useful so I usually specify basic = FALSE (to
get rid of these), but in the current context it is useful to override the default and specify
norm = TRUE so that we get statistics relating to the distribution of scores. Therefore, we
could execute:
[Link](dlf$day1, basic = FALSE, norm = TRUE)

Note that we have specified the variable day1 in the dlf dataframe, asked not to see the
basic statistics (basic = FALSE) but asked to see the normality statistics (norm = TRUE).
We can also use describe() and [Link]() with more than one variable at the same time,
using the cbind() function to combine two or more variables (see R’s Souls’ Tip 3.5).
describe(cbind(dlf$day1, dlf$day2, dlf$day3))
[Link](cbind(dlf$day1, dlf$day2, dlf$day3), basic = FALSE, norm = TRUE)

Note that in each case we have simply replaced a single variable with cbind(dlf$day1,
dlf$day2, dlf$day3) which combines the three variables day1, day2, and day3 into a single
object.
A second way to describe more than one variable is to select the variable names directly
from the data set (see section 3.9.1):
describe(dlf[,c("day1", "day2", "day3")])

[Link](dlf[, c("day1", "day2", "day3")], basic = FALSE, norm = TRUE)

4
There’s always a second way to do something with R. And often a third, fourth and fifth way. While writing this
book Jeremy and I would often look at each other’s bits (and sometimes what we’d written too) and then send an
email saying ‘oh, I didn’t know you could do that, I always use a different function in a different package’. People
can become quite attached to their ‘favourite’ way of doing things in R, but obviously we’re way too cool to have
favourite ways of doing stats, which is why I didn’t at all insist on adding reams of stuff on [Link]() because I
prefer it to Jeremy’s crappy old describe() function.
174 D I S C O VE R I N G STAT I ST I C S US I N G R

Remember that we can select rows and columns using [rows, columns], therefore, dlf[,
c(“day1”, “day2”, “day3”)] means from the dlf dataframe select all of the rows (because
nothing is specified before the comma) and select the columns labelled day1, day2, and
day3 (because we have specified c(“day1”, “day2”, “day3”)).

R ’ s S o ul s ’ T i p 5 . 1 Funny numbers 1

You might notice that R sometimes reports numbers with the letter ‘e’ placed in the mix just to confuse you. For
example, you might see a value such as 9.612 e−02 and many students find this notation confusing. Well, this
notation means 9.612 × 10−2 (which might be a more familiar notation, or could be even more confusing). OK,
some of you are still confused. Well think of e−02 as meaning ‘move the decimal place 2 places to the left’, so
9.612 e−02 becomes 0.09612. If the notation read 9.612 e−01, then that would be 0.9612, and if it read 9.612
e−03, that would be 0.009612. Likewise, think of e+02 (notice the minus sign has changed) as meaning ‘move
the decimal place 2 places to the right’. So 9.612 e+02 becomes 961.2.

The results of these commands are shown in Output 5.1 (describe) and Output 5.2
([Link]). These outputs basically contain the same values5 although they are presented
in a different notation in Output 5.2 (see R’s Souls’ Tip 5.1). We can see that, on average,
hygiene scores were 1.77 (out of 4) on day 1 of the festival, but went down to 0.96 and
0.98 on days 2 and 3, respectively. The other important measures for our purposes are the
skew and the kurtosis (see section 1.7.1). The values of skew and kurtosis should be zero
in a normal distribution. Positive values of skew indicate a pile-up of scores on the left of
the distribution, whereas negative values indicate a pile-up on the right. Positive values of
kurtosis indicate a pointy and heavy-tailed distribution, whereas negative values indicate
a flat and light-tailed distribution. The further the value is from zero, the more likely it is
that the data are not normally distributed. For day 1 the skew value is very close to zero
(which is good) and kurtosis is a little negative. For days 2 and 3, though, there is a skew
of around 1 (positive skew).
Although the values of skew and kurtosis are informative, we can convert these values
to z-scores. We saw in section 1.7.4 that a z-score is simply a score from a distribution
that has a mean of 0 and a standard deviation of 1. We also saw that this distribution has
known properties that we can use. Converting scores to a z-score can be useful (if treated
with suitable caution) because (1) we can compare skew and kurtosis values in different
samples that used different measures, and (2) we can see how likely our values of skew and
kurtosis are to occur. To transform any score to a z-score you simply subtract the mean of
the distribution (in this case zero) and then divide by the standard deviation of the distribu-
tion (in this case we use the standard error). Skew and kurtosis are converted to z-scores
in exactly this way.

S−0 K −0
zskewness = zkurtosis =
SEskewness SEkurtosis

5
The observant will notice that the values of kurtosis differ, this is because describe() produces an unbiased esti-
mate (DeCarlo, 1997) whereas [Link]() produces a biased one.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 175

In the above equations, the values of S (skew) and K (kurtosis) and their respective stan-
dard errors are produced by R. These z-scores can be compared against values that you
would expect to get by chance alone (i.e., known values for the normal distribution shown
in the Appendix). So, an absolute value greater than 1.96 is significant at p < .05, above
2.58 is significant at p < .01, and above 3.29 is significant at p < .001. Large samples
will give rise to small standard errors and so when sample sizes are big, significant values
arise from even small deviations from normality. In smallish samples it’s OK to look for
values above 1.96; however, in large samples this criterion should be increased to the 2.58
one and in very large samples, because of the problem of small standard errors that I’ve
described, no criterion should be applied. If you have a large sample (200 or more) it is
more important to look at the shape of the distribution visually and to look at the value of
the skew and kurtosis statistics rather than calculate their significance.
var n mean sd median trimmed mad min max range skew kurtosis se
1 809 1.77 0.69 1.79 1.77 0.70 0.02 3.69 3.67 0.00 -0.41 0.02
2 264 0.96 0.72 0.79 0.87 0.61 0.00 3.44 3.44 1.08 0.82 0.04
3 123 0.98 0.71 0.76 0.90 0.61 0.02 3.41 3.39 1.01 0.73 0.06
Output 5.1

day1 day2 day3

median 1.790000000 7.900000e-01 7.600000e-01
mean 1.770828183 9.609091e-01 9.765041e-01
[Link] 0.024396670 4.436095e-02 6.404352e-02
[Link].0.95 0.047888328 8.734781e-02 1.267805e-01
var 0.481514784 5.195239e-01 5.044934e-01
[Link] 0.693912663 7.207801e-01 7.102770e-01
[Link] 0.391857702 7.501022e-01 7.273672e-01
skewness -0.003155393 1.082811e+00 1.007813e+00
skew.2SE -0.018353763 3.611574e+00 2.309035e+00
kurtosis -0.423991408 7.554615e-01 5.945454e-01
kurt.2SE -1.234611514 1.264508e+00 6.862946e-01
normtest.W 0.995907247 9.083185e-01 9.077513e-01
normtest.p 0.031846386 1.281495e-11 3.804334e-07
Output 5.2

The [Link]() function produces skew.2SE and kurt.2SE, which are the skew and kur-
tosis value divided by 2 standard errors. Remember that z is significant if it is greater than
2 (well, 1.96), therefore this statistic is simply the equations above in a slightly different
format. We have said that if the skew divided by its standard error is greater than 2 then it
is significant (at p < .05), which is the same as saying that if the skew divided by 2 times
the standard error is greater than 1 then it is significant (at p < .05). In other words, if
skew.2SE or kurt.2SE are greater than 1 (ignoring the plus or minus sign) then you have
significant skew/kurtosis (at p < .05); values greater than 1.29 indicate significance at p
< .01, and above 1.65 indicate significance at p < .001. However, as I have just said, you
would only use this criterion in fairly small samples so you need to interpret these values
of skew.2SE or kurt.2SE cautiously.
For the hygiene scores, the values of skew.2SE are −0.018, 3.612, and 2.309 for days 1,
2 and 3 respectively, indicating significant skew on days 2 and 3; the values of kurt.2SE
are −1.235, 1.265, and 0.686, indicating significant kurtosis on days 1 and 2, but not day
3. However, bear in mind what I just said about large samples because our sample size is
pretty big so the histograms are better indicators of the shape of the distribution.
The output of [Link]() also gives us the Shapiro–Wilk test of normality, which we look
at in some detail in section 5.6. For the time being, just note that the test and its probability
value can be found in Output 5.2 labelled as normtest.W and normtest.p.
176 D I S C O VE R I N G STAT I ST I C S US I N G R

Changing how many decimal places are

R ’ s S o ul s ’ T i p 5 . 2 displayed in your output 1

Output 5.2 looks pretty horrible because of all of the decimal places and the scientific notation (i.e., 7.900000e–
01). Most of this precision is unnecessary for everyday purposes. However, we can easily convert our output
using the round() function. This function takes the general form:
round(object that we want to round, digits = x)

Therefore, we can stick an object into this function and then set digits to be the number of decimal places that we
want. For example, if we wanted Output 5.2 to be displayed to 3 decimal places we could execute:
round([Link](dlf[, c("day1", "day2", "day3")], basic = FALSE, norm = TRUE), digits
= 3)

Note that we have simply placed the original command ([Link](dlf[, c(“day1”, “day2”, “day3”)], basic = FALSE,
norm = TRUE)) within the round() function, and then set digits to be 3. The result is a more palatable output:
day1 day2 day3
median 1.790 0.790 0.760
mean 1.771 0.961 0.977
[Link] 0.024 0.044 0.064
[Link].0.95 0.048 0.087 0.127
var 0.482 0.520 0.504
[Link] 0.694 0.721 0.710
[Link] 0.392 0.750 0.727
skewness -0.003 1.083 1.008
skew.2SE -0.018 3.612 2.309
kurtosis -0.424 0.755 0.595
kurt.2SE -1.235 1.265 0.686
normtest.W 0.996 0.908 0.908
normtest.p 0.032 0.000 0.000

CRAMMING SAM’S TIPS Skew and kurtosis

v To check that the distribution of scores is approximately normal, we need to look at the values of skew and kurtosis in the
output.
v Positive values of skew indicate too many low scores in the distribution, whereas negative values indicate a build-up of high
scores.
v Positive values of kurtosis indicate a pointy and heavy-tailed distribution, whereas negative values indicate a flat and light-
tailed distribution.
v The further the value is from zero, the more likely it is that the data are not normally distributed.
v You can test the significance of these values of skew and kurtosis, but these tests should not be used in large samples
(because they are likely to be significant even when skew and kurtosis are not too different from normal).
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 177

5.5.3. Exploring groups of data 1

Sometimes we have data in which there are different groups of entities (cats
and dogs, different universities, people with depression and people without, Can I analyse
groups of data?
for example). There are several ways to produce basic descriptive statistics for
separate groups of people (and we will come across some of these methods in
section 5.6.1). However, I intend to use this opportunity to introduce you to
the by() function and reintroduce the subset() function from Chapter 3. These
functions allow you to specify a grouping variable which splits the data, or to
select a subset of cases.
You’re probably getting sick of the hygiene data from the Download Festival
so let’s use the data in the file [Link]. This file contains data regarding stu-
dents’ performance on an R exam. Four variables were measured: exam (first-
year R exam scores as a percentage), computer (measure of computer literacy
as a percentage), lecture (percentage of R lectures attended) and numeracy (a
measure of numerical ability out of 15). There is a variable called uni indicating whether
the student attended Sussex University (where I work) or Duncetown University. Let’s
begin by looking at the data as a whole.

[Link]. Running the analysis for all data 1

To begin with, open the file [Link] by executing:

rexam <- [Link]("[Link]", header=TRUE)

The variable uni will have loaded in as numbers rather than as text, because that was
how it was specified in the data file; therefore, we need to set the variable uni to be a factor
by executing (see section [Link]):
rexam$uni<-factor(rexam$uni, levels = c(0:1), labels = c("Duncetown
University", "Sussex University"))

Remember that this command takes the variable uni from the rexam dataframe (rexam$uni),
specifies the numbers used to code the two universities, 0 and 1 (levels = c(0:1)), and then
assigns labels to them so that 0 represents Duncetown University, and 1 represents Sussex
University (labels = c(“Duncetown University”, “Sussex University”)).

SELF-TEST

9 Using what you have learnt so far, obtain descriptive

statistics and draw histograms of first-year exam
scores, computer literacy, numeracy and lectures
attended.

Assuming you completed the self-test, you should see something similar to what’s in
Output 5.3 (I used [Link]()) and Figure 5.3. From Output 5.3, we can see that, on
average, students attended nearly 60% of lectures, obtained 58% in their R exam,
scored only 51% on the computer literacy test, and only 4.85 out of 15 on the numer-
acy test. In addition, the standard deviation for computer literacy was relatively small
178 D I S C O VE R I N G STAT I ST I C S US I N G R

FIGURE 5.3 0.025

Histograms
0.10
of the R exam
data 0.020
0.08

0.015
Density

Density
0.06

0.010
0.04

0.005 0.02

0.000 0.0
20 40 60 80 100 30 40 50 60 70
First Year Exam Score Computer Literacy

0.025

0.020 0.3

0.015
Density
Density

0.2

0.010

0.1
0.005

0.000 0.0
20 40 60 80 100 2 4 6 8 10 12 14
Percentage of Lectures Attended Numeracy

compared to that of the percentage of lectures attended and exam scores. The other
important measures are the skew and the kurtosis, and their associated tests of sig-
nificance. We came across these measures earlier on and found that we can interpret
absolute values of kurt.2SE and skew.2SE greater than 1, 1.29, and 1.65 as significant
p < .05, p < .01, and p < .001, respectively. We can see that for skew, numeracy scores
are significantly positively skewed (p < .001) indicating a pile-up of scores on the left
of the distribution (so most students got low scores). For kurtosis, prior exam scores
are significant (p < .05).
The histograms show us several things. The exam scores are very interesting because this
distribution is quite clearly not normal; in fact, it looks suspiciously bimodal (there are two
peaks, indicative of two modes). This observation corresponds with the earlier informa-
tion from the table of descriptive statistics. It looks as though computer literacy is fairly
normally distributed (a few people are very good with computers and a few are very bad,
but the majority of people have a similar degree of knowledge), as is the lecture attendance.
Finally, the numeracy test has produced very positively skewed data (i.e., the majority of
people did very badly on this test and only a few did well). This corresponds to what the
skew statistic indicated.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 179

exam computer lectures numeracy

median 60.000 51.500 62.000 4.000
mean 58.100 50.710 59.765 4.850
[Link] 2.132 0.826 2.168 0.271
[Link].0.95 4.229 1.639 4.303 0.537
var 454.354 68.228 470.230 7.321
[Link] 21.316 8.260 21.685 2.706
[Link] 0.367 0.163 0.363 0.558
skewness -0.104 -0.169 -0.410 0.933
skew.2SE -0.215 -0.350 -0.849 1.932
kurtosis -1.148 0.221 -0.285 0.763
kurt.2SE -1.200 0.231 -0.298 0.798
normtest.W 0.961 0.987 0.977 0.924
normtest.p 0.005 0.441 0.077 0.000
Output 5.3

Descriptive statistics and histograms are a good way of getting an instant picture of the
distribution of your data. This snapshot can be very useful: for example, the bimodal distri-
bution of R exam scores instantly indicates a trend that students are typically either very good
at statistics or struggle with it (there are relatively few who fall in between these extremes).
Intuitively, this finding fits with the nature of the subject: statistics is very easy once every-
thing falls into place, but before that enlightenment occurs it all seems hopelessly difficult.

[Link]. Running the analysis for different groups 1

If we want to obtain separate descriptive statistics for each of the universities, we can use
the by() function.6 The by() function takes the general form:
by(data = dataFrame, INDICES = grouping variable, FUN = a function that you
want to apply to the data)

In other words, we simply enter the name of our dataframe or variables that we’d like to anal-
yse, we specify a variable by which we want to split the output (in this case, it’s uni, because we
want separate statistics for each university), and we tell it which function we want to apply to
the data (in this case we could use describe or [Link]). Therefore, to get descriptive statistics
for the variable exam for each university separately using describe, we could execute:
by(data = rexam$exam, INDICES = rexam$uni, FUN = describe)

To do the same, but using [Link]() instead of describe() we could execute:

by(data = rexam$exam, INDICES = rexam$uni, FUN = [Link])

In both cases, we can get away with not explicitly using data, INDICES and FUN as long
as we order the variables in the order in the functions above; so, these commands have the
same effect as those above:
by(rexam$exam, rexam$uni, describe)
by(rexam$exam, rexam$uni, [Link])

Finally, you can include any options for the function you’re using by adding them in at the
end; for example, if you’re using [Link]() you can specify not to have basic statistics and
to have normality statistics by including those options:
by(rexam$exam, rexam$uni, [Link], basic = FALSE, norm = TRUE)

6
by() is what is known as a ‘wrapper’ function – that is, it takes a more complicated function and simplifies it
for people like me. by() is a wrapper for a very powerful and clever function, called tapply(), which can do all
sorts of things, but is harder to use, so we use by() instead, which just takes our commands and turns them into
commands for tapply().
180 D I S C O VE R I N G STAT I ST I C S US I N G R

If we want descriptive statistics for multiple variables, then we can use cbind() (see R’s
Souls’ Tip 3.5) to include them within the by() function. For example, to look at the
descriptive statistics of both the previous R exam and the numeracy test, we could execute:
by(cbind(data=rexam$exam,data=rexam$numeracy), rexam$uni, describe)

or
by(rexam[, c("exam", "numeracy")], rexam$uni, [Link], basic = FALSE,
norm = TRUE)

Note that the resulting Output 5.4 (which was created using describe rather than
[Link]) is split into two sections: first the results for students at Duncetown
University, then the results for those attending Sussex University. From these tables it
is clear that Sussex students scored higher on both their R exam (called V1 here) and
the numeracy test than their Duncetown counterparts. In fact, looking at the means
reveals that, on average, Sussex students scored an amazing 36% more on the R exam
than Duncetown students, and had higher numeracy scores too (what can I say, my
students are the best).
INDICES: Duncetown University
var n mean sd median trimmed mad min max range skew kurtosis se
V1 1 50 40.18 12.59 38 39.85 12.60 15 66 51 0.29 -0.57 1.78
V2 2 50 4.12 2.07 4 4.00 2.22 1 9 8 0.48 -0.48 0.29
---------------------------------------------------------------------
INDICES: Sussex University
var n mean sd median trimmed mad min max range skew kurtosis se
V1 1 50 76.02 10.21 75 75.70 8.90 56 99 43 0.26 -0.26 1.44
V2 2 50 5.58 3.07 5 5.28 2.97 1 14 13 0.75 0.26 0.43
Output 5.4

Next, we’ll look at the histograms. It might be possible to use by() with ggplot2() to draw
histograms, but if it is the command will be so complicated that no one will understand it.
A simple way, therefore, to create plots for different groups is to use the subset() function,
which we came across in Chapter 3 (section 3.9.2) to create an object containing only the
data in which we’re interested. For example, if we wanted to create separate histograms for
the Duncetown and Sussex Universities then we could create new dataframes that contain
data from only one of the two universities. For example, execute:
dunceData<-subset(rexam, rexam$uni=="Duncetown University")
sussexData<-subset(rexam, rexam$uni=="Sussex University")

These commands each create a new dataframe that is based on a subset of the rexam
dataframe; the subset is determined by the condition in the function. The first command
contains the condition rexam$uni==“Duncetown University”, which means that if the
value of the variable uni is exactly equal to the phrase “Duncetown University” then
the case is selected. In other words, it will retain all cases for which uni is Duncetown
University. Therefore, I’ve called the resulting dataframe dunceData. The second com-
mand does the same thing but this time specifies that uni must be exactly equal to the
phrase ‘Sussex University’. The resulting dataframe, sussexData, contains only the Sussex
University scores. This is a quick and easy way to split groups; however, you need to be
careful that the term you specify to select cases (e.g., ‘Duncetown University’) exactly
matches (including capital letters and spaces) the labelling in the data set otherwise you’ll
end up with an empty data set.
Having created our separate dataframes, we can generate histograms using the same
commands as before, but specifying the dataframe for the subset of data. For example, to
create a histogram of the numeracy scores for Duncetown University, we could execute:
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 181

[Link] <- ggplot(dunceData, aes(numeracy)) + opts(legend.

position = "none") + geom_histogram(aes(y = ..density..), fill = "white",
colour = "black", binwidth = 1) + labs(x = "Numeracy Score", y = "Density")
+ stat_function(fun=dnorm, args=list(mean = mean(dunceData$numeracy,
[Link] = TRUE), sd = sd(dunceData$numeracy, [Link] = TRUE)), colour = "blue",
size=1)
[Link]

Compare this code with that in section 5.5.1; note that it is exactly the same, but we have
used the dunceData dataframe instead of using the whole data set.7 We could create the
same plot for the Sussex University students by simply using sussexData in place of dunce-
Data in the command.
We could repeat these commands for the exam scores by replacing ‘numeracy’ with
‘exam’ throughout the commands above (this will have the effect of plotting exam scores
rather than numeracy scores). Figure 5.4 shows the histograms of exam scores and numer-
acy split according to the university attended. The first interesting thing to note is that for

Duncetown University Sussex University FIGURE 5.4

Distributions
0.005 of exam and
0.06 numeracy scores
0.004 for Duncetown
0.05 University and
0.003 Sussex University
0.04
Density

Density

students
0.03
0.002
0.02
0.001
0.01

0.000 0.00
20 30 40 50 60 60 70 80 90 100
First Year Exam Score First Year Exam Score

0.20 0.15

0.15
Density
Density

0.10

0.10

0.05
0.05

0.00 0.00
0 2 4 6 8 10 0 2 4 6 8 10 12 14
Numeracy Score Numeracy Score

7
Note that I have included ‘binwidth = 1’ (see Chapter 4) for the numeracy scores because it makes the result-
ing plot look better; for the other variables this option can be excluded because the default bin width produces
nice-looking plots.
182 D I S C O VE R I N G STAT I ST I C S US I N G R

exam marks the distributions are both fairly normal. This seems odd because the overall
distribution was bimodal. However, it starts to make sense when you consider that for
Duncetown the distribution is centred on a mark of about 40%, but for Sussex the distri-
bution is centred on a mark of about 76%. This illustrates how important it is to look at
distributions within groups. If we were interested in comparing Duncetown to Sussex it
wouldn’t matter that overall the distribution of scores was bimodal; all that’s important is
that each group comes from a normal distribution, and in this case it appears to be true.
When the two samples are combined, these two normal distributions create a bimodal
one (one of the modes being around the centre of the Duncetown distribution, and the
other being around the centre of the Sussex data). For numeracy scores, the distribution is
slightly positively skewed (there is a larger concentration at the lower end of scores) in both
the Duncetown and Sussex groups. Therefore, the overall positive skew observed before is
due to the mixture of universities.

SELF-TEST

9 Repeat these analyses for the computer literacy and

percentage of lectures attended and interpret the
results.

5.6. Testing whether a distribution is normal 1

Another way of looking at the problem is to see whether the distribution as

Is it possible to test a whole deviates from a comparable normal distribution. The Shapiro–Wilk
whether I am normal? test does just this: it compares the scores in the sample to a normally dis-
tributed set of scores with the same mean and standard deviation. If the test
is non-significant (p > .05) it tells us that the distribution of the sample is
not significantly different from a normal distribution. If, however, the test is
significant (p < .05) then the distribution in question is significantly differ-
ent from a normal distribution (i.e., it is non-normal). This test seems great:
in one easy procedure it tells us whether our scores are normally distributed
(nice!). However, it has limitations because with large sample sizes it is very
easy to get significant results from small deviations from normality, and so a
significant test doesn’t necessarily tell us whether the deviation from normality is enough
to bias any statistical procedures that we apply to the data. I guess the take-home message
is: by all means use these tests, but plot your data as well and try to make an informed
decision about the extent of non-normality.

5.6.1. Doing the Shapiro–Wilk test in R 1

We have already encountered the Shapiro–Wilk test as part of the output from the stat.
desc() function (see Output 5.2 and, for these data, Output 5.3). However, we can also use
the [Link]() function. This function takes the general form:
[Link](variable)
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 183

in which variable is the name of the variable that you’d like to test for normality. Therefore,
to test the exam and numeracy variables for normality we would execute:
[Link](rexam$exam)
[Link](rexam$numeracy)

The output is shown in Output 5.5. Note that the value of W corresponds to the value
of normtest.W, and the p-value corresponds to normtest.p from the [Link]() function
(Output 5.3). For each test we see the test statistic, labelled W, and the p-value. Remember
that a significant value (p-value less than .05) indicates a deviation from normality. For
both numeracy (p = .005) and R exam scores (p < .001), the Shapiro–Wilk test is highly
significant, indicating that both distributions are not normal. This result is likely to reflect
the bimodal distribution found for exam scores, and the positively skewed distribution
observed in the numeracy scores. However, these tests confirm that these deviations were
significant (but bear in mind that the sample is fairly big).
Shapiro-Wilk normality test

data: rexam$exam
W = 0.9613, p-value = 0.004991

Shapiro-Wilk normality test

data: rexam$numeracy
W = 0.9244, p-value = 2.424e-05
Output 5.5

As a final point, bear in mind that when we looked at the exam scores for separate
groups, the distributions seemed quite normal; now if we’d asked for separate Shapiro–
Wilk tests for the two universities we might have found non-significant results. In fact, let’s
try this out, using the by() function we came across earlier. We use [Link] as the FUN
instead of describe or [Link], which we have used before (although [Link] would also
give you the Shapiro–Wilk test as part of the output so you could use this function also):
by(rexam$exam, rexam$uni, [Link])
by(rexam$numeracy, rexam$uni, [Link])

You should get Output 5.6 for the exam scores, which shows that the percentages on the
R exam are indeed normal within the two groups (the p-values are greater than .05). This
is important because if our analysis involves comparing groups, then what’s important is
not the overall distribution but the distribution in each group.
rexam$uni: Duncetown University

Shapiro-Wilk normality test

data: dd[x, ]
W = 0.9722, p-value = 0.2829

-------------------------------------------------------------------
rexam$uni: Sussex University

Shapiro-Wilk normality test

data: dd[x, ]
W = 0.9837, p-value = 0.7151
Output 5.6
184 D I S C O VE R I N G STAT I ST I C S US I N G R

For numeracy scores (Output 5.7) the tests are still significant indicating non-normal
distributions both for Duncetown University (p = .015), and Sussex University (p = .007).
rexam$uni: Duncetown University

Shapiro-Wilk normality test

data: dd[x, ]
W = 0.9408, p-value = 0.01451

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

rexam$uni: Sussex University

Shapiro-Wilk normality test

data: dd[x, ]
W = 0.9323, p-value = 0.006787
Output 5.7

We can also draw Q-Q plots for the variables, to help us to interpret the results of the
Shapiro–Wilk test (see Figure 5.5).
qplot(sample = rexam$exam, stat="qq")
qplot(sample = rexam$numeracy, stat="qq")

The normal Q-Q chart plots the values you would expect to get if the distribution were
normal (theoretical values) against the values actually seen in the data set (sample values).
If the data are normally distributed, then the observed values (the dots on the chart) should
fall exactly along a straight line (meaning that the observed values are the same as you
would expect to get from a normally distributed data set). Any deviation of the dots from
the line represents a deviation from normality. So, if the Q-Q plot looks like a wiggly snake
then you have some deviation from normality. Specifically, when the line sags consistently
below the diagonal, or consistently rises above it, then this shows that the kurtosis differs
from a normal distribution, and when the curve is S-shaped, the problem is skewness.
In both of the variables analysed we already know that the data are not normal, and these
plots (see Figure 5.5) confirm this observation because the dots deviate substantially from
the line. It is noteworthy that the deviation is greater for the numeracy scores, and this is
consistent with the higher significance value of this variable on the Shapiro–Wilk test.

FIGURE 5.5 First Year Exam Scores Numeracy

Normal Q-Q plots 14
of numeracy and
R exam scores 12
80
10
sample

sample

60 8

6
40
4

20 2

−2 −1 0 1 2 −2 −1 0 1 2
theoretical theoretical
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 185

5.6.2. Reporting the Shapiro–Wilk test 1

The test statistic for the Shapiro–Wilk test is denoted by W; we can report the results in
Output 5.5 in the following way:
9
The percentage on the R exam, W = 0.96, p = .005, and the numeracy scores, W =
0.92, p < .001, were both significantly non-normal.

CRAMMING SAM’S TIPS Normality tests

v The Shapiro–Wilk test can be used to see if a distribution of scores significantly differs from a normal distribution.
v If the Shapiro–Wilk test is significant (p-value less than .05) then the scores are significantly different from a normal
distribution.
v Otherwise, scores are approximately normally distributed.
v Warning: In large samples this test can be significant even when the scores are only slightly different from a normal dis-
tribution. Therefore, they should always be interpreted in conjunction with histograms, or Q-Q plots, and the values of skew
and kurtosis.

5.7. Testing for homogeneity of variance 1

So far I’ve concentrated on the assumption of normally distributed data; however, at the
beginning of this chapter I mentioned another assumption: homogeneity of variance. This
assumption means that as you go through levels of one variable, the variance of the other
should not change. If you’ve collected groups of data then this means that the variance of
your outcome variable or variables should be the same in each of these groups. If you’ve
collected continuous data (such as in correlational designs), this assumption means that the
variance of one variable should be stable at all levels of the other variable. Let’s illustrate
this with an example. An audiologist was interested in the effects of loud concerts on peo-
ple’s hearing. So, she decided to send 10 people on tour with the loudest band she could
find, Motörhead. These people went to concerts in Brixton (London), Brighton, Bristol,
Edinburgh, Newcastle, Cardiff and Dublin and after each concert the audiologist measured
the number of hours after the concert that these people had ringing in their ears.
Figure 5.6 shows the number of hours that each person had ringing in his or her ears
after each concert (each person is represented by a circle). The horizontal lines represent
the average number of hours that there was ringing in the ears after each concert and these
means are connected by a line so that we can see the general trend of the data. Remember
that for each concert, the circles are the scores from which the mean is calculated. Now, we
can see in both graphs that the means increase as the people go to more concerts. So, after
the first concert their ears ring for about 12 hours, but after the second they ring for about
15–20 hours, and by the final night of the tour, they ring for about 45–50 hours (2 days).
So, there is a cumulative effect of the concerts on ringing in the ears. This pattern is found
in both graphs; the difference between the graphs is not in terms of the means (which are
roughly the same), but in terms of the spread of scores around the mean. If you look at the
left-hand graph, the spread of scores around the mean stays the same after each concert
186 D I S C O VE R I N G STAT I ST I C S US I N G R

(the scores are fairly tightly packed around the mean). Put it another way, if you measured
the vertical distance between the lowest score and the highest score after the Brixton con-
cert, and then did the same after the other concerts, all of these distances would be fairly
similar. Although the means increase, the spread of scores for hearing loss is the same at
each level of the concert variable (the spread of scores is the same after Brixton, Brighton,
Bristol, Edinburgh, Newcastle, Cardiff and Dublin). This is what we mean by homogeneity
of variance. The right-hand graph shows a different picture: if you look at the spread of
scores after the Brixton concert, they are quite tightly packed around the mean (the vertical
distance from the lowest score to the highest score is small), but after the Dublin show (for
example) the scores are very spread out around the mean (the vertical distance from the
lowest score to the highest score is large). This is an example of heterogeneity of variance:
that is, at some levels of the concert variable the variance of scores is different than other
levels (graphically, the vertical distance from the lowest to highest score is different after
different concerts).

FIGURE 5.6
70 70
Graphs illustrating
data with 60 60
homogeneous
50 50

Ringing (Hours)
Ringing (Hours)

(left) and
heterogeneous 40 40

(right) variances
30 30

20 20

10
10

0 0
y

on

ie

ff

in
to
y

on

ie

ff

in

em

rg

di
to

st

ub
em

rg

di

is
ht
st

ub

ar
bu
is

ca
ht

Br
ar

ad

D
bu

ig
ca

C
Br
ad

D
ig

in

ew
C

Br
Ac
in

ew
Br

Ed
Ac

Ed

N
N

n
n

to
to

ix
ix

Br
Br

Concert Concert

5.7.1. Levene’s test 1

Hopefully you’ve got a grip of what homogeneity of variance actually means. Now, how
do we test for it? Well, we could just look at the values of the variances and see whether
they are similar. However, this approach would be very subjective and probably prone to
academics thinking ‘Ooh look, the variance in one group is only 3000 times larger than the
variance in the other: that’s roughly equal’. Instead, in correlational analysis such as regres-
sion we tend to use graphs (see section 7.9.5) and for groups of data we tend to use a test
called Levene’s test (Levene, 1960). Levene’s test tests the null hypothesis that the variances
in different groups are equal (i.e., the difference between the variances is zero). It’s a very
simple and elegant test that works by doing a one-way ANOVA (see Chapter 10) conducted
on the deviation scores; that is, the absolute difference between each score and the mean of
the group from which it came (see Glass, 1966, for a very readable explanation).8 For now,
all we need to know is that if Levene’s test is significant at p ≤ .05 then we can conclude that
the null hypothesis is incorrect and that the variances are significantly different – therefore,
the assumption of homogeneity of variances has been violated. If, however, Levene’s test
8
We haven’t covered ANOVA yet, so this explanation won’t make much sense to you now, but in Chapter 10 we
will look in more detail at how Levene’s test works.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 187

is non-significant (i.e., p > .05) then the variances are roughly equal and the assumption
is tenable.

[Link]. Levene’s test with R Commander 1

First we’ll load the data into R Commander. Choose Data ⇒ Import data ⇒ from text file,
clipboard, or URL… and then select the file [Link] (see section 3.7.3). Before we can
conduct Levene’s test we need to convert uni to a factor because at the moment it is simply
0s and 1s so R doesn’t know that it’s a factor – see section 3.6.2 to remind yourself how to do
that. Once you have done this you should be able to select Statistics⇒Variances⇒Levene’s
test (you won’t be able to select it unless R can ‘see’ a factor in the dataframe). Choosing
this option in the menu opens the dialog box shown in Figure 5.7. You need to select a
grouping variable. R Commander has realized that you only have one variable that could
be the grouping variable – because it is the only factor – and that’s uni. Therefore, it has
already selected this variable.
Choose the variable on the right that you want to test for equality of variances across
the groups defined by uni. You can choose median or mean for the centring – the median
tends to be more accurate and is the default; I use this default throughout the book. Run
the analysis for both exam and numeracy. Output 5.8 shows the results.

FIGURE 5.7
Levene’s test in
R Commander

[Link]. Levene’s test with R 1

To use Levene’s test, we use the leveneTest() function from the car package. This function
takes the general form:
leveneTest(outcome variable, group, center = median/mean)
188 D I S C O VE R I N G STAT I ST I C S US I N G R

Therefore, we enter two variables into the function: first the outcome variable of which we
want to test the variances; and second, the grouping variable, which must be a factor. We
can just enter these variables and Levene’s test will centre the variables using the median
(which is slightly preferable), but if we want to override this default and centre using the
mean then we can add the option center = “mean”. Therefore, for the exam scores we
could execute:
leveneTest(rexam$exam, rexam$uni)
leveneTest(rexam$exam, rexam$uni, center = mean)

For the numeracy scores we would execute (note that all we have changed is the outcome
variable):
leveneTest(rexam$numeracy, rexam$uni)

[Link]. Levene’s test output 1

Output 5.8 shows the output for Levene’s test for exam scores (using the median), exam
scores (centring using the mean) and numeracy scores. The result is non-significant for the
R exam scores (the value in the Pr (>F) column is more than .05) regardless of whether
we centre with the median or mean. This indicates that the variances are not significantly
different (i.e., they are similar and the homogeneity of variance assumption is tenable).
However, for the numeracy scores, Levene’s test is significant (the value in the Pr (>F)
column is less than .05) indicating that the variances are significantly different (i.e., they
are not the same and the homogeneity of variance assumption has been violated).
> leveneTest(rexam$exam, rexam$uni)

Levene’s Test for Homogeneity of Variance (center = median)

Df F value Pr(>F)
group 1 2.0886 0.1516
98

> leveneTest(rexam$exam, rexam$uni, center = mean)

Levene’s Test for Homogeneity of Variance (center = mean)

Df F value Pr(>F)
group 1 2.5841 0.1112
98

> leveneTest(rexam$numeracy, rexam$uni)

Levene’s Test for Homogeneity of Variance (center = median)

Df F value Pr(>F)
group 1 5.366 0.02262 *
98
Output 5.8

5.7.2. Reporting Levene’s test 1

Levene’s test can be denoted with the letter F and there are two different degrees of free-
dom. As such you can report it, in general form, as F(df1, df2) = value, Pr (>F). So, for the
results in Output 5.8 we could say:
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 189

9
For the percentage on the R exam, the variances were similar for Duncetown and
Sussex University students, F(1, 98) = 2.09, ns, but for numeracy scores the variances
were significantly different in the two groups, F(1, 98) = 5.37, p = .023.

5.7.3. Hartley’s Fmax: the variance ratio 1

As with the Shapiro–Wilk test (and other tests of normality), when the sample size is
large, small differences in group variances can produce a Levene’s test that is significant
(because, as we saw in Chapter 1, the power of the test is improved). A useful double
check, therefore, is to look at Hartley’s Fmax – also known as the variance ratio (Pearson
& Hartley, 1954). This is the ratio of the variances between the group with the biggest
variance and the group with the smallest variance. This ratio was compared to critical
values in a table published by Hartley. Some of the critical values (for a .05 level of sig-
nificance) are shown in Figure 5.8 (see Oliver Twisted); as you can see, the critical values
depend on the number of cases per group (well, n − 1 actually), and the number of vari-
ances being compared. From this graph you can see that with sample sizes (n) of 10 per
group, an Fmax of less than 10 is more or less always going to be non-significant, with
15–20 per group the ratio needs to be less than about 5, and with samples of 30–60 the
ratio should be below about 2 or 3.

30 FIGURE 5.8
Number of Variances being Compared
2 Variances Selected critical
25 3 Variances
4 Variances values for
6 Variances
8 Variances Hartley’s Fmax test
20 10 Variances
Critical Value

15

10

0
5 6 7 8 9 10 12 15 20 30 60
n–1

OLIVER TWISTED Oliver thinks that my graph of critical values is stupid. ‘Look at that graph,’
he laughed, ‘it’s the most stupid thing I’ve ever seen since I was at Sussex
Please Sir, can I have Uni and I saw my statistics lecturer, Andy Fie…’. Well, go choke on your
some more … Hartley’s gruel you Dickensian bubo, because the full table of critical values is
Fmax? in the additional material for this chapter on the companion website.
190 D I S C O VE R I N G STAT I ST I C S US I N G R

CRAMMING SAM’S TIPS Homogeneity of variance

v Homogeneity of variance is the assumption that the spread of scores is roughly equal in different groups of cases, or more
generally that the spread of scores is roughly equal at different points on the predictor variable.
v When comparing groups, this assumption can be tested with Levene’s test.
v If Levene’s test is significant (Pr (>F) in the R output is less than .05) then the variances are significantly different in different
groups.
v Otherwise, homogeneity of variance can be assumed.
v The variance ratio is the largest group variance divided by the smallest. This value needs to be smaller than the critical values
in Figure 5.8.
v Warning: In large samples Levene’s test can be significant even when group variances are not very different. Therefore, it
should be interpreted in conjunction with the variance ratio.

5.8. Correcting problems in the data 2

The previous section showed us various ways to explore our data; we saw how to look for
problems with our distribution of scores and how to detect heterogeneity of variance. In
Chapter 4 we also discovered how to spot outliers in the data. The next question is what
to do about these problems.

5.8.1. Dealing with outliers 2

If you detect outliers in the data there are several options for reducing the impact of these
values. However, before you do any of these things, it’s worth checking that the data have
been entered correctly for the problem cases. If the data are correct then the three main
options you have are:
1 Remove the case: This entails deleting the data from the person who contributed the
outlier. However, this should be done only if you have good reason to believe that
this case is not from the population that you intended to sample. For example, if
you were investigating factors that affected how much cats purr and one cat didn’t
purr at all, this would likely be an outlier (all cats purr). Upon inspection, if you dis-
covered that this cat was actually a dog wearing a cat costume (hence why it didn’t
purr), then you’d have grounds to exclude this case because it comes from a different
population (dogs who like to dress as cats) than your target population (cats).
2 Transform the data: Outliers tend to skew the distribution and, as we will see in the
next section, this skew (and, therefore, the impact of the outliers) can sometimes be
reduced by applying transformations to the data.
3 Change the score: If transformation fails, then you can consider replacing the score.
This on the face of it may seem like cheating (you’re changing the data from what
was actually corrected); however, if the score you’re changing is very unrepresenta-
tive and biases your statistical model anyway then changing the score is the lesser of
two evils! There are several options for how to change the score:
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 191

a The next highest score plus one: Change the score to be one unit above the next
highest score in the data set.
b Convert back from a z-score: A z-score of 3.29 constitutes an outlier (see Jane
Superbrain Box 4.1), so we can calculate what score would give rise to a z-score
of 3.29 (or perhaps 3) by rearranging the z-score equation in section 1.7.4, which
– –
gives us X = (z × s) + X. All this means is that we calculate the mean (X) and stand-
ard deviation (s) of the data; we know that z is 3 (or 3.29 if you want to be exact)
so we just add three times the standard deviation to the mean, and replace our
outliers with that score.
c The mean plus two standard deviations: A variation on the above method is to use
the mean plus two times the standard deviation (rather than three times the stand-
ard deviation).

5.8.2. Dealing with non-normality and unequal variances 2

[Link]. Transforming data 2

This section is quite hair raising so don’t worry if it doesn’t make much sense – many
undergraduate courses won’t cover transforming data so feel free to ignore this section if
you want to.
We saw in the previous section that you can deal with outliers by transforming the data
and that these transformations are also useful for correcting problems with normality and
the assumption of homogeneity of variance. The idea behind transformations
is that you do something to every score to correct for distributional problems,
outliers or unequal variances. Although some students often (understandably) What do I do if
think that transforming data sounds dodgy (the phrase ‘fudging your results’ my data are
springs to some people’s minds!), in fact it isn’t because you do the same thing not normal?
9
to all of your scores. As such, transforming the data won’t change the relation-
ships between variables (the relative differences between people for a given
variable stay the same), but it does change the differences between different
variables (because it changes the units of measurement). Therefore, if you are
looking at relationships between variables (e.g., regression) it is alright just to
transform the problematic variable, but if you are looking at differences within
variables (e.g., change in a variable over time) then you need to transform all
levels of those variables.
Let’s return to our Download Festival data ([Link]) from earlier in the
chapter. These data were not normal on days 2 and 3 of the festival (section 5.4). Now, we
might want to look at how hygiene levels changed across the three days (i.e., compare the
mean on day 1 to the means on days 2 and 3 to see if people got smellier). The data for
days 2 and 3 were skewed and need to be transformed, but because we might later compare
the data to scores on day 1, we would also have to transform the day 1 data (even though
scores were not skewed). If we don’t change the day 1 data as well, then any differences in
hygiene scores we find from day 1 to day 2 or 3 will be due to us transforming one variable
and not the others.

9
Although there aren’t statistical consequences of transforming data, there may be empirical or scientific implica-
tions that outweigh the statistical benefits (see Jane Superbrain Box 5.1).
192 D I S C O VE R I N G STAT I ST I C S US I N G R

Table 5.1 Data transformations and their uses

Data Transformation Can Correct For
Log transformation (log(Xi)): Taking the logarithm of a set of numbers Positive skew,
squashes the right tail of the distribution. As such it’s a good way to reduce unequal variances
positive skew. However, you can’t take the log of zero or negative numbers,
so if your data tend to zero or produce negative numbers you need to add a
constant to all of the data before you do the transformation. For example, if you
have zeros in the data then do log(Xi + 1), or if you have negative numbers add
whatever value makes the smallest number in the data set positive.
Square root transformation (√Xi): Taking the square root of large Positive skew,
values has more of an effect than taking the square root of small values. unequal variances
Consequently, taking the square root of each of your scores will bring any
large scores closer to the centre – rather like the log transformation. As
such, this can be a useful way to reduce positive skew; however, you still
have the same problem with negative numbers (negative numbers don’t
have a square root).
Reciprocal transformation (1/Xi): Dividing 1 by each score also reduces Positive skew,
the impact of large scores. The transformed variable will have a lower unequal variances
limit of 0 (very large numbers will become close to 0). One thing to bear
in mind with this transformation is that it reverses the scores: scores that
were originally large in the data set become small (close to zero) after
the transformation, but scores that were originally small become big after
the transformation. For example, imagine two scores of 1 and 10; after
the transformation they become 1/1 = 1 and 1/10 = 0.1: the small score
becomes bigger than the large score after the transformation. However, you
can avoid this by reversing the scores before the transformation, by finding
the highest score and changing each score to the highest score minus the
score you’re looking at. So, you do a transformation 1/(XHighest−Xi).
Reverse score transformations: Any one of the above transformations Negative skew
can be used to correct negatively skewed data, but first you have to reverse
the scores. To do this, subtract each score from the highest score obtained,
or the highest score + 1 (depending on whether you want your lowest
score to be 0 or 1). If you do this, don’t forget to reverse the scores back
afterwards, or to remember that the interpretation of the variable is reversed:
big scores have become small and small scores have become big!

There are various transformations that you can do to the data that are helpful in correct-
ing various problems.10 However, whether these transformations are necessary or useful is
quite a complex issue (see Jane Superbrain Box 5.1). Nevertheless, because they are used by
researchers Table 5.1 shows some common transformations and their uses.

[Link]. Choosing a transformation 2

Given that there are many transformations that you can do, how can you decide which one
is best? The simple answer is trial and error: try one out and see if it helps and if it doesn’t

10
You’ll notice in this section that I keep writing Xi. We saw in Chapter 1 that this refers to the observed score for
the ith person (so, the i could be replaced with the name of a particular person, thus for Graham, Xi = XGraham =
Graham’s score, and for Carol, Xi = XCarol = Carol’s score).
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 193

then try a different one. If you are looking at differences between variables you must apply
the same transformation to all variables (you cannot, for example, apply a log transforma-
tion to one variable and a square root transformation to another). This can be quite time
consuming.

that their conclusion was incorrect, which Levine and

Dunlap (1983) contested in a response to the response.
Finally, in a response to the response to the response,
Games (1984) pointed out several important questions
to consider:

1. The central limit theorem (section 2.5.1) tells us

that in big samples the sampling distribution will
JANE SUPERBRAIN 5.1 be normal regardless, and this is what’s actually
important, so the debate is academic in anything
To transform or not to transform, that is the
other than small samples. Lots of early research
question 3
did indeed show that with samples of 40 the nor-
mality of the sampling distribution was, as pre-
Not everyone agrees that transforming data is a good
dicted, normal. However, this research focused
idea; for example, Glass, Peckham and Sanders (1972),
on distributions with light tails and subsequent
in a very extensive review, commented that ‘the payoff of
work has shown that with heavy-tailed distributions
normalizing transformations in terms of more valid prob-
larger samples would be necessary to invoke the
ability statements is low, and they are seldom considered
central limit theorem (Wilcox, 2005). This research
to be worth the effort’ (p. 241). In which case, should we
suggests that transformations might be useful for
bother?
such distributions.
The issue is quite complicated (especially for this
2. By transforming the data you change the hypoth-
early in the book), but essentially we need to know
esis being tested (when using a log transformation
whether the statistical models we apply perform better
and comparing means you change from comparing
on transformed data than they do when applied to data
arithmetic means to comparing geometric means).
that violate the assumption that the transformation cor-
Transformation also means that you’re now address-
rects. If a statistical model is still accurate even when
ing a different construct than the one originally
its assumptions are broken it is said to be a robust test
measured, and this has obvious implications for
(section 5.8.4). I’m not going to discuss whether particu-
interpreting that data (Gelman & Hill, 2007; Grayson,
lar tests are robust here, but I will discuss the issue for
2004).
particular tests in their respective chapters. The question
3. In small samples it is tricky to determine normality one
of whether to transform is linked to this issue of robust-
way or another (tests such as Shapiro–Wilk will have
ness (which in turn is linked to what test you are perform-
low power to detect deviations from normality and
ing on your data).
graphs will be hard to interpret with so few data points).
A good case in point is the F-test in ANOVA (see
4. The consequences for the statistical model of apply-
Chapter 10), which is often claimed to be robust (Glass
ing the ‘wrong’ transformation could be worse than the
et al., 1972). Early findings suggested that F performed
consequences of analysing the untransformed scores.
as it should in skewed distributions and that transform-
ing the data helped as often as it hindered the accuracy As we will see later in the book, there is an exten-
of F (Games & Lucas, 1966). However, in a lively but sive library of robust tests that can be used and which
informative exchange, Levine and Dunlap (1982) showed have considerable benefits over transforming data. The
that transformations of skew did improve the perform- definitive guide to these is Wilcox’s (2005) outstanding
ance of F; however, in a response, Games (1983) argued book.
194 D I S C O VE R I N G STAT I ST I C S US I N G R

5.8.3. Transforming the data using R 2

[Link]. Computing new variables 2

Transformations are very easy using R. We use one of two general commands:

newVariable <- function(oldVariable)

in which function is the function we will use to transform the variable. Or possibly:

newVariable <- arithmetic with oldVariable(s)

Let’s first look at some of the simple arithmetic functions:

+ Addition: We can add two variables together, or add a constant to our variables. For
example, with our hygiene data, ‘day1 + day2’ creates a column in which each row
contains the hygiene score from the column labelled day1 added to the score from the
column labelled day2 (e.g., for participant 1: 2.65 + 1.35 = 4). In R we would execute:
dlf$day1PlusDay2 <- dlf$day1 + dlf$day2
which creates a new variable day1PlusDay2 in the dlt dataframe based on adding the
variables day1 and day2.

− Subtraction: We can subtract one variable from another. For example, we could
subtract the day 1 hygiene score from the day 2 hygiene score. This creates a new
variable in our dataframe in which each row contains the score from the column
labelled day1 subtracted from the score from the column labelled day2 (e.g., for
participant 1: 1.35 − 2.65 = −1.30). Therefore, this person’s hygiene went down by
1.30 (on our 5-point scale) from day 1 to day 2 of the festival. In R we would execute:
dlf$day2MinusDay1 <- dlf$day2 - dlf$day1
which creates a new variable day2MinusDay1 in the dlf dataframe based on subtracting
the variable day1 from day2.

∗ Multiply: We can multiply two variables together, or we can multiply a variable by any
number. In R, we would execute:
dlf$day2Times5 <- dlf$day1 * 5
which creates a new variable day2Times5 in the dlf dataframe based on multiplying
day1 by 5.

∗∗ Exponentiation: Exponentiation is used to raise the preceding term by the power

OR^ of the succeeding term. So ‘day1**2’ or ‘day1^2’ (it doesn’t matter which you use)
creates a column that contains the scores in the day1 column raised to the power of
2 (i.e., the square of each number in the day1 column: for participant 1, 2.652 =7.02).
Likewise, ‘day1**3’ creates a column with values of day1 cubed. In R, we would
execute either:
dlf$day2Squared <- dlf$day2 ** 2
or
dlf$day2Squared <- dlf$day2 ^ 2
both of which create a new variable day2Squared in the dlf dataframe based on
squaring values of day2.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 195

< Less than: This is a logical operator – that means it gives the answer TRUE (or 1)
or FALSE (or 0). If you typed ‘day1 < 1’, R would give the answer TRUE to those
participants whose hygiene score on day 1 of the festival was less than 1 (i.e., if day1 was
0.9999 or less). So, we might use this if we wanted to look only at the people who were
already smelly on the first day of the festival. In R we would execute:
dlf$day1LessThanOne <- dlf$day1 < 1

to create a new variable day1LessThanOne in the dlf dataframe for which the values
are TRUE (or 1) if the value of day1 is less than 1, but FALSE (or 0) if the value of day1 is
greater than 1.
<= Less than or equal to: This is the same as above but returns a response of TRUE (or 1) if the
value of the original variable is equal to or less than the value specified. In R we would execute:
dlf$day1LessThanOrEqualOne <- dlf$day1 <= 1
to create a new variable day1LessThanOrEqualOne in the dlf dataframe for which the
values are TRUE (or 1) if the value of day1 is less than or equal to 1, but FALSE (or 0) if
the value of day1 is greater than 1.
> Greater than: This is the opposite of the less than operator above. It returns a response
of TRUE (or 1) if the value of the original variable is greater than the value specified. In R
we would execute:
dlf$day1GreaterThanOne <- dlf$day1 > 1
to create a new variable day1GreaterThanOne in the dlf dataframe for which the values
are TRUE (or 1) if the value of day1 is greater than 1, but FALSE (or 0) if the value of
day1 is less than 1.
>= Greater than or equal to: This is the same as above but returns a response of TRUE (or
1) if the value of the original variable is equal to or greater than the value specified. In R
we would execute:
dlf$day1GreaterThanOrEqualOne <- dlf$day1 >= 1
to create a new variable day1GreaterThanOrEqualOne in the dlf dataframe for which
the values are TRUE (or 1) if the value of day1 is greater than or equal to 1, but FALSE
(or 0) if the value of day1 is less than 1.
== Double equals means ‘is equal to?’ It’s a question, rather than an assignment, like a single
equals (=). Therefore, if we write something like dlf$gender == “Male” we are asking ‘is the
value of the variable gender in the dlf dataframe equal to the word ‘Male’? In R, if we executed:
dlf$male <- dlf$gender == "Male"
we would create a variable male in the dlf dataframe that contains the value TRUE if the
variable gender was the word ‘Male’ (spelt as it is specified, including capital letters) and
FALSE in all other cases.
!= Not equal to. The opposite of ==. In R, if we executed:
dlf$notMale <- dlf$gender != "Male"
we would create a variable notMale in the dlf dataframe that contains the value TRUE
if the variable gender was not the word ‘Male’ (spelt as it is specified including capital
letters) and FALSE otherwise.

Some of the most useful functions are listed in Table 5.2, which shows the standard form
of the function, the name of the function, an example of how the function can be used
and what R would output if that example were used. There are several basic functions for
calculating means, standard deviations and sums of columns. There are also functions such
as the square root and logarithm that are useful for transforming data that are skewed, and
we will use these functions now.
196 D I S C O VE R I N G STAT I ST I C S US I N G R

Table 5.2 Some useful functions

Function Name Input example Output
rowMeans() Mean for a rowMeans(cbind For each row, R calculates the mean
row (dlf$day1, dlf$day2, hygiene score across the three days of the
dlf$day3), [Link] = festival. [Link] tells R whether to exclude
TRUE) missing values from the calculation (see R’s
Souls’ Tip 5.3).
rowSums() Sums for rowSums(cbind For each row, R calculates the sum of the
a row (dlf$day1, dlf$day2, hygiene scores across the three days of the
dlf$day3), [Link] = festival. [Link] tells R whether to exclude
TRUE) missing values from the calculation (see R’s
Souls’ Tip 5.4).
sqrt() Square sqrt(dlf$day2) Produces a column containing the square
root root of each value in the column labelled day2
abs() Absolute abs(dlf$day1) Produces a variable that contains the
value absolute value of the values in the column
labelled day1 (absolute values are ones
where the signs are ignored: so −5
becomes +5 and +5 stays as +5)
log10() Base 10 log10(dlf$day1) Produces a variable that contains the logarithm
logarithm (to base 10) values of the variable day1.
log() Natural log10(dlf$day1) Produces a variable that contains the natural
logarithm logarithm values of the variable day1.
[Link]() Is [Link](dlf$day1) This is used to determine if a variable is
missing? missing or not. If the variable is missing,
the case will be assigned TRUE (or 1); if
the case is not missing, the case will be
assigned FALSE (or 0).

R ’ s S o ul s ’ T i p 5 . 3 The [Link]() function and missing data 3

If we want to count missing data, we can use [Link](). For example, if we want to know whether a person is missing
for their day 2 hygiene score, we use:
dlf$missingDay2 <- [Link](dlf$day2)

But we can then use that variable in some clever ways. How many people were missing on day 2? Well, we know
that the variable we just created is TRUE (or 1) if they are missing, so we can just add them up:
sum(dlf$missingDay2)

If we want to be lazy, we can embed those functions in each other, and not bother to create a variable:
(sum([Link](dlf$day2))

which tells us that 546 scores are missing. What proportion of scores is that? Well, we have a 1 if they are missing,
and a zero if not. So the mean of that variable will be the proportion which are missing:
mean([Link](dlf$day2))

This tells us that the mean is 0.674, so 67.4% of people are missing a hygiene score on day 2.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 197

[Link]. The log transformation in R 2

Now we’ve found out some basic information about the how to compute variables, let’s
use it to transform our data. To transform the variable day1, and create a new variable
logday1, we execute this command:
dlf$logday1 <- log(dlf$day1)

This command creates a variable called logday1 in the dlf dataframe, which contains values
that are the natural log of the values in the variable day1.
For the day 2 hygiene scores there is a value of 0 in the original data, and there is no
logarithm of the value 0. To overcome this we should add a constant to our original
scores before we take the log of those scores. Any constant will do, provided that it
makes all of the scores greater than 0. In this case our lowest score is 0 in the data set so
we can simply add 1 to all of the scores and that will ensure that all scores are greater
than zero.
The advantage of adding 1 is that the logarithm of 1 is equal to 0, so people who scored
a zero before the transformation score a zero after the transformation. To do this transfor-
mation we would execute:
dlf$logday1 <- log(dlf$day1 + 1)

This command creates a variable called logday1 in the dlf dataframe, which contains values
that are the natural log of the values in the variable day1 after 1 has been added to them.

SELF-TEST

9 Have a go at creating similar variables logday2

and logday3 for the day2 and day3 variables. Plot
histograms of the transformed scores for all three
days.

[Link]. The square root transformation in R 2

To do a square root transformation, we run through the same process, by using a name
such as sqrtday1. Therefore, to create a variable called sqrtday1 that contains the square
root of the values in the variable day1, we would execute:
dlf$sqrtday1 <- sqrt(day1)

SELF-TEST

9 Repeat this process for day2 and day3 to create

variables called sqrtday2 and sqrtday3. Plot
histograms of the transformed scores for all three
days.

[Link]. The reciprocal transformation in R 2

To do a reciprocal transformation on the data from day 1, we don’t use a function, we use
an arithmetic expression: 1/variable. However, the day 2 data contain a zero value and if
198 D I S C O VE R I N G STAT I ST I C S US I N G R

we try to divide 1 by 0 then we’ll get an error message (you can’t divide by 0). As such
we need to add a constant to our variable just as we did for the log transformation. Any
constant will do, but 1 is a convenient number for these data. We could use a name such as
recday1, and to create this variable we would execute:
dlf$recday1 <- 1/(dlf$day1 + 1)

SELF-TEST

9 Repeat this process for day2 and day3. Plot

histograms of the transformed scores for all three
days.

[Link]. The ifelse() function in R 2

The ifelse() function is used to create a new variable, or change an old variable, depending
on some other values. This function takes the general form:
ifelse(a conditional argument, what happens if the argument is TRUE, what
happens if the argument if FALSE)

This function needs three arguments: a conditional argument to test, what to do if the test
is true, and what to do if the test is false. Let’s use the original data where there was an
outlier in the day1 hygiene score. We can detect this outlier because we know that the high-
est score possible on the scale was 4. Therefore, we could set our conditional argument to
be dlf$day1 > 4, which means we’re saying ‘if the value of day1 is greater than 4 then …’.
The rest of the function tells it what to do, for example, we might want to set it to missing
(NA) if the score is over 4, but keep it as the old score if the score is not over 4. In which
case we could execute this command:
dlf$day1NoOutlier <- ifelse(dlf$day1 > 4, NA, dlf$day1)

This command creates a new variable called day1NoOutlier which takes the value NA if
day1 is greater than 4, but is the value of day1 if day1 is less than 4:

If yes, then the

new variable is set
to NA (Missing)

dlf$day1NoOutlier <–ifelse(dlf$day1 > 5, NA, dlf$day1)

If no, then the new

Test is day1 variable is set to be
greater than 5? the value of the
old variable
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 199

R ’ s S o u l s ’ T i p 5 .4 Careful with missing data 3

If you have any missing data in your variables, you need to be careful when using functions such as rowMeans(),
to get the answer that you want. The problem is what you do when you have some missing values. Here’s a prob-
lem: I have 2 oranges and 3 apples. How many fruits do I have? Obviously, I have a total of 5 fruits.
You have 2 oranges, and we don’t know how many apples – this value is missing. How many fruits do you
have? We could say that you have 2. Or we could say that we don’t know: the answer is missing. If you add apples
and oranges in R, most functions will tell you that the answer is NA (unknown).
apples <- 2
oranges <- NA
apples + oranges

[1] NA

The rowSums and rowMeans functions will allow you to choose what to do with missing data, by using the [Link]
option, which asks ‘should missing values (na) be removed (rm)?’
To obtain the mean hygiene score across three days, removing anyone with any missing values, we would use:
dlf$meanHygiene <- rowMeans(cbind(dlf$day1, dlf$day2, dlf$day3))

But a lot of people would be missing. If we wanted to use everyone who had at least one score for the three days,
we would add [Link]=TRUE:
dlf$meanHygiene <- rowMeans(cbind(dlf$day1, dlf$day2, dlf$day3), [Link] = TRUE)

But what would we do if we had 100 days of hygiene scores? And if we didn’t mind if people were missing one or
two scores, but we didn’t want to calculate a mean for people who only had one score? Well, we’d use the [Link]()
function first, to count the number of missing variables.
dlf$daysMissing <- rowSums (cbind ([Link](dlf$day1),
[Link](dlf$day2),
[Link](dlf$day3)))

(It’s OK to break a command across rows like that, and sometimes it makes it easier to see that you didn’t make
a mistake.) Then we can use the ifelse() function to calculate values only for those people who have a score on
at least two days:
dlf$meanHygiene <- ifelse(dlf$daysMissing < 2, NA,
rowMeans(cbind( dlf$day1,
dlf$day2,
dlf$day3),
[Link]=TRUE))

Notice how I’ve used spacing so it’s clear which arguments go with which function? That makes it (slightly) easier
to avoid making mistakes.11

[Link]. The effect of transformations 2

Figure 5.9 shows the distributions for days 1 and 2 of the festival after the three different
transformations. Compare these to the untransformed distributions in Figure 5.2. Now,
you can see that all three transformations have cleaned up the hygiene scores for day 2:

11
It still took me three tries to get this right.
200 D I S C O VE R I N G STAT I ST I C S US I N G R

the positive skew is reduced (the square root transformation in particular has been useful).
However, because our hygiene scores on day 1 were more or less symmetrical to begin
with, they have now become slightly negatively skewed for the log and square root trans-
formation, and positively skewed for the reciprocal transformation!12 If we’re using scores

FIGURE 5.9 Day 1 of Download Day 2 of Download

Distributions of
2.0 1.4
the hygiene data
on day 1 and day 1.2
2 after various
1.5 1.0
transformations
Density

Density
0.8
1.0
0.6

0.4
0.5
0.2

0.0 0.0
0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5
Log Transformed Hygiene Score in Day 1 Log Transformed Hygiene Score in Day 2

1.5
1.5

1.0
Density

Density

1.0

0.5
0.5

0.0 0.0
0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5
Square Root of Hygiene Score on Day 1 Square Root of Hygiene Score in Day 2

5
2.5

4
2.0

3
Density
Density

1.5

2 1.0

1 0.5

0 0.0
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0
Reciprocal of Hygiene Score in Day 1 Reciprocal of Hygiene Score in Day 2

12
The reversal of the skew for the reciprocal transformation is because, as I mentioned earlier, the reciprocal has
the effect of reversing the scores.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 201

from day 2 alone then we could use the transformed scores; however, if we wanted to look
at the change in scores then we’d have to weigh up whether the benefits of the transforma-
tion for the day 2 scores outweigh the problems it creates in the day 1 scores – data analysis
can be frustrating sometimes!

5.8.4. When it all goes horribly wrong 3

It’s very easy to think that transformations are the answers to all of your broken assump-
tion prayers. However, as we have seen, there are reasons to think that transformations
are not necessarily a good idea (see Jane Superbrain Box 5.1), and even if you think that
they are they do not always solve the problem, and even when they do solve
the problem they often create different problems in the process. This happens
more frequently than you might imagine (messy data are the norm). What do I do if my
If you find yourself in the unenviable position of having irksome data then transformation
there are some other options available to you (other than sticking a big samu- doesn’t work?
rai sword through your head). The first is to use a test that does not rely on the
assumption of normally distributed data, and as you go through the various
chapters of this book I’ll point out these tests – there is also a whole chapter
dedicated to them later on.13 One thing that you will quickly discover about
non-parametric tests is that they have been developed for only a fairly limited
range of situations. So, happy days if you want to compare two means, but sad
and lonely days listening to Joy Division if you have a complex experimental
design.
A much more promising approach is to use robust methods (which I mentioned in Jane
Superbrain Box 5.1). These tests have developed as computers have got more sophisticated
(doing these tests without computers would be only marginally less painful than ripping
off your skin and diving into a bath of salt). How these tests work is beyond the scope of
this book (and my brain), but two simple concepts will give you the general idea. Some
of these procedures use a trimmed mean. A trimmed mean is simply a mean based on the
distribution of scores after some percentage of scores has been removed from each extreme
of the distribution. So, a 10% trimmed mean will remove 10% of scores from the top and
bottom before the mean is calculated. With trimmed means you have to specify the amount
of trimming that you want; for example, you must decide to trim 5%, 10% or perhaps even
20% of scores. A similar robust measure of location is an M-estimator, which differs from a
trimmed mean in that the amount of trimming is determined empirically. In other words,
rather than the researcher deciding before the analysis how much of the data to trim, an
M-estimator determines the optimal amount of trimming necessary to give a robust esti-
mate of, say, the mean. This has the obvious advantage that you never over- or under-trim
your data; however, the disadvantage is that it is not always possible to reach a solution. In
other words, robust tests based on M-estimators don’t always give you an answer.
We saw in Chapter 2 that the accuracy of the mean depends on a symmetrical distribu-
tion, but a trimmed mean (or M-estimator) produces accurate results even when the dis-
tribution is not symmetrical, because by trimming the ends of the distribution we remove
outliers and skew that bias the mean. Some robust methods work by taking advantage of
the properties of the trimmed mean and M-estimator.

13
For convenience a lot of textbooks refer to these tests as non-parametric tests or assumption-free tests and
stick them in a separate chapter. Actually neither of these terms are particularly accurate (none of these tests is
assumption-free) but in keeping with tradition I’ve put them in a chapter on their own (Chapter 15), ostracized
from their ‘parametric’ counterparts and feeling lonely.
202 D I S C O VE R I N G STAT I ST I C S US I N G R

The second general procedure is the bootstrap (Efron & Tibshirani, 1993). The idea of
the bootstrap is really very simple and elegant. The problem that we have is that we don’t
know the shape of the sampling distribution, but normality in our data allows us to infer
that the sampling distribution is normal (and hence we can know the probability of a par-
ticular test statistic occurring). Lack of normality prevents us from knowing the shape of
the sampling distribution unless we have big samples (but see Jane Superbrain Box 5.1).
Bootstrapping gets around this problem by estimating the properties of the sampling dis-
tribution from the sample data. In effect, the sample data are treated as a population from
which smaller samples (called bootstrap samples) are taken (putting the data back before
a new case is drawn). The statistic of interest (e.g., the mean) is calculated in each
sample, and by taking many samples the sampling distribution can be estimated (rather
like in Figure 2.7). The standard error of the statistic is estimated from the standard devia-
tion of this sampling distribution created from the bootstrap samples. From this standard
error, confidence intervals and significance tests can be computed. This is a very neat way
of getting around the problem of not knowing the shape of the sampling distribution. The
bootstrap can be used in conjunction with trimmed means and M-estimators. For a fairly
gentle introduction to the concept of bootstrapping see Wright, London, and Field (2011).
There are numerous robust tests based on trimmed means, bootstrapping and
M-estimators described by Rand Wilcox (Figure 5.10) in his definitive text (Wilcox, 2005).
He has also written functions in R to do these tests (which, when you consider the number
of tests in his book, is a feat worthy of anyone’s respect and admiration). We cover quite a
few of these tests in this book.
There are two ways to access these functions: from a package, and direct from Wilcox’s
website. The package version of the tests is called WRS (although it is what’s known as a
beta version, which means it is not complete).14 To access this package in R we need to
execute:
[Link]("WRS", repos="[Link]
library(WRS)

This is a standard install procedure, but note that we have to include repos=[Link]
[Link] because it is not a full package and this instruction tells R where to
find the package. This package is not always implemented in the most recent versions of R
(because it is only a beta) and it is not kept as up to date as Wilcox’s webpage, so although
we tend to refer to the package, to be consistent with the general ethos of downloading
packages, you should also consider sourcing the functions from Wilcox’s website. One
advantage of the website is that he keeps the functions very up to date. To source the func-
tions from his website, execute:
source("[Link]

This command uses the source() function to access the webpage where Wilcox stores
the functions (as a text file). Rallfun-v14 is the name of the file (short for ‘R all functions –
version 14’). Without wishing to state the obvious, you need to be connected to the Internet
for this command to work. Depending on this book’s shelf-life, it is possible that the name
of the file might change (most likely to Rallfun-v15 or Rallfun-v16), so if you get an error
try replacing the v14 at the end with v15 and so on. It’s also possible that Rand might move
his webpage ([Link] in which case Google him, locate the lat-
est Rallfun file and replace the URL in the source function above with the new one. Having
either loaded the package or sources the file from the web, you now have access to all of
the functions in Wilcox’s book.

14
Actually, all of the functions are there, but there is very little documentation about what they do, which is why
it is only at the ‘beta’ stage rather than being a full release.
CHAPTE R 5 E X PLOR I NG A SSU MPTI ONS 203

FIGURE 5.10
The absolute
legend that is Rand
Wilcox, who is the
man you almost
certainly ought to
thank if you want
to do a robust test
in R

What have I discovered about statistics? 1

‘You promised us swans,’ I hear you cry, ‘and all we got was normality this, homosome-
thingorother that, transform this, it’s all a waste of time that. Where were the bloody
swans?!’ Well, the Queen owns them all so I wasn’t allowed to have them. Nevertheless,
this chapter did negotiate Dante’s eighth circle of hell (Malebolge), where data of deliber-
ate and knowing evil dwell. That is, data that don’t conform to all of those pesky assump-
tions that make statistical tests work properly. We began by seeing what assumptions
need to be met for parametric tests to work, but we mainly focused on the assumptions
of normality and homogeneity of variance. To look for normality we rediscovered the
joys of frequency distributions, but also encountered some other graphs that tell us about
deviations from normality (Q-Q plots). We saw how we can use skew and kurtosis values
to assess normality and that there are statistical tests that we can use (the Shapiro–Wilk
test). While negotiating these evildoers, we discovered what homogeneity of variance is,
and how to test it with Levene’s test and Hartley’s Fmax. Finally, we discovered redemp-
tion for our data. We saw we can cure their sins, make them good, with transformations
(and on the way we discovered some of the uses of the by() function and the transforma-
tion functions). Sadly, we also saw that some data are destined always to be evil.
We also discovered that I had started to read. However, reading was not my true pas-
sion; it was music. One of my earliest memories is of listening to my dad’s rock and soul
records (back in the days of vinyl) while waiting for my older brother to come home
from school, so I must have been about 3 at the time. The first record I asked my parents
to buy me was ‘Take on the World’ by Judas Priest, which I’d heard on Top of the Pops (a
now defunct UK TV show) and liked. This record came out in 1978 when I was 5. Some
people think that this sort of music corrupts young minds. Let’s see if it did …
204 D I S C O VE R I N G STAT I ST I C S US I N G R

R packages used in this chapter

car psych
ggplot2 Rcmdr
pastecs

R functions used in this chapter

abs() qplot()
by() rowMeans()
cbind() rowSums()
describe() round()
dnorm() [Link]()
ifelse() source()
[Link]() sqrt()
leveneTest() [Link]()
log() stat_function()
log10() tapply()

Key terms that I’ve discovered

Bootstrap Normally distributed data
Hartley’s Fmax Parametric test
Heterogeneity of variance Q-Q plot
Homogeneity of variance Quantile
Independence Robust test
Interval data Shapiro–Wilk test
Levene’s test Transformation
Log Trimmed mean
M-estimator Variance ratio

Smart Alex’s tasks

G Task 1: Using the [Link] data from Chapter 4, check the assumptions of
normality and homogeneity of variance for the two films (ignore gender): are the
assumptions met? 1
G Task 2: Remember that the numeracy scores were positively skewed in the RExam.
dat data (see Figure 5.5)? Transform these data using one of the transformations
described in this chapter: do the data become normal? 2
Answers can be found on the companion website.

Further reading
Tabachnick, B. G., & Fidell, L. S. (2007). Using multivariate statistics (5th ed.). Boston: Allyn &
Bacon. (Chapter 4 is the definitive guide to screening data!)
Wilcox, R. R. (2005). Introduction to robust estimation and hypothesis testing (2nd ed.). Burlington,
MA: Elsevier. (Quite technical, but this is the definitive book on robust methods.)
Wright, D. B., London, K., & Field, A. P. (2011). Using bootstrap estimation and the plug-in principle
for clinical psychology data. Journal of Experimental Psychopathology, 2(2), 252–270. (A fairly
gentle introduction to bootstrapping in R.)