Hi to all!
I’m trying to fit a hierarchical model to (processed) EEG data from multiple patients with many timepoints. The predicted variable peak_alpha is the frequency (in Hz) of the peak of Alpha waves, predictor ce_mac is the MAC of the used volatile anaesthetic, time is the timestamp in seconds, different for every patient. Our hypothesis is, that peak_alpha decreases with increasing ce_mac. We want to account for differences between individuals with group-level effects, so we chose a mixed model approach.
Data looks like this:
head(data, 10)
#> # A tibble: 10 x 4
#> # Groups: pid [1]
#> pid time ce_mac peak_alpha
#> <fct> <int> <dbl> <dbl>
#> 1 19 22 1 8.17
#> 2 19 37 1 8.17
#> 3 19 52 1 8.01
#> 4 19 67 1 8.1
#> 5 19 82 1 8.47
#> 6 19 97 1 8.22
#> 7 19 112 1 7.74
#> 8 19 127 1 7.88
#> 9 19 142 1 7.96
#> 10 19 157 1 8.17
tail(data, 10)
#> # A tibble: 10 x 4
#> # Groups: pid [1]
#> pid time ce_mac peak_alpha
#> <fct> <int> <dbl> <dbl>
#> 1 1071 3825 0.931 8.28
#> 2 1071 3840 0.935 7.91
#> 3 1071 3855 0.935 8.34
#> 4 1071 3870 0.937 8.37
#> 5 1071 3885 0.942 8.03
#> 6 1071 3900 0.940 8.2
#> 7 1071 3915 0.941 8.26
#> 8 1071 3930 0.941 7.85
#> 9 1071 3945 0.945 8.2
#> 10 1071 3960 0.946 8.17
Created on 2020-07-27 by the reprex package (v0.3.0)
A sample of 9 patients looks like this:
ggplot(data = sample, aes(x = ce_mac, y = peak_alpha)) +
geom_point() +
facet_wrap(vars(pid))
ggplot(data = sample, aes(x = time, y = ce_mac)) +
geom_point() +
facet_wrap(vars(pid))
ggplot(data = sample, aes(x = time, y = peak_alpha)) +
geom_point() +
facet_wrap(vars(pid))
Created on 2020-07-27 by the reprex package (v0.3.0)
The outcome variable is lognormally distributed and truncated with a lower bound of 6, so we ran the following model:
priors <-
prior(normal(0, 1), class = b) +
prior(normal(0, 1), class = sd) +
prior(normal(0, 1), class = sigma) +
prior(student_t(30, 2.3, 0.2), class = Intercept)
model <- brm(
bf(peak_alpha | trunc(lb = 6) ~ ce_mac + (1 + ce_mac || pid)),
family = lognormal(),
prior = priors,
iter = 4000,
warmup = 1000,
control = list(adapt_delta = 0.999, max_treedepth = 15)
)
check_hmc_diagnostics(model$fit)
#>
#> Divergences:
#> 0 of 12000 iterations ended with a divergence.
#>
#> Tree depth:
#> 0 of 12000 iterations saturated the maximum tree depth of 15.
#>
#> Energy:
#> E-BFMI indicated no pathological behavior.
summary(model)
#> Family: lognormal
#> Links: mu = identity; sigma = identity
#> Formula: peak_alpha | trunc(lb = 6) ~ ce_mac + (1 + ce_mac || pid)
#> Data: pluck(helper.model_list, "data", helper.model_num) (Number of observations: 14326)
#> Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup samples = 12000
#>
#> Group-Level Effects:
#> ~pid (Number of levels: 74)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.23 0.02 0.19 0.27 1.00 2094 4152
#> sd(ce_mac) 0.23 0.02 0.20 0.28 1.00 2628 5195
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 2.39 0.03 2.34 2.45 1.01 1136 2071
#> ce_mac -0.30 0.03 -0.35 -0.24 1.00 1220 2471
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.08 0.00 0.08 0.08 1.00 23269 8752
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
Created on 2020-07-27 by the reprex package (v0.3.0)
mcmc_plot(model.actual, type = "trace")
#> No divergences to plot.
mcmc_plot(model.actual, type = "pairs")
Created on 2020-07-27 by the reprex package (v0.3.0)
What I’m concerned about, is the fact, that we need >10000 iterations to get a reasonable ESS and the ACF plot presenting like this:
mcmc_plot(model.actual, type = "acf")
Created on 2020-07-27 by the reprex package (v0.3.0)
My question:
- Do we have to be concerned about this? Can we trust the posterior estimates?
- Do you need more information?
- Operating System: CentOS Cluster
- brms Version: 2.13.0