Stochastic Data: Spaced Bimodal (5Ma and 20Ma)

When the peaks are well spaced, binomfit should do a good job fitting the peak ages.

library(FTanalysis)

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

b <- 0.037725003875

# Upper and lower age limits for the random uniform age model
firstAge = 5
secondAge = 20
peaksAt = c(firstAge, secondAge)

ageLabels = c(0, 2, 5, 10, 20, 40, 80)

Create synthetic dataset

Induced track data

For the synthetic sample, we need a set of nI to base it on. We choose to set nI using one of the alpine samples.

input = readAlpsFile(i = 2)

## Warning: number of items read is not a multiple of the number of columns
## Warning: number of items read is not a multiple of the number of columns
## Warning: number of items read is not a multiple of the number of columns
## Warning: number of items read is not a multiple of the number of columns
## Warning: number of items read is not a multiple of the number of columns

nI = c(input$nI, input$nI)  # Create list of induced tracks using the Alpine dataset tracks
nGrains = length(nI)

Expected number of spontaneous tracks given the bimodal age model

synthetic_LambdaS <- sampleLambdaS_ClosureAgePeaks(nI, peaksAt, b, plot = FALSE)

## [1] "2 peak synthetic"

Generate a sample of nS

nS <- sampleNs(synthetic_LambdaS$lambdaS)

The values of nS are: 1, 3, 2, 1, 1, 2, 3, 7, 16, 0, 0, 1, 20, 8, 0, 17, 0, 1, 0, 7, 11, 0, 0, 0, 5, 1, 1, 0, 2, 9, 19, 3, 2, 0, 1, 19, 0, 7, 2, 20, 8, 2, 4, 23, 8, 0, 19, 26, 18, 10, 1, 6, 2, 6, 1, 1, 4, 1, 2, 5, 0, 1, 2, 21, 6, 0, 2, 0, 2, 0, 1, 6, 3, 3, 1, 1, 5, 2, 0, 9, 2, 17, 4, 4, 0, 0, 4, 1, 9, 4, 28, 0, 0, 1, 3, 12, 3, 13, 10, 10, 2, 0

The values of nI are: 44, 46, 58, 11, 16, 146, 23, 55, 192, 5, 106, 12, 258, 54, 2, 171, 44, 129, 9, 111, 63, 8, 41, 49, 86, 67, 50, 3, 48, 90, 210, 170, 150, 13, 21, 172, 3, 86, 132, 354, 80, 53, 21, 334, 101, 15, 130, 446, 133, 90, 18, 44, 46, 58, 11, 16, 146, 23, 55, 192, 5, 106, 12, 258, 54, 2, 171, 44, 129, 9, 111, 63, 8, 41, 49, 86, 67, 50, 3, 48, 90, 210, 170, 150, 13, 21, 172, 3, 86, 132, 354, 80, 53, 21, 334, 101, 15, 130, 446, 133, 90, 18

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)

## Warning: position_stack requires constant width: output may be incorrect

Calculate the central age

BINOMFIT_CentralAge(FTdataset)

## $centralAge
## [1] 13.86
## 
## $centralAgeStError
## [1] 1.209
## 
## $ChiSquared
## [1] 285.9
## 
## $degFreedom
## [1] 101
## 
## $eta
##  [1] 0.05316 0.05350 0.05369 0.05383 0.05391 0.05394 0.05395 0.05396
##  [9] 0.05396 0.05396 0.05396 0.05396 0.05396 0.05396 0.05396 0.05396
## [17] 0.05396 0.05396 0.05396 0.05396
## 
## $dispersion
## [1] 63.71

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_4AgeModels(FTdataset, resultsOutput1, resultsOutput2, 
    resultsOutput3, resultsOutput4, 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 82 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: position_stack requires constant width: output may be incorrect
## Warning: Removed 82 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: position_stack requires constant width: output may be incorrect
## Warning: Removed 82 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: position_stack requires constant width: output may be incorrect
## Warning: Removed 82 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).

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

Compare FTanalysis results with the Bimodal Ages (in blue)

Bimodal Ages: The input ages are at 5, 20 FTanalysis: 2 peaks fitted at 4.2986, 20.7275

comparisonPlot <- compareModelSolutions(resultsList, favouredBIC_nPeaks)
comparisonPlot <- comparisonPlot + geom_vline(xintercept = peaksAt, colour = "blue") + 
    xlim(0, 60)
comparisonPlot

## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).