GAM_SD_prior_simulation_split.Rmd

---
title: "GAM_BETA_SD_Prior_Simulation_split_slope"
author: "Adam C. Smith"
date: "30/03/2022"
output: pdf_document
bibliography: references.bib
---

```{r setup, message=FALSE,echo=FALSE}
options(scipen=99999)
library(tidyverse)
library(cmdstanr)
library(mgcv) 
library(patchwork)
library(kableExtra)

```

# Prior Simulation of the SD of GAM parameters

Running a prior simulation for alternative prior distributions on the SD of the GAM BETA parameters, which are the hyperparameters that control the shape and wiggliness of the survey-wide mean smoothed population trajectory. This simulation compares the effects of alternative priors on derived estimates of the 1,2, 5, 10 and 20 year population trend estimates that would result from a population trajectory estimated with the same spline basis function used in this paper.

This application of the thin-plate regression spline basis (pg 215, \[@wood2017a\]), treats the linear component of the basis as an unpenalized part of the model "null space". Therefore, in this simulation, the linear parameter is fixed at 0 (stable population), and only the penalized parameters that affect the non-linear component of the trajectory are simulated.

The SD parameter in this prior simulation controls the scale and variation of the parameters that link the spline basis function to non-linear component of the estimated smooth population trajectory. The posterior estimate of this SD parameter is what determines the complexity (wiggliness) and magnitude of the population change (\[@crainiceanu2005; @bürkner2017; @goodrich2020\]).

The semi-parametric nature of smooths and the variation among alternative basis functions complicates the use of default or informative priors ([@lemoine2019; @banner2020]). We've used a prior simulation to translate these alternative priors into intuitive values of population trends that are directly interpretable as biological parameters. We then compare the prior distribution of trends that would result from these alternative priors to:

1.  the collection of realised trend estimates from a different statistical model applied to the North American Breeding Bird Survey ([@link2020]).

2.  our own prior knowledge about probable rates of change in wild bird populations at continental scales.

We compared half normal and half t-distributions for the priors on the Standard Deviation on the GAM parameters.

1.  normal

2.  t-distribution with 3 degrees of freedom

And, for each of the half-normal and half-t-distributions, we compared 5 different values to set the prior-scale: (0.5, 1, 2, 3, and 4). Given the log-link in the trend model and the scaling of the low-rank thin-plate regression spline with identifiability constraints ([@wood2020]), these 5-values of prior-scales should cover the range of plausible parameter values.

# Selected prior

We suggest a half t-distribution, with a scale parameter = 2, fits the realised distributions of 1, 5, 10, and 20 year trends for most bird species surveyed by the BBS (i.e., most bird species with the best information on population trends at continental scales and for \~50 years). This half-t prior results in prior distributions of trends that encompass the realized distributions and includes long tails that cover the range of plausible trend estimates without including large amounts of prior probability mass at implausibly extreme values.

```{r load saved plots, echo=FALSE}
load("data/prior_t_sel_split.RData")

```

```{r plot spoiler,echo=FALSE, fig.show='asis', fig.dim = c(6,6), fig.cap= "Observed distributions of the absolute values of 1, 5, 10, and 20-year population trends from the BBS data using non-GAM models (in black), and the simulated prior distribution (in green) of trends from the spline smooth basis used in this paper, with a half-t (df = 3 and scale parameter = 1) prior distribution on the standard deviation of the spline parameters"}

 print(overp_t2)
```

\newpage

------------------------------------------------------------------------

# Details on the prior simulations

Here's the code that runs the simulations in Stan.

