1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(HDInterval) # for HPD intervals
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(posterior) # for posterior draws
library(ggeffects) # for partial effects plots
library(patchwork) # for multi-panel figures
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(geoR) # framework for stats, modelling and visualisation
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")

2 Scenario

Freitas et al. (2016) investigated the effect of sea temperature on habitat selection and habitat use of acoustically tagged Atlantic cod (Gadus morhua) a the Norwegian Skagerrak coast. More specifically, they explored the relationship between sea temperature and how deep individually tagged Atlantic cod dove during the day and night. These data were collected between late June and early December 2013, although the duration of the data varies across the different individual cod.

The data are in the file freitas.csv in the data folder.

Table 1: Format of the freitas.csv data file

FISH	DATE	DEPTH_MEAN_DAY	DEPTH_MEAN_NIGHT	TEMPERATURE
Cod_6755	2013-06-28	19.64416667	19.59565217	18.00541667
Cod_6755	2013-06-29	19.98275862	19.72678571	17.75595833
Cod_6755	2013-06-30	19.94709302	20.040625	17.181125
Cod_6755	2013-07-01	19.96827957	19.94153846	17.49952174
Cod_6755	2013-07-02	17.74377907	19.93	17.601375
Cod_6755	2013-07-03	14.47237037	19.12307692	17.6055
...	...

Table 2: Description of the variables in the freitas data file

FISH:	Individual Atlantic code fish ID - Random variable
DATE:	Date of observation - Covariate variable
DEPTH_MEAN_DAY:	Mean dive depth that day - Response variable
DEPTH_MEAN_NIGHT:	Mean dive depth that night - Response variable
TEMPERATUREd:	Sea surface temperature - Predictor variable

3 Read in the data

freitas <- read_csv("../data/freitas.csv", trim_ws = TRUE)

Rows: 5083 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): FISH
dbl  (3): DEPTH_MEAN_DAY, DEPTH_MEAN_NIGHT, TEMPERATURE
date (1): DATE

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(freitas)

Rows: 5,083
Columns: 5
$ FISH             <chr> "Cod_6755", "Cod_6755", "Cod_6755", "Cod_6755", "Cod_…
$ DATE             <date> 2013-06-28, 2013-06-29, 2013-06-30, 2013-07-01, 2013…
$ DEPTH_MEAN_DAY   <dbl> 19.64417, 19.98276, 19.94709, 19.96828, 17.74378, 14.…
$ DEPTH_MEAN_NIGHT <dbl> 19.59565, 19.72679, 20.04062, 19.94154, 19.93000, 19.…
$ TEMPERATURE      <dbl> 18.00542, 17.75596, 17.18113, 17.49952, 17.60138, 17.…

## Explore the first 6 rows of the data
head(freitas)

# A tibble: 6 × 5
  FISH     DATE       DEPTH_MEAN_DAY DEPTH_MEAN_NIGHT TEMPERATURE
  <chr>    <date>              <dbl>            <dbl>       <dbl>
1 Cod_6755 2013-06-28           19.6             19.6        18.0
2 Cod_6755 2013-06-29           20.0             19.7        17.8
3 Cod_6755 2013-06-30           19.9             20.0        17.2
4 Cod_6755 2013-07-01           20.0             19.9        17.5
5 Cod_6755 2013-07-02           17.7             19.9        17.6
6 Cod_6755 2013-07-03           14.5             19.1        17.6

str(freitas)

