Plotting CAR(1) trend smooths

[npatt01]

Hello,

I have a long running dataset of pelagic surveys recording species counts. The surveys run ~monthly, but many months are missing, or may have multiple trips undertaken. Sometimes multiple trips can even run on the same day. The trips run from select ports to the same survey location. I am modelling single species abundance through time with:

mvgam(formula = count ~ - 1 +
                     s(ws, k = 6, bs = 'cr') +
                     s(port, bs = 're'),
                  trend_formula = ~ 1 +
                    s(time, k = 6, bs = 'cr') +
                    s(month, bs = 'cc', k = 12) +
                    te(wd, ws, bs = c('cr', 'cr'), k = c(8, 4)) +
                    s(sst, k = 4, bs = 'cr') +
                    s(sst_anom, k = 4, bs = 'cr') +
                    s(soi, k = 4, bs = 'cr') +
                    s(siod, k = 4, bs = 'cr') +
                    s(sam, k = 4, bs = 'cr'),
                  trend_model = CAR(1),
                  data = data,
                  family = poisson(),
                  noncentred = T,
                  samples = 1000,
                  control = list(adapt_delta = 0.99, max_treedepth = 20),
                  parallel = T)

The issue I'm running into is that when using conditional_effects() to produce plots of time:

conditional_effects(m1, effects = 'time', type = 'expected')

I get odd returns:

I'm guessing this is due to the CAR(1) trend model, as I wasn't encountering it when running a discrete time model on a subsetted version of the same dataset. In that dataset, time was evenly spaced and surveys were undertaken from a single port. I wasn't including a trend_model or re term in those models. In that case conditional_effects() plots looked like:

I am aiming to produce hindcast() and conditional_effects() plots to present the species trends through time — any help would be much appreciated :)

FYI model residuals:

and summary:

GAM observation formula:
count ~ s(ws, k = 6, bs = "cr") + s(port, bs = "re") - 1

GAM process formula:
~1 + s(time, k = 6, bs = "cr") + s(month, bs = "cc", k = 12) + 
    te(wd, ws, bs = c("cr", "cr"), k = c(8, 4)) + s(sst, k = 4, 
    bs = "cr") + s(sst_anom, k = 4, bs = "cr") + s(soi, k = 4, 
    bs = "cr") + s(siod, k = 4, bs = "cr") + s(sam, k = 4, bs = "cr")

Family:
poisson

Link function:
log

Trend model:
CAR(p = 1)


N process models:
1 

N series:
1 

N timepoints:
455 

Status:
Fitted using Stan 
4 chains, each with iter = 1500; warmup = 500; thin = 1 
Total post-warmup draws = 4000


GAM observation model coefficient (beta) estimates:
           2.5%    50% 97.5% Rhat n_eff
s(ws).1   -1.10 -0.067  0.78 1.00   793
s(ws).2   -0.92 -0.058  0.64 1.01   670
s(ws).3   -0.77 -0.050  0.61 1.01   620
s(ws).4   -1.60 -0.053  1.30 1.01   668
s(ws).5   -0.82  0.070  1.10 1.00   774
s(port).1 -2.60 -0.360  1.80 1.00  2956
s(port).2 -3.20 -0.790  1.40 1.00  2514
s(port).3 -2.50 -0.200  2.10 1.00  2895
s(port).4 -3.90 -0.940  1.30 1.00  1923

GAM observation model group-level estimates:
               2.5%   50% 97.5% Rhat n_eff
mean(s(port)) -2.50 -0.49   1.5    1  3884
sd(s(port))    0.11  0.46   2.2    1  1155

Approximate significance of GAM observation smooths:
        edf Ref.df Chi.sq p-value
s(ws)     5      5   0.07    1.00
s(port)   4      4   1.14    0.89

Process model AR parameter estimates:
        2.5%  50% 97.5% Rhat n_eff
ar1[1] 0.017 0.33  0.87    1  1181

Process error parameter estimates:
         2.5% 50% 97.5% Rhat n_eff
sigma[1]  2.8 3.5   4.3    1   949

GAM process model coefficient (beta) estimates:
                      2.5%      50%  97.5% Rhat n_eff
