Tag Archives: ggplot

p-values and null hypotheses

Almost all statistical tests that archaeologists and historians use are based on the concept of Null Hypothesis Significance Testing or NHST. If someone told you that you should be using some stats in this paper this is probably what they mean for it. The comment is usually accompanied by a reference to a thick statistics manual where you repeatedly find a concept called p-value and other obscure terms.

A good understanding of p-values and NHST is essential for the current scientist; otherwise you won’t understand the results of a large percentage of papers you should be reading. The approach itself is rather simple, but current teaching of statistics bury it into a myriad of complex tests that you are supposed to learn in an almost ritualistic way, thus making everything more complicated than it really is (T-test, Chi-square test, Kruskal-Wallis test…seriously, who chooses these names??).

Let’s take a step back and see what NHST is really about.

the concept of the null hypothesis

Imagine that you manage a lab and you need to test if a newly developed antibiotic drug kills bacteria. The trial seems straightforward: you put some of the new drug into a Petri dish with the bacteria and you wait to see if it works (also, apologies for my ignorance on lab trials to any lab technician that could be listening…). Your working hypothesis is that the drug is effective and kills bacteria, so what are the possible outcomes?

  1. it works! bacteria are dead and the new drug is effective
  2. it doesn’t work and bacteria are still alive.

In the first case we would accept the working hypothesis while in the second one we would reject it. However, it is not as simple as it seems because you could have lots of grey results between the 2 extremes: some bacteria may still be alive, bacteria could be dying for different causes than the drug, or even by pure chance.

Someone came up with a simpler approach to this problem that can be summed up as follows: what are the chances that you get the same result by pure chance? If these chances (known as the p-value of the test) are low enough then we could reject the idea of pure chance, thus confirming that in the trial something was killing bacteria. This is NHST: if we can reject the null hypothesis of pure chance with a good degree of significance then the alternate hypothesis of something going on is selected as a better explanation of the evidence.

This degree of significance is in reality defined as a threshold for the p-value; if the chances are typically lower than 5% (p-value < 0.05) then we will reject the null hypothesis.

when to apply and when not to apply NHST

