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_BengalData)
input = BengalData$data
benchmarkSolutions = BengalData$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: 2, 1, 1, 24, 2, 2, 14, 4, 2, 16, 4, 0, 3, 1, 2, 0, 0, 23, 1, 0, 8, 10, 5, 0, 0, 15, 5, 10, 0, 11, 0, 5, 5, 1, 3, 2, 18, 2
The values of nI are: 40, 89, 15, 578, 36, 28, 286, 93, 162, 196, 114, 105, 47, 74, 21, 15, 19, 373, 14, 8, 135, 218, 55, 12, 127, 328, 61, 341, 11, 207, 39, 276, 295, 57, 65, 19, 402, 115
– 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)
## Warning: position_stack requires constant width: output may be incorrect
BINOMFIT_CentralAge(FTdataset)
## $centralAge
## [1] 9.532
##
## $centralAgeStError
## [1] 1.007
##
## $ChiSquared
## [1] 56.63
##
## $degFreedom
## [1] 37
##
## $eta
## [1] 0.03827 0.03788 0.03781 0.03775 0.03771 0.03768 0.03767 0.03766
## [9] 0.03765 0.03764 0.03764 0.03764 0.03764 0.03764 0.03764 0.03764
## [17] 0.03764 0.03764 0.03764 0.03764
##
## $dispersion
## [1] 35.22
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)
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 | -83.54 | 170.71 | 2.97 |
| 2 | 2.00 | -78.42 | 167.75 | 0.00 |
| 3 | 3.00 | -78.42 | 175.03 | 7.29 |
| 4 | 4.00 | -78.42 | 182.31 | 14.56 |
| 5 | 5.00 | -78.42 | 189.58 | 21.84 |
| 6 | 6.00 | -78.42 | 196.85 | 29.11 |
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 | 6.06 |
| 3 | 3 | 4 | 6.42 |
| 4 | 4 | 5 | 6.90 |
| 5 | 5 | 6 | 7.41 |
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: position_stack requires constant width: output may be incorrect
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
## Warning: position_stack requires constant width: output may be incorrect
## Warning: Removed 87 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_segment).
## Warning: Removed 3 rows containing missing values (geom_text).
Benchmark analysis: 2 peaks fitted at 3.65, 11.72
FTanalysis: 2 peaks fitted at 3.6701, 11.7292
Below, we show a summary of the ages for each of the fitted models and the benchmark data from the literature.
comparisonPlot <- compareModelSolutions(benchmarkSolutions, resultsList, favouredBIC_nPeaks)
comparisonPlot
