timcdlucas / INLAutils

R package providing utilities for INLA
Other
26 stars 10 forks source link

inlasloo fails if lat or long is called 'y' #52

Closed timcdlucas closed 5 years ago

timcdlucas commented 5 years ago

library(sp)
data(meuse)

coords <- meuse[, c('x', 'y')] %>% scale
dataf1 <- sp::SpatialPointsDataFrame(coords = coords, data = meuse[, -c(1:2)])

mesh <- inla.mesh.2d(loc = sp::coordinates(dataf1), max.edge = c(0.2, 0.5), cutoff = 0.1)
spde <- inla.spde2.matern(mesh, alpha=2) # SPDE model is defined
A <- inla.spde.make.A(mesh, loc = sp::coordinates(dataf1)) # projector matrix
dataframe <- data.frame(dataf1) # generate dataframe with response and covariate
modform <- cadmium ~ -1 + y.intercept + ffreq + om + soil + lime + f(spatial.field, model = spde)
modform2 <- cadmium ~ -1 + y.intercept + ffreq + om + soil + lime

# make index for spatial field
s.index <- inla.spde.make.index(name="spatial.field",n.spde=spde$n.spde)

## Prepare the data
stk <- inla.stack(data=list(cadmium=dataframe$cadmium),
                      A=list(A,1), 
                      effects=list(c(s.index,list(y.intercept=1)),
                                   list(dataframe[, 7:10])),
                      tag='est')

out <- inla(modform, family = 'normal', Ntrials = 1,
            data = inla.stack.data(stk, spde = spde),
            control.predictor = list(A = inla.stack.A(stk), link = 1),
            control.compute = list(config = TRUE), 
            control.inla = list(int.strategy = 'eb'))
out.field <- inla.spde2.result(out,'spatial.field', spde, do.transf = TRUE)
range.out <- inla.emarginal(function(x) x, out.field$marginals.range.nominal[[1]])

# parameters for the SLOO process
ss <- 20 # sample size to process (number of SLOO runs)
# define the radius of the spatial buffer surrounding the removed point. 
rad <- min(range.out, max(dist(coords)) / 4) 
# Make sure it isn't bigger than 25% of the study area (see Le Rest et al.(2014))
alpha <- 0.05 # rmse and mae confidence intervals (1-alpha)

# run the function to compare both models
cv <- inlasloo(dataframe = dataframe, 
               long = 'x', lat = 'y',
               y = 'cadmium', ss = ss, 
               rad = rad, 
               modform = list(modform, modform2),
               mesh = mesh, family = 'normal',
               mae = TRUE)

because of the funky ordering of these lines.

    colnames(dataframe)[colnames(dataframe) == y] <- "y"
    colnames(dataframe)[colnames(dataframe) == long] <- "long"
    colnames(dataframe)[colnames(dataframe) == lat] <- "lat"
timcdlucas commented 5 years ago

There was a lot of mess with the tests. Don't know why the meuse tests failed on travis.