R-inla output

[Subramanian Swaminathan]

Dear all,

I am new to R-inla. I am running spatio-temporal model with interactions for disease mapping. The model gives and output with a (1) DIC value and also (2) "Deviance information criterian (DIC, saturate): NAN".  I am confused, why I am getting "NaN' (2) , despite geting a DIC value in (1). For your kind information the following are the codes I used for fitting a model covariates with spatio-temporal interactions (all 4 types). I would be grateful, if any suggesions from the group.

Thanks advance

ssubra 

---

rm(list=ls())

library(lme4)

library(INLA)

library(reshape2)

library(tidyverse)

require(splancs)

require(sp)

require(fields)

require(maptools)

require(lattice)

require(abind)

library(rgdal)

library(spdep)

library(dplyr)

library(RColorBrewer)

library(ggplot2)

library(ggregplot)

# Priors

prec_period <- list(c(1, 0.00040145), c(1, 0.00160580), c(1,0.04014488), c(1, 0.00005))

prec_space <- list(c(1, 0.00048684), c(1, 0.00194738), c(1, 0.04868439), c(1, 0.00005))

prec_interact <- list(c(1, 0.00078138), c(1, 0.00312552), c(1, 0.07813805), c(1, 0.00005))

prec_interact2 <- list(c(1,0.00040145), c(1, 0.00160580), c(1, 0.04014488), c(1, 0.00005))

##############import shapefile#########################

#shp <- readOGR(dsn = "H:/blockmapVl", layer = "blockmapVL") #491 BLOCKS

shp <- readOGR(dsn = "D:/VL model/502 blocks/shp", layer = "merged_bihar") #502 BLOCKS

# define adjacency matrix for spatial model

nb <- poly2nb(shp)

# translate to INLA format and save

nb2INLA("D:/VL model/502 blocks/map/map.adj", nb)

# read in the INLA graph we just created

g <- inla.read.graph(filename = "D:/VL model/502 blocks/map/map.adj")

data<-read.csv("D:/VL model/502 blocks/data/EM_DATA 2013 to 2016_training_3E.csv")

###### Spatio-temporal Models WITH interactions ########

#----- Model 7 Type I interaction ----#

formula.intI1<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym", graph=g,hyper = prec_space) +f(period_id,model="rw1") + 

  #f(period_id2,model="iid",hyper = prec_period,initial = 1) 

  +f(block.period,model="iid",hyper = prec_interact)

mod.intI1 <- inla(formula.intI1,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                  control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                  control.predictor=list(compute=TRUE), verbose=T,

                  control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intI1)

modelintI1<-capture.output(summary(mod.intI1))

write.table(modelintI1,file="D:/VL model/502 blocks/sum_2/modelintI1.txt",sep=",",row.names = F,quote = F)

fit=mod.intI1$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintI1.csv")

formula.intI2<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym", graph=g,hyper = prec_space) +f(period_id,model="rw2",hyper = prec_period) + 

  #f(period_id2,model="iid",hyper = prec_period,initial = 1) 

  +f(block.period,model="iid",hyper = prec_interact)

mod.intI2 <- inla(formula.intI2,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                  control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                  control.predictor=list(compute=TRUE), verbose=T,

                  control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intI2)

modelintI2<-capture.output(summary(mod.intI2))

write.table(modelintI2,file="D:/VL model/502 blocks/sum_2/modelintI2.txt",sep=",",row.names = F,quote = F)

fit=mod.intI2$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintI2.csv")

#--- Model 8 Type II interaction ---#

formula.intII1<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw1",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_period) +f(block_id2,model="iid", group=period_id3,control.group=list(model="rw1"),hyper = prec_interact) 

mod.intII1 <- inla(formula.intII1,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                   control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                   control.predictor=list(compute=TRUE), verbose=T,

                   control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intII1)

modelintII1<-capture.output(summary(mod.intII1))

write.table(modelintII1,file="D:/VL model/502 blocks/sum_2/modelintII1.txt",sep=",",row.names = F,quote = F)

fit=mod.intII1$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintII1.csv")

formula.intII2<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw2",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_period) +f(block_id2,model="iid", group=period_id3,control.group=list(model="rw2"),hyper = prec_interact) 

mod.intII2 <- inla(formula.intII2,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                   control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                   control.predictor=list(compute=TRUE), verbose=T,

                   control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intII2)

