INLA model with binary response, one predictor (continuous on 0-1 scale) and spatial term - issue with the model outcome

[Anna Szyniszewska]

Hi,

As described in the subject, I built an INLA model with a binary response (`cb`), one predictor (`oldr`) (continuous on 0-1 scale) and a spatial term to account for strong spatial autocorrelation in the data.

My points are distributed in linear clusters (along the roads) and although I am new to INLA I hope I set the model right. 

My questions are:

1) Are my priors set up correctly?

2) Is the model (result1) written correctly?

3) Why in the following lines of the code (at the end of the full code and also pasted below for a reference) I am getting NULL making it impossible to look further at model results? Are my model outcomes rubbish or it does not run fully?

result1.res <- inla.spde2.result(result1, 'spatial.field', spde, do.transf=TRUE)

result1.res$summary.log.range.nominal

Thanks!

Full code is here:

# Data 

n <- 125 # number of observations 

lon <- c(22.58,22.5,22.37,22.21,22.02,21.98,21.94,20.97,21.08,21.15,21.37,21.51,21.58,21.66,20.89,20.87,20.79,20.77,20.82,20.88,21.03,21.15,21.34,25.24,25.31,25.34,25.38,25.39,25.36,25.33,25.3,25.29,25.24,25.17,25.16,25.16,25.23,25.26,25.22,25.39,25.43,25.45,25.51,25.57,25.63,25.68,25.76,25.44,25.52,25.62,25.72,25.81,25.89,25.97,26.05,22.76,27.9,27.87,27.84,27.81,27.78,27.75,27.72,27.69,27.66,27.59,27.63,27.71,27.7,27.69,27.72,27.77,27.75,27.72,27.67,27.71,27.73,27.82,27.89,27.63,27.67,27.76,27.92,28.03,28.12,28.23,28.39,28.7,28.77,28.93,29.28,29.33,29.22,29.13,29.21,29.24,25.6,25.6,25.61,25.62,25.56,25.57,25.59,25.58,25.6,25.69,25.8,25.97,26.1,26.17,26.17,26.21,26.22,26.22,26.29,26.4,26.16,26.15,26.12,26.15,25.95,25.83,25.7,25.65,25.59)

lat <- c(-18.71,-18.67,-18.62,-18.6,-18.43,-18.37,-18.3,-18.04,-18.03,-18.04,-18.03,-18.07,-18.1,-18.15,-18.13,-18.26,-18.54,-19.06,-19.8,-19.69,-19.54,-19.51,-19.48,-18.06,-17.74,-17.65,-17.57,-17.45,-17.37,-17.27,-17.19,-17.1,-16.98,-16.91,-16.78,-16.69,-15.82,-16.22,-16.29,-17.41,-17.32,-17.24,-17.16,-17.07,-16.98,-16.83,-16.73,-17.51,-17.56,-17.57,-17.57,-17.58,-17.62,-17.66,-17.64,-18.66,-18.13,-18.03,-17.93,-17.83,-17.73,-17.63,-17.53,-17.43,-17.41,-17.29,-17.06,-16.59,-16.75,-16.87,-16.46,-16.38,-16.28,-16.16,-16,-15.84,-15.68,-15.52,-15.38,-15.1,-15.17,-15.17,-15.31,-15.41,-15.41,-15.41,-15.36,-15.42,-15.46,-15.58,-15.64,-15.7,-15.7,-15.83,-15.99,-16.14,-14.87,-14.87,-14.86,-14.86,-14.85,-14.83,-14.8,-14.78,-14.75,-14.87,-14.76,-14.74,-14.73,-14.75,-14.71,-14.71,-14.7,-14.75,-14.74,-14.66,-14.81,-14.86,-14.83,-14.8,-15.22,-15.36,-15.41,-15.38,-15.34)

cb <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,0,1,1,1,1,1,0,1,1,1)

oldr <- c(0,0.01,0,0,0,0,0,0.2,0.13,0.05,0.03,0.03,0.03,0.01,0.21,0.03,0.01,0,0,0.01,0,0,0,0.02,0,0.01,0.06,0.07,0.04,0.04,0.11,0.15,0.29,0.4,0.59,0.64,0.67,0.63,0.63,0.07,0.07,0.08,0.06,0.07,0.2,0.14,0.12,0.17,0.02,0.03,0.01,0.02,0.01,0.12,0.1,0.01,0.71,0.39,0.12,0.13,0.08,0.06,0.06,0.04,0.04,0.13,0.02,0.63,0.1,0.08,0.62,0.36,0.19,0.16,0.06,0.09,0.04,0.09,0.19,0.26,0.16,0.35,0.19,0.06,0.17,0.13,0.13,0.19,0.19,0.11,0.99,0.99,0.99,0.32,0.28,0.35,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.74,0.74,0.72,0.74,0.78,0.84,0.85,0.86,0.86,0.7,0.78,0.78,0.69,0.78,0.78,0.78,0.78,0.3,0.27,0.75,0.62,0.8)

 

# Libraries

library(R2jags) 

library(INLA) 

library(ggplot2) 

library(reshape)

loc <- as.matrix(data.frame(lon,lat)) 

mesh = inla.mesh.create.helper(points.domain=loc,

                               max.edge=c(0.2, 0.5))

proj.obs = inla.mesh.projector(mesh, loc=loc) 

proj.pred = inla.mesh.projector(mesh, loc=mesh$loc) 

spde = inla.spde2.matern(mesh,

                         B.tau=cbind(log(1), 1, 0),

                         B.kappa=matrix(c(log(sqrt(8)/0.2), 0, 1), 1, 3))

y = cb 

covar = oldr 

field <- inla.qsample(n=1, Q=inla.spde2.precision(spde, theta=c(0,0)))[,1] 

A.obs = inla.spde.make.A(mesh, loc=loc)

A.pred = inla.spde.make.A(mesh, loc=proj.pred$loc)

stack.obs =

  inla.stack(data=list(y=y),

             A=list(A.obs, 1),

             effects=list(c(list(intercept=rep(1, mesh$n)),

                            inla.spde.make.index("spatial", spde$n.spde)),

                          covar=covar),

             tag="obs")

formula = y ~ -1 + intercept + covar + f(spatial, model=spde)

result1 = inla(formula,

               data=inla.stack.data(stack.obs, spde=spde),

               family="binomial",

               Ntrials = n,

               #control.family=list(link='logit'),

               control.predictor=list(A=inla.stack.A(stack.obs), compute=TRUE),

               verbose=TRUE)

plot(y,result1$summary.fitted.values[inla.stack.index(stack.obs,"obs")$data,"mean"])

result1$summary.fix

result1$summary.hyperpar

## Problem here - is my model wrongly drafted or the results are rubbish?

result1.res <- inla.spde2.result(result1, 'spatial.field', spde, do.transf=TRUE)

result1.res$summary.log.range.nominal