spc_tbl_ [5,083 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ FISH            : chr [1:5083] "Cod_6755" "Cod_6755" "Cod_6755" "Cod_6755" ...
 $ DATE            : Date[1:5083], format: "2013-06-28" "2013-06-29" ...
 $ DEPTH_MEAN_DAY  : num [1:5083] 19.6 20 19.9 20 17.7 ...
 $ DEPTH_MEAN_NIGHT: num [1:5083] 19.6 19.7 20 19.9 19.9 ...
 $ TEMPERATURE     : num [1:5083] 18 17.8 17.2 17.5 17.6 ...
 - attr(*, "spec")=
  .. cols(
  ..   FISH = col_character(),
  ..   DATE = col_date(format = ""),
  ..   DEPTH_MEAN_DAY = col_double(),
  ..   DEPTH_MEAN_NIGHT = col_double(),
  ..   TEMPERATURE = col_double()
  .. )
 - attr(*, "problems")=<externalptr>

freitas |> datawizard::data_codebook()

freitas (5083 rows and 5 variables, 5 shown)

ID | Name             | Type      | Missings |        Values |          N
---+------------------+-----------+----------+---------------+-----------
1  | FISH             | character | 0 (0.0%) |      Cod_6755 | 101 (2.0%)
   |                  |           |          |      Cod_6761 | 104 (2.0%)
   |                  |           |          |      Cod_6765 | 105 (2.1%)
   |                  |           |          |      Cod_6766 |  93 (1.8%)
   |                  |           |          |      Cod_6773 | 107 (2.1%)
   |                  |           |          |      Cod_6776 |   6 (0.1%)
   |                  |           |          |      Cod_6778 |  53 (1.0%)
   |                  |           |          |      Cod_6779 | 156 (3.1%)
   |                  |           |          |      Cod_6780 | 107 (2.1%)
   |                  |           |          |      Cod_6783 | 107 (2.1%)
   |                  |           |          |         (...) |           
---+------------------+-----------+----------+---------------+-----------
2  | DATE             | date      | 0 (0.0%) |    2013-06-28 |  44 (0.9%)
   |                  |           |          |    2013-06-29 |  45 (0.9%)
   |                  |           |          |    2013-06-30 |  44 (0.9%)
   |                  |           |          |    2013-07-01 |  44 (0.9%)
   |                  |           |          |    2013-07-02 |  45 (0.9%)
   |                  |           |          |    2013-07-03 |  45 (0.9%)
   |                  |           |          |    2013-07-04 |  42 (0.8%)
   |                  |           |          |    2013-07-05 |  42 (0.8%)
   |                  |           |          |    2013-07-06 |  42 (0.8%)
   |                  |           |          |    2013-07-07 |  42 (0.8%)
   |                  |           |          |         (...) |           
---+------------------+-----------+----------+---------------+-----------
3  | DEPTH_MEAN_DAY   | numeric   | 0 (0.0%) | [1.68, 41.41] |       5083
---+------------------+-----------+----------+---------------+-----------
4  | DEPTH_MEAN_NIGHT | numeric   | 0 (0.0%) | [0.05, 37.85] |       5083
---+------------------+-----------+----------+---------------+-----------
5  | TEMPERATURE      | numeric   | 0 (0.0%) |  [6.72, 21.1] |       5083
-------------------------------------------------------------------------

freitas |> modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
DEPTH_MEAN_DAY	5053	0	17.7	6.4	1.7	19.5	41.4
DEPTH_MEAN_NIGHT	5053	0	15.3	7.3	0.0	16.7	37.9
TEMPERATURE	159	0	16.5	3.5	6.7	17.8	21.1
FISH	N	%
Cod_6755	101	2.0
Cod_6761	104	2.0
Cod_6765	105	2.1
Cod_6766	93	1.8
Cod_6773	107	2.1
Cod_6776	6	0.1
Cod_6778	53	1.0
Cod_6779	156	3.1
Cod_6780	107	2.1
Cod_6783	107	2.1
Cod_6784	62	1.2
Cod_6785	111	2.2
Cod_6787	92	1.8
Cod_6791	110	2.2
Cod_6793	111	2.2
Cod_6795	111	2.2
Cod_7266	109	2.1
Cod_7267	112	2.2
Cod_7270	114	2.2
Cod_7272	56	1.1
Cod_7274	26	0.5
Cod_7279	102	2.0
Cod_7280	159	3.1
Cod_7282	57	1.1
Cod_7283	112	2.2
Cod_7285	59	1.2
Cod_7286	46	0.9
Cod_8981	107	2.1
Cod_8982	54	1.1
Cod_8984	158	3.1
Cod_8985	156	3.1
Cod_8986	159	3.1
Cod_8987	48	0.9
Cod_8988	159	3.1
Cod_8999	27	0.5
Cod_9000	159	3.1
Cod_9002	65	1.3
Cod_9003	150	3.0
Cod_9004	159	3.1
Cod_9005	159	3.1
Cod_9031	159	3.1
Cod_9032	151	3.0
Cod_9034	159	3.1
Cod_9035	83	1.6
Cod_9043	159	3.1
Cod_9044	159	3.1
Cod_9045	102	2.0
Cod_9046	63	1.2

4 Data preparation

For the current worksheet, we are going to restrict ourselfs to a single individual fish (Cod_9044) so that we can concentrate on dependency issues without the extra distraction of dependency issues introduced by random effects.

We will also restrict ourselves to exploring the mean diving depth during daylight hours as the response (i.e. we will ignore the MEAN_DEPTH_NIGHT variable).

filter to Cod_9044
declare categorical variables as factors
create a numerical version of DATE

freitas <- freitas |>
  filter(FISH == "Cod_9044") |>
  droplevels() |>
  mutate(
    DAY = as.numeric(as.factor(DATE)),
    DEPTH = DEPTH_MEAN_DAY
  )

5 Exploratory data analysis

As the response represents a strictly positive real number (mean depth dived during the day), it would make sense to start by considering a Gamma model.

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{Gamma}(\mu_i, \phi)\\ ln(\mu_i) &=\boldsymbol{\beta} \bf{X_i} \end{align} \]

