# This script illustrates use of the TimeFRAME package in different model configurations on an example dataset
# The models are run using default values, the TimeFRAME documentation describes in detail
# the parameters that can be changed and edited.
rm(list = ls())
require(lubridate)
library(TimeFRAME)
library(tibble)
library(dplyr)
library(tidyr)
library(ggplot2)
set.seed(1)
# Load the data
data = read.csv("data3.csv")
sd = data[,c("d15N","d15NSP","d18O")]*0+2
# Basic plot of the data
par(mfrow=c(3,1),mar=c(2,4,1,1))
plot(data$t,data$d15N,col="blue",ylab="d15N")
lines(data$t,data$d15N,col="blue")
plot(data$t,data$d15NSP,col="blue",ylab="d15NSP")
lines(data$t,data$d15NSP,col="blue")
plot(data$t,data$d18O,col="blue",ylab="d18O")
lines(data$t,data$d18O,col="blue")
## Independent fits
model <- frame_model(n2o_sources, frac = n2o_frac) # Load the model using default sources and fractionation
fit.independent <- fit_frame(model, data[,c("d15N","d15NSP","d18O")], sd=sd, t = data$t)
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
autoplot(fit.independent) # autoplot to look
## Independent fits, no fD, no d18O
# To run this, we just create the model without these rows/columns in the sources and fractionation
model <- frame_model(n2o_sources[1:2,c(1,2,4,5)], frac = n2o_frac[c(1,2,4,5)])
fit.independent.nofD <- fit_frame(model, data[,c("d15N","d15NSP")], sd=sd[,1:2], t = data$t)
## Warning: There were 25 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
autoplot(fit.independent.nofD) # autoplot to look
## GP prior
model <- frame_model(n2o_sources, frac = n2o_frac)
fit.gp <- fit_gp(model, data[,c("d15N","d15NSP","d18O")], sd=sd, t = data$t)
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
autoplot(fit.gp) # autoplot to look
## DGP
model <- frame_model(n2o_sources, frac = n2o_frac)
fit.dgp <- fit_dgp(model, data[,c("d15N","d15NSP","d18O")], sd=sd, t = data$t)
autoplot(fit.dgp) # autoplot to look
## spline
model <- frame_model(n2o_sources, frac = n2o_frac)
fit.spline <- fit_glm(model, data[,c("d15N","d15NSP","d18O")], sd=sd, t = data$t)
## Warning: The largest R-hat is 1.06, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
autoplot(fit.spline) # autoplot to look
