calculating confidence intervals for surface grain size samples

Introduction

This document presents a step-by-step demonstration of how the GSDtools package created by Brett Eaton and Dan Moore can be applied to your bed surface grain size data. The package implements the use of binomial sampling theory to estimate percentile confidence ranges, as described by Eaton et. al (in review). The package contains a number of functions to easily visualize the data, but this example focuses on how to generate estimates of the upper and lower confidence bounds, and how to export them for use in a spreadsheet.

Setting up your working environment

We assume that you have already installed the R language on your computer, and that you have a graphical interface in which to run the code we present in this document. To install R, start at the R Project website. We recommend using the R Studio graphical interface to run R code (which can be downloaded here), but you may use whatever interface you are most comfortable with.

Installing the GSDtools package

The GSDtools package is hosted on Brett Eaton’s GitHub page, and is freely accessible to download. To do this, we will use commands in the devtools package. If you have not installed devtools yet, you can simply run the following command.

if("devtools" %in% rownames(installed.packages()) == FALSE) 
  {install.packages("devtools")}

Now we can install the GSDtools package from GitHub.

if("GSDtools" %in% rownames(installed.packages()) == FALSE){
  devtools::install_github("bceaton/GSDtools")}

Finally, we need to tell R that we wish to use the functions in the GSDtools package by running the following command.

At this point, we should be ready to go with our analysis.

Entering the data

Users familiar with R will know that there are various ways of importing data and creating the necessary variables for analysis. To keep things simple, we will enter the grain size data as a string of numbers, the cumulative percent finer data as another string of data, and then combine them into a single table of data called a data.frame in R. The following line of code creates a list of grain size data in the standard Wentworth size classes.

Now lets enter the cumulative proportion finer data that corresponds to each size.

Finally, we combine these two strings of number into a single data frame for analysis.

You can verify that the data match up properly by looking at the table.

print(my.data)

We can also plot the data to ensure that it all looks right.

plot(my.data, type = "b", log = "x", xlim = c(1, 1000))

Estimating the confidence bounds

Once we have a data.frame with the necessary data that has variables called size and probs, we can calculate any percentile of interest (e.g. \(D_{50}\)), and its associated confidence range, but only if we know the number of stones that were measured to define the grain size distribution. Let’s assume that our distribution was based on a Wolman sample of 103 stones.

This can be done for a single percentile as follows:

This produces a data.frame that records the percentile being estimated, the estimate of the percentile, and the upper and lower bounds for that estimate.

print(my.D50)

A more useful example would be to estimate the percentiles between \(D_{10}\) and \(D_{90}\) in increments of 5.

This produces a data.frame with several lines of data, each representing an estimated percentile, along with its confidence interval.

head(my.Dvalues)  #this command only shows the first few lines of data

Exporting the data

To export the data as a .csv file, we simply specify a file name, and run the following command.

The file created by this command will be saved to whatever your working director is. If you don’t know what your working directory is, you can find out by running the command.

getwd()

You can use setwd or the user interface menus to navigate to the directory where you want to save your file. In R Studio, you go first to the Session menu, then to Set Working Directory, and finally you choose Choose Directory.

The file you just generated can be imported into a spreadsheet, and you can generate whatever kind of plots you would like using the data. If you

Visualizing the data

The GSDtools package also includes functions to make it easier to plot the confidence bounds. The key one is PolyCI. Below, we plot the grain size data as before, then add a polygon the represents the confidence interval using the same inputs as were used above to generate my.Dvalues.

{plot(my.data, type = "b", log = "x", xlim = c(1, 1000))
my.polygon = PolyCI(cfd = my.data, n = my.n, P = my.percentile, alpha = my.alpha)
polygon(my.polygon, lty = 0, col = rgb(0,0,1,0.25))  #makes a blue polygon 
abline(h = c(0.1, 0.9), lty = 2)}  #add lines indicating range of confidence interval

Comparing two distributions

There are two tools for comparing samples to determine if they are statistically similar or not. One tool uses raw data comprising measurements of b axis diameters for individual stones. The other tool uses data binned into size classes. We will use a random number generator to produce two samples of raw data to demonstrate the use of the first tool (CompareRAWs), then we will bin that data to produce a cumulative frequency distribution to demonstrate the use of the second tool (CompareCFDs).

We can create two samples as follows:

The CompareRAWs tool uses a re-sampling (with replacement) approach to produce bootstrap confidence intervals for the differences in user-specified percentiles.

first.comparison = CompareRAWs(Raw1,Raw2, P = 50)  #compare the 50th percentile
print(first.comparison)

Hint: You can play around with the sample size (n1 and n2) to figure out when the differences in the 50th percentile are statistically significant, and when they are not.

We can also compare a number of percentiles at the same time.

scnd.comparison = CompareRAWs(Raw1,Raw2, P = seq(10,90,10))  #compare the 50th percentile
print(scnd.comparison)

hint: try running the sample generation code (i.e., Raw1 = 2^rnorm(n1, mean = 4.1, sd = 1)) a few times with the same values and see how the results change.

Usually, we do not have measurements of individual grain sizes, we have binned data. We can transform our samples into cumulative frequency distributions as follows.

GSD1 = MakeCFD(Raw1)
GSD2 = MakeCFD(Raw2)
head(GSD1)  #lets look at first few lines of the results for our first sample

We can use the inverse transform sampling approach to get a bootstrap estimate of the confidence interval for the differences between percentiles using the CompareCFDs tool.

thrd.comparison = CompareCFDs(GSD1, GSD2, n1, n2, P = seq(10,90,10))
print(thrd.comparison)

Because we are using binned data, the CompareCFDs tool will sometimes produce slightly different results than the CompareRAWs tool, but our tests of the tools shows that they are generally produce relatively similar results. The differences are presumably introduced by the act of binning the raw data into 0.5\(\phi\) classes.