bremla: Bayesian regression modeling of layered archives

In eirikmn/bremla: Bayesian Regression Modeling of Layer-Counted Archives

Bayesian regression modeling of layered archives

Description

Fits a regression model to the data and produces chronology simulations.

Usage

bremla(
  formula,
  data,
  reference.label = NULL,
  x.label = NULL,
  y.label = NULL,
  nsims = 10000,
  events = NULL,
  synchronization = NULL,
  control.fit = NULL,
  control.sim = NULL,
  control.linramp = NULL,
  control.transition_dating = NULL,
  control.bias = NULL,
  print.progress = FALSE
)

Arguments

formula	formula describing the predictor used in fitting the data. Partial covariates (psi) can be filled in by specifying the degree and events in events argument.
data	data.frame containing response and covariates. Must include 'age' and 'depth'. The first row is used for initializing 'age' and 'depth', hence the rest can be set to NA. Covariates related to the events (psi) can be filled in automatically by specifying the degree and events in events argument.
reference.label	Character label of reference timescale. Used in plot,summary.
x.label	Character label for the x-axis (depth).
y.label	Character label for the y-axis (age).
nsims	Number of chronologies to be simulated.
events	List object describing locations, degree and other informations about the specific climatic periods used in the predictor. If used must include at least events\$locations. See events.default for more details.
synchronization	List containing locations and details for sampling tie-points used to synchronize our time scale.
control.fit	List containing specifications for the fitting procedure. See control.fit.default for more details.
control.sim	List containing specifications for sampling chronologies (synchronous or otherwise). See control.sim.default for more details.
control.linramp	List containing specifications for fitting a linear ramp model to given data. If used, must include control.linramp\$proxy and control.linramp\$interval. See control.linramp.default for details.
control.transition_dating	List containing specifications for estimating a given transition. See control.transition_dating.default for more details.
control.bias	List containing specifications for adding a potential unknown stochastic counting bias. See control.bias.default for more details.
print.progress	Boolean. If TRUE progress will be printed to screen

Value

Returns an S3 object of class "bremla". This includes output from all functions nested within the bremla function. Including fitted marginals and summary statistics, simulated chronologies, time spent on each step.

Author(s)

Eirik Myrvoll-Nilsen, [e-mail]

See Also

bremla_modelfitter,bremla_chronology_simulation

Examples

if(inlaloader()){
## Simulation example
require(stats)
set.seed(1)
n <- 1000
phi <- 0.8
sigma <- 1.2
a_lintrend <- 0.3; a_proxy = 0.8
dy_noise <- as.numeric(arima.sim(model=list(ar=c(phi)),n=n,sd=sqrt(1-phi^2)))
lintrend <- seq(from=10,to=15,length.out=n)

proxy <- as.numeric(arima.sim(model=list(ar=c(0.9)),n=n,sd=sqrt(1-0.9^2)))*5
dy <- a_lintrend*lintrend + a_proxy*proxy + sigma*dy_noise

y0 = 11700;z0=1200
age = y0+cumsum(dy)
depth = 1200 + 1:n*0.05
depth2 = depth^2/depth[1]^2 #normalize for stability


formula = dy~-1+depth2 + proxy
data = data.frame(age=age,dy=dy,proxy=proxy,depth=depth,depth2=depth2)
data = rbind(c(y0,NA,NA,z0,NA),data) #First row is only used to extract y0 and z0.

events=list(locations=c(1210,1220,1240))
control.fit = list(noise="ar1", hyperprior = list(
                                        prec = list(prior="loggamma", param = c(1,0.01)),
                                        rho = list(prior="normal", param = c(0,0.15)))
                                        )
synchronization=list(locations=depth[c(100,400,700)],method="gauss",
                           params=list(mean=c(age[c(100,400,700)]+c(30,-100,50)),
                                       sd=c(50,20,100)
                                       )
                     )
control.sim=list(synchronized=TRUE,
                 summary=list(compute=TRUE))


object = bremla(formula,data,nsims=5000,reference.label="simulated timescale",
                         events=events,
                         synchronization=synchronization,
                         control.fit=control.fit,
                         control.sim=control.sim,
                         print.progress=TRUE)
summary(object)
plot(object)

}

if(inlaloader()){
### Real data example ###
require(stringr)
data("event_intervals")
data("events_rasmussen")
data("NGRIP_5cm")


age = NGRIP_5cm$age
depth = NGRIP_5cm$depth
d18O = NGRIP_5cm$d18O
proxy=d18O
formula = dy~-1+depth2 + proxy
depth2 = depth^2/depth[1]^2 #normalize for stability

data = data.frame(age=age,dy=c(NA,diff(age)),depth=depth,depth2=depth2,proxy=proxy)

lowerints = which.index(event_intervals$depth_int_lower.m, depth[2:length(depth)])
upperints = which.index(event_intervals$depth_int_upper.m, depth[2:length(depth)])

eventnumber=13 #number between 1 and 29. specifies which transition to consider
transitionlabel = str_sub(event_intervals$onsetlabel[eventnumber],
                      str_locate(event_intervals$onsetlabel[eventnumber],"GI")[1])
depth.reference = event_intervals$NGRIP_depth_m[eventnumber]
age.reference = event_intervals$GICC_age.yb2k[eventnumber]
interval = lowerints[eventnumber]:upperints[eventnumber]

nsims=5000
events=list(locations = events_rasmussen$depth,
            locations_unit="depth",degree=1)
synchronization = list(method="adolphi", locations = c(FALSE,TRUE,TRUE,TRUE,TRUE))
control.fit=list(method="inla")
control.sim=list(synchronized = 2)

control.linramp=list(proxy=proxy,interval=interval,interval.unit="index",
                     depth.reference=depth.reference,
                     label=transitionlabel,depth.label="d18O (permil)")
control.transition_dating=list(age.ref=age.reference,age.label="years (yb2k)")
object = bremla(formula,data,reference.label="GICC05",
                   nsims=nsims,
                   events=events,
                   synchronization=synchronization,
                   control.fit=control.fit,
                   control.sim=control.sim,
                   control.linramp=control.linramp,
                   control.transition_dating=control.transition_dating,
                   print.progress=TRUE )
  summary(object)
  plot(object)
  }

eirikmn/bremla documentation built on Jan. 25, 2025, 4:41 a.m.