timcdlucas / INLAutils

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

inlasloo error: Intercept not found #64

Open cperk opened 1 year ago

cperk commented 1 year ago

Hello,

I'm trying to use inlasloo for my INLA model but get the error above. In your example with the meuse dataset I see you call the intercept "y.intercept". I call mine "Intercept".

Is that the issue? Do i need to re-run my model after renaming "Intercept" to "y.intercept" so that inlasloo recognizes it?

thx!!

-Carrie

timcdlucas commented 1 year ago

Hi,

Yeah this part of the function is not coded very sensibly. You don't need to re-run your model, you just need to name the intercept "y.intercept" in the formula you put into modform.

So in the meuse example you can run the original model with intercept and it'll still work as long as you then use y.intercept in modform.

dataframe <- data.frame(dataf1) #generate dataframe with response and covariate
modform<-stats::as.formula(paste('y ~ -1+ x1 + x2 + intercept + f(spatial.field, model=spde)'))
stk <- INLA::inla.stack(data=list(y=dataframe$y),A=list(A, 1),
                        effects=list(list(spatial.field=1:spde$n.spde),
                                     list(intercept=rep(1,length(dataframe$y)),covariate=dataframe[c(-1)])),tag='est')
out <- INLA::inla(modform, family='normal',Ntrials = 1, data=INLA::inla.stack.data(stk, spde=spde),
                  control.predictor = list(A=INLA::inla.stack.A(stk),link=1),
                  control.compute = list( config=TRUE),control.inla = list(int.strategy='eb'))
out.field <- INLA::inla.spde2.result(out,'spatial.field', spde, do.transf=TRUE)
range.out <- INLA::inla.emarginal(function(x) x, out.field$marginals.range.nominal[[1]])

# parameters for the SLOO process
ss <- 1#sample size to process (number of SLOO runs)
rad <- range.out*0.15#define the radius of the spatial buffer surrounding the removed point
mesh <- mesh#use the mesh of the model
dataframe <- dataframe#dataframe with response 'y' and covariates 'x1', 'x2'
dataframe$y <- round(runif(length(dataframe$y), 1, 12))#create positive discrete response
modform <- stats::as.formula(paste('y ~ -1+ y.intercept + x1 + f(spatial.field, model=spde)'))
family <- list('gamma')#one model
ntrials <- rep(round(max(dataframe$y,na.rm=TRUE)*2),length(dataframe$y))# create ntrials for Binomial family
alpha <- 0.05#rmse and mae confidence intervals (1-alpha)

# run the function
inlasloo(dataframe=dataframe, long='long', lat='lat',y= 'y', ss=ss, rad=rad,
         modform=modform, mesh=mesh, family=family,alpha=0.05,mae=TRUE,ds=TRUE,sqroot=FALSE)

We should fix this so it just uses whatever you used!

Tim

cperk commented 1 year ago

Thanks so much for the quick reply, Tim!

I had actually already started re-running the model using "y.intercept" when I sent you my first message. But I must still be doing something wrong, because inlasloo keeps crashing. Here is my code. At first I thought it might be that I have "f(w, model = spde)" in my model formula rather than "f(spatial.field, model=spde)," but even when I put "spatial.field" into the formula I use for modform, it doesn't work.

coords <- as.matrix(datafile.final2[c('x', 'y')])/1000
spde <- inla.spde2.pcmatern(
  mesh=mesh, alpha=2,
  prior.range=c(max.edge/10, 0.01), ### P(range < max.edge/10) = 0.01
  prior.sigma=c(1, 0.01), ### P(sigma > 1) = 0.01
  constr=TRUE)
A.s <- inla.spde.make.A(
  mesh=mesh, loc=coords)
N <- nrow(datafile.final2)

X <- data.frame(y.intercept=rep(1,N), datafile.final2$salinity.std, datafile.final2$Nit.std,datafile.final2$distance.between.beds.std,datafile.final2$dIICconnector_4000.std,datafile.final2$dIICflux_6000.std)

colnames(X) <- c("y.intercept","salinity.std","Nit.std", "distance.between.beds.std","dIICconnector_4000.std","dIICflux_6000.std")

stack3s <- inla.stack(
  data=list(y=datafile.final2$pres, link=1, Ntrials=1),
  A=list(A.s,  1),
  effects=list(w=1:spde$n.spde, 
               X = as.data.frame(X)), 
  tag="est")

formula3s <-
  y ~ -1 + y.intercept + salinity.std + Nit.std  + distance.between.beds.std  + dIICconnector_4000.std  + dIICflux_6000.std +f(w, model = spde)

spat_only_allyrs_sept27Questforsal2022_constrained_utils <- inla(
  formula3s, num.threads='6:-4',
  family='binomial', ### the zeros are spatially structured so do not need zero inflation (if try, the additional probability is small)
  Ntrials=Ntrials,
  data=inla.stack.data(stack3s),
  control.compute=ctrl.compute,
  control.predictor=list(
    A=inla.stack.A(stack3s), compute=T, link=1,
    quantiles=c(0.025, 0.5, 0.975)),
  control.inla=ctrl.inla,
  control.fixed=ctrl.fixed,
  verbose=T)

out.field <- inla.spde2.result(spat_only_allyrs_sept27Questforsal2022_constrained_utils,'w', spde, do.transf = TRUE)
range.out <- inla.emarginal(function(x) x,
                            out.field$marginals.range.nominal[[1]])

ss <- 1 # sample size to process (number of SLOO runs)

# define the radius of the spatial buffer surrounding the removed point.
# Make sure it isn't bigger than 25% of the study area (Le Rest et al.(2014))
rad <- range.out*0.15
alpha <- 0.05 # RMSE and MAE confidence intervals (1-alpha)
set.seed(199)

modform <- y ~ -1 + y.intercept + salinity.std + Nit.std  + distance.between.beds.std  + dIICconnector_4000.std  + dIICflux_6000.std + f(spatial.field, model=spde)

cv <- inlasloo(dataframe = datafile.final2,
               long = 'x', lat = 'y',
               y = 'pres', ss = ss,
               rad = rad,
               modform = modform,
               mesh = mesh, family = 'binomial',
               mae = F,ntrials = 1)

Thank you for your help!