How can I calculate the difference between the marginals.fitted.values of two INLA models?

Hello,

The marginals.fitted.values of an INLA model gives you a list of marginal posterior distributions, not point estimates, which is why simple subtraction isn't working.

As you've correctly understood, taking the difference between two marginal distributions isn't as straightforward as subtracting means or modes. Instead, you can perform the operation on the whole distributions.

One possible approach to your problem could be to sample from these distributions, subtract these samples, and then look at the distribution of these differences.

Below is a basic sketch of how you might do this in R:

# Number of samples to draw from each marginal posterior
num_samples <- 1000

# Function to draw samples from a given marginal posterior
draw_samples <- function(marginal, n) {  q <- seq(0, 1, length.out = n)  quantile <- sapply(marginal$quantile, function(x) approx(x$x, x$y, q)$y)  sample(quantile, n, replace = TRUE)
}

# Draw samples from each marginal fitted value for both models
samples1 <- sapply(model1$marginals.fitted.values, draw_samples, n = num_samples)
samples2 <- sapply(model2$marginals.fitted.values, draw_samples, n = num_samples)

# Subtract the samples
differences <- samples1 - samples2

# Now, each column of 'differences' represents the distribution of differences
# for a given location.
# You can summarize these distributions (e.g., calculate the mean or median difference)
# or visualize them to understand the change in prevalence.

This script draws samples from each marginal posterior of the fitted values, then subtracts these samples. This gives you a distribution of differences for each location. You could then summarize these distributions (e.g., by calculating the mean or median) or visualize them to understand how the disease prevalence has changed.

Remember, each marginal distribution in the marginals.fitted.values list is a different location, so this will give you a separate distribution of differences for each location.

This is a somewhat advanced topic and might seem difficult if you're new to Bayesian statistics or spatial modeling. I'd recommend discussing this approach with a mentor or statistician if you can, and cross-validating the results with other forms of analysis if possible.