where $\boldsymbol{\beta}$ is a vector of the fixed parameters and $\bf{X}$ is the model matrix representing the overall intercept and effects of temperature.

freitas |> ggplot(aes(y = DEPTH)) +
  geom_boxplot()

## freitas |> ggplot(aes(y = Dist)) +
##   geom_boxplot()
freitas |> ggplot(aes(y = DEPTH)) +
  geom_boxplot() +
  scale_y_log10()

freitas |> ggplot(aes(y = TEMPERATURE)) +
  geom_boxplot() +
  scale_y_log10()

freitas |> ggplot(aes(y = TEMPERATURE)) +
  geom_boxplot()

## freitas |> ggplot(aes(y = Dist, x = Day)) +
##   geom_point() +
##   geom_line(aes(group = ID))
freitas |> ggplot(aes(y = DEPTH, x = TEMPERATURE)) +
  geom_point() +
  geom_line(aes(colour = DATE))

6 Fit the model

brms

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

freitas_form <- bf(I(DEPTH) ~ scale(TEMPERATURE),
                family=gaussian(link = 'identity'))
options(width=150)
freitas_form |> get_prior(data = freitas)

                   prior     class             coef group resp dpar nlpar lb ub       source
                  (flat)         b                                                   default
                  (flat)         b scaleTEMPERATURE                             (vectorized)
 student_t(3, 14.1, 4.3) Intercept                                                   default
    student_t(3, 0, 4.3)     sigma                                         0         default

options(width=80)

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

$\beta_0$: normal centred at 134 with a standard deviation of 65
- mean of 134: since median(log(freitas$Dist_moved/1000))
- sd pf 0.7: since mad(log(freitas$Dist_moved/1000))
$\beta_1$: normal centred at 0 with a standard deviation of 2
- sd of 1: since mad(log(Dist_moved/1000))/mad(scale(Day))
$\sigma$: half-t centred at 0 with a standard deviation of 65 OR
- sd pf 65: since mad(log(freitas$Dist_moved/1000))
$\sigma$: gamma with shape parameters of 2 and 1
$\beta_0$: normal centred at 1.5 with a standard deviation of 1.5
- mean of 1.5: since mean(log(freitas$Dist_moved/1000)) or mean(asinh(freitas$Dist_moved/(2*1))/log(exp(1)))
- sd of 1.5: since sd(log(freitas$Dist_moved/1000)) or sd(asinh(freitas$Dist_moved/(2*1))/log(exp(1)))

freitas |>
  summarise(
    Int_mu = median(DEPTH),
    Int_sd = mad(DEPTH),
    b_sd = mad(DEPTH)
  )

# A tibble: 1 × 3
  Int_mu Int_sd  b_sd
   <dbl>  <dbl> <dbl>
1   14.1   4.32  4.32

freitas |>
  mutate(DEPTH = log(DEPTH)) |>
  summarise(
    Int_mu = median(DEPTH),
    Int_sd = mad(DEPTH),
    b_sd = mad(DEPTH)
  )

# A tibble: 1 × 3
  Int_mu Int_sd  b_sd
   <dbl>  <dbl> <dbl>
1   2.65  0.300 0.300

priors <- prior(normal(14, 4.5), class = "Intercept") +
  prior(normal(0, 4.5), class = "b") +
  prior(student_t(3, 0, 4.5), class = "sigma")
