[R-sig-ME] Predicting the overall population trend from a GLMM with many sites (MCMCglmm)

[Jarrod Hadfield]

Hi,

Thinking about it more, the expected proportional change in abundance per year is perhaps the most readily digestible number for an ecologist. To obtain it (with assumptions):

The slope for site_i is

slope_i = b1+b2*latitude_i+b_i

where b1 is the main year effect, b2 is the year by latitude interaction and b_i is the random slope deviation for site i.

The b_i are assumed normal and independent of latitude. If we assume latitude is also normal with mean m and variance V, then the slopes are normal with mean b1+b2*m (since E[b_i]=0) and variance Vs+V*b2^2 (since b1 and b2 are constants). Here, Vs=VAR(b_i) [i.e. the slope variance from the model].

The exponentiated slopes are log normal with mean equal to exp(b1+b2*m+(Vs+V*b2^2)/2).

You can obtain the posterior for this as b1, b2 and Vs are parameters in the model. You would have to use point estimates for m and V (you haven't modelled latitude) so any uncertainty in the expected proportional change in abundance per year does not include uncertainty about the mean and variance of site's latitudes.

m<-mean(df$latitude)
V<-var(df$latitude)


posterior_mean<-exp(model$Sol["year"]+ model$Sol["year:latitude"]*m+(model$VCV[,"year:year.site"]+V*model$Sol["year:latitude"]^2)/2)

Cheers,

Jarrod


From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of Jarrod Hadfield <j.hadfield using ed.ac.uk>
Date: Thursday, 14 September 2023 at 16:29
To: Chris Oosthuizen <w.chris.oosthuizen using gmail.com>, r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
Subject: Re: [R-sig-ME] Predicting the overall population trend from a GLMM with many sites (MCMCglmm)
HI,

Using marginal= ~us(1 + year):site in the predict function does not marginalise sites it marginalises *site effects*: the prediction will stary vary across sites because the sites will be at varying latitudes and latitude has not been marginalised.

On the data scale, the relationship between the expected count and year is not linear and so I think it is easier to work on the link (log) scale where the relationship is linear. On this scale the mean slope is simply

b1+b2*mean(latitude)

where b1 is the main effect of year, and b2 the year by latitude interaction.

If you want to work with the average of the exponentiated slopes (i.e. the expected proportional change in abundance per year) this would require more work, but I�m not sure that is what you want?

Cheers,

Jarrod




From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of Chris Oosthuizen <w.chris.oosthuizen using gmail.com>
Date: Thursday, 14 September 2023 at 16:06
To: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
Subject: [R-sig-ME] Predicting the overall population trend from a GLMM with many sites (MCMCglmm)
This email was sent to you by someone outside the University.
You should only click on links or attachments if you are certain that the email is genuine and the content is safe.

Dear all,

I have intermittent population counts of a single species at many sites
that are situated at different latitudes. I want to estimate the overall
population trend across all sites (how much did the overall population
decline in 30 years?). I would appreciate clarification on precisely how I
should calculate a single, overall population trend estimate from a
random-slope mixed model.

The model we fit is:
model <- MCMCglmm(count ~ year * latitude, random = ~us(1 + year):site,
                rcov=~units, family="poisson", pr = T, data = df�)

If I were interested in site-level effects I would predict with random
effects (set marginal to NULL):
pred <- data.frame(predict(model, newdata=df, type="response",
marginal=NULL, interval="confidence", posterior="all"))

We can also marginalise the random effects � this prediction is �at the
population level� and it may have poor site-level predictions if the random
effect is important (e.g. some sites are increasing and others are
decreasing):
pred_margin <- data.frame(predict(model, newdata=df, type="response",
marginal=model$Random$formula, interval="confidence", posterior="all"))

Both predictions above are made using a �newdata� data frame that contains
�site� � this seems to be required by the �predict� function. Thus, we
always end up with site-level data as output in the prediction data frame,
regardless of whether the random effects were marginalised or not.

I don�t know the best statistical approach to obtain a single population
trend estimate. My approach would be to predict with random effects, and
sum the counts per site to get an overall estimate. This is the approach
that a colleague took (although he disregarded the random effects). His
calculations do not use the �predict� function, but directly calculates
population sizes from the model coefficients, for every posterior draw (to
propagate uncertainty). He then estimates population sizes at two points
(in 1990 and 2020) and asks whether the population decreased between these
two points.
This approach and reproducible code is here:
https://drive.google.com/drive/folders/1_thQQW9Vy3SxGi6WBkk-1CuqSpE6XaYy?usp=sharing

It is important that I get this right, so I hope someone can give me
guidance. Thanks in advance!
Chris

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. Is e buidheann carthannais a th� ann an Oilthigh Dh�n �ideann, cl�raichte an Alba, �ireamh cl�raidh SC005336.

        [[alternative HTML version deleted]]

	[[alternative HTML version deleted]]