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)