Uncertainty analysis

Uncertainty visualization

Cartographers, and more recently those working in Information Visualization, have been concerned for some time with visual variables, or adopting Munzner (2014)’s vocabulary visual channels, that might be used to represent uncertainty information. Figure 1 is adapted from Kinkeldey, MacEachren, and Schiewe (2014) and displays visual variables that could be used to encode uncertainty information. Ideally, visual variables for representing uncertainty should be intuitive, logically related to notions of precision and accuracy, whilst also allowing sufficient discriminative power when deployed in data dense visualizations.

Kinkeldey, MacEachren, and Schiewe (2014) provides an overview of empirical research (controlled experiments) that explore the effectiveness of proposed visual variables against these criteria. As intuitive signifiers of uncertainty, or lack of precision, fuzziness and location have been shown to work well. Slightly less intuitive, but nevertheless successful in terms of discrimination are size, transparency and colour value. Sketchiness is another intuitive signifier, proposed in Boukhelifa et al. (2012) and Wood et al. (2012). As with many visual variables, sketchiness is probably best considered as an ordinal visual variable to the extent that there is a limited range of sketchiness levels that can be discriminated. An additional feature of sketchiness is that, different from the other visual variables, it is related to informality – this may be desirable in certain contexts, less so in others (see Wood et al. 2012).

Figure 1: Visual variables that can be used to represent levels of uncertainty information.

From within the Uncertainty Visualization literature, a key maxim that consistently appears for successful uncertainty visualization is that:

Things that are not precise should not be encoded with symbols that look precise.

The go-to example, which appears in most primers on uncertainty visualization, is the National Hurricane Center’s Cone of Uncertainty graphic (Figure 2) used to represent major storm events. The cone starts at the storm’s current location and spreads out to represent the projected path of the storm as determined by National Hurricane Center’s modelling. The main problem is that the cone implies that the storm is growing as we move away from its current location, when in fact this is not the case. Instead there is more uncertainty in the areas that could be affected by the storm the further away those areas are from the storm’s current location. The second problem is that the cone uses strong lines that imply precision. The temptation is to think that anything contained by the cone is unsafe and anything outside of it is safe. This is of course not what is suggested by the model; rather that areas beyond the cone fall outside some chosen threshold probability.

In Jo Wood’s redesign, colour value is used to represent four binned categories of storm probability suggested by the model. In representing the cone of uncertainty on a map, and incorporating geographic context – administrative boundaries – it is also natural to make quite specific and binary judgements about response based on this context. For example, if the cone is close to but not overlapping a state boundary, a decision might be made not to mount a response/prepare for the hurricane. Jo’s approach is to use curve schematisation (see Van Goethem et al. 2014) to symbolise states using a single line. This provides context but in a way that may discourage binary thinking; precise inferences of location are not possible as the area and state borders are very obviously not exact.

Figure 2: Cone of Uncertainty produced by National Hurricane Center and example re-design by Jo Wood.

Frequency framing

The visual variables in Figure 1 could be used to represent different categories of uncertainty information, not just those associated with parameters – for example, locational or temporal uncertainty at data collection. The rest of the session addresses techniques for quantifying, and charting idioms (Munzner 2014) for representing, parameter uncertainty: estimated injury severity rates for pedestrian-vehicle crashes in our STATS19 road crash dataset.

Often parameters are represented as probabilities, or relative frequencies – ratios or percentages describing the probability of some event happening. It is notoriously difficult to develop intuition around these sorts of relative frequencies. In the STATS19 dataset, we might wish to compare the injury severity rate of pedestrain-vehicle road crashes – the proportion of all crashes that resulted in a serious injury or fatality (KSI) – taking place between two local authority areas, say Bristol and Sheffield. There is in fact quite a difference in the injury severity rate between these two local authority areas in 2019: 15% for Bristol (35 out of 228 reported crashes were KSI) versus 50% for Sheffield (124 out of 248 reported crashes were KSI).

This feels like quite a large difference, but it is difficult to imagine or experience these differences in probabilities when written down or encoded visually using say bar length. Icon arrays are used extensively in health communication and have been demonstrated to be effective at communicating probabilities of event outcomes. They offload the thinking that happens when comparing ratios – the internal weighing up of numerators and denominators. In the example in Figure 3, icon arrays are used to compare the two injury severity rates for Bristol and Sheffield. Each crash is a square and crashes are coloured according to whether they resulted in a serious injury or fatality (KSI, dark red) or slight injury (light red).