freitas_form <- bf(DEPTH ~ scale(TEMPERATURE),
  family = gaussian()
)
## freitas_form <- bf(Dist ~ scale(Day) + scale(TOD),
##                 family=Gamma(link = "log"))
get_prior(freitas_form, data = freitas)

                   prior     class             coef group resp dpar nlpar lb ub
                  (flat)         b                                             
                  (flat)         b scaleTEMPERATURE                            
 student_t(3, 14.1, 4.3) Intercept                                             
    student_t(3, 0, 4.3)     sigma                                         0   
       source
      default
 (vectorized)
      default
      default

## priors <- prior(normal(4.2, 0.2), class = 'Intercept') +
##     prior(normal(0, 0.2), class = 'b') +
##   prior(gamma(0.01, 0.01), class = 'shape')
freitas_brm2 <- brm(freitas_form,
  data = freitas,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 100,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

freitas_brm2 |>
  conditional_effects("TEMPERATURE") |>
  plot(points = TRUE)

The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.

freitas_brm3 <- update(freitas_brm2,
  sample_prior = "yes",
  chains = 3, cores = 3,
  refresh = 100
)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

freitas_brm3 |>
  conditional_effects("TEMPERATURE") |>
  plot(points = TRUE)

freitas_brm3 |> SUYR_prior_and_posterior()

7 MCMC sampling diagnostics

brms

stan plots

The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

freitas_brm3$fit |> stan_trace()

freitas_brm3$fit |> stan_trace(inc_warmup = TRUE)

The chains appear well mixed and very similar

stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

freitas_brm3$fit |> stan_ac()

There is no evidence of auto-correlation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

freitas_brm3$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

freitas_brm3$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

freitas_brm3$fit |> stan_dens(separate_chains = TRUE)

8 Model validation

brms

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_dots
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

freitas_brm3 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to represent the shape of the observed data reasonably well

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

# freitas_brm3 |> pp_check(type = 'error_scatter_avg')

The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments.

intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

freitas_brm3 |> pp_check(type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

## freitas_brm3 |> pp_check(group='DENSITY', type='intervals')

Whilst the modelled predictions do a reasonable job of representing the observed data, the observed data do appear to be more varied than the model is representing.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(freitas_brm3)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

freitas_resids <- make_brms_dharma_res(freitas_brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(freitas_resids)) +
  wrap_elements(~ plotResiduals(freitas_resids, form = factor(rep(1, nrow(freitas))))) +
  wrap_elements(~ plotResiduals(freitas_resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(freitas_resids))

freitas_resids |> testTemporalAutocorrelation(time = freitas$DAY)


    Durbin-Watson test

data:  simulationOutput$scaledResiduals ~ 1
DW = 0.61291, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0

## freitas_resid1 <- freitas_resids |>
##   recalculateResiduals(group = freitas$Day, aggregateBy = mean)
##   ## recalculateResiduals(group = interaction(freitas$Day,  freitas$ID),  aggregateBy = mean)
## freitas_resid1 |> testTemporalAutocorrelation(time=unique(freitas$Day))
autocor_check(freitas, freitas_brm3, variable = "DAY", n.sim = 250)

variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram

Using  as id variables

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

`geom_smooth()` using formula = 'y ~ x'

residuals(freitas_brm3)[, "Estimate"] |> acf()

Conclusions:

the simulated residuals do suggest that there might be a dispersion issue
it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.

9 Refit model

freitas_form <- bf(
  DEPTH ~ scale(TEMPERATURE) +
    ar(time = DAY, p = 1, cov = TRUE),
  family = gaussian()
)
get_prior(freitas_form, data = freitas)

                   prior     class             coef group resp dpar nlpar lb ub
                  (flat)        ar                                        -1  1
                  (flat)         b                                             
                  (flat)         b scaleTEMPERATURE                            
 student_t(3, 14.1, 4.3) Intercept                                             
    student_t(3, 0, 4.3)     sigma                                         0   
       source
      default
      default
 (vectorized)
      default
      default

priors <- prior(normal(14, 4.5), class = "Intercept") +
  prior(normal(0, 4.5), class = "b") +
  prior(student_t(3, 0, 4.5), class = "sigma") +
  prior(normal(0, 1), class = "ar", lb = -1, ub = 1)

freitas_brm4 <- brm(freitas_form,
  data = freitas,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 1000,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

freitas_brm4 |>
  conditional_effects("TEMPERATURE") |>
  plot(points = TRUE) |>
  _[[1]]

## freitas_brm4 |> SUYR_prior_and_posterior()

9.1 stan plots