modelintII2<-capture.output(summary(mod.intII2))

write.table(modelintII2,file="D:/VL model/502 blocks/sum_2/modelintII2.txt",sep=",",row.names = F,quote = F)

fit=mod.intII2$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintII2.csv")

#--- Model 9 Type III interaction ---#

formula.intIII1<- count ~  1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw1",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_period) +f(block_id2,model="iid",   group=period_id3,control.group=list(model="besag", graph=g),hyper = prec_uspace)

mod.intIII1 <- inla(formula.intIII1,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                    control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                    control.predictor=list(compute=TRUE), verbose=T,

                    control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intIII1)

modelintIII1<-capture.output(summary(mod.intIII1))

write.table(modelintIII1,file="D:/VL model/502 blocks/sum_2/modelintIII1.txt",sep=",",row.names = F,quote = F)

fit=mod.intIII1$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintIII1.csv")

formula.intIII2<- count ~  1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw2",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_period) +f(block_id2,model="iid",     group=period_id3,control.group=list(model="besag", graph=g),hyper = prec_uspace)

mod.intIII2 <- inla(formula.intIII2,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                    control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                    control.predictor=list(compute=TRUE), verbose=T,

                    control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

modelintIII2<-capture.output(summary(mod.intIII2))

write.table(modelintIII2,file="D:/VL model/502 blocks/sum_2/modelintIII2.txt",sep=",",row.names = F,quote = F)

fit=mod.intIII2$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintIII2.csv")

#--- Model 10 Type IV interaction ---#

formula.intIV1<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw1",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_interact2) + f(block_id2,model="besag", graph=shp.adj,group=period_id3,control.group=list(model="rw2"),hyper = prec_uspace)

mod.intIV1 <- inla(formula.intIV1,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                   control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                   control.predictor=list(compute=TRUE), verbose=T,

                   control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intIV1)

modelintIV1<-capture.output(summary(mod.intIV1))

write.table(modelintIV1,file="D:/VL model/502 blocks/sum_2/modelintIV1.txt",sep=",",row.names = F,quote = F)

fit=mod.intIV1$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintIV1.csv")

formula.intIV2<- count ~ 1+Maximum_temp+EVI+Minimum_temp+BIO.1+BIO.3+BIO.12+popdens+Soil_moisture +

  f(block_id,model="bym",graph=g,hyper = prec_space) +f(period_id,model="rw2",hyper = prec_period) + 

  f(period_id2,model="iid",hyper = prec_interact2) + f(block_id2,model="besag", graph=shp.adj,group=period_id3,control.group=list(model="rw2"),hyper = prec_uspace)

mod.intIV2 <- inla(formula.intIV2,family="nbinomial",data=data,offset=log(E2),# changed from 'pop' to 'E'

                   control.inla = list(tolerance = 1e-9, diagonal=10, strategy="laplace"),

                   control.predictor=list(compute=TRUE), verbose=T,

                   control.compute=list(waic=TRUE,dic=TRUE,cpo=TRUE))

summary(mod.intIV2)

modelintIV2<-capture.output(summary(mod.intIV1))

write.table(modelintIV2,file="D:/VL model/502 blocks/sum_2/modelintIV2.txt",sep=",",row.names = F,quote = F)

fit=mod.intIV2$summary.fitted.values

write.csv(fit,"D:/VL model/502 blocks/out 2/modelintIV2.csv")

[Helpdesk]

that could happens with large uncertainty, as it integrates over the
posterior for the linear predictors.

make sure you're using a recent testing version and R-4.0, as there were
some changes in the code to improve this sometime last year.

If it still fail, you can look for the local.dic output to pinpoint
which data point(s) that cause the isse, as it can all be due to one
observation

I need to run it here to see what happens. If so, please provide script
and data so I can rerun. Also, using relative paths to filename, so I
can run it without making tons of changes...

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/381cf805-93af-45fa-b93e-d4923d1ca218n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Subramanian Swaminathan]

Many thanks for your quick response. Please find attached the script and the required data files for your kind perusal. I am running the model in RStudio ver. 1.3.1056 and R version 4.0.2 (2020-06-22) under windows ver 10. I am sorry that I do not know how to set relative paths to file name. Could you kindly let me know how this can be done, so that I can prepare the script with relative paths and send it again.

I am unable send with attachment, "getting error posting message"

Please help me, how shold I send the data files

Thanks in advance

