Thank you so much to Cameron Patrick for reading over my answer below and making several important corrections. I will edit my answer to highlight some of those corrections - @samuelcyk please read this answer once more to make sure you don’t miss the corrections.
Ok, so your model looks at trends in Diversity over time at each of a number of sites. Not sure how many sites you have but ideally at least 5, as per a commonly used rule of thumb. You would need these 5 sites to be “broadly similar” to each other and representative of a larger number of sites about which you wish to make inferences.
The model as you have it assumes the site-specific trends have identical shapes over time for each of your observed sites (as well as for each of the unobserved sites). The only thing that changes from one site to another is the “height” of the trend curves.
It doesn’t really matter how you define your Year variable unless you want to explicitly force the smooth of Year to be 0 in the first year of study. If Year = 0 in the first year of study, Year = 1 in the second year, etc., then it’s easier to interpret the global intercept in the “classical” sense provided you force the smooth of Year to be 0 when Year = 0. The global intercept would then represent the expected response value in the first year of study for a “typical” site (i.e., a site with a random intercept of 0). How do you force the smooth of Year to be 0 when Year = 0? You can use the option pc = 0 inside the smooth of Year: s(Year, pc = 0). But this model parametrization comes with some disadvantages too.
The default model parametrization in GAM or GAMM models assumes the smooth of Year, s(Year), is subject to a zero-sum constraint. This means that the average of the smooth across all study years will be equal to 0. The interpretation of the global intercept will be quite different now: it represents the average response value (i.e., Diversity) over the years of study for a “typical” site (i.e., a site with a random intercept equal to 0). So assume from now on we are dealing with the default model parametrization.
For that same “typical” site, the smooth effect of year would represent the shape of the time trend but “stripped of its height”. That is because the smooth is “centered about 0” to make it identifiable. (For a gaussian family, brms identifies this intercept as a population-level intercept.)
So to get the “complete” time trend for the “typical” site, you need to get its proper height (via the global intercept) PLUS its shape (via the smooth of Year). Of course, if all you cared about was the shape of the trend, the smooth alone will suffice.
It is my understanding that you will get the “complete” time trend for the “typical” site using posterior_epred() with re.form set to NA (as this will set the random effect of Site to 0). Additionally, in the newdata argument, you will need to feed a grid of values for Year. You can also set Site to NA in the newdata. (Ultimately, newdata needs to be a dataframe.)
Sites that are included in your model (i.e., they are present in the data) will still get the exact shape of time trend in your model as that for the “typical” site but the “height” of that trend will be made up of two components: 1. the global intercept; 2. a Site-specific deviation from the global intercept (aka the random effect of Site).
Once the model is fitted, these two components will be estimated and added up together to give the overall height of the time trend for a specific site included in the data.
That overall height will represent the expected response value at that site in the first year of study (i.e., when Year = 0).
It is my understanding that you will get the “complete” time trend for any site present in the data by using posterior_epred() with re.form set to NULL. Additionally, in the newdata argument, you will need to feed the specific Site you have in mind, as well as a grid of values for Year. (Again, newdata needs to be a dataframe.)
I know @AWoodward referred to the smooth in a model without random effects of Site as a “marginal effect”, but I think all we can say is that it is a “partial effect” (i.e., a “piece” we are adding to the intercept to get the expected response value)? I am hoping he can clarify the meaning of that terminology? If we take the difference between two consecutive years, Year and Year + 1, the difference in expected response values at the “typical” site would be s(Year + 1) - s(Year). So it seems to me that a “marginal” effect of Year would actually be the first derivative of the smooth of Year, not the actual smooth of Year? The actual smooth of Year is just a “partial” effect. This hurts my brain!
It is also my understanding that if you wanted to predict Diversity for the year following the last year included in your study (i.e., forecast Diversity one year ahead), then you would use posterior_predict().
Even with posterior_predict(), you have (at least) two situations: a) forecast Diversity one-year ahead at a “typical” site; b) forecast Diversity one-year ahead at one of the sites included in the data. (Of course, there is also the possibility of forecasting at a new site from your target universe of sites, that was not included in the data.)
For a) and b), I think you can specify re.form and newdata exactly as suggested above but be sure to now use posterior_predict() instead of posterior_epred().
It seems to me that posterior_smooths() would only give you the shape of the time trend at the sites you specify in its formula, but not the height. Because your model assumes the same shape at all sites (which is a very strong assumption, hopefully supported by your data), you should get the same shape whether you are looking at the “typical” site or a specific site from your data. While I haven’t tried it myself, I would imagine pisterior_smooths() should allow you to specify newdata and re.form the same way I suggested above.
I hope someone smarter than me can confirm that everything I said here makes sense. 🤭