Benchmark Example: FT Analysis Tarim W2 Dataset

This is a basic example showing how to analyse a set of spontaneous and induced tracks using my implementation of binomfit. This basically tries to emulate Brandon's program BINOMFIT V1.8 in terms of functionality and analyses.

This example is demonstrated within an R Workbook and the analysis is performed using the code you can read below. This html page is the output of the analysis so it is easy for you to follow how the code works. In this example, we are loading the MtTom data which comes in the FTanalysis package.

Your own scipt would differ in that it would be loaded from a separate file.

library(FTanalysis)

## Loading required package: ggplot2
## Loading required package: scales
## Loading required package: zipfR
## Loading required package: xtable
## Loading required package: grid

data(FT_TarimW2Data)
input = TarimW2Data$data
benchmarkSolutions = TarimW2Data$benchmarkResults

ageLabels = c(10, 20, 50, 100, 200)

nS = input$nS  # List of spontaneous tracks
nI = input$nI  # List of induced tracks
dPar = input$SqNum

The values of nS are: 8, 17, 19, 9, 7, 25, 2, 14, 6, 7, 7, 6, 16, 6, 16, 12, 18, 32, 5, 17, 16, 29, 39, 26, 18, 20, 41, 26, 13, 25, 48, 19

The values of nI are: 133, 136, 210, 146, 58, 257, 35, 160, 109, 133, 35, 137, 209, 141, 190, 310, 175, 197, 131, 304, 374, 405, 482, 369, 445, 240, 107, 146, 191, 364, 123, 471

Create the FTdataset object and the data summary plots

– This is what FTanalysis used to hold all the input data and various processed forms of the data.

FTdataset <- makeFTdataset(nS = nS, nI = nI, rhoD = input$RhoD, relErrRhoD = input$RERhoD, 
    c = 1, K = 1, Zeta = input$Zeta, relErrZeta = input$REZeta, SqSize = NULL, 
    geomFactor = 0.5)

trackCountSummaryPlot <- plotTrackCountSummary(FTdataset)
PDplotNoOverlay <- PDplot(FTdataset, resultsOutput, plotType = 1, zeroNsOffset = 0.5)
plot2 <- PDplot(FTdataset, resultsOutput, plotType = 6, zeroNsOffset = 0.5)

layout <- matrix(c(1, 2, 3), nrow = 1, byrow = TRUE)
multiplot(trackCountSummaryPlot, PDplotNoOverlay, plot2, layout = layout)

Calculate the central age

BINOMFIT_CentralAge(FTdataset)

## $centralAge
## [1] 24.15
## 
## $centralAgeStError
## [1] 3.567
## 
## $ChiSquared
## [1] 273.8
## 
## $degFreedom
## [1] 31
## 
## $eta
##  [1] 0.07595 0.07882 0.08261 0.08425 0.08452 0.08453 0.08452 0.08452
##  [9] 0.08452 0.08452 0.08452 0.08452 0.08452 0.08452 0.08452 0.08452
## [17] 0.08452 0.08452 0.08452 0.08452
## 
## $dispersion
## [1] 70.48

Perform the binomfit analysis

Run Binomfit analysis for different numbers of peaks and hold these results in a list.

resultsOutput1 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 1, K = input$K, 
    details = input$details, verbose = TRUE)
resultsOutput2 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 2, K = input$K, 
    details = input$details, verbose = FALSE)
resultsOutput3 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 3, K = input$K, 
    details = input$details, verbose = FALSE)
resultsOutput4 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 4, K = input$K, 
    details = input$details, verbose = FALSE)
resultsOutput5 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 5, K = input$K, 
    details = input$details, verbose = FALSE)
resultsOutput6 <- BINOMFIT(FTdataset, peakAgeModel = 2, PkNum = 6, K = input$K, 
    details = input$details, verbose = FALSE)

resultsList <- list(resultsOutput1, resultsOutput2, resultsOutput3, resultsOutput4, 
    resultsOutput5, resultsOutput6)

Model comparison

BIC model comparison

logLikeArray = c()
for (i in 1:6) {
    logLikeArray <- c(logLikeArray, resultsList[[i]]$logLike)
}

nAgeArray <- c(1, 2, 3, 4, 5, 6)
BICarray <- -2 * logLikeArray + (2 * nAgeArray - 1) * log(FTdataset$nGrain)
deltaBIC <- BICarray - min(BICarray)

BICdf <- data.frame(nAges = nAgeArray, logLike = logLikeArray, BIC = BICarray, 
    deltaBIC = deltaBIC)

xt <- xtable(BICdf)
BICmodelComparisonPlot <- ggplot(BICdf, aes(x = nAges, y = deltaBIC)) + geom_point(color = "blue") + 
    geom_line()

print(xt, type = "html")

BICmodelComparisonPlot

n = which(BICarray == min(BICarray))
favouredBIC_nPeaks <- resultsList[[n]]$PkNum
favouredBIC_ages <- resultsList[[n]]$PeakAgeResults

Chi2 model comparison (needs a little tidying up)

ModelCompProb <- c()
F <- c()

for (i in 1:5) {
    tmp = getChi2Comp(resultsList[[i]]$ChiSq, resultsList[[i]]$degFreedom, resultsList[[i + 
        1]]$ChiSq, resultsList[[i + 1]]$degFreedom)
    ModelCompProb[i] <- tmp$P
}

df = data.frame(Model1 = seq(1, 5), Model2 = seq(2, 6), Prob_F_ByChanceAlone = ModelCompProb * 
    100)

print(xtable(df), type = "html")

ggplot(df, aes(x = Model2, y = Prob_F_ByChanceAlone)) + geom_point(color = "blue") + 
    geom_line()

Summary Plots of the binomfit results

This plot summarises the data in the first 3 figures, and the models fitted by FTanalysis in the remainder. Pairs of PD plots and Radial Plots are shown for models with upto 6 peaks. The prefered model is indicated using the minimum BIC.

makeBinomfitSummaryPlot_6AgeModels(FTdataset, resultsOutput1, resultsOutput2, 
    resultsOutput3, resultsOutput4, resultsOutput5, resultsOutput6, ageLabels, 
    dataTrasformStyle = "arcsinTransformation")

## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 79 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).

Compare FTanalysis results with Benchmark Solutions

Benchmark analysis: 3 peaks with the minimum at ~17Ma

FTanalysis: 3 peaks fitted at 16.8995, 39.8364, 99.8927

The green line below shows the age quoted in the literature as the minimum peak age consistent with the data

comparisonPlot <- compareModelSolutions(resultsList, favouredBIC_nPeaks)
comparisonPlot <- comparisonPlot + geom_vline(xintercept = c(17), colour = "green")
comparisonPlot