ssubra

[Subramanian Swaminathan]

Using INLA_20.11.16-3 built under version 4.0.3.

ssubra

[Subramanian Swaminathan]

Hi Rue,

I am unable to send the data files as attachement. Could you suggest an alternative..

ssubra

[Subramanian Swaminathan]

Hi Rue,

I have sent the script and data files via he...@r-inla.org..

Please help me out to resolve the above-mentioned issue.

Thanks

ssubra

[Subramanian Swaminathan]

Hi Rue,

Kindly let me know whether you were able to access the data files and are intact.

Thanks and regards

S. Subramanian 

[Helpdesk]

I think one issue is here

f(block_id2,model="iid",
group=period_id3,control.group=list(model="besag", graph=g),hyper =
prec_uspace)

its idd kronecker besag. with this, the sum to zero constraints are not
added, then you need to write this as besag kronecker iid (which is the
same as besag replicated

On Fri, 2021-05-07 at 09:58 -0700, Subramanian Swaminathan wrote:
> Hi Rue,
>
> Kindly let me know whether you were able to access the data files and
> are intact.
>
> Thanks and regards
> S. Subramanian 
>
> On Friday, May 7, 2021 at 1:48:14 AM UTC+5:30 Subramanian Swaminathan
> wrote:
> > Hi Rue,
> >
> > I have sent the script and data files via he...@r-inla.org..
> >
> > Please help me out to resolve the above-mentioned issue.
> >
> > Thanks
> > ssubra
> >
> > On Thursday, May 6, 2021 at 11:59:24 PM UTC+5:30 Subramanian
> > Swaminathan wrote:
> > >
> > > Hi Rue,
> > >
> > > I am unable to send the data files as attachement. Could you
> > > suggest an alternative..
> > >
> > > ssubra
> > >

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit

> https://groups.google.com/d/msgid/r-inla-discussion-group/809ef9bd-c426-4b25-9bea-a8a48b2994f7n%40googlegroups.com

[Subramanian Swaminathan]

Many thanks for the clarifications via helpdesk. I tried running the type I and II interactions models with your modified script, It stopped with the following error: Kindly help me how should I takle them:

Error in inla.cygwin.run.command(paste("cygpath -u ", inla.cygwin.protect.special.chars(filename)),  :

 Cannot find the CYGWIN installation: C:/cygwin

Tried installing CYGWIN, which resulted in the following error:

> install.packages("CYGWIN")
Installing package into ‘C:/Users/Subra/Documents/R/win-library/4.0’
(as ‘lib’ is unspecified)
Warning in install.packages :
  package ‘CYGWIN’ is not available (for R version 4.0.2)

Also, could you update me on the status of scripting type III and IV interaction models by the PhD student.

Thanks, I look forward to hearing from you

 

ssubra

[Helpdesk]

you have to remove the

inla.call="remote"

argument, as it was for me only

[Subramanian Swaminathan]

Thanks for your quick response. Now it works.

ssubra

[Subramanian Swaminathan]

Hi Rue,

Using your modified script I could run the type I and type II interactions models along with covariates. This is just a gentle reminder about the script for type III and IV interactions.  I would appreciate if you could kindly share the script for the type III and type IV interactions, if it has already been done by your PhD student.

Thanks

ssubra

[Helpdesk]

this is work in progress...

On Wed, 2021-05-12 at 08:51 -0700, Subramanian Swaminathan wrote:
> Hi Rue,
>

[Subramanian Swaminathan]

Thanks a lot, will look forward to hear from you

Regards

ssubra

[Subramanian Swaminathan]

Hi Rue,

This is just a gentle remainder related to modifying the script for type III and type IV interactions. I understand that the script for these models are in progress with one of your PhD student. I would be grateful if you could kindly share the modified script if the work is already completed.

Thanks in advance

ssubra 

[Helpdesk]

it will take time still

[Subramanian Swaminathan]

Hi Rue,

I am just refering to my earlier mail i.r. of modifying script for running type II-IV interaction models for my high dimension data, which I shared with you. May I request you to kindly let me know the progress in this respect by your PhD student.

Thanks,

ssubra

[Subramanian Swaminathan]

Hi Rue,

This is just gentle remainder to know whether your team could do the r-script running type I-IV interaction models for hhigh dimensional data which shared with you long back. I am eager to receive the r-script.

Thanks and regards

ssubra