Figure 3: Icon array displaying injury severity rates for Pedestrian-Vehicle crashes.

There are compelling examples of icon arrays being used in data journalism, most obviously to support communication of probabilities in political polling. You might remember that at the time of 2016 US Presidential Election, there was much criticism levelled at pollsters, even the excellent FiveThirtyEight, for not correctly calling the result. Huffpost gave Trump a 2% chance of winning the election, The New York Times 15% and FiveThirtyEight 28%. Clearly, the Huffpost estimate was really quite off, but thinking about FiveThirtyEight’s prediction, how surprised should we be if an outcome does in fact occur, which is predicted to happen with a probability of almost a third?

The risk theatre (Figure 4) is a variant of an icon array. In this case it represents polling probabilities as seats of a theatre – a dark seat represents a Trump victory. If you imagine buying a theatre ticket and being randomly allocated to a seat, how confident would you be about not sitting in a “Trump” seat in the FiveThirtyEight image? The distribution of dark seats suggests that the 28% risk of a Trump victory according to the model is not negligible.

Figure 4: Risk theatre of different election eve forecasts by Justin Gross as it appears in Washington Post

Quantifying uncertainty in frequencies

In the icon arrays above I made little of the fact that the sample size varies between the two recorded crash rates. Partly this was because the differences were reasonably small. When looking at injury severity rates across all local authorities in 2019, however, there is substantial variation in the rates and the sample size. Bromsgrove has a very low injury severity rate based on a small sample size (4%, or one out of 27 crashes resulting in KSI); Cotswold has a very high injury severity rate based on a small sample size (75%, or 14 out of 19 crashes resulting in KSI). With some prior knowledge of these areas, one might expect the difference in KSI rates to be in this direction, but would we expect the difference to be of this order of magnitude? Just eight more KSIs recorded in Bromsgrove makes its KSI rate equivalent to that of Bristol’s.

Although STATS19 is a population dataset to the extent that it contains data on every crash recorded by the Police, it makes sense that the more data on which our KSI rates are based, the more certainty we have in them being reliable estimates of injury severity – ones that might be used to predict injury severity in future years. So we can treat our observed injury rates as being derived from samples of an (unobtainable) population. Our calculated injury severity rates are parameters that try to represent, or estimate, this population.

Although this formulation might seem unnecessary, from here we can apply some statistical concepts to quantify uncertainty around our parameter estimates. We assume:

1. The variable of interest, KSI rate, has an unobtainable population mean and standard deviation.
2. That a sample will consist of a set of observations from this unobtainable population and that these samples could vary in size.
3. From any sample we can calculate a mean and standard deviation, which will differ from the population mean and standard deviation.
4. That we can derive a sampling distribution and generate an array of estimates that could be obtained from repeating the sampling process many, many times.
5. This range of the sampling distribution could then be used to quantify how precise are the estimates. Generally the larger the sampling distribution, the more precise – the less uncertain – the estimate.

In session 6 and session 7, we used Confidence Intervals to estimate uncertainty around regression coefficients. These described the range of the sampling distribution – the range of values that coefficients estimated from a large number of resamples could take. From early stats courses, you might have learnt how these can be calculated using statistical theory. Better still, we can derive these empirically via bootstrapping. A bootstrap resample involves taking a random sample with replacement from the original sample and of the same size as the original sample. From this resample, a parameter estimate can be derived, in this case the KSI rate. And this process can be repeated many times to generate an empirical sampling distribution for the parameter. The standard error can be calculated from the standard deviation of the sampling distribution. There are several reasons why bootstrapping is a really useful procedure: it can be applied to almost any sample statistic, makes no distributional assumptions and can work on quite complicated sampling designs.

Presented in Figure 5 are KSI rates with error bars used to display 95% Confidence Intervals generated from a bootstrap procedure in which 1000 resamples were taken with replacement. Upper and lower limits were lifted from .025 and .975 percentile positions of the bootstrap sampling distribution. Assuming that the observed data are drawn from a wider (unobtainable) population, the 95% Confidence Intervals demonstrate that whilst Cotswold recorded a very large KSI rate, sampling variation means that this figure could be much lower (or higher), whereas for Bristol and Sheffield, where our KSI rate is derived from more data, the range of plausible values that the KSI rate might take due to sampling variation is much smaller – there is less uncertainty associated with their KSI rates.

Figure 5: KSI rates for pedestrian-vehicle crashes in selected local authorities with 95% CIs (derived from 1000 resample bootstraps).

