Stochastic Data: Random Uniform

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

b <- 0.037725003875

minAge = 5
maxAge = 20

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 give the step function age model

synthetic_LambdaS <- sampleLambdaS_StepFunction(nI, minAge, maxAge, b)

## [1] "Step function synthetic"

Generate a sample of nS

nS <- sampleNs(synthetic_LambdaS$lambdaS)

The values of nS are: 1, 3, 1, 1, 0, 8, 0, 4, 6, 0, 5, 0, 14, 3, 1, 9, 1, 7, 0, 8, 4, 0, 2, 2, 6, 5, 1, 0, 1, 5, 8, 15, 9, 1, 1, 4, 0, 8, 18, 11, 4, 0, 1, 23, 8, 1, 3, 23, 13, 9, 1, 2, 0, 2, 0, 2, 11, 0, 1, 12, 0, 1, 1, 16, 2, 0, 16, 4, 6, 0, 9, 2, 0, 1, 3, 2, 2, 3, 0, 5, 3, 15, 15, 9, 0, 0, 6, 0, 1, 15, 36, 1, 4, 1, 27, 3, 0, 6, 9, 7, 3, 3

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)

Calculate the central age

BINOMFIT_CentralAge(FTdataset)

## $centralAge
## [1] 13.14
## 
## $centralAgeStError
## [1] 0.8345
## 
## $ChiSquared
## [1] 139.8
## 
## $degFreedom
## [1] 101
## 
## $eta
##  [1] 0.05268 0.05227 0.05202 0.05179 0.05161 0.05149 0.05142 0.05137
##  [9] 0.05135 0.05133 0.05132 0.05132 0.05132 0.05132 0.05131 0.05131
## [17] 0.05131 0.05131 0.05131 0.05131
## 
## $dispersion
## [1] 33.19

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: Removed 85 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 85 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 85 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 85 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 Random Uniform Model Age Limits (in blue)

FTanalysis: 2 peaks fitted at 9.2722, 17.923

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