library(tracerer)
log_file <- get_tracerer_path("beast2_example_output.log")
estimates <- parse_beast_log(log_file)
knitr::kable(estimates)
Sample posterior likelihood prior treeLikelihood TreeHeight BirthDeath birthRate2 relativeDeathRate2
0 -76.40152 -66.56126 -9.840263 -66.56126 0.6383792 -2.932507 1.0000000 0.5000000
1000 -68.68523 -60.11853 -8.566701 -60.11853 0.7003758 -1.658946 2.8041208 0.4952765
2000 -70.06839 -58.93807 -11.130319 -58.93807 1.4085182 -4.222564 0.7901276 0.3674313
3000 -69.03990 -58.90834 -10.131561 -58.90834 0.9220334 -3.223806 1.8637476 0.2496224
4000 -74.15268 -59.98232 -14.170365 -59.98232 1.8159958 -7.262610 1.7823070 0.3519155
5000 -71.17163 -61.56657 -9.605062 -61.56657 0.7408032 -2.697307 1.7849504 0.7085961
6000 -69.64176 -58.73713 -10.904631 -58.73713 1.0967406 -3.996875 1.1182466 0.3702605
7000 -71.92546 -61.43600 -10.489461 -61.43600 1.0088078 -3.581706 1.2610806 0.4008574
8000 -69.59441 -58.89381 -10.700593 -58.89381 0.8291479 -3.792838 0.3909076 0.5845426
9000 -69.69113 -62.40904 -7.282093 -62.40904 0.4529637 -0.374338 1.1123238 0.6983194
10000 -71.86884 -60.73520 -11.133636 -60.73520 0.7693617 -4.225880 1.5626756 0.7107459

Display the ESSes:

esses <- rep(NA, ncol(estimates))
burn_in_fraction <- 0.1
for (i in seq_along(estimates)) {
  # Trace with the burn-in still present
  trace_raw <- as.numeric(t(estimates[i]))

  # Trace with the burn-in removed
  trace <- remove_burn_in(trace = trace_raw, burn_in_fraction = 0.1)

  # Store the effectice sample size
  esses[i] <- calc_ess(trace, sample_interval = 1000)
}

# Note that the first value of three is nonsense:
# it is the index of the sample. I keep it in
# for simplicity of writing this code
expected_esses <- c(3, 10, 10, 10, 10, 7, 10, 9, 6)
testit::assert(all(expected_esses - esses < 0.5))

df_esses <- data.frame(esses)
rownames(df_esses) <- names(estimates)
knitr::kable(df_esses)
esses
Sample 2.668464
posterior 10.000000
likelihood 10.000000
prior 10.000000
treeLikelihood 10.000000
TreeHeight 6.657254
BirthDeath 10.000000
birthRate2 8.905181
relativeDeathRate2 6.217762

I will melt these, to plot all at once:

melted <- reshape2::melt(estimates, id.vars = c("Sample")) # nolint use uppercase variable name just like BEAST2

ggplot2::ggplot(
  melted, 
  ggplot2::aes(
    x = melted$Sample, 
    y = melted$value,
    color = as.factor(melted$variable)
  )
) + ggplot2::geom_line() +
  ggplot2::ggtitle("Parameter estimates")

plot of chunk unnamed-chunk-4