This is the type of scenarios where NHST was originally developed because in a controlled environment you can only design experiments with 2 results (it works or it doesn't work). However, this binary outcome is not necessarily found in other disciplines and specially it is problematic in archaeology.

Imagine that you are studying settlement patterns in late prehistoric Europe and you want to identify the geographical properties of these settlements. You would test if the gradient of terrain or slope is relevant here by testing the null hypothesis that it is not relevant.

Would you?

Really? If you think about it this null hypothesis is meaningless because we know that you won't be living at a slope of 40 degrees and hanging from a rope. You will obviously reject the null hypothesis! Moreover, this result doesn't mean that slope was particularly relevant during late prehistory because slope is always relevant to human settlements for obvious functional reasons: we don't have wings to fly.

The choice of a meaningless null hypothesis is a classic example of bad statistics; you should always evaluate if the approach you are applying makes sense or you need to tackle the research question with a different method.

What you could certainly do is to assess if the geographical features exhibited by settlements changed over time. Let's see how NHST works in practice.

late prehistory in the NE Iberian Peninsula

The transition from Late Bronze Age (LBA) to Early Iron Age (EIA) in current Catalonia has been extensively studied during the last decades because of its interesting dynamics. Population seemed to shift from high-elevation pastures towards low-level agricultural regions; at the same time you find an increasing amount of artefacts brought to this area though long-range trade. Both processes seemed to increase the level of regional connectivity so…can we test this?

load the dataset

We will use here a dataset collected by Maria Yubero-Gómez who studied the period for her PhD and was published in:
* Yubero-Gómez, M., Rubio-Campillo, X., López-Cachero, F.J. and Esteve-Gràcia, X., 2015. Mapping changes in late prehistoric landscapes: a case study in the Northeastern Iberian Peninsula. Journal of Anthropological Archaeology, 40, pp.123-134. https://doi.org/10.1016/j.jaa.2015.07.002

To load the simplified dataset into R:

[generic linenumbers="False"] sites <- read.csv("https://git.io/v9XNs")

You have here 5 variables:

  • the name of the archaeological site
  • distance to the closest natural route (i.e. main routes you would take to cross the region)
  • slope of the site in degrees
  • height in meters
  • period of the site (LBA or EIA; if a site was active during both periods it will have 2 entries)

perform Exploratory Data Analysis

Let’s first focus on height. The first thing you need to do in any quantitative analysis is to explore the dataset with some plots:
[generic linenumbers=”False”] library(ggplot2)
ggplot(sites, aes(x=height)) + geom_density(col=”grey70″, aes(fill=period, group=period), alpha=0.9) + facet_wrap(~period) + scale_fill_manual(values=c(“indianred2”, “skyblue3″)) + theme_bw() + theme(legend.position=”none”)

Empirical density estimates of height (m.) per period

These are the density estimates for the evidence we collected. The distribution looks different for both periods, don’t they? Let’s take a look at the other 2 variables (slope and distance to natural routes):

[generic linenumbers=”False”] ggplot(sites, aes(x=slope)) + geom_density(col=”grey70″, aes(fill=period, group=period), alpha=0.9) + facet_wrap(~period) + scale_fill_manual(values=c(“indianred2”, “skyblue3″)) + theme_bw() + theme(legend.position=”none”)

Empirical density estimates of slope (degrees) per period

and the last one:

[generic linenumbers=”False”] ggplot(sites, aes(x=dist)) + geom_density(col=”grey70″, aes(fill=period, group=period), alpha=0.9) + facet_wrap(~period) + scale_fill_manual(values=c(“indianred2”, “skyblue3″)) + theme_bw() + theme(legend.position=”none”)

Empirical density estimates of distance to natural route (m,) per period

They don’t look as different as the slope. Specifically, the distribution of distances seems almost identical! However, the human eye is too good at identifying patterns so we may be finding similarities where there are not; we need to statistically test if they are different.

test the null hypothesis

Our null hypothesis can be defined as follows: the distribution of slopes is the same for both periods. If we cannot reject this null hypothesis then it means that the differences we see in the plots are caused by pure chance; on the other hand, if we can reject it then we can say that there was a change on the relevance of slope over the transition LBA-EIA.

The perfect test to evaluate this null hypothesis is known as Kolmogorov-Smirnov test. The input of the test are the 2 distributions and the output is the probability that they were created by the same process. To summarize:

  • If p-value > 0.05 then it means that there was a change between the 2 periods
  • If p-value < 0.05 then it means that there was no significant change from LBA to EIA.

Let's create 2 new data frames (one for each period) and apply the Kolmogorov-Smirnov test to the 3 variables:

[generic linenumbers="False"] lbaSites <- subset(sites, period=="LBA")
eiaSites <- subset(sites, period=="EIA")

ks.test(eiaSites$slope, lbaSites$slope)
ks.test(eiaSites$height, lbaSites$height)
ks.test(eiaSites$dist, lbaSites$dist)

The results for slope should be something like:
[generic linenumbers="False"] Two-sample Kolmogorov-Smirnov test

data: eiaSites$slope and lbaSites$slope
D = 0.27647, p-value = 0.006207
alternative hypothesis: two-sided

Ignore the warning as it is not relevant right now. The p-value you see is extremely low so we can safely reject the null hypothesis that the distribution of slope does not vary over the transition. Nice! The same can be said for elevation as the value is also low; however, the p-value of distances is really, really (really) high. The chances are over 90% so we cannot reject in this case the null hypothesis that the typical distance to natural route changed from LBA to EIA.

final remarks

Some last comments on this analysis:

  • Attention! A low p-value does not mean that the difference between distributions is large. You can have a very significant yet small difference between distributions and it would give you a very low p-value (e.g. all sites at an elevation of 200m against all sites at an elevation of 220m.).
  • The choice of 0.05 as the significance threshold is a convention. Yes, we could choose 0.041 or 0.1 but you will see 0.05 as the limit in 90% of papers (also, if you see a paper using 0.1…don't trust the results…seriously)
  • there are like 2023 statistical tests besides K-S for any type of variable: categorical, discrete, multivariate…in any case I think KS is one of the most useful for archaeologists

The approaches based on Null-Hypothesis Significance Testing have a lot of issues and there is a current debate on science about the appropriateness of the method. Alternatives such as Bayesian inference and model selection are becoming increasingly popular. However, NHST is still the standard statistical approach in archaeology so I hope this post was useful to understand the philosophy behind the method.

Loading and visualizing Shapefiles in R

Geographical Information Systems such as QGIS are common in the typical archaeologist’s toolbox. However, the complexities of our datasets mean that we need additional tools to identify patterns in time and space.

Imagine that you are working with a set of points in a vector format such as a Shapefile. It can be a set of settlements, C14 dates or shipwrecks…it doesn’t matter; you will want to combine spatial analysis with other methods. This is particularly true when you are trying to get a feeling of the dataset by performing Exploratory Data Analysis, so you need to compile summary statistics and create some meaningful visualisations.

My personal approach is to combine QGIS with R and exploit the awesome ggplot package to get some insight into the dataset beyond spatial patterns. R has some nice spatial functionality, but the cool thing is to integrate these methods with other plots to get a general picture of the case study you are studying.

Example: Aircraft crashes in the Orkney islands

The Orkney islands were a key location for the British war effort during the Two World Wars because it hosted Scapa Flow: the main naval base of the Royal Navy. For this reason it saw an intense aerial activity of all sorts: squadrons defending the base, German bombers attacking it, aircraft carrier operations…you name it. Inevitably this activity generated aircraft crashes due to accidents and combat and these events have been compiled by different initiatives such as the Project Adair-Whitaker or the ARGOS group.

I created a nice Shapefile of aircraft crash sites based on the information provided by Canmore. How can I explore the dataset beyond the spatial dynamics? How many German aircrafts were shot down? What types and squadrons sustained more losses? Let’s take a look!

Download Data

The first thing to do is to load the Shapefile. This format requires of several files so I created a zip file with all the set. First of all we need to download it to a temporary directory and extract its contents:

[generic linenumbers=”False”] tmp <- tempdir()
url <- "https://github.com/xrubio/pastByNumbers/raw/master/data/aircraft_orkney.zip"
file <- basename(url)
download.file(url, file)
unzip(file, exdir = tmp)

Load the spatial data

We will use the rgdal library to load the contents of the Shapefile into R:
[generic linenumbers=”False”] library(rgdal)
aircraftShp <- readOGR(dsn = tmp, layer="aircraft_orkney")

aircraftShp is a variable of type SpatialPointsDataFrame from package sp. It is a container of some interesting information such as the bounding box delimiting our spatial entities or the Coordinate Reference System that it is being used. Ideally we would like to get the attributes of the spatial entities as a data frame so we can use it with most R packages:

[generic linenumbers=”False”] aircraft <- as.data.frame(aircraftShp)

Maps and planes

First of all let’s plot the spatial coordinates of the aircraft crash sites against a nice background. I created a rectangle location based on the bounding box of aircraftShp:

[generic linenumbers=”False”] aircraftShp@bbox
location <- c(-3.85, 58.7, -1.95, 59.4)

We can use ggmap to download the background:
[generic linenumbers=”False”] library(ggplot2)
bgMap <- get_map(location=location, source="stamen", maptype="watercolor")
ggmap(bgMap) + geom_point(data=aircraft, aes(x=coords.x1, y=coords.x2))

Aircraft crash sites around Scapa Flow

Exploring other dimensions

This is ok but it is something we can already do with any GIS so…what’s the deal? The thing is that aircraft is an R Data Frame so we can plot any variable using a graphic library, something that R does way better than any GIS. For example, let’s take a look at the temporal dimension:

[generic linenumbers=”False”] ggplot(aircraft, aes(x=year, fill=force)) + geom_histogram(binwidth=1) + facet_wrap(~war, scales=”free_x”)

Losses by force and year

We can also create a table of planes used by the different British squadrons that suffered losses:

[generic linenumbers=”False”] ggplot(na.omit(aircraft), aes(y=type, col=force, x=as.factor(sqdn))) + geom_raster()

Aircraft types per squadron

Combining everything

We can finally create an infographic combining space, time and other stuff with the library gridExtra:

First, we create nice version of the plots we just saw:

[generic linenumbers=”False”] g1 <- ggmap(bgMap, extent="panel", darken = c(.4,"white")) + geom_point(data=aircraft, aes(x=coords.x1, y=coords.x2, col=war), size=2) + geom_text(data=aircraft, aes(x=coords.x1, y=coords.x2,label=type), family = "Trebuchet MS", color="grey40", check_overlap=T, nudge_y=0.01) + theme(legend.position="bottom") + ggtitle("aircraft losses") + xlab("") + ylab("") + scale_color_manual(values=c("indianred2", "goldenrod1"))

g2 <- ggplot(na.omit(aircraft), aes(y=type, col=force, x=as.factor(sqdn))) + geom_point(size=3) + theme_bw() + theme(legend.position="bottom") + xlab("squadron") + ylab("") + ggtitle("Britain – aircraft types per squadron")

g3 <- ggplot(aircraft, aes(x=year, fill=force)) + geom_histogram(binwidth=1, col="black") + facet_wrap(~war, scales="free_x", ncol=1) + theme_bw() + theme(legend.position="bottom") + ggtitle("losses per year") + xlab("") + ylab("")

And we can then combine the 3 plots into a nice visualization:
[generic linenumbers=”False”] library(gridExtra)
grid.arrange(arrangeGrob(arrangeGrob(g1,g2, heights=c(2/3,1/3), ncol=1), g3, widths=c(2/3, 1/3), ncol=2))

Aicraft Losses in Orkney islands during the two World Wars


These visualizations allow us to identify several interesting things going on:

  • The Second World War had way more accidents than the First one.
  • It would seem that Luftwaffe stopped flying over Scapa Flow after the first years of WW2
  • The Fleet Air Arm had a very large increase on losses after 1941.
  • Most squadrons only flew 1 type of plane being the exception the 771 Naval Air Squadron. It makes total sense given the diversity of missions flown by this squadron.