(Intercept)_trend   -8.700 -5.70000 -3.000 1.00  1202
s(time).1_trend     -0.440  0.25000  1.100 1.00  1338
s(time).2_trend     -0.500  0.49000  1.600 1.01   930
s(time).3_trend     -0.560  0.36000  1.400 1.00   921
s(time).4_trend     -0.620  0.10000  0.810 1.00  1359
s(time).5_trend     -1.700 -0.46000  0.570 1.00  1035
s(month).1_trend    -6.500 -3.50000 -1.800 1.00   715
s(month).2_trend    -3.000 -0.96000  0.850 1.00  1963
s(month).3_trend     1.100  2.80000  5.500 1.00   645
s(month).4_trend     4.100  6.10000  9.400 1.01   484
s(month).5_trend     4.100  6.00000  8.600 1.01   480
s(month).6_trend     3.100  4.90000  7.400 1.01   404
s(month).7_trend     1.300  3.10000  5.500 1.01   541
s(month).8_trend    -0.210  1.70000  4.500 1.01   699
s(month).9_trend    -2.600 -0.69000  1.700 1.00  1507
s(month).10_trend   -5.100 -2.90000 -1.100 1.00  1456
te(wd,ws).1_trend   -0.740  0.04200  1.100 1.01   391
te(wd,ws).2_trend   -0.840  0.02000  0.920 1.01   747
te(wd,ws).3_trend   -2.700 -0.11000  1.300 1.02   176
te(wd,ws).4_trend   -0.530  0.00021  0.530 1.00   482
te(wd,ws).5_trend   -1.200 -0.07700  0.500 1.01   226
te(wd,ws).6_trend   -0.510 -0.03300  0.280 1.01   382
te(wd,ws).7_trend   -0.790 -0.00520  0.670 1.02   208
te(wd,ws).8_trend   -0.660 -0.03600  0.470 1.00   420
te(wd,ws).9_trend   -0.680 -0.04800  0.330 1.01   299
te(wd,ws).10_trend  -0.530 -0.03600  0.310 1.00   427
te(wd,ws).11_trend  -0.630  0.03400  1.000 1.00   378
te(wd,ws).12_trend  -0.910 -0.07000  0.460 1.01   227
te(wd,ws).13_trend  -0.910 -0.06400  0.430 1.01   267
te(wd,ws).14_trend  -0.570 -0.03300  0.370 1.01   488
te(wd,ws).15_trend  -0.670  0.07400  1.600 1.01   197
te(wd,ws).16_trend  -1.100 -0.06900  0.480 1.01   145
te(wd,ws).17_trend  -0.790 -0.05300  0.380 1.01   192
te(wd,ws).18_trend  -0.760 -0.04200  0.470 1.01   366
te(wd,ws).19_trend  -0.710  0.09000  2.000 1.02   119
te(wd,ws).20_trend  -0.910 -0.03400  0.490 1.01   207
te(wd,ws).21_trend  -0.610 -0.02400  0.360 1.01   273
te(wd,ws).22_trend  -0.540 -0.01600  0.410 1.01   471
te(wd,ws).23_trend  -0.660  0.07700  2.000 1.02   102
te(wd,ws).24_trend  -0.620  0.01300  0.650 1.00   761
te(wd,ws).25_trend  -0.270  0.03300  0.440 1.00   738
te(wd,ws).26_trend  -1.100 -0.04000  0.680 1.01   305
te(wd,ws).27_trend  -0.640  0.02500  1.200 1.01   187
te(wd,ws).28_trend  -0.870  0.14000  1.900 1.01   182
te(wd,ws).29_trend  -0.640  0.08800  1.100 1.01   443
te(wd,ws).30_trend  -0.710  0.03500  0.830 1.01   780
te(wd,ws).31_trend  -3.200 -0.11000  1.300 1.03   107
s(sst).1_trend      -0.690 -0.00680  0.660 1.00  1869
s(sst).2_trend      -1.100  0.00340  1.200 1.00  1854
s(sst).3_trend      -1.000 -0.02400  0.910 1.00  1774
s(sst_anom).1_trend -1.000  0.06300  1.400 1.00  1084
s(sst_anom).2_trend -1.000  0.07100  1.500 1.00  1210
s(sst_anom).3_trend -1.700 -0.08200  1.300 1.00  1193
s(soi).1_trend      -0.044  2.10000  5.700 1.01   513
s(soi).2_trend      -0.012  2.20000  5.700 1.01   475
s(soi).3_trend      -6.800 -2.60000  0.082 1.01   481
s(siod).1_trend     -2.000 -0.18000  0.980 1.00   734
s(siod).2_trend     -1.400 -0.07200  1.000 1.00   773
s(siod).3_trend     -1.000  0.26000  2.300 1.00   801
s(sam).1_trend      -0.860  0.04500  1.000 1.00   937
s(sam).2_trend      -0.620  0.07200  0.940 1.00   992
s(sam).3_trend      -1.000  0.05800  1.300 1.00  1198

