library(sf)
library(sp)
library(ggplot2)

First just loading the stuff we need:

source("../R/Model_visualization_functions.R")
model_0 <- readRDS("../R/output/model_0.RDS")
model_1 <- readRDS("../R/output/model_1.RDS")
model_area <- readRDS("../R/output/model_area.RDS")
#model_logarea <- readRDS("../R/output/model_logarea.RDS")

stk.pred <- readRDS("../R/output/stkpred.RDS")
Projection <- CRS("+proj=longlat +ellps=WGS84")

I ran three models where the only change was in the E used. All code for fitting the models is on my Github-page, here is the main script used. I also wanted to do a model with E = log_area and one where I didn’t specify E at all (so that it would just use the default, in order to check if this gave the exact same model as the one with E = rep(1, ...)), but both these models caused INLA to crash, and I decided to not try to figure out the issue now.

Here are the means of the parameter estimates for all the models:

mean_df <- data.frame(model_0 = model_0$model$summary.fixed$mean,
            model_1 = model_1$model$summary.fixed$mean,
            model_area = model_area$model$summary.fixed$mean)
row.names(mean_df) <- row.names(model_0$model$summary.fixed['mean'])
mean_df
##                        model_0       model_1    model_area
## int.survey        1.860150e-01  6.239805e+00  9.570329e+00
## int.artsobs       8.102272e-01  5.448119e+00  9.439464e+00
## decimalLongitude -6.878653e-04 -1.940457e-01 -4.710899e-02
## decimalLatitude   1.474794e-04 -5.056418e-02 -1.486342e-01
## log_area          1.994093e-03  1.509153e-02 -6.800193e-01
## perimeter_m       7.303207e-08 -2.619230e-08 -4.732566e-06
## eurolst_bio10    -5.472810e-05 -1.661580e-03  5.426807e-03
## SCI               1.681673e-03  2.885602e-02  1.797914e-01
## distance_to_road  3.732053e-07  9.276972e-07  9.789153e-05
## HFP               4.876571e-05 -2.741610e-04 -6.101042e-03

And here are the standard deviations of the parameter estimates:

sd_df <- data.frame(model_0 = model_0$model$summary.fixed$sd,
            model_1 = model_1$model$summary.fixed$sd,
            model_area = model_area$model$summary.fixed$sd)
row.names(sd_df) <- row.names(model_0$model$summary.fixed['sd'])
sd_df
##                       model_0      model_1   model_area
## int.survey       4.973408e-02 1.019556e+01 2.310702e+00
## int.artsobs      6.013037e-02 1.179417e+01 2.402444e+00
## decimalLongitude 1.964126e-03 1.143198e-01 7.438028e-02
## decimalLatitude  9.563549e-04 2.987900e-02 4.409913e-02
## log_area         7.732439e-03 2.057746e-02 2.902495e-02
## perimeter_m      8.402543e-07 1.023414e-06 1.435191e-06
## eurolst_bio10    4.549567e-04 3.893602e-03 8.029512e-03
## SCI              3.022715e-02 1.144494e-01 1.478734e-01
## distance_to_road 3.653346e-06 1.950357e-05 2.779279e-05
## HFP              2.030865e-03 5.120807e-03 7.977416e-03

And here are the log-marginal likelihoods:

mlik_df <- data.frame(model_0 = model_0$model$mlik,
            model_1 = model_1$model$mlik,
            model_area = model_area$model$mlik)
mlik_df
##                                        model_0   model_1 model_area
## log marginal-likelihood (integration) 183.6986 -961.4921  -1271.853
## log marginal-likelihood (Gaussian)    183.6511 -961.9546  -1271.775

Here are the plotted predictions for each of the models:

Finally, here are the model summaries, just to see what the full model was: (observation that I don’t know how to interpret: only for model_0 has the dic been calculated.)

