carstm_model_inla = (
O, # p parameter list
DS = ,
sppoly =,
=,
fit = ,
space_id=, time_id=, cyclic_id=,
fn_fit=(="fit_", fileext=".RDS"),
redo_fit = ,
redo_posterior_simulations = ,
summarize_simulations=,
theta=,
compress="qs-preset",
compression_level=3,
qs_preset="high",
toget = ("modelinfo", "summary", "random_spatial", "predictions"),
nposteriors=0,
posterior_simulations_to_retain=( "summary", "predictions" ),
exceedance_threshold=,
deceedance_threshold=,
exceedance_threshold_predictions=,
deceedance_threshold_predictions=,
=,
ndiscretization = 1024L,
eps = 1e-12,
) {
if (!(DS)) {
if (DS=="modelled_fit") {
if (!(fn_fit)) {
if ((fn_fit)) {
if (("\\.RDS$", fn_fit)) {
fit = aegis::read_write_fast(fn_fit)
} {
( fn_fit )
}
}
if ((fit)) ("Modelled fit was not found.")
}
( fit )
}
if (DS=="carstm_modelled_summary") { # used for look up tables
O =
fn_modelinfo = ( "_fit~", "_modelinfo~", fn_fit, fixed= )
if ((fn_modelinfo)) O = (O, read_write_fast( =fn_modelinfo ) )
fn_summary = ( "_fit~", "_summary~", fn_fit, fixed= )
if ((fn_summary)) O = (O, =read_write_fast( =fn_summary ) )
fn_randomeffects = ( "_fit~", "_randomeffects~", fn_fit, fixed= )
if ((fn_randomeffects)) O = (O, randomeffects=read_write_fast( =fn_randomeffects ) )
fn_preds = ( "_fit~", "_predictions~", fn_fit, fixed= )
if ((fn_preds)) O = (O, predictions=read_write_fast( =fn_preds ) )
fn_samples = ( "_fit~", "_samples~", fn_fit, fixed= )
if ((fn_samples)) O = (O, samples=read_write_fast( =fn_samples ) )
( O )
}
if (DS=="carstm_modelinfo") {
fn_modelinfo = ( "_fit~", "_modelinfo~", fn_fit, fixed= )
O = read_write_fast( =fn_modelinfo )
( O )
}
if (DS=="carstm_summary") {
fn_summary = ( "_fit~", "_summary~", fn_fit, fixed= )
Osummary = read_write_fast( =fn_summary )
( Osummary )
}
if (DS=="carstm_randomeffects") {
fn_randomeffects = ( "_fit~", "_randomeffects~", fn_fit, fixed= )
Orandom = read_write_fast( =fn_randomeffects )
( Orandom )
}
if (DS=="carstm_predictions") {
fn_preds = ( "_fit~", "_predictions~", fn_fit, fixed= )
Opredictions = read_write_fast( =fn_preds )
( Opredictions)
}
if (DS=="carstm_samples") {
fn_samples = ( "_fit~", "_samples~", fn_fit, fixed= )
Osamples = read_write_fast( =fn_samples )
( Osamples )
}
}
if (("debug")) if (()) if () ()
run_start = ()
### 1. Prepare inla_args and vars used in both modelling and extraction
inla_args = () # INLA options to be passed directly to it
# some checks
if (!("verbose", inla_args)) inla_args"verbose"]]=
be_verbose = inla_args"verbose"]]
if (("num.threads", inla_args)) {
num.threads = inla_args"num.threads"]]
} {
num.threads = "1:1"
}
inla.setOption(num.threads=num.threads)
# n cores to use for posterior marginals in mcapply .. can be memory intensive so make it a bit less than "num.threads" ..
if (("mc.cores", inla_args)) {
mc.cores = inla_args"mc.cores"]]
} if ( !(num.threads) ) {
mc.cores = ( ((num.threads, ":") )[1] )
} {
mc.cores = 1
}
# Formula parsing
if ( !("formula", inla_args ) ) {
if (("formula", O ) ) {
inla_args"formula"]] = O"formula"]]
if ( be_verbose ) ( "Formula found in options, O" )
} {
( "Model formula was not provided")
}
}
O"fm"]] = parse_formula( inla_args"formula"]] )
if ( be_verbose ) {
( "Formula parsed as follows. Check in case of errors:" )
(O"fm"]])
}
# labels .. not all used but defining the here makes things simpler below
# shorten var names
fmre = O"fm"]]$random_effects
vO = O"fm"]]$offset_variable
vY = O"fm"]]$dependent_variable
vS = O"fm"]]$vn$S
vT = O"fm"]]$vn$T
vU = O"fm"]]$vn$U
vS2 = O"fm"]]$vn$S2
vT2 = O"fm"]]$vn$T2
vU2 = O"fm"]]$vn$U2
vS3 = O"fm"]]$vn$
# family related O
if ( !("family", inla_args ) ) {
if ( ("family", O) ) {
if ( be_verbose ) ( "Family found in options, O" )
inla_args"family"]] = O"family"]]
} {
inla_args"family"]] = "gaussian"
}
}
if (!("family", inla_args)) ( "Model family was not provided")
O"family"]] = inla_args"family"]] # copy to O in case it was provided as inla_args
if ( inla_args"family"]] == "gaussian" ) {
lnk_function = inla.link.identity
lnk_function_betas = inla.link.identity
} if ( inla_args"family"]] == "lognormal" ) {
# inla uses identity link and simply transforms Y :: log(X)~normaLl()
# however to back-transform properly to user scale, we override with a log-exp link
# need to manually intervene in a few calcs below (with identity link)
lnk_function = inla.link.log
lnk_function_betas = inla.link.log
} if ( ( ".*poisson", inla_args"family"]])) {
lnk_function = inla.link.log
lnk_function_betas = inla.link.log
} if ( ( ".*nbinomial", inla_args"family"]])) {
lnk_function = inla.link.log
lnk_function_betas = inla.link.log
invlink_id = "exp"
} if ( ( ".*binomial", inla_args"family"]])) {
lnk_function = inla.link.logit
lnk_function_betas = inla.link.log
}
O"invlink"]] = invlink = (x) lnk_function( x, inverse= )
O"invlink_betas"]] = invlink_betas = (x) lnk_function_betas( x, inverse= )
# changes in INLA data structures noted in 2023:
# posterior samples is on link scale ((for both experimental and classical))
# marginals.fitted.values are on response scale (for both experimental and classical)
# summary.fitted.values on response scale ((for both experimental and classical))
# 2024: no more experimental mode .. just use default
# --------------------------------------------
# --------------------------------------------
# --------------------------------------------
### START Modelling
if (redo_fit) {
# permit passing a function rather than data directly .. less RAM usage in parent call
if (!("data", inla_args)) {
if (("data", O)) {
if ( be_verbose) ( "Data not passed as an argument, using data found in options, O")
inla_args"data"]] = O"data"]] # priority to inla_args
O"data"]] = # reduce mem usage
}
}
if (!("data", inla_args)) {
if (!()) {
inla_args"data"]] =
=
}
}
if (((inla_args"data"]])=="character")) {
if ( be_verbose) ( "Data is a function. Getting data ...")
Odata = ( ((=inla_args"data"]]) ) )
inla_args"data"]] = Odata
Odata =
}
if ((inla_args"data"]], "try-error")) {
inla_args"data"]] =
}
if ((inla_args"data"]])) ("Data not found")
setDT(inla_args"data"]])
if ( !("carstm_model_label", O) ) {
if ( ("model_label", O) ) {
O"carstm_model_label"]] = O"model_label"]]
} {
("carstm_model_label should be specified")
}
}
DM = ( fmre$dimensionality )
dims = ("s", "t", "c")
RR = ( dims, ( dims, DM ) )
SS = (
("space", "time", "cyclic")[match(RR, dims)],
collapse="-"
)
if ( ("dimensionality", O) ) {
if (O"dimensionality"]] != SS ) {
("Dimensionality parameters specified (left) do not match those in formula (right): " )
( O"dimensionality"]], " vs ", SS )
( "This is OK when predicting to one time slice, but verify that this is your intent ...")
}
} {
O"dimensionality"]] = SS
("Check dimensionality guessed from formula: ", O"dimensionality"]])
}
if ( !(O"fm"]]"vn"]]"S"]]) ) {
# this is about model data inputs not dimensionality of outputs
if ((sppoly)) if (("sppoly", O)) sppoly = O$sppoly
if ((sppoly)) if (("areal_units")) {
sppoly = areal_units( O )
if (!(sppoly)) sppoly =
}
if ((sppoly)) ( "sppoly is required")
# test to see if sppoly has been altered:
sp_nb_auid =
sp_nb_auid = ( ((sppoly)$NB_graph)$region.id, silent=) # order of neighbourhood graph
if (!("try-error", sp_nb_auid)) {
o = ( sppoly$AUID, sp_nb_auid)
if ((((o))) !=1) {
( "Neighbourhood indices and polygons AUID order are mismatched .. this will not end well" )
}
o =
O"space_id"]] = sp_nb_auid
sp_nb_auid =
}
O"sppoly"]] = sppoly # copy in case mapping directly from O
# the master / correct sequence of the AU's and neighbourhood matrix index values
if (!(space_id)) {
O"space_id"]] = space_id
space_id =
} {
if (("space_id", (sppoly)) ) {
O"space_id"]] = ( (sppoly)$space_id )
} if (("region.id", (sppoly)) ) {
O"space_id"]] = ( (sppoly)$region.id )
} if (("AUID", sppoly) ) {
O"space_id"]] = ( sppoly"AUID"]] )
}
}
if (!("space_id", O)) ( "space_id could not be determined from data")
if (!("space_id", (sppoly)) ) (sppoly)$space_id = O"space_id"]] # copy as attribute in case
O"space_n"]] = ( O "space_id" ]] )
# better formatted areal unit names (for reporting or plotting) (outside of carstm))
if (!("space_name", O)) {
if (("space_name", (sppoly))) {
O"space_name"]] = (sppoly)$space_name
} {
O"space_name"]] = O"space_id"]]
}
}
if (!("space_name", (sppoly)) ) (sppoly)$space_name = O"space_name"]] # copy as attribute in case
if (("AUID_label", sppoly)) {
# long name useful in some cases (for plotting)
O"space_label"]] = sppoly"AUID_label"]]
if (!("space_label", (sppoly)) ) (sppoly)$space_label = O"space_label"]] # copy as attribute in case
}
if ( (O"fm"]]"vn"]]"T"]]) ) {
# force to carry a "time" to make dimensions of predictions simpler to manipulate
inla_args"data"]]"time"]] = -1
}
missingS = ((inla_args"data"]][vS]] ))
if (( missingS ) > 0 ) {
( "Data areal units and space_id (from sppoly) do not match ... this should not happen:")
( (missingS) )
( "Dropping them from analysis and assuming neighbourhood relations are OK even though this is unlikely!")
inla_args"data"]] = inla_args"data"]] -missingS , ]
}
missingS =
missingS = ( ( O"space_id"]], ( inla_args"data"]][vS]] ) ) )
if ((missingS) > 0) {
( "No. of areal unique units in data do not match those in sppoly:", ( (missingS), sep=", "))
}
missingS =
}
if ( !(O"fm"]]"vn"]]"T"]]) | !(O"fm"]]"vn"]]"U"]]) ) {
# carstm internal ID
if (!(time_id)) {
O"time_id"]] = time_id
time_id =
} {
if (("yrs", O)) {
if ((O"yrs"]])) {
O"time_id"]] = (O"yrs"]])
} {
O"time_id"]] = (O"yrs"]]) # redundant but make explicit
}
}
}
if (!("time_id", O)) ( "time_id or yrs need to be provided")
O"time_n"]] = ( O "time_id" ]] )
# for plot labels, etc .. time (year) gets swapped out for time index (outside of carstm)
if (!("time_name", O)) {
if (("yrs", O)) {
O"time_name"]] = (O "yrs" ]] )
} {
O"time_name"]] = (O "time_id" ]] )
}
}
# sub-annual time
if (!(O"fm"]]"vn"]]"U"]]) ) {
# carstm internal ID
if (!(cyclic_id)) {
O"cyclic_id"]] = cyclic_id
cyclic_id =
} {
if (("cyclic_levels", O)) {
if ((O"cyclic_levels"]])) {
O"cyclic_id"]] = (O"cyclic_levels"]])
} {
O"cyclic_id"]] = O"cyclic_levels"]]
}
}
}
if (!("cyclic_id", O)) ( "cyclic_id or cyclic_levels need to be provided")
O"cyclic_n"]] = ( O "cyclic_id" ]] )
# for plot labels, etc .. season (time of year) gets swapped out for cyclic index (outside of carstm)
if (!("cyclic_name", O)) {
if (("cyclic_levels", O)) {
O"cyclic_name"]] = (O "cyclic_levels" ]] )
} {
O"cyclic_name"]] = (O "cyclic_id" ]] )
}
}
}
}
ii = ((inla_args"data"]][vY]]))
if ( be_verbose ) ("Prepping data ")
# dev.new()
# hist( inla_args[["data"]][[vY]][ii] , main="Histogram of input variable to model" )
mq = ( inla_args"data"]][vY]][ii] , probs=(0.025, 0.5, 0.975) )
O"data_range"]] = (
=(inla_args"data"]][vY]][ii] ),
=(inla_args"data"]][vY]][ii] ),
=(inla_args"data"]][vY]][ii] ),
=(inla_args"data"]][vY]][ii] ),
lb=mq[1],
=mq[2],
ub=mq[3]
) # on data /user scale not internal link
mq =
# prefilter/transformation (e.g. translation to make all positive)
if (("data_transformation", O)) inla_args"data"]][vY]] = O$data_transformation$forward( inla_args"data"]][vY]] )
# get hyper param scalings
# temp Y var on link scale:
# on user scale
yl = inla_args"data"]][vY]]
if ( ( ".*binomial", inla_args"family"]])) {
# for binomial, prob=0,1, becomes infinite so minor fix for hyper parameter prior approximation
tweak = 0.05 # tail truncation prob
if (("habitat.threshold.quantile", O)) tweak = O"habitat.threshold.quantile"]]
yl yl==1 ] = 1 - tweak
yl yl==0 ] = tweak
}
yl = lnk_function( yl ) # necessary in case of log(0)
# offsets need to be close to 1 in user scale ( that is log(1)==0 in internal scale ) in experimental mode .. rescale
if ( !(vO) ) {
# link function applied to offsets here .. do not need to send log()
if (("log[[:space:]]*[(]", vO)) ("Probably do not want to transform the offset .. it is done internally in , unlike glm, inla, etc")
obs = 1:(inla_args"data"]])
if (("tag", inla_args"data"]])) {
obso = (inla_args"data"]]"tag"]]=="observations")
if ((obso) > 3) obs = obso
obso =
}
inla_args"data"]][vO]] = lnk_function( inla_args"data"]][vO]])
obs =
yl = yl - inla_args"data"]][vO]]
}
ll = ((yl))
mqi = ( yl[ll], probs=(0.025, 0.5, 0.975) )
O"data_range_internal"]] = (
=(yl[ll]),
=(yl[ll]),
=(yl[ll]),
=(yl[ll]),
lb=mqi[1],
=mqi[2],
ub=mqi[3]
) # on data /user scale not internal link
mqi =
O"predictions_range"]] = ( ( inla_args"data"]]"tag"]] == "predictions" ) )
O"priors"]] = H = inla_hyperparameters(
reference_sd = O"data_range_internal"]]"sd"]],
alpha=0.5,
O"data_range_internal"]]"median"]]
)
# sd biased (high) due to 0's being dropped .. but we use pc.priors which shrink it to 0, so is ok as more robust ... :)
m = yl = ii = ll = fy = ol =
()
outputdir = (fn_fit)
if ( !(outputdir)) ( outputdir, recursive=, showWarnings= )
if (("carstm_prediction_surface_parameters", O)) O"carstm_prediction_surface_parameters"]] =
# make these temporary indices here to drop inla_args early and reduce RAM usage and make things easier later
O"ipred"]] = ( inla_args"data"]]"tag"]]=="predictions" )
if ( O"dimensionality"]] == "space" ) {
# filter by S and T in case additional data in other areas and times are used in the input data
O"matchfrom"]] = (
space=O"space_id"]][inla_args"data"]][vS]][O"ipred"]]]]
)
}
if (O"dimensionality"]] == "space-time" ) {
O"matchfrom"]] = (
space=O"space_id"]][inla_args"data"]][vS]][O"ipred"]]]] ,
=O"time_id"]][inla_args"data"]][vT]][O"ipred"]]] ]
)
}
if ( O"dimensionality"]] == "space-time-cyclic" ) {
O"matchfrom"]] = (
space=O"space_id"]][inla_args"data"]][vS]][O"ipred"]]]] ,
=O"time_id"]][inla_args"data"]][vT]][O"ipred"]]]],
cyclic=O"cyclic_id"]][inla_args"data"]][vU]][O"ipred"]]]]
)
}
fn_modelinfo = ( "_fit~", "_modelinfo~", fn_fit, fixed= )
read_write_fast( =O, =fn_modelinfo, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
( "Model info saved as: \n", fn_modelinfo )
()
if (("debug")) if (()) if ( =="fit") ()
# access H, cyclic_values, etc
inla_args".parent.frame"]]=()
# check INLA options
if (!("control.inla", inla_args)) inla_args"control.inla"]] = (strategy = "auto", optimiser="default")
if (!("control.predictor", inla_args)) inla_args"control.predictor"]] = ()
# these are required to know the scale of predictions
inla_args"control.predictor"]]$compute = # force
inla_args"control.predictor"]]$link = 1 # force
if (!("control.mode", inla_args ) ) inla_args"control.mode"]] = ( = )
if (!(theta) ) inla_args"control.mode"]]$theta= theta
if (!("control.compute", inla_args)) inla_args"control.compute"]] = (dic=, waic=, cpo=, config=, return.marginals.predictor= )
# if (!exists("control.fixed", inla_args)) inla_args[["control.fixed"]] = list(mean.intercept=0.0, prec.intercept=0.001, mean=0, prec=0.001)
if (!("control.fixed", inla_args)) inla_args"control.fixed"]] = H$fixed
setDF(inla_args"data"]]) # in case .. INLA requires this ?
fit = ( ( inla, inla_args ) )
if ((fit, "try-error" )) {
inla_args"control.inla"]]$strategy="laplace"
inla_args"control.inla"]]$optimiser="gsl"
inla_args"control.inla"]]$=1
inla_args"control.inla"]]$cmin=0
inla_args"control.inla"]]$diagonal=1e-8
fit = ( ( inla, inla_args ) )
}
if ((fit, "try-error" )) {
inla_args"control.inla"]]$strategy="laplace"
inla_args"control.inla"]]$optimiser="gsl"
inla_args"control.inla"]]$=1
inla_args"control.inla"]]$h = 0.00075
inla_args"control.inla"]]$step.factor = 0.0075
inla_args"control.inla"]]$cmin=0
inla_args"control.inla"]]$diagonal=1e-8
fit = ( ( inla, inla_args ) )
}
if ((fit, "try-error" )) {
inla_args"control.inla"]]$strategy="laplace"
inla_args"control.inla"]]$optimiser="gsl"
inla_args"control.inla"]]$h = 0.00125
inla_args"control.inla"]]$step.factor = 0.0125
inla_args"control.inla"]]$cmin=0
inla_args"control.inla"]]$diagonal=1e-6
fit = ( ( inla, inla_args ) )
}
if ((fit, "try-error" )) {
inla_args"safe"]] =
inla_args"control.inla"]]$h =
inla_args"control.inla"]]$step.factor =
inla_args"control.inla"]]$cmin = 0
inla_args"control.inla"]]$diagonal = 1e-6
fit = ( ( inla, inla_args ) )
}
if ((fit, "try-error" )) {
inla_args"safe"]] =
inla_args"control.inla"]]$h =
inla_args"control.inla"]]$step.factor =
inla_args"control.inla"]]$cmin = 0
inla_args"control.inla"]]$diagonal = 1e-6
inla_args"control.inla"]]$int.strategy='eb'
fit = ( ( inla, inla_args ) )
}
if ((fit, "try-error" )) ( "Model fit error" )
# to improve hyper param estimates..
if (("improve.hyperparam.estimates", O)) {
if (O"improve.hyperparam.estimates"]]) {
fit = inla.hyperpar(fit, dz=0.6, diff.logdens=9 ) # get improved estimates for the hyperparameters
}
}
if ((fit)) ("model fit error")
if ("try-error" %in% (fit) ) ("model fit error")
setDT(inla_args"data"]]) # revert to DT for faster / efficient operations
fit$.args = # reduce size
inla_args= ;
()
read_write_fast( =fit, =fn_fit, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
if (be_verbose) ( "\nModel fit saved as: \n", fn_fit )
} # END redo fit
run_fit_end = ()
### END Modelling
# --------------------------------------------
# --------------------------------------------
# --------------------------------------------
if ((fit)) {
if (be_verbose) ( "\nLoading fit as: redo_fit=FALSE ...\n", fn_fit )
fit =aegis::read_write_fast( fn_fit )
}
if ((fit)) {
( "Fit file not found, set option: redo_fit=TRUE" )
()
}
fn_modelinfo = ( "_fit~", "_modelinfo~", fn_fit, fixed= )
if ((fn_modelinfo)) {
O = aegis::read_write_fast( fn_modelinfo )
} {
("modelinfo not found")
}
if (("debug")) if (()) if (=="extract") ()
if (("data_transformation", O)) {
backtransform = ( b ) {
b,1] = O$data_transformation$backward( b,1] )
( b )
}
}
list_simplify = (x) ( ( ( x )))
exceedance_prob = (x, threshold) {1 - inla.pmarginal( = threshold, x)}
deceedance_prob = (x, threshold) { inla.pmarginal( = threshold, x)}
# local functions
list_to_dataframe = (Y ){
if (!(Y)) {
Y = (Y)
}
Z = ((Y, (x) (, x)))
(Z) = (Y)
(Z) = (Y)
(Z)
}
apply_generic = () (, mc.cores=mc.cores ) # drop-in for lapply
apply_simplify = () ((, mc.cores=mc.cores ), higher = ) # drop in for sapply
# not used, for reference
apply_generic_serial = () ( ) # drop-in for lapply .. serial
apply_simplify_serial = () (( )) # drop in for sapply .. serial
sqrt_safe = ( a, eps=eps ) ( ( ( a, eps ), 1/eps ) )
marginal_clean = ( w) {
i <- (!((w)) )
if ( (i) > 0) w = w-i,]
w = w[order(w,1]),]
(w)
}
test_for_error = ( Z ) {
if ( "try-error" %in% (Z) ) ("error")
if ((Z)) {
if (( ((Z, (o) (o, "try-error"))) )) ("error")
} if ((Z) ){
if (( (Z, "try-error"))) {
("error")
} if ((("error", m, ignore.case =))) {
("error")
}
}
( "good" )
}
marginal_summary = (Z, invlink= ) {
if (!(invlink)) {
Z = ( apply_generic( Z, inla.tmarginal, fun=invlink ) )
if (test_for_error(Z) =="error") {
Z = ( apply_generic( Z, inla.tmarginal, fun=invlink, n=ndiscretization, method = "linear" ) )
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
(m) = "try-error"
(m)
}
Z = ( apply_generic( Z, marginal_clean ) )
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
Z = ( apply_generic( Z, inla.zmarginal, silent= ) )
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
Z = ( ( Z ), silent=)
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
Z = ( list_simplify( Z), silent= )
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
(Z)
} {
Z = ( apply_generic( Z, inla.zmarginal, silent= ), silent=)
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
Z = ( ( Z ), silent=)
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
Z = ( list_simplify( Z ), silent=)
if (test_for_error(Z) =="error") {
(m) = "try-error"
(m)
}
(Z)
}
}
inla_tokeep = ("mean", "sd", "0.025quant", "0.5quant", "0.975quant")
tokeep = ("mean", "sd", "quant0.025", "quant0.5", "quant0.975")
# some info required for posterior simulations
O"names.fixed"]] = fit$names.fixed
if (!("nposteriors", O)) O"nposteriors"]] = nposteriors
if (("debug")) if (()) if (=="summary") ()
if ( "summary" %in% toget) { # add marginals to summary
Osummary = ()
Osummary"all.hypers"]] = INLA:::inla.all.hyper.postprocess(fit$all.hyper)
Osummary"hypers"]] = fit$marginals.hyperpar
Osummary"fixed"]] = fit$marginals.fixed
Osummary"dic"]] = fit$dic[c("dic", "p.eff", "dic.sat", "mean.deviance")]
Osummary"waic"]] = fit$waic[c("waic", "p.eff")]
Osummary"mlik"]] = fit$mlik[2]
Osummary"direct"]] = (fit)
# print(Osummary[["direct"]])
# remove a few unused but dense data objects
Osummary"direct"]]$linear.predictor =
Osummary"direct"]]$waic$local.waic =
Osummary"direct"]]$waic$local.p.eff =
Osummary"direct"]]$dic$local.dic =
Osummary"direct"]]$dic$local.p.eff =
Osummary"direct"]]$dic$local.dic.sat =
Osummary"direct"]]$dic$ =
if (be_verbose) ("Extracting parameter summary" )
# extract and back transform where possible
if (( "marginals.fixed", fit)) {
m = fit$marginals.fixed # make a copy to do transformations upon
fi = ("Intercept", (m) )
m = ( apply_generic( m, (x) marginal_clean( inla.tmarginal( invlink_betas, x)) ), silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, (x) marginal_clean( inla.tmarginal( invlink_betas, x, n=ndiscretization, method="linear" )) ), silent= )
}
if ( ("data_transformation", O)) {
m[fi]] = inla.tmarginal( O$data_transformation$backward, m[fi]] , n=ndiscretization ) # on user scale
}
m = (apply_simplify( m, FUN=inla.zmarginal, silent= ), silent=)
if (test_for_error(m) =="error") {
("Problem with marginals for intercept (probably too variable), doing a simple inverse transform, SD is blanked out")
W = fit$summary.fixed "(Intercept)", 1:5]
(W) = tokeep
W,(1,3:5)] = invlink_betas( W,(1,3:5)] )
W,2] =
} {
W = cbind ( t (m) ) #
W = list_to_dataframe( W , tokeep, =] )
}
W$ID = (W)
Osummary"fixed_effects"]] = W
m =
W =
()
}
if (( "marginals.hyperpar", fit)) {
# hyperpar (variance components)
hyps = (fit$summary.hyperpar)
prcs = ( "^Precision.*", hyps, value= )
if ((prcs) > 0) {
m = fit$marginals.hyperpar[prcs]
m = ( apply_generic( m, inla.tmarginal, fun=(y) 1/sqrt_safe( y, eps ) ), silent=)
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=(y) 1/sqrt_safe( y, eps ), n = ndiscretization, method="linear" ), silent= )
}
m = ( apply_generic( m, inla.zmarginal, silent= ), silent=)
m = ( ( m ), silent=)
if (test_for_error(m) =="error") {
if (be_verbose) {
( "NAN or Inf values encountered in marginals. ")
( "Try an alternate parameterization as model may be stuck in local optima or boundary conditions. ")
( "Copying fit summaries directly rather than from marginals ... ")
}
m = fit$summary.hyperpar[prcs,1:5]
m,(1,3:5)] = 1/sqrt_safe( m,(1,3:5)], eps )
(m) = tokeep
}
m = ( list_simplify( m), silent=)
(m) = ("Precision for", "SD", (m) )
(m) = (" for", "", (m) )
(m) = (" the", "", (m) )
Osummary"random_effects"]] = m, tokeep, =]
m =
()
}
# update phi's, lambda's (used in besagproper2 -- Leroux model) .. etc
rhos = ( "^Rho.*|^GroupRho.*", hyps, value= )
phis = ( "^Phi.*", hyps, value= )
other = ( "^Lambda.*|^Diagonal.*|zero-probability.*", hyps, value= )
known = ( rhos, phis, other )
unknown = ( hyps, (prcs, known) )
hyps = prcs = other = known =
if ((rhos) > 0) {
m = marginal_summary( fit$marginals.hyperpar rhos ], invlink= )
if (test_for_error(m) =="error") {
m = fit$summary.hyperpar[rhos, 1:5]
(m) = tokeep
}
Osummary"random_effects"]] = ( Osummary"random_effects"]], m, tokeep, =] )
m =
}
rhos =
if ((phis) > 0) {
m = marginal_summary( fit$marginals.hyperpar phis ], invlink= )
if (test_for_error(m) =="error") {
m = fit$summary.hyperpar[phis, 1:5]
(m) = tokeep
}
Osummary"random_effects"]] = ( Osummary"random_effects"]], m, tokeep, =] )
m =
}
phis =
if ((unknown) > 0) {
m = marginal_summary( fit$marginals.hyperpar unknown ], invlink= ) # default to no transform
if ((m) && m=="error") {
# alternatively: m[,"mode"] = apply_simplify( fit$marginals.hyperpar[ unknown ], FUN=function(x) inla.mmarginal( x ))
m = fit$summary.hyperpar[unknown, 1:5]
(m) = tokeep
}
Osummary"random_effects"]] = ( Osummary"random_effects"]], m, tokeep, =])
m =
}
unknown =
Osummary"random_effects"]] = list_to_dataframe( Osummary"random_effects"]] )
Osummary"random_effects"]]$ID = ( Osummary"random_effects"]] )
()
if ((fit"marginals.random"]]) > 0) {
if (("debug")) if (()) if ( =="random_covariates") ()
raneff = ( ( fit$marginals.random ), (vS, vS2, vS3 ) )
for (rnef raneff) {
m = marginal_summary( fit$marginals.random[rnef]], invlink=invlink_betas )
if (test_for_error(m) =="error") {
( "failed to transform marginals .. copying directly from INLA summary instead: ", rnef)
m = fit$summary.random[rnef]], inla_tokeep, = ]
(m) = tokeep
Osummary[rnef]] = m
} {
Osummary[rnef]] = m, tokeep, =]
}
Osummary[rnef]]$ID = fit$summary.random[rnef]]$ID
Osummary[rnef]] = list_to_dataframe( Osummary[rnef]] )
}
m = raneff =
()
}
}
if (be_verbose) {
# message( "")
# message( "Random effects:")
# print( Osummary[["random_effects"]] )
# message( "\n--- NOTE --- 'SD *' from marginal summaries are on link scale")
# message( "--- NOTE --- SD * from posteriors simulations are on user scale")
# message( "")
}
# save summary
fn_summary = ( "_fit~", "_summary~", fn_fit, fixed= )
read_write_fast( =Osummary, =fn_summary, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
Osummary =
()
} # end parameters
# random effects (random spatial and random spatiotemporal) from INLA's marginals
Orandom = ()
# copy a few items to keep each output autonomous
Orandom$sppoly = sppoly
Orandom$fm = O$fm
if (("time_name", O)) Orandom$time_name = O$time_name
if (("space_name", O)) Orandom$space_name = O$space_name
if (("cyclic_name", O)) Orandom$cyclic_name = O$cyclic_name
# separate out random spatial and randomm spatiotemporal (as they can be large arrays)
# NOTE ::: marginals and summary of fitted.values are on the response scale
# ... BUT .. random effects are not and require inverse transform (check with inst/scripts/test_inla.R)
if ("random_spatial" %in% toget) {
# space only
Z =
iSP = ( fmre$dimensionality=="s" & fmre$level=="main")
if ((iSP) > 0 ) {
if (be_verbose) ("Extracting random spatial effects" )
if (("debug")) if (()) if ( =="random_spatial") ()
matchto = ( space=O"space_id"]] )
W = ( , =( O"space_n"]], (tokeep) ), =( space=O"space_name"]], stat=tokeep ) )
((W))[1] = vS # need to do this in a separate step ..
Orandom[vS]] = () # space as a main effect vS==vS
if ((iSP) == 1) {
# vS = fmre$vn[ iSP ] # == vS as it is a single spatial effect
model_name = fmre$model iSP ] # iid (re_unstructured)
m = fit$marginals.random[vS]]
# NOTE ::: marginals and summary of fitted.values are on the response scale
# ... BUT .. random effects are not and require inverse transform (check with inst/scripts/test_inla.R)
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas ) , silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas, n=ndiscretization, method="linear" ) , silent= )
}
m = ( apply_generic( m, inla.zmarginal, silent= ), silent= )
m = ( ( m ), silent=)
m = ( list_simplify( m ), silent= )
# single spatial effect (eg in conjuction with besag) .. indexing not needed but here in case more complex models ..
if (test_for_error(m) =="error") {
( "failed to transform random_spatial marginals .. copying directly from INLA summary instead .. sd's will not be correct")
m = fit$summary.random[vS]], inla_tokeep ]
(m) = tokeep
}
if ( model_name %in% ("bym", "bym2") ) {
# bym2 effect is coded by INLA as a double length vector: re_total and re_neighbourhood
Z = expand.grid( space=O"space_id"]], = ("re_total", "re_neighbourhood"), stringsAsFactors = )
# extract re_total main effects
ire = (Z$=="re_total")
matchfrom = ( space=Z"space"]][ire] )
for (k 1:(tokeep)) {
W,k] = reformat_to_array( input = (m[ire, tokeep[k]]), matchfrom=matchfrom, matchto=matchto )
}
Orandom[vS]] "re_total"]] = W , tokeep, =]
} {
# single spatial effect that is not bym2
# this is redundant with iSP being a single factor, but approach is generalizable for higher dims
Z = expand.grid( space=O"space_id"]], ="re_neighbourhood", stringsAsFactors = )
}
ine = (Z$=="re_neighbourhood")
matchfrom = ( space=Z"space"]][ine] )
W] =
for (k 1:(tokeep)) {
W,k] = reformat_to_array( input = (m[ine, tokeep[k]]), matchfrom=matchfrom, matchto=matchto )
}
Orandom[vS]] "re_neighbourhood"]] = W , tokeep, =]
}
if ((iSP) == 2) {
for (j 1:(iSP)) {
vnS = fmre$vn iSP[j] ]
model_name = fmre$model iSP[j] ]
m = fit$marginals.random[vnS]]
# NOTE ::: marginals and summary of fitted.values are on the response scale
# ... BUT .. random effects are not and require inverse transform (check with inst/scripts/test_inla.R)
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas ), silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas, n=ndiscretization, method="linear"), silent= )
}
m = ( apply_generic( m, inla.zmarginal, silent=) , silent= )
m = ( ( m ), silent=)
m = ( list_simplify(m), silent= )
if (test_for_error(m) =="error") {
( "failed to transform random_spatial marginals .. copying directly from INLA summary instead.. sd's will not be correct")
m = fit$summary.random[vnS]], inla_tokeep ]
(m) = ("ID", tokeep)
}
# single spatial effect (eg besag, etc)
matchfrom = ( space=O"space_id"]] )
W] =
for (k 1:(tokeep)) {
W,k] = reformat_to_array( input = (m, tokeep[k]]), matchfrom=matchfrom, matchto=matchto )
}
Orandom[vnS]] [model_name]] = ( W , tokeep, =], ID=(W) )
}
}
Z = m = matchfrom =
()
}
} # end random spatial effects
matchfrom = i1 = i2=
Z = W = m = space = space1 = space2 = skk1 = skk2 = iSP =
()
if ("random_spatiotemporal" %in% toget ) {
# space-time
iST = ( fmre$dimensionality=="st" & fmre$level=="main")
if ((iST) > 0 ) {
if (be_verbose) ("Extracting random spatiotemporal effects" )
if (("debug")) if (()) if ( =="random_spatiotemporal") ()
matchto = ( space=O"space_id"]], =O"time_id"]] )
W = ( ,
=( O"space_n"]], O"time_n"]], (tokeep) ),
=( space=O"space_name"]], =O"time_name"]], stat=tokeep ) )
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
if ((iST) == 1) {
vnST = fmre$vn iST ]
model_name = fmre$model iST ]
if ((vnST, fit$marginals.random )) {
m = fit$marginals.random[vnST]]
# NOTE ::: marginals and summary of fitted.values are on the response scale
# ... BUT .. random effects are not and require inverse transform (check with inst/scripts/test_inla.R)
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas ) , silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas, n=ndiscretization, method="linear" ) , silent= )
}
m = ( apply_generic( m, inla.zmarginal, silent=), silent= )
m = ( ( m ), silent=)
m = ( list_simplify( m ), silent= )
if (test_for_error(m) =="error") {
( "!!! Failed to transform random_spatial marginals .. copying directly from INLA summary instead.. sd's will not be correct")
m = fit$summary.random[vnST]], inla_tokeep ]
(m) = tokeep
}
if ( model_name %in% ("bym", "bym2") ) {
# bym2 effect: re_total and re_neighbourhood with annual results
Z = expand.grid( space=O"space_id"]], = ("re_total", "re_neighbourhood"), =O"time_id"]], stringsAsFactors = )
# spatiotemporal interaction effects
ire = (Z$=="re_total")
matchfrom = ( space=Z"space"]][ire], =Z"time"]][ire] )
for (k 1:(tokeep)) {
W,,k] = reformat_to_array( input = (m[ire, tokeep[k]]), matchfrom = matchfrom, matchto = matchto )
}
Orandom[vnST]] "re_total"]] = W ,, tokeep, =]
} {
# besag effect: with annual results
Z = expand.grid( space=O"space_id"]], ="re_neighbourhood", =O"time_id"]], stringsAsFactors = )
}
# spatiotemporal interaction effects
ine = (Z$=="re_neighbourhood")
matchfrom = ( space=Z"space"]][ine], =Z"time"]][ine] )
for (k 1:(tokeep)) {
W,,k] = reformat_to_array( input = (m[ine,tokeep[k]]), matchfrom=matchfrom, matchto=matchto )
}
Orandom[vnST]] "re_neighbourhood"]] = W ,, tokeep, =]
}
}
if ((iST) == 2) {
for (j 1:(iST)) {
vnST = fmre$vn iST[j] ]
model_name = fmre$model iST[j] ]
m = fit$marginals.random[vnST]]
# NOTE ::: marginals and summary of fitted.values are on the response scale
# ... BUT .. random effects are not and require inverse transform (check with inst/scripts/test_inla.R)
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas ), silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=invlink_betas, n=ndiscretization, method="linear"), silent= )
}
m = ( apply_generic( m, inla.zmarginal, silent= ), silent= )
m = ( ( m ), silent=)
m = ( list_simplify( m ), silent= )
if (test_for_error(m) =="error") {
( "!!! Failed to transform random_spatiotemporal marginals .. copying directly from INLA summary instead.. sd's will not be correct")
m = fit$summary.random[vnST]], inla_tokeep ]
(m) = tokeep
}
Z = expand.grid( space=O"space_id"]], =model_name, =O"time_id"]], stringsAsFactors = )
jst = (Z$==model_name)
matchfrom = ( space=Z"space"]][jst], =Z"time"]][jst] )
for (k 1:(tokeep)) {
W,,k] = reformat_to_array( input = (m[jst, tokeep[k] ]), matchfrom = matchfrom, matchto = matchto )
}
Orandom[vnST]] [model_name]] = W ,, tokeep, =]
m =
}
}
Z = W = m = space_time = space_time1 = space_time2 = skk1 = skk2 =
}
} # end random spatio-temporal effects
# save random effects to file
if ((Orandom) > 0) {
fn_randomeffects = ( "_fit~", "_randomeffects~", fn_fit, fixed= )
read_write_fast( =Orandom, =fn_randomeffects, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
Orandom =
()
}
Opredictions = ()
# copy a few items to keep each output autonomous
Opredictions$sppoly = sppoly
Opredictions$fm = O$fm
if (("time_name", O)) Opredictions$time_name = O$time_name
if (("space_name", O)) Opredictions$space_name = O$space_name
if (("cyclic_name", O)) Opredictions$cyclic_name = O$cyclic_name
if ("predictions" %in% toget ) {
# **************** README
# marginals.fitted.values are on **response scale** and already incorporates offsets
# summary.fitted.values on **response scale**
# posterior samples are on **link scale** and already incorporates offsets
# **************** README
if (be_verbose) ("Extracting posterior predictions" )
if (("debug")) if (()) if ( =="predictions") ()
if (("marginals.fitted.values", fit)) {
if ((fit"marginals.fitted.values"]]) > 0 ) {
m = fit$marginals.fitted.values[O"ipred"]]] # on **response scale** & already incorporates offsets
if (("data_transformation", O)) m = apply_generic( m, backtransform )
if (O"family"]] == "lognormal") {
# only lognormal needs inverse link due to how inla treats it
m = ( apply_generic( m, inla.tmarginal, fun=invlink ), silent= )
if (test_for_error(m) =="error") {
m = ( apply_generic( m, inla.tmarginal, fun=invlink, n=ndiscretization, method="linear"), silent= )
}
}
m = ( apply_generic( m, inla.zmarginal, silent= ), silent=)
m = ( list_simplify( ( m ) ), silent=)
if (test_for_error(m) =="error") {
if (O"family"]] == "lognormal") {
( "The results are on log-scale and will need to be back-transformed")
}
m = fit$summary.fitted.values[O"ipred"]], inla_tokeep ] # on **response scale** & already incorporates offsets
(m) = tokeep
}
if ( O"dimensionality"]] == "space" ) {
W = ( ,
=( O"space_n"]], ((m)) ),
=( space=O"space_name"]], stat=(m) ) )
((W))[1] = vS # need to do this in a separate step ..
matchto = ( space=O"space_id"]] ) # matchfrom already created higher up
for (k 1:((m))) {
W,k] = reformat_to_array( input=(m,k]), matchfrom=O"matchfrom"]], matchto=matchto )
}
Opredictions"predictions"]] = W, tokeep, =]
m = W =
}
if (O"dimensionality"]] == "space-time" ) {
W = ( ,
=( O"space_n"]], O"time_n"]], ((m)) ),
=( space=O"space_name"]], =O"time_name"]], stat=(m) )
)
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
# matchfrom already created higher up
matchto = ( space=O"space_id"]], =O"time_id"]] )
for (k 1:((m))) {
W,,k] = reformat_to_array( input=(m,k]), matchfrom=O"matchfrom"]], matchto=matchto)
}
Opredictions"predictions"]] = W,, tokeep, =]
m = W =
}
if ( O"dimensionality"]] == "space-time-cyclic" ) {
W = ( ,
=( O"space_n"]], O"time_n"]], O"cyclic_n"]], ((m)) ),
=( space=O"space_name"]], =O"time_name"]], cyclic=O"cyclic_name"]], stat=(m) ) )
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
((W))[3] = vU # need to do this in a separate step ..
matchto = ( space=O"space_id"]], =O"time_id"]], cyclic=O"cyclic_id"]] )
for (k 1:((m))) {
W,,,k] = reformat_to_array( input=(m,k]), matchfrom=O"matchfrom"]], matchto=matchto )
}
Opredictions"predictions"]] = W,,, tokeep, =]
m = W =
}
}
}
# save summary
fn_preds = ( "_fit~", "_predictions~", fn_fit, fixed= )
read_write_fast( =Opredictions, =fn_preds, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
Opredictions =
()
}
O"ipred"]] =
# --------------------------------------------
# --------------------------------------------
# --------------------------------------------
run_post_samples = ()
if (redo_posterior_simulations) {
if (("debug")) if (()) if ( =="posterior_samples") ()
if ((nposteriors)) {
if (("nposteriors", O)) {
nposteriors = O$nposteriors
} {
("nposteriors not found, defaulting to 0")
nposteriors = 0
}
}
( "Sampling from joint posteriors: n = ", nposteriors )
if (nposteriors > 0) {
S = ( inla.posterior.sample( nposteriors, fit, add.names=, num.threads=mc.cores ), silent=)
if ( "try-error" %in% (S) ) {
S = ( inla.posterior.sample( nposteriors, fit, add.names=, use.improved.mean= ), silent=)
}
if ( "try-error" %in% (S) ) {
S = ( inla.posterior.sample( nposteriors, fit, add.names=, skew.corr = ), silent=)
}
if ( "try-error" %in% (S) ) {
S = ( inla.posterior.sample( nposteriors, fit, add.names=, use.improved.mean=, skew.corr = ), silent=)
}
if ( "try-error" %in% (S) ) ("posterior sampling error")
if (0) {
# not used at the moment ..
nlogdens = (S[1]]$logdens)
logdens = (, =( nlogdens, nposteriors ) )
for (i 1:nposteriors) {
logdens,i] = (S[i]]$logdens)
}
logdens_names = (S[1]]$logdens)
logdens = format_results( logdens, =logdens_names ) # same seq as space_id ( == attributes(space)$row.names )
}
if (be_verbose) ( "\nPosterior simulations complete. Extracting required components ... " )
fit = ; () # no longer needed .. remove from memory
Osamples = ()
# copy a few items to keep each output autonomous
Osamples$sppoly = sppoly
Osamples$fm = O$fm
if (("time_name", O)) Osamples$time_name = O$time_name
if (("space_name", O)) Osamples$space_name = O$space_name
if (("cyclic_name", O)) Osamples$cyclic_name = O$cyclic_name
if (summarize_simulations) Osamples"summary"]] = ()
# extraction here as no need for fit object anymore (reduces RAM requirements)
if ((exceedance_threshold)) if (("exceedance_threshold", O)) exceedance_threshold = O"exceedance_threshold"]]
if ((deceedance_threshold)) if (("deceedance_threshold", O)) deceedance_threshold = O"deceedance_threshold"]]
if ((deceedance_threshold_predictions)) if (("deceedance_threshold_predictions", O)) deceedance_threshold_predictions = O"deceedance_threshold_predictions"]]
if ((exceedance_threshold_predictions)) if (("exceedance_threshold_predictions", O)) exceedance_threshold_predictions = O"exceedance_threshold_predictions"]]
for (z ("tag", "start", "length") ) (z, (S)".contents"]][z]] ) # index info
if ( "summary" %in% posterior_simulations_to_retain ) {
# useful to check representativeness of samples
if (("debug")) if (()) if ( =="posterior_samples_summary") ()
# -- check variable names here
# posteriors
flabels= O"names.fixed"]]
nfixed = (O"names.fixed"]])
Osamples"fixed_effects"]] = (, =( nfixed, O"nposteriors"]] ) )
fkk = inla_get_indices(O"names.fixed"]], tag=tag, =, len=, model="direct_match")
fkk = (fkk)
for (i 1:O"nposteriors"]]) Osamples"fixed_effects"]],i] = S[i]]$latent[fkk,]
Osamples"fixed_effects"]] = invlink_betas(Osamples"fixed_effects"]])
(Osamples"fixed_effects"]]) = flabels
if (summarize_simulations) {
Osamples"summary"]]"fixed_effects"]] = posterior_summary(format_results( fixed, =flabels))
}
fkk = fixed =
other_random = ( (O$random), (vS, vS2 ) )
if ((other_random) > 0 ) {
nrandom = ((O$random[other_random], ))
Osamples"random_effects"]] = (, =(nrandom, O"nposteriors"]] ) )
rkk = inla_get_indices(other_random, tag=tag, =, len=, model="direct_match" ) # if bym2, must be decomposed
rkk = ( rkk )
for (i 1:O"nposteriors"]]) random,i] = S[i]]$latent[rkk,]
rlabels = (S[1]]$latent[rkk,])
(Osamples"random_effects"]]) = rlabels
if (summarize_simulations) {
Osamples"summary"]]"random_effects"]] = posterior_summary( format_results( Osamples"random_effects"]], =rlabels ) )
}
rkk = nrandom= rlabels=
}
hyper_names = (S[1]]$hyperpar)
nhyperpar = (hyper_names)
if (nhyperpar == 0) ("No hyperparameters? ... control.mode(restart=TRUE) might help" )
Osamples"hyperpars"]] = ( ,=(nhyperpar, O"nposteriors"]] ) )
for (i 1:O"nposteriors"]]) Osamples"hyperpars"]],i] = S[i]]$hyperpar
k = ("Precision", hyper_names)
if ((k) > 0) {
Osamples"hyperpars"]][k,] = 1 / (Osamples"hyperpars"]][k,])
hyper_names[k] = ("Precision for", "SD", hyper_names[k] )
}
hyper_names = ("for ", "", hyper_names )
hyper_names = ("the ", "", hyper_names )
( Osamples"hyperpars"]] ) = hyper_names
if (summarize_simulations) {
Osamples"summary"]]"hyperpars"]] = posterior_summary( format_results( Osamples"hyperpars"]], =hyper_names) )
}
hyper_names = nhyperpar = k =
} # end summary
k = known = unknown = m =
()
if ( "random_spatial" %in% posterior_simulations_to_retain ) {
if (("debug")) if (()) if ( =="posterior_samples_spatial") ()
# POSTERIOR samples (on link scale)
space = (, =( O$space_n, O"nposteriors"]] ) )
(space) = O"space_name"]]
space1 = space2 =
fmre = O"fm"]]$random_effects
iSP = ( fmre$dimensionality=="s" & fmre$level=="main")
matchto = ( space=O"space_id"]] )
if ((iSP) == 1) {
stx1 = ("^", fmre$vn[iSP], "$", sep="")
if (fmre$model[iSP] %in% ("bym2", "bym") ) {
# special case bym2 has two factors rolled together
space1 = (, =( O$space_n, O"nposteriors"]] ) )
space2 = (, =( O$space_n, O"nposteriors"]] ) )
(space1) = O"space_name"]]
(space2) = O"space_name"]]
skk1 = inla_get_indices(stx1, tag=tag, =, len=, model="bym2" ) # if bym2, must be decomposed
for (i 1:O"nposteriors"]]) {
space,i] = S[i]]$latent[skk1"re_total"]],] # ICAR + IID
space2,i] = S[i]]$latent[skk1"re_neighbourhood"]],] # ICAR
}
space1 = space - space2 # re_unstructured (IID)
} {
# single spatial effect of some kind
skk1 = inla_get_indices(stx1, tag=tag, =, len= )
skk1 = (skk1)
for (i 1:O"nposteriors"]]) {
space,i] = S[i]]$latent[skk1,] # ICAR or IID
}
}
}
if ((iSP) == 2) {
# Assume additive
space1 = (, =( O$space_n, O"nposteriors"]] ) )
space2 = (, =( O$space_n, O"nposteriors"]] ) )
(space1) = O"space_name"]]
(space2) = O"space_name"]]
model_name1 = fmre$model iSP[1] ]
model_name2 = fmre$model iSP[2] ]
stx1 = ("^", fmre$vn[iSP[1]], "$", sep="")
stx2 = ("^", fmre$vn[iSP[i]], "$", sep="")
skk1 = inla_get_indices(stx1, tag=tag, =, len= ) # if bym2, must be decomposed
skk2 = inla_get_indices(stx2, tag=tag, =, len= ) # if bym2, must be decomposed
for (i 1:O"nposteriors"]]) {
space1,i] = S[i]]$latent[skk1[1]],]
space2,i] = S[i]]$latent[skk2[1]],]
}
space = space1 + space2
("Two spatial effects .. assuming they are are additive")
}
space = invlink_betas(space) # the overall spatial random effect
(space) = O"space_name"]]
matchfrom = ( space=O"space_id"]] )
Osamples[vS]] = ()
if (!(exceedance_threshold)) {
if (be_verbose) ("Extracting random spatial effects exceedence" )
for ( b 1:(exceedance_threshold)) {
m = apply ( space, 1, FUN=(x) ( (x > exceedance_threshold[b]) ) ) / O"nposteriors"]]
m = reformat_to_array( input = m, matchfrom=matchfrom, matchto = matchto )
((m))[1] = vS
( m )[vS]] = O"space_name"]]
Osamples[vS]] "exceedance"]] [as.character(exceedance_threshold[b])]] = ( m , tokeep, =], ID=(m) )
m =
}
}
if (!(deceedance_threshold)) {
if (be_verbose) ("Extracting random spatial effects deceedance" )
# redundant but generalizable to higher dims
for ( b 1:(deceedance_threshold)) {
m = apply ( space, 1, FUN=(x) ( (x < deceedance_threshold[b]) ) ) / O"nposteriors"]]
m = reformat_to_array( input = m, matchfrom=matchfrom, matchto=matchto )
((m))[1] = vS
( m )[vS]] = O"space_name"]]
Osamples[vS]] "deceedance"]] [as.character(deceedance_threshold[b])]] = ( m , tokeep, =], ID=(m) )
m =
}
}
if (summarize_simulations) {
Osamples"summary"]] [vS]] = ()
m = posterior_summary( format_results( space, =O"space_name"]] ) )
W = ( , =( O"space_n"]], (tokeep) ), =( space=O"space_name"]], stat=tokeep ) )
((W))[1] = vS # need to do this in a separate step ..
for (k 1:(tokeep)) {
W,k] = reformat_to_array( input = m, tokeep[k]], matchfrom=matchfrom, matchto=matchto )
}
Osamples"summary"]] [vS]] "re_total"]] = W, tokeep, =]
}
Osamples[vS]]"re_total"]] = space # already inverse link scale (line 1517)
if ( "random_spatial12" %in% posterior_simulations_to_retain ) {
if (!(space1)) Osamples [vS]] [model_name1]] = invlink_betas(space1)
if (!(space2)) Osamples [vS]] [model_name2]] = invlink_betas(space2)
}
} # END random spatial post samples
matchfrom = i1 = i2=
Z = W = m = space = space1 = space2 = skk1 = skk2 = iSP =
()
if ( "random_spatiotemporal" %in% posterior_simulations_to_retain ) {
# posterior simulations (on link scale)
if (("debug")) if (()) if ( =="posterior_samples_spatiotemporal") ()
fmre = O"fm"]]$random_effects
iST = ( fmre$dimensionality=="st" & fmre$level=="main")
matchto = ( space=O"space_id"]], =O"time_id"]] )
space_time1 = space_time2 = space_time = ( ,
=( O"space_n"]] * O"time_n"]] , O"nposteriors"]] ) )
L = CJ( =O"time_id"]], space=O"space_id"]] ) # note:: CJ has reverse order vs expand.grid
stlabels = (L"space"]], L"time"]], sep="_")
if ((iST) == 1) {
stx1 = ("^", fmre$vn[iST], "$", sep="")
if (fmre$model[iST] %in% ("bym", "bym2") ) {
# special case bym2 has two factors rolled together
skk1 = inla_get_indices(stx1, tag=tag, =, len=, model="bym2" ) # if bym2, must be decomposed
for (i 1:O"nposteriors"]]) {
space_time,i] = S[i]]$latent[skk1"re_total"]],]
space_time2,i] = S[i]]$latent[skk1"re_neighbourhood"]],]
}
space_time1 = space_time - space_time2 # re_unstructured
(space_time1) = stlabels
(space_time2) = stlabels
(space_time) = stlabels
} {
# single spatial effect of some kind
skk1 = inla_get_indices(stx1, tag=tag, =, len= ) # if bym2, must be decomposed
skk1 = (skk1)
for (i 1:O"nposteriors"]]) {
space_time,i] = S[i]]$latent[skk1,]
}
}
}
if ((iST) == 2) {
model_name1 = fmre$model iST[1] ]
model_name2 = fmre$model iST[2] ]
stx1 = ("^", fmre$vn[iST[1]], "$", sep="")
stx2 = ("^", fmre$vn[iST[2]], "$", sep="")
skk1 = inla_get_indices(stx1, tag=tag, =, len= ) # if bym2, must be decomposed
skk2 = inla_get_indices(stx2, tag=tag, =, len= ) # if bym2, must be decomposed
for (i 1:O"nposteriors"]]) {
space_time1,i] = S[i]]$latent[skk1[1]],]
space_time2,i] = S[i]]$latent[skk2[1]],]
}
space_time = space_time1 + space_time2 # assume additive
(space_time1) = stlabels
(space_time2) = stlabels
}
space_time = invlink_betas(space_time)
(space_time) = stlabels
Z = expand.grid( space=O"space_id"]], ="re_total", =O"time_id"]], stringsAsFactors = )
jst = (Z$=="re_total")
matchfrom = ( space=Z"space"]][jst], =Z"time"]][jst] )
Osamples[vnST]] = ()
if (!(exceedance_threshold)) {
if (be_verbose) ("Extracting random spatiotemporal errors exceedence" )
Osamples[vnST]] "exceedance"]] = ()
for ( b 1:(exceedance_threshold)) {
m = apply ( space_time, 1, FUN=(x) ( (x > exceedance_threshold[b] ) ) ) / O"nposteriors"]]
m = reformat_to_array( input=m, matchfrom=matchfrom, matchto=matchto )
((m))[1] = vS
((m))[2] = vT
( m )[vS]] = O"space_name"]]
( m )[vT]] = O"time_name"]]
Osamples[vnST]] "exceedance"]] [as.character(exceedance_threshold[b])]] = m
}
m =
}
if (!(deceedance_threshold)) {
if (be_verbose) ("Extracting random spatiotemporal errors deceedance" )
Osamples[vnST]] "deceedance"]] = ()
for ( b 1:(deceedance_threshold)) {
m = apply ( space_time, 1, FUN=(x) ( (x < deceedance_threshold) ) ) / O"nposteriors"]]
m = reformat_to_array( input = m, matchfrom=matchfrom, matchto=matchto )
((m))[1] = vS
((m))[2] = O"time_name"]]
( m )[vS]] = O"space_name"]]
( m )[vT]] = vT
Osamples[vnST]] "deceedance"]] [as.character(deceedance_threshold[b])]] = m
}
m =
}
if (summarize_simulations) {
W = ( ,
=( O"space_n"]], O"time_n"]], (tokeep) ),
=( space=O"space_name"]], =O"time_name"]], stat=tokeep ) )
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
m = posterior_summary(format_results( space_time, =stlabels ))
for (k 1:(tokeep)) {
W,,k] = reformat_to_array( input = m, tokeep[k]], matchfrom=matchfrom, matchto=matchto )
}
Osamples"summary"]][vnST]] = ()
Osamples"summary"]][vnST]] "re_total"]] = W,, tokeep, =]
}
if ( "random_spatiotemporal" %in% posterior_simulations_to_retain ) {
Osamples [vnST]] "re_total"]] = space_time # already invlinked .. line 1626
}
if ( "random_spatiotemporal12" %in% posterior_simulations_to_retain ) {
if (!(space_time1)) Osamples [vnST]] [model_name1]] = invlink_betas(space_time1)
if (!(space_time2)) Osamples [vnST]] [model_name2]] = invlink_betas(space_time2)
}
} # END sp-temp post samples
Z = W = m = space_time = space_time1 = space_time2 = skk1 = skk2 =
()
if ( "predictions" %in% posterior_simulations_to_retain ) {
# marginals.fitted.values are on response scale (for both experimental and classical) and already incorporates offsets
# summary.fitted.values on response scale ((for both experimental and classical))
# posteriors are on link scale and already incorporates offsets
if (("debug")) if (()) if ( =="posterior_samples_predictions") ()
if (be_verbose) ("Extracting posterior predictions" )
# prepare prediction simulations (from joint posteriors)
ptx = "^Predictor"
preds = O"predictions_range"]][1]:O"predictions_range"]][2]
npredictions = (preds)
pkk = inla_get_indices(ptx, tag=tag, =, len=)
pkk = (pkk)[preds]
predictions = (, =( npredictions, O"nposteriors"]] ) )
for (i 1:O"nposteriors"]]) predictions,i] = S[i]]$latent[pkk, ]
predictions = invlink(predictions) # was on link scale .. now on response scale
pkk = preds = S =
()
if ( ("data_transformation", O)) predictions = O$data_transformation$backward( predictions )
if ( O"dimensionality"]] == "space" ) {
matchto = ( space=O"space_id"]] )
W = ( ,
=( O"space_n"]], O"nposteriors"]] ),
=( space=O"space_name"]], sim=1:O"nposteriors"]] ) )
((W))[1] = vS # need to do this in a separate step ..
for (k 1:O"nposteriors"]] ) {
W,k] = reformat_to_array( input=(predictions,k]), matchfrom=O"matchfrom"]], matchto=matchto )
}
Osamples"predictions"]] = W, =] # on response scale
W =
}
if (O"dimensionality"]] == "space-time" ) {
matchto = ( space=O"space_id"]], =O"time_id"]] )
W = ( ,
=( O"space_n"]], O"time_n"]], O"nposteriors"]] ),
=( space=O"space_name"]], =O"time_name"]], sim=1:O"nposteriors"]] ) )
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
for (k 1:O"nposteriors"]] ) {
W,,k] = reformat_to_array( input=(predictions,k]), matchfrom=O"matchfrom"]], matchto=matchto )
}
Osamples"predictions"]] = W,,, =] # on response scale
W =
}
if ( O"dimensionality"]] == "space-time-cyclic" ) {
matchto = ( space=O"space_id"]], =O"time_id"]], cyclic=O"cyclic_id"]] )
W = ( ,
=( O"space_n"]], O"time_n"]], O"cyclic_n"]], O"nposteriors"]] ),
=( space=O"space_name"]], =O"time_name"]], cyclic=O"cyclic_name"]], sim=1:O"nposteriors"]] ) )
((W))[1] = vS # need to do this in a separate step ..
((W))[2] = vT # need to do this in a separate step ..
((W))[3] = vU # need to do this in a separate step ..
for (k 1:O"nposteriors"]] ) {
W,,,k] = reformat_to_array( input=(predictions,k]), matchfrom=O"matchfrom"]], matchto=matchto )
}
Osamples"predictions"]] = W,,, , =] # on response scale
W =
}
if (!(exceedance_threshold_predictions)) {
if (be_verbose) ("Extracting de/exceedance of predictions" )
Osamples"predictions_exceedance"]] = ()
Osamples"predictions_exceedance"]] "exceedance"]] = ()
for ( b 1:(exceedance_threshold_predictions)) {
m = apply ( predictions, 1, FUN=(x) ( (x > exceedance_threshold_predictions[b] ) ) ) / O"nposteriors"]]
W = reformat_to_array( input=m, matchfrom=O"matchfrom"]], matchto=matchto )
((W))[1] = vS
( W )[vS]] = O"space_name"]]
m =
if ( (O"matchfrom"]]) > 1 ) {
((W))[2] = vT
( W )[vT]] = O"time_name"]]
}
if ((O"matchfrom"]]) > 2 ) {
((W))[3] = vU
( W )[vU]] = O"cyclic_name"]]
}
Osamples"predictions_exceedance"]] "exceedance"]] (exceedance_threshold_predictions[b]) ]]= W
W =
}
}
if (!(deceedance_threshold_predictions)) {
Osamples"predictions_deceedance"]] = ()
Osamples"predictions_deceedance"]] "deceedance"]] = ()
for ( b 1:(deceedance_threshold_predictions)) {
m = apply ( predictions, 1, FUN=(x) ( (x < deceedance_threshold_predictions[b] ) ) ) / O"nposteriors"]]
W = reformat_to_array( input=m, matchfrom=O"matchfrom"]], matchto=matchto )
((W))[1] = vS
( W )[vS]] = O"space_name"]]
m =
if ( (O"matchfrom"]]) > 1 ) {
((W))[2] = vT
( W )[vT]] = O"time_name"]]
}
if ( (O"matchfrom"]]) > 2 ) {
((W))[3] = vU
( W )[vU]] = O"cyclic_name"]]
}
Osamples"predictions_deceedance"]] "deceedance"]] (deceedance_threshold_predictions[b]) ]]= W
W =
}
}
predictions =
# save summary
fn_samples = ( "_fit~", "_samples~", fn_fit, fixed= )
read_write_fast( =Osamples, =fn_samples, compress=compress, compression_level=compression_level, qs_preset=qs_preset )
Osamples =
()
} # END predictions posteriors loop
} # nposteriors > 0
} # END redo posterior samples
end_post_samples = ()
# --------------------------------------------
# --------------------------------------------
# --------------------------------------------
run_end = ()
# --------------------------------------------
# --------------------------------------------
# --------------------------------------------
if (be_verbose) {
("---------------------------------------")
("Run times:")
("Fit: ", ((run_fit_end, run_start)))
("Summarize results from marginals: ", ((run_post_samples, run_fit_end)))
("Posterior simulations: ", ((end_post_samples, run_post_samples)))
("Total: ", ((run_end, run_start)))
("---------------------------------------")
}
( fn_fit )
}
