Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

inlasloo error: Intercept not found #64

Open
cperk opened this issue Jan 12, 2023 · 2 comments
Open

inlasloo error: Intercept not found #64

cperk opened this issue Jan 12, 2023 · 2 comments

Comments

@cperk
Copy link

cperk commented Jan 12, 2023

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
Copy link
Owner

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
Copy link
Author

cperk commented Jan 13, 2023

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants