Estimating uncertainty in bed surface grain sizes from binomial theory

1 Brett’s GitHub Page

For other R Notebooks on our evolving ideas about geomorphology, see Brett’s GitHub Page

2 Overview

One of the basic measurements using in the study of gravel bed rivers is the bed surface grain size distribution. This distribution is used to compare one site to another, and to infer relative sediment supply by calculating the armour ratio. It is also often used to track changes in stream conditions following a change in sediment supply or flow regime. However, while geomorphologists often discuss bed surface texture variability, this is almost never done in the context of uncertainty of the estimates from a particular sample. Simply put, we often do not consider the statistical power of our samples, which depends on the size of the sample taken (i.e., how many stones were measured), and on the nature of the distribution being sampled.

However, binomial theory can be used to estimate the uncertainty of any estimate of the surface percentile.

3 Grain size data

let’s first load some data, and turn it into a data frame containing all measurements from all locations along the stream table. Since most bed surface samples are taken in the field, we will apply a 1:25 scale ratio to the data in order to translate our model grain sizes into prototype equivalent grain sizes.

rm(list = ls())
library(reshape2)
Lr = 25  #specify the scale ratio for the model: prototype relation
grain.sizes = read.csv('Data/GSD1HIGHQ_4_grainmeasurements.csv')
all.data = melt(grain.sizes, variable.name = "location", value.name = "D")

No id variables; using all as measure variables

all.data$D = all.data$D * Lr

Now lets perform a cumulative frequency distribution analysis on the data.

max.phi = ceiling(log2(max(all.data$D)))
min.phi = floor(log2(min(all.data$D)))
x = seq(from = min.phi, to = max.phi, by = 0.5)
results = hist(all.data$D, 
               breaks = 2^x,
               plot = FALSE)
cfd = data.frame(x)
cfd$y = c(0, cumsum(results$counts))/sum(results$counts)
plot(2^cfd$x, cfd$y,
     type = "b",
     col = "blue",
     log = "x",
     xlab = "grain size (mm)",
     ylab = "cum. prop. finer")

We can estimate the percentiles from these data as follows:

phi.est = approx(x = cfd$y, y = cfd$x, xout = 0.50)
phi.est$di = 2^phi.est$y

Note that we have interpolated based on the phi sizes, and then converted those values to grain sizes in mm. A comparison shows that there are small but consistent biases introduced by applying a linear interpolation to the grain size data rather than the phi data. The fact that we typically display the grain size using log2(D) as the x axis and then often estimate grain sizes graphically implies that we believe any interpolation should linearly relate the cumulative proportion finer to log2(D).

size.est = approx(x = cfd$y, y = 2^cfd$x, xout = 0.50)
phi.est$di

[1] 32.86458

size.est$y

[1] 33.0196

These errors are small compared to the uncertainties associated with the estimates of each percentile, based on binomial theory.

4 Estimating confidence intervals from binomial theory

We can use binomial theory to estimate the uncertainty of each percentile (in terms of percentile values), and then translate these uncertainties into ranges of grain sizes using the sample data.

For example, we know that for the median size, the expected number of observations less than the median value ought to be \(np\), where \(n\) is the total number of observations and \(p\) is the percentile of interest (i.e. 0.50 in this case). However, binomial theory implies that the actual number of observations less than this value is a normally distributed quantity with a mean value of \(np\) and a standard deviation of \(\sqrt{np(1-p)}\). It is quite easily to calculate the confidence interval for any percentile in terms of the percentile values themselves.

alpha = 0.05
percentile = seq(0.1,0.9,0.1)
n = c(100,200,400)
tcrit = qt(1 - alpha/2, n)
ci.upper = matrix(data = NA, nrow = length(percentile), ncol = length(n))
ci.lower = matrix(data = NA, nrow = length(percentile), ncol = length(n))
for (i in seq_along(percentile)){
  sigma = sqrt(n*percentile[i]*(1-percentile[i]))
  ci.upper[i,] = round((n*percentile[i] + tcrit*sigma)/n, digits = 2)
  ci.lower[i,] = round((n*percentile[i] - tcrit*sigma)/n, digits = 2)
}
ci.percentile = data.frame(percentile,
                           ci.lower[,1],
                           ci.upper[,1],
                           ci.lower[,2],
                           ci.upper[,2],
                           ci.lower[,3],
                           ci.upper[,3] )
names(ci.percentile) =  c('Percentile', 'L100', 'U100', 'L200', 'U200','L400', 'U400' )
ci.percentile

For example, the 50th percentile of the distribution estimated from a sample of 100 observations could fall between the 40th percentile and the 60th percentile of the sample data. For a sample of 400 observations, the 50th percentile would be expected to fall between the 45th and the 55th percentile of the sample.

In order to estimate the confidence interval in mm, it is necessary to know something about the underlying distribution. A function that estimates the percentile size and confidence bounds for any percentile and confidence limit, \(\alpha\), can be written as follows:

ErrorBounds = function(n, p, cfd, alpha = 0.05){
  tcrit = qt(1 - alpha/2, n)
  sigma = sqrt(n*p*(1-p))
  p.upper = (n*p + tcrit*sigma)/n
  p.lower = (n*p - tcrit*sigma)/n
  bounds = c(p.lower, p.upper, p)
  estimates = approx(x = cfd$y, y = cfd$x, xout = bounds)
  estimates$y = 2^estimates$y
return(estimates)
}

The error bounds around our sample of 500 measurements for various percentiles are easily calculated:

conf.50 = ErrorBounds(500, 0.5, cfd)
conf.84 = ErrorBounds(500, 0.84, cfd)
conf.90 = ErrorBounds(500, 0.90, cfd)
conf.95 = ErrorBounds(500, 0.95, cfd)

For such a large sample, the uncertainty range for the \(D_{50}\) is 30 – 36 mm. The error bounds around the \(D_{84}\) are larger, at about 59 – 74 mm. It is even bigger for the largest grain sizes: the uncertainty of the \(D_{90}\) is about 74 – 98 mm.

The size of the sample has a great effect on the uncertainty range. For example, consider the confidence intervals about the \(D_{50}\), calculated and plotted below.

N = seq(from = 100, to = 500, by = 5)
D50.civn = matrix(data = NA, nrow = length(N), ncol = 3)
for (i in seq_along(N)){
  D50.civn[i,] = ErrorBounds(N[i], 0.50, cfd)[[2]]
}
matplot(N, D50.civn,
     type = "l",
     log = "y",
     ylim = 2^c(min.phi, max.phi),
     lty = 1)
abline(h = 2^x,
       lwd = 0.5,
       lty = 2)

abline(v = seq(100,500,50),
       lwd = 0.5,
       lty = 2)

For a sample of 100 stones, the confidence interval for \(D_{50}\) is 27 – 40 mm, but for a sample of 200 stones the confidence interval is 29 – 38 mm, and for a 400 stone sample is 30 – 36 mm. (The estimate of \(D_{50}\) in all cases is 33 mm). Generally, the confidence interval is approximately symmetrical about the estimate, but this is not always true, and depends on the underlying distribution. Thus, while the confidence interval represents about 37% of the estimate value for \(N = 100\), it represents only 18% for \(N = 400\). Clearly, the desired level of precision must be considered when designing a bed surface sampling strategy.

The issue is rather more important when estimating larger-than-average grain sizes. Consider how the uncertainty for the \(D_{84}\) varies with the size of the sample.

N = seq(from = 100, to = 500, by = 5)
D84.civn = matrix(data = NA, nrow = length(N), ncol = 3)
for (i in seq_along(N)){
  D84.civn[i,] = ErrorBounds(N[i], 0.84, cfd)[[2]]
}
matplot(N, D84.civn,
     type = "l",
     log = "y",
     ylim = 2^c(min.phi, max.phi),
     lty = 1)
abline(h = 2^x,
       lwd = 0.5,
       lty = 2)

abline(v = seq(100,500,50),
       lwd = 0.5,
       lty = 2)

For a sample of 100 stones, the confidence interval for \(D_{84}\) is 55 – 92 mm, but for a 400 stone sample is 59 – 75 mm. (The estimate of \(D_{84}\) is 63 mm). Because of the underlying distribution, the confidence interval is not symmetrical about the estimate, particularly for small sample sizes. In this case the confidence interval represents about 58% of the estimate value for \(N = 100\), and 26% for \(N = 400\).

For the \(D_{95}\), the analysis looks like this:

N = seq(from = 100, to = 500, by = 5)
D95.civn = matrix(data = NA, nrow = length(N), ncol = 3)
for (i in seq_along(N)){
  D95.civn[i,] = ErrorBounds(N[i], 0.95, cfd)[[2]]
}
matplot(N, D95.civn,
     type = "l",
     log = "y",
     ylim = 2^c(min.phi, max.phi),
     lty = 1)
abline(h = 2^x,
       lwd = 0.5,
       lty = 2)

abline(v = seq(100,500,50),
       lwd = 0.5,
       lty = 2)

For a sample of 100 stones, the confidence interval for \(D_{95}\) is 89 – 173 mm, but for a 400 stone sample is 99 – 121 mm. (The estimate of \(D_{95}\) is 109 mm). For small values of \(N\), the confidence interval is not symmetrical about the estimate. In this case the confidence interval represents about 77% of the estimate value for \(N = 100\), and 20% for \(N = 400\).

5 Choosing a sample size

In order to choose an appropriate sample size, it is first necessary to have a sense of the shape of the distribution, and to know what percentiles are of interest. For example, if we assume the data used herein is reasonably representative of the bed surface of interest (in terms of distribution shape, at least), then we can easily model the relation between the confidence interval width and the number of measurements taken.

ci.width50 = 100*(D50.civn[,2] - D50.civn[,1])/D50.civn[,3]
plot(N, ci.width50, 
     type = "l",
     ylab = "C.I. Width (% of est.)",
     col = "blue")
abline(h = seq(5,95,5),
       lwd = 0.5,
       lty = 2)

abline(v = seq(100,500,50),
       lwd = 0.5,
       lty = 2)

If we wish to have a confidence interval no wider than about 30% of the estimate of the \(D_{50}\), then we would need a sample size of 155.

However, if we wish to have equivalent confidence in an estimate of \(D_{90}\), our sample size would have to be larger.

D90.civn = matrix(data = NA, nrow = length(N), ncol = 3)
for (i in seq_along(N)){
  D90.civn[i,] = ErrorBounds(N[i], 0.90, cfd)[[2]]
}
ci.width90 = 100*(D90.civn[,2] - D90.civn[,1])/D90.civn[,3]
plot(N, ci.width90, 
     type = "l",
     ylab = "C.I. Width (% of est.)",
     col = "blue")
abline(h = seq(5,95,5),
       lwd = 0.5,
       lty = 2)

abline(v = seq(100,500,50),
       lwd = 0.5,
       lty = 2)

In fact, we would need a sample of size \(N =\) 415 have a confidence interval no greater than 30% of the estimate.

6 Actual Sample Variability

We can use our sample of 500 measurements to estimate this sampling effect by subsetting the data, and using the subset to calculate the median grain size. In this case, we assume that our sample of 500 measurements represents the underlying population of grains, and that the true median value is 33 mm. Since the sample of 500 measurements is in fact a sample itself of the underlying population, this assumption is not strictly true, and could introduce some small distortions in the comparision.

k = 100  #set the sample size
m = 1000  #set the number of times to sample the data
sample.50 = vector(mode = "numeric", length = m)
for(i in seq_along(sample.50)){
  sample.data = sample(all.data$D, size = k)
  sample.hist = hist(sample.data, 
               breaks = 2^x,
               plot = FALSE)
  sample.cfd = c(0, cumsum(sample.hist$counts))/sum(sample.hist$counts)
  sample.50[i] = 2^approx(x = sample.cfd, y = x, xout = 0.50)[[2]]
}
boxplot(sample.50,
        ylim = D50.civn[which(N == k),1:2]*c(0.8, 1.2))
abline(h = D50.civn[which(N == k),],
       col = "red",
       lty = 2)

On the plot, the confidence bounds estimated from binomial theory and the median based on all data are shown as dashed red lines, and the boxplot presents the distribution of estimates based on sampling the original distribution with a sample size of 100. There are two points to be made: (1) the effect of sample size does introduce significant uncertainty into the estimate of the median size; and (2) the binonmial theory captures this uncertainty very well.

For this analysis, we find 13 observations that are outside the 95% confidence intervals predicted by binomial theory. Based on the sample size, we would expect 50 observations to fall outside of the confidence intervals.

For fun, lets see what happens with a larger sample size.

k = 200  #set the sample size
m = 1000  #set the number of times to sample the data
sample.50 = vector(mode = "numeric", length = m)
for(i in seq_along(sample.50)){
  sample.data = sample(all.data$D, size = k)
  sample.hist = hist(sample.data, 
               breaks = 2^x,
               plot = FALSE)
  sample.cfd = c(0, cumsum(sample.hist$counts))/sum(sample.hist$counts)
  sample.50[i] = 2^approx(x = sample.cfd, y = x, xout = 0.50)[[2]]
}
boxplot(sample.50,
        ylim = D50.civn[which(N == k),1:2]*c(0.8, 1.2))
abline(h = D50.civn[which(N == k),1:2],
       col = "red",
       lty = 2)

We find that 13 observations fall outside the 95% confidence intervals predicted by binomial theory; this is presumably related to the fact that samples of 200 observations taken from a population of 500 observations will naturally have a smaller variability than if they were taken from the true grain size population (i.e., the samples can hardly be considered to be independent). Nevertheless, there is still substantial variability in the estimate of the median, and it is similar to the range predicted by binomial theory.

And what about larger grain sizes, like the \(D_{84}\)?

k = 100  #set the sample size
m = 1000  #set the number of times to sample the data
sample.84 = vector(mode = "numeric", length = m)
for(i in seq_along(sample.84)){
  sample.data = sample(all.data$D, size = k)
  sample.hist = hist(sample.data, 
               breaks = 2^x,
               plot = FALSE)
  sample.cfd = c(0, cumsum(sample.hist$counts))/sum(sample.hist$counts)
  sample.84[i] = 2^approx(x = sample.cfd, y = x, xout = 0.84)[[2]]
}
boxplot(sample.84,
        ylim = D84.civn[which(N == k),1:2]*c(0.8, 1.2))
abline(h = D84.civn[which(N == k),],
       col = "red",
       lty = 2)

For this grain size, We find that 13 observations fall outside the 95% confidence intervals.