Visualizing uncertainty in frequencies

Error bars are a space-efficient way of conveying parameter uncertainty. However, remembering our maxim for uncertainty visualization – that things that are not precise should not be encoded with symbols that look precise – they do have problems. The hard borders can lead to a binary/categorical thinking (see Correll and Gleicher 2014). Certain values within a Confidence Interval are more probable than others and so we should endeavour to use a visual encoding that reflects this. Matt Kay’s excellent ggdist package extends ggplot2 with a range of chart idioms for representing these sorts of intervals. In Figure 6 I have replaced the error bars with half eye plots and interval bars, which give greater visual saliency to parameter estimates that are more likely.

Figure 6: KSI rates for pedestrian-vehicle crashes in selected local authorities with uncertainty estimates.

The STATS19 dataset is released annually, and given the wide uncertainty bands for certain local authorities, it might be instructive to explore the stability of local authority KSI rates year-on-year. In Figure 5 these KSI rates are represented with the bold line and the faint lines are superimposed bootstrap resamples. The overall line probably most clearly demonstrates volatility in the KSI rates for Cotswold and Bromsgrove due to small sample sizes. The observed increase in KSI rates for Sheffield since 2015 does appear to be a genuine one, although may also be affected by uncertainty around data collection – changes to how police record injury severity.

Figure 7: Year-on-yar KSI rates for pedestrian-vehicle crashes in selected local authorities with bootstrap resamples superimposed.

The superimposed lines in the figure above are a form of ensemble visualization. An alternative approach might have been to animate over the bootstrap resamples to generate a Hypoethetical Outcome Plot (HOP) (Hullman, Resnick, and Adar 2015). HOPs convey a sense of uncertainty by animating over random draws of a distribution. As there is no single outcome to anchor to (although I have actually included one below), HOPs force viewers to account for uncertainty, recognising that less probable outcomes may be possible – essentially to think distributionally (e.g. Figure 8).

Figure 8: HOP of year-on-yar KSI rates for pedestrian-vehicle crashes in selected local authorities.

Multiple comparisons

In road safety monitoring, a common ambition is to compare crash rates across local authorities on a map. We might represent injury severity rates as Risk Ratios (RR) comparing the observed injury severity in each local authority to a benchmark, say the injury severity rate we would expect to see nationally. RRs are an intuitive measure of effect size: \(RRs >1\) indicate that the injury severity rate is greater than the national average; \(RRs <1\) that it is less than the national average. As they are a ratio of ratios, and therefore agnostic to sample size, RRs can nevertheless be unreliable. Two ratios might be compared that have very different sample sizes and no compensation is made for the one that contains more data. We know this to be a problem in our dataset and so we will want to generate confidence intervals for our RRs.

In Figure 9, RRs for each local authority is represented with a | icon, which is then encoded according RR values by continuously varying icon angles: \(RRs >1\) are angled to the right /, \(RRs <1\) to the left \. You might remember this approach being used to represent electoral Swing in session 3.

Figure 9: Risk Ratios comparing local authority KSI rates to the GB average in 2019. Significance based on 90% bootstrap CI.

Whilst this is likely successful at conveying these effect sizes, we will want to de-emphasise those for which there is less certainty. One way of doing this would be to vary the lightness of the icons according to estimated standard errors. In Epidemiology and other domains, more common practice is to use null hypothesis testing (NHST). Here we set up an assumption (null hypothesis) that the difference in the injury severity rate for any local authority and the national average equals zero (\(RR = 1.0\)). Our bootstrap Confidence Interval is then used to derive probabilities that random sampling generate data similar to what was observed if this null hypothesis were true. The confidence interval can be used to derive this probability. If a 95% confidence interval excludes the expected parameter value under the null hypothesis (i.e. \(RR \not= 1.0\)), then the null hypothesis is rejected and the observed RR is deemed statistically significant. In Figure 9 RRs \(>1\) that are statistically significant coloured red /, RRs \(<1\) that are statistically significant coloured red are coloured blue \ and those that are not statistically significant are coloured grey (/ | \).

A problem with this approach is one familiar to Statisticians, but that is rarely addressed in visual data analyses: the multiple comparison problem. Whenever a NHST is made, there is a chance that the result observed is in fact a false alarm. In the plot above which uses a 90% confidence level, the “false positive rate” is expected to be 10% or 1/10. When many tests are considered simultaneously, as in Figure 9, the number of these false alarms begins to accumulate. There are corrections that can be applied for this: test statistics can be adjusted and made more conservative. But these corrections have consequences. Too severe a correction can result in statistical tests that are underpowered and result in an elevated false negative rate, where a statistical test fails to detect an effect that truly exists. See Brunsdon and Charlton (2011) for an interesting discussion of this in the context of mapping crime rates.

