Partial multivariate normal prior for fixed effects

[dsjohnson]

Hi all,
I noticed that there is an option to place a multivariate normal prior on the fixed effects, e.g.,

sdmTMB(density ~ 0 + depth_scaled + depth_scaled2 + as.factor(year),
    data = d, mesh = pcod_spde, family = tweedie(link = "log"),
    priors = sdmTMBpriors(b = mvnormal(rep(0, 6), diag(1, 6))),
    spatial = "off", spatiotemporal = "off"
  )

Is it also possible to to have a partial prior, say b[1:3] ~ mvn(), but b[-c(1:3)] has no (flat) prior? If so, I assume the mean is specified with NAs like the univariate specification, but how is the covariate matrix specified?
Best wishes,
Devin

[ericward-noaa]

Yes, see this example using the 'pcod' dataset. There's 4 years (2011, 2013, 2015, 2017) and I'm leaving one with no prior (2017) while increasing the prior information for the others. I should preface this too by saying that these are informative with respect to the variances -- not covariances. You can use informative covariance matrices -- in which case you'd replace any of the diag() below with a scale matrix (used as Sigma here:

sdmTMB/R/priors.R

Line 125 in a00e32d

mvnormal <- function(location = 0, scale = diag(length(location))) {

)

library(ggplot2)
library(sdmTMB)
d <- pcod_2011
d$fyear <- as.factor(d$year)

pcod_spde <- pcod_mesh_2011

no_prior <- sdmTMB(density ~ 0 + fyear,
                data = d, mesh = pcod_spde, family = tweedie(link = "log"),
                spatial = "off", spatiotemporal = "off"
)

prior_1 <- sdmTMB(density ~ 0 + fyear,
                data = d, mesh = pcod_spde, family = tweedie(link = "log"),
                priors = sdmTMBpriors(b = mvnormal(c(3.7,3.7,3.7,NA), diag(c(0.1,0.1,0.1,NA), 4))),
                spatial = "off", spatiotemporal = "off"
)

prior_01 <- sdmTMB(density ~ 0 + fyear,
                   data = d, mesh = pcod_spde, family = tweedie(link = "log"),
                   priors = sdmTMBpriors(b = mvnormal(c(3.7,3.7,3.7,NA), diag(c(0.01,0.01,0.01,NA), 4))),
                   spatial = "off", spatiotemporal = "off"
)

prior_001 <- sdmTMB(density ~ 0 + fyear,
                   data = d, mesh = pcod_spde, family = tweedie(link = "log"),
                   priors = sdmTMBpriors(b = mvnormal(c(3.7,3.7,3.7,NA), diag(c(0.001,0.001,0.001,NA), 4))),
                   spatial = "off", spatiotemporal = "off"
)

df <- rbind(tidy(no_prior),
            tidy(prior_1),
            tidy(prior_01),
            tidy(prior_001))
df$prior <- c(rep("NA",4), rep("0.1",4), rep("0.01",4), rep("0.001",4))

df$prior <- factor(df$prior,levels = c("NA","0.1","0.01","0.001"))
ggplot(df, aes(prior, estimate)) + 
  geom_hline(aes(yintercept=3.7,col="red")) + 
  geom_pointrange(aes(ymin=conf.low, ymax = conf.high)) + 
  facet_wrap(~ term, scale="free_y") 

I think the other way to see these changes is by inspecting the covariance matrices. E.g.

cov2cor(no_prior$sd_report$cov.fixed[1:4,1:4])
cov2cor(prior_001$sd_report$cov.fixed[1:4,1:4])

[dsjohnson]

Great, thanks @ericward-noaa !