```{r run_simulations,eval=FALSE}

# this is all set-up
library(tidyverse)
library(cmdstanr)
library(mgcv) 



for(pp in c("norm","t")){
  for(prior_scale in c(0.5,1,2,3,4)){
    
    tp = paste0("GAM_split_",pp,prior_scale,"_rate")
    
    #STRATA_True <- log(2)
    output_dir <- "output/"
    out_base <- tp
    csv_files <- paste0(out_base,"-",1:3,".csv")
    
    
    if(pp == "norm"){
      pnorm <- 1
    }
    if(pp == "t"){
      pnorm <- 0
    }
    
    #if(!file.exists(paste0(output_dir,csv_files[1]))){
      
      nyears = 54 #to match the time-scales of BBS and CBC analyses
      dat = data.frame(year = 1:nyears)
      nknots = 13  
      nknots_realised = nknots-2 #removes the two non-penalized components (mean and linear) 
      lin_component = nknots-1 # the mgcv function below orders the basis such that
                # the linear component is hte final column in the prediction matrix
      M = mgcv::smoothCon(s(year,k = nknots, bs = "tp"),data = dat,
                          absorb.cons=TRUE,#this drops the constant
                          diagonal.penalty=TRUE) ## If TRUE then the smooth is reparameterized to turn the penalty into an identity matrix, with the final diagonal elements zeroed (corresponding to the penalty nullspace). This fits with the Bayesian interpretation of the complexity penalty as the inverse of the variance of the i.i.d. collection of parameters.
      
      
      year_basis = M[[1]]$X
      
      
      
      
      stan_data = list(#scalar indicators
        nyears = nyears,
        
        
        #GAM structure
        nknots_year = nknots_realised,
        year_basis = year_basis,
        
        prior_scale = prior_scale,
        pnorm = pnorm,
        lin_component = lin_component
      )
      
      
      
      
      # Fit model ---------------------------------------------------------------
      
      print(paste("beginning",tp,Sys.time()))
      
      mod.file = "models/GAM_split_prior_sim.stan"
      
      ## compile model
      model <- cmdstan_model(mod.file)
      
      
      # Initial Values ----------------------------------------------------------
      
      
      init_def <- function(){ list(sdbeta = runif(1,0.01,0.1),
                                   BETA_raw = rnorm(nknots_realised,0,0.01))}
      
      stanfit <- model$sample(
        data=stan_data,
        refresh=100,
        chains=2, iter_sampling=1000,
        iter_warmup=500,
        parallel_chains = 2,
        #pars = parms,
        adapt_delta = 0.8,
        max_treedepth = 14,
        seed = 123,
        init = init_def,
        output_dir = output_dir,
        output_basename = out_base)
      
      
      #stanfit1 <- as_cmdstan_fit(files = paste0(output_dir,csv_files))
      
      
      save(list = c("stanfit","stan_data","csv_files",
                  "out_base"),
         file = paste0(output_dir,"/",out_base,"_gamye_iCAR.RData"))
    
    
    
  
  
}#end prior_scale loop
}#end pp loop



```

## The simulation model