So there is no canonical solution to this multiple comparison problem, which I see a lot in Geography where health and other outcomes are mapped visually. Presenting the RRs in their spatial context, and providing full information around RRs that are not significant, helps with making uncertainty type judgements. For example, we can attach more certainty to RRs that are labelled statistically significant and whose direction is consistent with its neighbours than those that are exceptional from their neighbours. But this is certainly an interesting area for those working in visual data analysis.

Import

The template file lists the required packages – tidyverse, sf, tidymodels (for working with the bootstraps), ggdist and distributional for generating plots of parameter uncertainty and gganimate for the hypothetical outcome plot. You may not have saved the dataset of STATS19 pedestrian-vehicle road crashes generated in session 4 (ped_veh) and so I have placed this on a repo for download. Code for loading these data into your session is in the 08-template.Rmd file.

Plot icon arrays

Figure 10: Icon array displaying injury severity rates for Pedestrian-Vehicle crashes.

I quickly explored extensions to ggplot2 for generating icon arrays, but found no stand-out package. They can be generated reasonably easily using geom_tile() and some standard data generation functions. In the example in Figure 3 the icon arrays varied with sample size of the two Local Authorities being compared. A more generalisable approach would be for icon arrays to be regularly sized and used to represent probabilities of, in this case injury severity, if a stated rate of crashes occurred (e.g. Figure 10).

We need to generate the array data: a data frame of array locations (candidate crashes) with values representing whether the crash is slight or KSI depending on the observed KSI rate. In the code below, we set up a 10x10 grid of row and column locations and populate these with values for the selected local authorities (Oxford and Fareham) using base R’s sample() function.

array_data <- tibble(
  row=rep(1:10, times=1, each=10),
  col=rep(1:10, times=10, each=1),
  Oxford=
    sample(
    c(rep(TRUE, times=1, each=17), rep(FALSE, times=1, each=83)),
    size=100, replace=FALSE),
  Fareham=
    sample(
    c(rep(TRUE, times=1, each=41), rep(FALSE, times=1, each=59)),
    size=100, replace=FALSE)
)

Plot code:

array_data %>%
  pivot_longer(cols=c(Oxford,Fareham), names_to="la", values_to="is_ksi") %>%
  ggplot(aes(x=row,y=col, fill=is_ksi)) +
  geom_tile(colour="#ffffff", size=1) +
  scale_fill_manual(values=c("#fee0d2","#de2d26"), guide=FALSE)+
  facet_wrap(~la)

Plot specification:

1. Data: The array data, with pivot_longer() so that we can facet by local authority.
2. Encoding: x- and y-position according to the array locations and filled on whether the sampled crash is KSI or slight.
3. Marks: geom_tile() for drawing square icons.
4. Scale: scale_fill_manual() is supplied with values that are dark (KSI) and light (slight) red.
5. Facets: facet_wrap() for faceting on local authority.
6. Setting: Tiles are given large, white borders (geom_tile(colour="#ffffff", size=1)).

Generate bootstrap estimates of parameter uncertainty

The code for generating bootstrap resamples, stored in rate_boots, initially looks formidable. It is a template that is nevertheless quite generalisable, and so once learnt can be extended and applied to suit different use cases.

rate_boots <- ped_veh %>%
  mutate(
    is_ksi=accident_severity!="Slight",
    year=lubridate::year(date)
  ) %>%
  filter(year==2019,
         local_authority_district %in% c("Bristol, City of", "Sheffield", "Bromsgrove", "Cotswold")
         ) %>%
  select(local_authority_district, is_ksi) %>%
  nest(-local_authority_district) %>%
  mutate(la_boot=map(data, bootstraps, times=1000, apparent=TRUE)) %>%
  select(-data) %>%
  unnest(la_boot) %>%
  mutate(
    is_ksi=map(splits, ~ analysis(.) %>% pull(is_ksi)),
    ksi_rate=map_dbl(is_ksi, ~mean(.x)),
    sample_size=map_dbl(is_ksi, ~length(.x))
  ) %>%
  select(-c(splits, is_ksi))

Code description:

