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_TarimW3Data)
input = TarimW3Data$data
benchmarkSolutions = TarimW3Data$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: 13, 42, 2, 80, 11, 7, 12, 53, 14, 121, 12, 2, 91, 19, 5, 10, 48, 55, 4, 4, 5, 3, 9, 16, 26, 6, 13, 29, 39, 65, 229, 133, 34, 32, 19, 7, 36, 9, 17, 5, 9, 51, 41, 16, 51, 35, 28, 71, 45, 19
The values of nI are: 220, 114, 7, 229, 104, 54, 176, 586, 182, 168, 188, 101, 458, 185, 12, 95, 595, 110, 64, 24, 45, 152, 59, 117, 79, 69, 108, 328, 87, 79, 568, 323, 346, 195, 181, 34, 280, 65, 131, 57, 22, 648, 482, 183, 720, 436, 273, 123, 120, 241
– 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)
BINOMFIT_CentralAge(FTdataset)
## $centralAge
## [1] 46.2
##
## $centralAgeStError
## [1] 6.255
##
## $ChiSquared
## [1] 1053
##
## $degFreedom
## [1] 49
##
## $eta
## [1] 0.1428 0.1472 0.1498 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503
## [11] 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503 0.1503
##
## $dispersion
## [1] 80.66
Run Binomfit analysis for different numbers of peaks and hold these results in a list.
resultsOutput1 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 1, K = input$K,
details = input$details, verbose = TRUE)
resultsOutput2 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 2, K = input$K,
details = input$details, verbose = FALSE)
resultsOutput3 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 3, K = input$K,
details = input$details, verbose = FALSE)
resultsOutput4 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 4, K = input$K,
details = input$details, verbose = FALSE)
resultsOutput5 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 5, K = input$K,
details = input$details, verbose = FALSE)
resultsOutput6 <- BINOMFIT(FTdataset, peakAgeModel = 1, PkNum = 6, K = input$K,
details = input$details, verbose = FALSE)
resultsList <- list(resultsOutput1, resultsOutput2, resultsOutput3, resultsOutput4,
resultsOutput5, resultsOutput6)
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")
| nAges | logLike | BIC | deltaBIC | |
|---|---|---|---|---|
| 1 | 1.00 | -449.84 | 903.59 | 508.04 |
| 2 | 2.00 | -214.94 | 441.61 | 46.06 |
| 3 | 3.00 | -199.80 | 419.16 | 23.61 |
| 4 | 4.00 | -184.08 | 395.55 | 0.00 |
| 5 | 5.00 | -184.11 | 403.43 | 7.87 |
| 6 | 6.00 | -184.09 | 411.21 | 15.66 |
BICmodelComparisonPlot
n = which(BICarray == min(BICarray))
favouredBIC_nPeaks <- resultsList[[n]]$PkNum
favouredBIC_ages <- resultsList[[n]]$PeakAgeResults
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")
| Model1 | Model2 | Prob_F_ByChanceAlone | |
|---|---|---|---|
| 1 | 1 | 2 | 0.00 |
| 2 | 2 | 3 | 0.04 |
| 3 | 3 | 4 | 0.61 |
| 4 | 4 | 5 | 4.65 |
| 5 | 5 | 6 | 5.13 |
ggplot(df, aes(x = Model2, y = Prob_F_ByChanceAlone)) + geom_point(color = "blue") +
geom_line()
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 76 rows containing missing values (geom_path).
## Warning: Removed 76 rows containing missing values (geom_path).
## Warning: Removed 76 rows containing missing values (geom_path).
## Warning: Removed 76 rows containing missing values (geom_path).
## Warning: Removed 76 rows containing missing values (geom_path).
## Warning: Removed 76 rows containing missing values (geom_path).
Benchmark analysis: 3 peaks with the minimum at ~22Ma with Lmax=36.31
FTanalysis: 4 peaks fitted at 22.6168, 46.8645, 103.9773, 179.6844
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(22), colour = "green")
comparisonPlot