Approximate significance of GAM process smooths:
            edf Ref.df  Chi.sq p-value    
s(time)       5      5   26.84    0.75    
s(month)     10     10 4365.80 8.1e-06 ***
te(wd,ws)    31     31    0.71    1.00    
s(sst)        3      3    0.02    1.00    
s(sst_anom)   3      3    0.14    1.00    
s(soi)        3      3  261.94    0.18    
s(siod)       3      3    3.99    0.93    
s(sam)        3      3    0.73    0.99    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Stan MCMC diagnostics:
✔ No issues with effective samples per iteration
✔ Rhat looks good for all parameters
✔ No issues with divergences
✔ No issues with maximum tree depth

[nicholasjclark]

Hi @npatt01, thanks for the post. I can't say exactly what is happening, but it could be that the smooth of time in your trend_formula is competing a lot with the CAR() process. Sometimes it can be very hard to estimate multiple effects of time in the same model. Have you tried using gratia::draw() to look at the smooth and see how it is behaving? The very large value of sigma could be an indication of this as well. Do you feel like you need the CAR(1)? An alternative could be to use a gp() of time in place of the smooth and not include the CAR() process, just to see what that looks like. Or if you do want to keep the autoregressive component, you may do well to impose a more regularizing prior on sigma to reign it in and not let it jump all over the place.

[npatt01]

Thanks @nicholasjclark! I do not think the CAR process is needed, but have been including it to allow for uneven sampling (i.e. to resolve Error: One or more series in data is missing observations for one or more timepoints). gratia::draw() produces reasonable looking smooths (I think):

I attempted removing the CAR(1) process and substituting the s() on time for gp(), but had the same missing time points error message. Is there another way of including uneven sampling without instituting the CAR process?

[nicholasjclark]

Hi @npatt01. At present this can't be done if you are using a trend_formula, but you could put all your effects into the main formula and not use a State Space model. This would allow you to use gp() effects of time and would probably get you closer to what you want I think

[npatt01]

Thank you for your help @nicholasjclark! Sorry for my delay in following up from your comment and apologies in advance for my wandering response!

I restructured the model so everything was in the main_formula with a gp() on time. I still found I needed to have evenly spaced samples though (i.e. I still recieved Error: One or more series in data is missing observations for one or more timepoints), so I included the CAR(1) trend_model. The gp() on time dramatically improved model hindcasts, but the gp() smooth of time appeared to be misbehaving given the shape of the data:

plot_mvgam_series(data = data, y = "count")

Residuals:

Gratia::draw(mod_gp) (not sure why time is shown as s(time) with a TPRS basis as it is a gp())

conditional_effects(mod_gp, effects = 'time', type = 'expected')

What I am hoping to achieve is this:

Where the top plot shows a response scale smooth on time with the first derivative included to highlight periods of significant change. The bottom plot shows the median of the models hindcast to give better insight into the seasonal patterns observed. My idea in presenting both outputs (i.e. the smooth and the hindcast) is that a reader is given both the overarching trend observed in a given species abundance (smooth), while also allowed a more nuanced insight into the species seasonal highs and lows (hindcast).

The model represented in those plots was produced used a subset of my survey data (i.e. removing instances where more than one survey took place in a given month), with NA's introduced into the response during missing time periods. Now that I am attempting to include all surveys I am running the issues with producing the response scale smooths when running conditional_effects(m1, effects = 'time', type = 'expected').

If there is another way of integrating the unevenly sampled data without relying on a CAR trend_model that would amazing. Otherwise, it might be best to return to subsetting the data to be evenly sampled.