1. Setup: The first mutate is straightforward – we generate a binary is_ksi variable identifying whether the crash resulted in a KSI and extract the year from the crash date. We then filter on crashes that took place in 2019 and the four comparator local authorities. In order to generate bootstrap resamples for each local authority, we need to nest() on local authority. You will remember that nest() creates a special type of column (a list-column) in which the values of the column is a list of data frames – in this case the crash data for each local authority. So running the code up to and including the nest(), you will notice that a data frame is returned which contains four rows, corresponding to the filtered local authorities, and a list-column called data, each element of which contains a data frame of varying dimensions (lengths), according to the number of crashes recorded in each local authority.
2. Generate bootstraps resamples: In the mutate() that follows purrr’s map() function is used to iterate over the list of datasets and the bootstraps() function is used to generate 1000 bootstrap resamples for each nested dataset. The new column la_boot is a list-column this time containing a list of bootstrap datasets.
3. Calculate sample estimates: We unnest() the la_boot column to return a dataset with a row for each boostrap resample and a list-column named splits which contains the bootstrap data. Again we map() over each element of splits to calculate the ksi_rate for each of the bootstrap datasets. The first call to map() extracts the is_ksi variable; the second is just a convenient way of calculating a rate from this (remembering that is_ksi is a binary variable); the third collects the sample size for each of the bootstraps which of course is the same as the number of crashes recorded for each local authority.

Plot parameter estimates with uncertainty information

With ggdist, the code for generating the plots of KSI rates with estimates of parameter uncertainty is straightforward and very similar to the error bar plots in the previous two sessions.

Plot code:

rate_boots %>%
  group_by(local_authority_district) %>%
  mutate(std.error=sd(ksi_rate)) %>%
  filter(id=="Apparent") %>%
  ggplot(aes(x=reorder(local_authority_district, ksi_rate), y=ksi_rate)) +
  stat_dist_gradientinterval(
    aes(dist = dist_normal(mu=ksi_rate, sigma=std.error)),
    point_size = 1.5
    ) +
  coord_flip()

Plot specification:

1. Data: The rate_boots data frame is grouped by local authority and in the mutate() we calculate an estimate of bootstrap standard error, the standard deviation of the sampling distribution, and filter all rows where id=="Apparent" – this contains the KSI rate for the observed (unsampled) data.
2. Encoding: x- position varies according to local authority and y-position according to KSI rate. The estimated KSI rate and bootsrtap standard error is also passed to stat_dist_gradientinterval(), the ggdist function for producing gradient plots.
3. Marks: stat_dist_gradientinterval() for drawing the gradients and point estimates.
4. Setting: coord_flip() for easy reading of local authority names.

Ensemble plots and hypothetical outcome plots

To generate bootstrap resamples on local authority and year, necessary for the year-on-year analysis, we can use the same template for rate_boots, the only difference is that we select() and nest() on the year as well as local_authority_district column. You may also not wish to resample 1000 times for each year and local authority as this tends to slow down the computation (for obvious reasons).

rate_boots_temporal <- ped_veh %>%
  ...
  ... %>%
  select(local_authority_district, is_ksi, year) %>%
  nest(-c(local_authority_district, year)) %>%
  ...
  ...
  ..

The ensemble plot is again reasonably straightforward:

rate_boots_temporal %>%
  ggplot(aes(x=year, y=ksi_rate)) +
  geom_line(data=. %>%  filter(id=="Apparent"), aes(group=id), size=.7) +
  geom_line(data=. %>%  filter(id!="Apparent"), aes(group=id), alpha=.05, size=.2) +
  facet_wrap(~local_authority_district)

Plot specification:

1. Data: The rate_boots_temporal data frame. Note that we include two line layers, one with the observed data (data=. %>% filter(id=="Apparent"), one with the bootsrap data (data=. %>% filter(id!="Apparent").
2. Encoding: x- position varies according to year, y-position according to KSI rate.
3. Marks: geom_line for drawing lines.
4. Facets: facet_wrap() for faceting on local authority.
5. Setting: The bootstrap lines are de-emphasised by making the alpha and size channels very small.

The Hypothetical Outcome Plot (HOP) can be easily created using the gganimate package, simply by adding a call to transition_states() at the end of the plot specification:

rate_boots_temporal %>%
  filter(id!="Apparent") %>%
  ggplot(aes(x=year, y=ksi_rate)) +
  geom_line(aes(group=id), size=.6) +
  facet_wrap(~local_authority_district)+
  transition_states(id, 0,1)