bremla_modelfitter: Bremla model fitting

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

Bremla model fitting

Description

Fits the bremla linear regression model to observations.

Usage

bremla_modelfitter(object, control.fit, print.progress = FALSE)

Arguments

object	List object which is the output of function bremla_prepare
control.fit	List containing specifications for fitting procedure. See control.fit.default for details.
print.progress	Boolean. If TRUE progress will be printed to screen

Value

Returns the same object list from the input, but appends information from the fitting procedure in object\$fitting. Including fitted values, residuals, posterior distributions and summary statistics for both fixed and random effects.

Author(s)

Eirik Myrvoll-Nilsen, [e-mail]

See Also

bremla_prepare,bremla_chronology_simulation

Examples

if(inlaloader()){
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)))
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(ncores=2,noise="ar1")
object = bremla_prepare(formula,data,nsims=5000,reference.label="simulated timescale",
                        events = events,
                        control.fit=control.fit)
object = bremla_modelfitter(object, print.progress=TRUE)
summary(object)
plot(object)
}

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