summary(model_0$model)
## 
## Call:
##    c("INLA::inla(formula = Formula, family = c(\"poisson\", \"binomial\"), 
##    ", " data = INLA::inla.stack.data(stck), E = 
##    INLA::inla.stack.data(stck)$e, ", " offset = offset.formula, Ntrials = 
##    INLA::inla.stack.data(stck)$Ntrials, ", " verbose = verbose, 
##    control.compute = list(waic = waic, dic = dic), ", " control.predictor 
##    = list(A = INLA::inla.stack.A(stck), link = NULL, ", " compute = TRUE), 
##    control.family = list(list(link = \"log\"), ", " list(link = 
##    \"cloglog\")), control.results = list(return.marginals.random = FALSE, 
##    ", " return.marginals.predictor = FALSE), control.fixed = 
##    control.fixed)" ) 
## Time used:
##     Pre = 4.48, Running = 513, Post = 0.475, Total = 518 
## Fixed effects:
##                    mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## int.survey        0.186 0.050      0.088    0.186      0.283  0.186   0
## int.artsobs       0.810 0.060      0.692    0.810      0.928  0.810   0
## decimalLongitude -0.001 0.002     -0.005   -0.001      0.003 -0.001   0
## decimalLatitude   0.000 0.001     -0.002    0.000      0.002  0.000   0
## log_area          0.002 0.008     -0.013    0.002      0.017  0.002   0
## perimeter_m       0.000 0.000      0.000    0.000      0.000  0.000   0
## eurolst_bio10     0.000 0.000     -0.001    0.000      0.001  0.000   0
## SCI               0.002 0.030     -0.058    0.002      0.061  0.002   0
## distance_to_road  0.000 0.000      0.000    0.000      0.000  0.000   0
## HFP               0.000 0.002     -0.004    0.000      0.004  0.000   0
## 
## Random effects:
##   Name     Model
##     shared_field SPDE2 model
##    bias_field SPDE2 model
## 
## Model hyperparameters:
##                           mean    sd 0.025quant 0.5quant 0.975quant   mode
## Theta1 for shared_field  1.512 0.370      0.785    1.512      2.241  1.513
## Theta2 for shared_field -1.298 0.355     -2.008   -1.294     -0.611 -1.279
## Theta1 for bias_field    0.455 0.285     -0.100    0.452      1.020  0.444
## Theta2 for bias_field   -1.393 0.310     -2.017   -1.386     -0.800 -1.360
## 
## Expected number of effective parameters(stdev): 35.11(3.94)
## Number of equivalent replicates : 105.84 
## 
## Deviance Information Criterion (DIC) ...............: -2248.43
## Deviance Information Criterion (DIC, saturated) ....: NaN
## Effective number of parameters .....................: 5.05
## 
## Marginal log-Likelihood:  183.65 
## Posterior marginals for the linear predictor and
##  the fitted values are computed
summary(model_1$model)
## 
## Call:
##    c("INLA::inla(formula = Formula, family = c(\"poisson\", \"binomial\"), 
##    ", " data = INLA::inla.stack.data(stck), E = 
##    INLA::inla.stack.data(stck)$e, ", " offset = offset.formula, Ntrials = 
##    INLA::inla.stack.data(stck)$Ntrials, ", " verbose = verbose, 
##    control.compute = list(waic = waic, dic = dic), ", " control.predictor 
##    = list(A = INLA::inla.stack.A(stck), link = NULL, ", " compute = TRUE), 
##    control.family = list(list(link = \"log\"), ", " list(link = 
##    \"cloglog\")), control.results = list(return.marginals.random = FALSE, 
##    ", " return.marginals.predictor = FALSE), control.fixed = 
##    control.fixed)" ) 
## Time used:
##     Pre = 3.91, Running = 301, Post = 0.545, Total = 305 
## Fixed effects:
##                    mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## int.survey        6.240 10.196    -20.415    6.632     27.383  5.939   0
## int.artsobs       5.448 11.794    -23.716    6.028     29.526  5.395   0
## decimalLongitude -0.194  0.114     -0.447   -0.183     -0.002 -0.156   0
## decimalLatitude  -0.051  0.030     -0.119   -0.047     -0.001 -0.041   0
## log_area          0.015  0.021     -0.025    0.015      0.055  0.015   0
## perimeter_m       0.000  0.000      0.000    0.000      0.000  0.000   0
## eurolst_bio10    -0.002  0.004     -0.009   -0.002      0.006 -0.002   0
## SCI               0.029  0.114     -0.197    0.029      0.252  0.031   0
## distance_to_road  0.000  0.000      0.000    0.000      0.000  0.000   0
## HFP               0.000  0.005     -0.010    0.000      0.010  0.000   0
## 
## Random effects:
##   Name     Model
##     shared_field SPDE2 model
##    bias_field SPDE2 model
## 
## Model hyperparameters:
##                          mean    sd 0.025quant 0.5quant 0.975quant  mode
## Theta1 for shared_field  1.16 0.287      0.608     1.16       1.74  1.14
## Theta2 for shared_field -2.87 0.759     -4.442    -2.83      -1.45 -2.72
## Theta1 for bias_field    1.39 0.321      0.830     1.37       2.08  1.27
## Theta2 for bias_field   -2.50 0.713     -4.108    -2.41      -1.37 -2.04
## 
## Expected number of effective parameters(stdev): 28.37(3.57)
## Number of equivalent replicates : 130.97 
## 
## Marginal log-Likelihood:  -961.96 
## Posterior marginals for the linear predictor and
##  the fitted values are computed
summary(model_area$model)
## 
## Call:
##    c("INLA::inla(formula = Formula, family = c(\"poisson\", \"binomial\"), 
##    ", " data = INLA::inla.stack.data(stck), E = 
##    INLA::inla.stack.data(stck)$e, ", " offset = offset.formula, Ntrials = 
##    INLA::inla.stack.data(stck)$Ntrials, ", " verbose = verbose, 
##    control.compute = list(waic = waic, dic = dic), ", " control.predictor 
##    = list(A = INLA::inla.stack.A(stck), link = NULL, ", " compute = TRUE), 
##    control.family = list(list(link = \"log\"), ", " list(link = 
##    \"cloglog\")), control.results = list(return.marginals.random = FALSE, 
##    ", " return.marginals.predictor = FALSE), control.fixed = 
##    control.fixed)" ) 
## Time used:
##     Pre = 4.28, Running = 293, Post = 0.784, Total = 298 
## Fixed effects:
##                    mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## int.survey        9.570 2.311      6.047    9.208     15.040  8.465   0
## int.artsobs       9.439 2.402      5.732    9.079     15.089  8.341   0
## decimalLongitude -0.047 0.074     -0.192   -0.049      0.106 -0.052   0
## decimalLatitude  -0.149 0.044     -0.251   -0.142     -0.079 -0.130   0
## log_area         -0.680 0.029     -0.737   -0.680     -0.623 -0.680   0
## perimeter_m       0.000 0.000      0.000    0.000      0.000  0.000   0
## eurolst_bio10     0.005 0.008     -0.010    0.005      0.021  0.005   0
## SCI               0.180 0.148     -0.113    0.181      0.467  0.183   0
## distance_to_road  0.000 0.000      0.000    0.000      0.000  0.000   0
## HFP              -0.006 0.008     -0.022   -0.006      0.009 -0.006   0
## 
## Random effects:
##   Name     Model
##     shared_field SPDE2 model
##    bias_field SPDE2 model
## 
## Model hyperparameters:
##                           mean    sd 0.025quant 0.5quant 0.975quant   mode
## Theta1 for shared_field -1.894 0.172     -2.225   -1.897     -1.547 -1.909
## Theta2 for shared_field  0.128 0.230     -0.335    0.133      0.567  0.151
## Theta1 for bias_field   -1.063 0.287     -1.596   -1.075     -0.472 -1.119
## Theta2 for bias_field   -0.471 0.324     -1.144   -0.455      0.126 -0.396
## 
## Expected number of effective parameters(stdev): 140.30(11.65)
## Number of equivalent replicates : 26.48 
## 
## Marginal log-Likelihood:  -1271.78 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Notice how the longitude and latitude are more significant in the non-zero-E models!

dots_whiskers_inla(model_0)

dots_whiskers_inla(model_1)

dots_whiskers_inla(model_area)