The Stan simulation model is very simple. We implemented it in Stan to match the implementation in the full model (although this simulation doesn't require MCMC sampling).

```{stan, output.var = "GAM_split_prior_sim.stan", eval = FALSE}
// simple GAM prior simulation

data {
  int<lower=1> nyears;
  real<lower=0>  prior_scale; //scale of the prior distribution
  int<lower=0,upper=1> pnorm; // indicator for the prior distribution 0 = t, 1 = normal
  // data for spline s(year)
  int<lower=1> nknots_year;  // number of knots in the penalized components of the basis function for year
  int<lower=1> lin_component;  // column of the basis that represents the linear component
                              // is also the final column of the basis matrix for the thin-plate regression spline used here
  matrix[nyears, lin_component] year_basis; // basis function matrix
}

parameters {
  real<lower=0> sdbeta;    // sd of spline coefficients  
  vector[nknots_year] BETA_raw;         // unscaled spline coefficients
}
 
transformed parameters { 
  vector[lin_component] BETA;         // spatial effect slopes (0-centered deviation from continental mean slope B)
  vector[nyears] smooth_pred;

    BETA[1:nknots_year] = sdbeta * BETA_raw; //scaling the spline parameters
    BETA[lin_component] = 0; //ensures that the linear component == 0
     smooth_pred = year_basis * BETA; //log-scale smooth trajectory
  }
  
model {

//Conditional statements to select the prior distribution
if(pnorm == 1){
 sdbeta ~ normal(0,prior_scale); //prior on sd of GAM parameter variation
}
if(pnorm == 0){
  sdbeta ~ student_t(3,0,prior_scale); //prior on sd of GAM parameter variation
}

   BETA_raw ~ normal(0,1); //non-centered parameterisation

}

 generated quantities {
  //estimated smooth on a count-scale
   vector[nyears] nsmooth = exp(smooth_pred);
    
  }



```

We then summarized the estimated trajectories as well as all possible 1-year, 5-, 10-, and 20-year, trends from the alternative priors.

```{r summarising,eval=FALSE}
source("Functions/posterior_summary_functions.R")

nsmooth_out <- NULL
trends_out <- NULL
summ_out <- NULL

for(pp in c("norm","t")){
  for(prior_scale in c(0.5,1,2,3,4)){
    
    tp = paste0("GAM_split_",pp,prior_scale,"_rate")
    
    #STRATA_True <- log(2)
    output_dir <- "output/"
    out_base <- tp
    csv_files <- paste0(out_base,"-",1:3,".csv")
    

load(paste0(output_dir,"/",out_base,"_gamye_iCAR.RData"))

summ = stanfit$summary()

summ <- summ %>% 
  mutate(prior_scale = prior_scale,
         distribution = pp)

nsmooth_samples <- posterior_samples(stanfit,
                                 parm = "nsmooth",
                                 dims = c("Year_Index"))



BETA_samples <- posterior_samples(stanfit,
                                  parm = "BETA",
                                  dims = c("k"))
BETA_wide <- BETA_samples %>% 
  pivot_wider(.,id_cols = .draw,
              names_from = k,
              names_prefix = "BETA",
              values_from = .value)


nsmooth_samples <- nsmooth_samples %>% 
  left_join(., BETA_wide,by = ".draw") %>% 
  mutate(prior_scale = prior_scale,
         distribution = pp)



nyears = max(nsmooth_samples$Year_Index)
# function to calculate a %/year trend from a count-scale trajectory
trs <- function(y1,y2,ny){
  tt <- (((y2/y1)^(1/ny))-1)*100
}

for(tl in c(2,6,11,21)){ #estimating all possible 1-year, 10-year, and full trends
  ny = tl-1
  yrs1 <- seq(1,(nyears-ny),by = ny)
  yrs2 <- yrs1+ny
  for(j in 1:length(yrs1)){
    y2 <- yrs2[j]
    y1 <- yrs1[j]
    
nyh2 <- paste0("Y",y2)
nyh1 <- paste0("Y",y1)
trends <- nsmooth_samples %>% 
  filter(Year_Index %in% c(y1,y2)) %>% 
  select(.draw,.value,Year_Index) %>% 
  pivot_wider(.,names_from = Year_Index,
              values_from = .value,
              names_prefix = "Y") %>%
  rename_with(.,~gsub(pattern = nyh2,replacement = "YE", .x)) %>% 
  rename_with(.,~gsub(pattern = nyh1,replacement = "YS", .x)) %>% 
  group_by(.draw) %>% 
  summarise(trend = trs(YS,YE,ny))%>% 
  mutate(prior_scale = prior_scale,
         distribution = pp,
         first_year = y1,
         last_year = y2,
         nyears = ny)
trends_out <- bind_rows(trends_out,trends)
}
}
nsmooth_out <- bind_rows(nsmooth_out,nsmooth_samples)
summ_out <- bind_rows(summ_out,summ)
print(paste(pp,prior_scale))

  }#prior_scale
}# pp

save(file = "output/GAM_split_prior_sim_summary.RData",
     list = c("nsmooth_out",
              "trends_out",
              "summ_out"))



```

# Comparing simulation priors to realised data

## Realised trend estimates

First, here is the distribution of long-term trends from a different model for the BBS data from 1966-2019, for 426 species.

```{r realised bbs trends, fig.show='asis', fig.dim = c(8, 8)}
bbs_indices_usgs <- read.csv("data/Index_best_1966-2019_core_best.csv",
                             colClasses = c("integer",
                                            "character",
                                            "integer",
                                            "numeric",
                                            "numeric",
                                            "numeric"))

bbs_continental_inds <- bbs_indices_usgs %>% 
  filter(Region == "SU1") #just hte continental estimates

# function to calculate a %/year trend from a count-scale trajectory
trs <- function(y1,y2,ny){
  tt <- (((y2/y1)^(1/ny))-1)*100
}

miny = min(bbs_continental_inds$Year)
maxy = max(bbs_continental_inds$Year)
bbs_continental_trends <- NULL

for(tl in c(2,6,11,21)){ #estimating all possible 1-year, 2-year, 5-year, 10-year, and 20-year trends, with no uncertainty, just the point estimates based on the comparison of posterior means fo annual indices
  ny = tl-1
  yrs1 <- seq(miny,(maxy-ny),by = 1)
  yrs2 <- yrs1+ny
  for(j in 1:length(yrs1)){
    y2 <- yrs2[j]
    y1 <- yrs1[j]
    
nyh2 <- paste0("Y",y2)
nyh1 <- paste0("Y",y1)

tmp <- bbs_continental_inds %>% 
  filter(Year %in% c(y1,y2)) %>% 
  select(AOU,Index,Year) %>% 
  pivot_wider(.,names_from = Year,
              values_from = Index,
              names_prefix = "Y") %>%
  rename_with(.,~gsub(pattern = nyh2,replacement = "YE", .x)) %>% 
  rename_with(.,~gsub(pattern = nyh1,replacement = "YS", .x)) %>% 
  drop_na() %>% 
  group_by(AOU) %>% 
  summarise(trend = trs(YS,YE,ny))%>% 
  mutate(first_year = y1,
         last_year = y2,
         nyears = ny,
         abs_trend = abs(trend),
         t_years = paste0(ny,"-year trends"))

bbs_continental_trends <- bind_rows(bbs_continental_trends,tmp)
}
}

t_quants <- bbs_continental_trends %>% 
  group_by(t_years) %>% 
  summarise(x99 = quantile(abs_trend,0.99),
            x995 = quantile(abs_trend,0.995))

bbs_continental_trends <- bbs_continental_trends %>% 
  mutate(t_years = factor(t_years,
                          levels = c("1-year trends",
                                     "5-year trends",
                                     "10-year trends",
                                     "20-year trends"),
                          ordered = TRUE))

realised_long_bbs_hist <- ggplot(data = bbs_continental_trends,
                            aes(abs_trend,after_stat(density),
                                colour = t_years))+
  geom_freqpoly(breaks = seq(0,20,0.5),center = 0)+
  xlab("Absolute value of long-term BBS trends USGS models (1966-2019)")+
  theme_bw()+
  scale_y_continuous(limits = c(0,1))+
scale_colour_viridis_d(begin = 0.1,end = 0.9)
print(realised_long_bbs_hist)

```

## Summarize the long-term trends from the prior simulations

```{r prior_trend_distributions}
 
realised_long_bbs_hist <- ggplot(data = bbs_continental_trends,
                            aes(abs_trend,after_stat(density)))+
  geom_freqpoly(breaks = seq(0,20,0.5),center = 0,
                alpha = 1,size = 1)+
  xlab("Absolute value of all possible trends from USGS models (1966-2019)")+
  theme_bw()+
  scale_y_continuous(limits = c(0,0.75))+
facet_wrap(vars(t_years))


load("output/GAM_split_prior_sim_summary.RData")

trends_out <- trends_out %>% 
  mutate(abs_trend = abs(trend),# absolute values of trends
         scale_factor = factor(prior_scale,ordered = TRUE))%>% 
  mutate(t_years = factor(paste0(nyears,"-year trends"),
                          levels = c("1-year trends",
                                     "5-year trends",
                                     "10-year trends",
                                     "20-year trends"),
                          ordered = TRUE)) #just for plotting


trends_normal <- trends_out %>% 
  filter(distribution == "norm")

trends_t <- trends_out %>% 
  filter(distribution == "t")

```

The distributions for the normal priors do a reasonable job of covering the range of possible long-term trend values, but the heavier-tailed t-distributions seem to better fit the shape of the distributions.

```{r overplots long-term}

  overp_normal <- realised_long_bbs_hist +
    geom_freqpoly(data = trends_normal,
                aes(abs_trend,after_stat(density),
                    colour = scale_factor),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\nSD half-normal")+
  coord_cartesian(xlim = c(0,20))
    
  #print(overp_normal)
 

  overp_t <- realised_long_bbs_hist +
    geom_freqpoly(data = trends_t,
                aes(abs_trend,after_stat(density),
                    colour = scale_factor),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\nSD half-t df=3")+
  coord_cartesian(xlim = c(0,20))
 
  #print(overp_t)
 
  
```

```{r save plots for start long, echo=FALSE}
trends_t2 <- trends_t %>% 
  filter(prior_scale == 1)

  overp_t2 <- realised_long_bbs_hist +
    geom_freqpoly(data = trends_t2,
                aes(abs_trend,after_stat(density)),
                colour = "darkgreen",
                breaks = seq(0,100,1),center = 0,
                alpha = 0.8)+
  coord_cartesian(xlim = c(0,20))
 
  save(list = "overp_t2",file = "data/prior_t_sel_split.Rdata")
```

The half-t-distributions with a scale value of 1 or 2 fit the shape of the realised trend distribution reasonably well, and the long-tail includes significant prior mass at and beyond the observed maximum absolute values of trends.

```{r plot_prior_trend_distributions,echo=FALSE, fig.show='asis', fig.dim = c(8, 10)}
 print(overp_normal / overp_t)
```

In addition, each of the prior distributions include some prior mass at relatively extreme values, so that these priors are unlikely to overwhelm data that support a steep rate of change. The broader priors (scales of 3 or 4), include prior mass at trend values that are extremely unlikely for a wild population to sustain over more than a few years. For example, the 99th percentiles extend to values that are truly extreme (40 - 50%/year).

```{r prior_table}
quant_long_tends <- trends_out %>% 
  group_by(distribution,prior_scale) %>% 
  summarise(median_abs_t = median(abs_trend),
            U80 = quantile(abs_trend,0.80),
            U90 = quantile(abs_trend,0.90),
            U99 = quantile(abs_trend,0.99),
            pGTmax = length(which(abs_trend > 30))/length(abs_trend))

kable(quant_long_tends, booktabs = TRUE,
      digits = 3,
      format.args = list(width = 7), 
      col.names = c("Distribution",
                    "Prior Scale",
                    "Median Prior Predicted Distribution",
                    "80th percentile",
                    "90th percentile",
                    "99th percentile",
                    "Proportion of prior distribution > 30%/year"),
      caption = "Prior simulated distribution quantiles for all possible, 1-, 5-, 10-, and 20-year trends. Prior simulations for two prior distributions with 5 different scales (normal and t)") %>%
  kable_styling(font_size = 8)%>%
column_spec(column = 1:7,width = "2cm")

 
```

# Priors on the linear component of the trajectory

The overall mean rate of change for a continental landbird population in North America, is something for which we have reasonably good prior information. For a 30-50 year time-period, sustained steep rates of population change are relatively rare. For example, the realised long-term population trends from a different BBS model suggest that long-term trends beyond 10%/year are very rare.

```{r realised bbs long-term trends, fig.show='asis', fig.dim = c(8, 8)}

bbs_trends_usgs <- read.csv("data/BBS_1966-2019_core_best_trend.csv")


## selecting survey-wide trend estimates
bbs_trends_usgs_long <-bbs_trends_usgs %>% 
  filter(Region == "SU1")%>%  
  select(Trend,Species.Name) %>% 
  mutate(abs_trend = abs(Trend)) %>% #calculating absolute values of the trends
  arrange(-abs_trend)

G_long_usgs <- max(bbs_trends_usgs_long$abs_trend)
  
realised_long_bbs_freq1 <- ggplot(data = bbs_trends_usgs_long,
                            aes(abs_trend,after_stat(density)))+
  geom_freqpoly(breaks = seq(0,13,0.5),center = 0)+
  xlab("Absolute value of long-term BBS trends USGS models (1966-2019)")+
  theme_bw()+
  scale_y_continuous(limits = c(0,1))
print(realised_long_bbs_freq1)

```

The maximum absolute value of an observed long-term trends are for Cave Swallow and Eurasian Collared Dove, which have increased at an annual rate of `r round(G_long_usgs,1)` %/year. This annual rate of change implies an approximate `r signif(((((G_long_usgs/100)+1)^53)-1)*100,2)` % overall increase in the populations since 1966. As such, we feel this represents example of an "extreme" long-term trend that is unlikely to be observed in most of the BBS dataset. For example, the next largest values of trend in the data is a `r round((bbs_trends_usgs_long$abs_trend[3]),1)` %/year increase in Canada Goose populations.

```{r realised trend header,echo=FALSE}

kable(bbs_trends_usgs_long[1:3,], booktabs = TRUE,
      digits = 3,
      caption = "Most extreme long-term BBS trends") %>%
  kable_styling(font_size = 8)

```

And the largest absolute value of trends for a declining species is \< 4%/year for King Rail, Bank Swallow, and Lark Bunting.

```{r realised declines header,echo=FALSE}
bbs_declines_usgs_long <-bbs_trends_usgs_long %>% 
  filter(Trend < 0)%>%  
  arrange(-abs_trend)

kable(bbs_declines_usgs_long[1:3,], booktabs = TRUE,
      digits = 3,
      caption = "Most extreme long-term BBS declines") %>%
  kable_styling(font_size = 8)

```

We could reasonably use this realised distribution of trends to set a somewhat informative prior on the hyperparameter for the linear component of the model.

For example, a simple normal or t-prior on the linear parameter, with a standard deviation between 0.025 and 0.05 would fit the realised distribution reasonably well, while excluding extremely unlikely values. Even priors with a SD value of 0.05 or 0.1, which may seem relatively informative compared to a "non-informative" or "flat" prior on the slope such as a standard normal distribution, will include some prior mass in regions that appear to be extremely unlikely (> 20%/year over 50 years).

```{r linear_priors, echo=FALSE}
sds <- c(0.025,0.04,0.05,0.1)
lin_priors_n <- NULL
lin_priors_t <- NULL
N = 10000
for(ss in sds){
  tn <- abs((exp(rnorm(N,0,1)*ss)-1)*100)
  tmp <- data.frame(prior_scale = paste0(ss,"SD of normal"),
                    distribution = "normal",
                    abs_trend = tn)
  lin_priors_n <- bind_rows(lin_priors_n,tmp)
  
    tt <- abs((exp(rt(N,3)*ss)-1)*100)
  tmpt <- data.frame(prior_scale = paste0(ss,"SD of t(df = 3)"),
                    distribution = "t",
                    abs_trend = tt)
  lin_priors_t <- bind_rows(lin_priors_t,tmpt)
  
}



overp_lin_t <- realised_long_bbs_freq1 +
    geom_freqpoly(data = lin_priors_t,
                aes(abs_trend,after_stat(density),
                    colour = prior_scale),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\nSD half-t df=3")+
  coord_cartesian(xlim = c(0,20))



overp_lin_n <- realised_long_bbs_freq1 +
    geom_freqpoly(data = lin_priors_n,
                aes(abs_trend,after_stat(density),
                    colour = prior_scale),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\n normal")+
  coord_cartesian(xlim = c(0,20))
```

```{r linear_priors_plot }
print(overp_lin_n/overp_lin_t)


```



```{r linear_priors_show, eval=FALSE}
sds <- c(0.025,0.04,0.05,0.1)
lin_priors_n <- NULL
lin_priors_t <- NULL
N = 10000
for(ss in sds){
  tn <- abs((exp(rnorm(N,0,1)*ss)-1)*100)
  tmp <- data.frame(prior_scale = paste0(ss,"SD of normal"),
                    distribution = "normal",
                    abs_trend = tn)
  lin_priors_n <- bind_rows(lin_priors_n,tmp)
  
    tt <- abs((exp(rt(N,3)*ss)-1)*100)
  tmpt <- data.frame(prior_scale = paste0(ss,"SD of t(df = 3)"),
                    distribution = "t",
                    abs_trend = tt)
  lin_priors_t <- bind_rows(lin_priors_t,tmpt)
  
}



overp_lin_t <- realised_long_bbs_freq1 +
    geom_freqpoly(data = lin_priors_t,
                aes(abs_trend,after_stat(density),
                    colour = prior_scale),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\nSD half-t df=3")+
  coord_cartesian(xlim = c(0,20))



overp_lin_n <- realised_long_bbs_freq1 +
    geom_freqpoly(data = lin_priors_n,
                aes(abs_trend,after_stat(density),
                    colour = prior_scale),
                breaks = seq(0,100,1),center = 0,
                alpha = 0.5)+
    scale_colour_viridis_d(begin = 0.5,alpha = 0.8,
                           "Prior Scale\n normal")+
  coord_cartesian(xlim = c(0,20))

print(overp_lin_n/overp_lin_t)


```



# References