BYM Model & Space – time interactions, got Hessian WARNING，and Different results on different computers?

271 views
Skip to first unread message

wqdon...@gmail.com

unread,
Sep 29, 2018, 10:03:08 AM9/29/18
to R-inla discussion group
Dear all,

I tried to run the R code in Chapter 7 of this book，the name is《Spatial and spatio-temporal Bayesian models with R-INLA》。

and I got Hessian WARNING，and Different results on different computers. I do not know how to deal with it.

I try to section 7.1.2 Space– time interactions，there are four models here,

I wanna say the model 1&4 seem is fine,but model2&3 seem to have some problems.

the book code is following:

## the book code
> require(INLA)
> inla.setOption(scale.model.default=FALSE)
> require(splancs)
> require(sp)
> require(fields)
> require(maptools)
> require(lattice)
> require(abind)

...

###################################################
### Code for Section 7.1.1
###################################################
# You neeed a folder called "LBW in Georgia" inside your working directory (my.dir)
# with the data downloaded from

#import Data
#and Data preprocessing

....

###################################################
### Code for Section 7.1.2 (run the code for Section 7.1.1 first) 
###################################################
#--- Type I interaction ---#
formula.intI<- y ~ + f(ID.area,model="bym", graph=Georgia.adj) +
    f(ID.year,model="rw2") + f(ID.year1,model="iid") + f(ID.area.year,model="iid")

mod.intI <- inla(formula.intI,family="poisson",data=data,E=E, 
                 control.predictor=list(compute=TRUE), 
                 control.compute=list(dic=TRUE,cpo=TRUE))

#--- Type II interaction ---#
ID.area.int <- data$ID.area
ID.year.int <- data$ID.year
formula.intII<- y ~ f(ID.area,model="bym",graph=Georgia.adj) +
  f(ID.year,model="rw2") + f(ID.year1,model="iid") +
  f(ID.area.int,model="iid", group=ID.year.int,control.group=list(model="rw2")) 

mod.intII <- inla(formula.intII,family="poisson",data=data,E=E, 
                  control.predictor=list(compute=TRUE), 
                  control.compute=list(dic=TRUE,cpo=TRUE))

#--- Type III interaction ---#
formula.intIII<- y ~  f(ID.area,model="bym",graph=Georgia.adj) +
  f(ID.year,model="rw2") + 
  f(ID.year1,model="iid") +
  f(ID.year.int,model="iid", group=ID.area.int,control.group=list(model="besag", graph=Georgia.adj))

mod.intIII <- inla(formula.intIII,family="poisson",data=data,E=E, 
                   control.predictor=list(compute=TRUE), 
                   control.compute=list(dic=TRUE,cpo=TRUE))

#--- Type IV interaction ---#
formula.intIV<- y ~ f(ID.area,model="bym",graph=Georgia.adj) +
  f(ID.year,model="rw2") +
  f(ID.year1,model="iid") + 
  f(ID.area.int,model="besag", graph=Georgia.adj,group=ID.year.int,control.group=list(model="rw2"))

mod.intIV <- inla(formula.intIV,family="poisson",data=data,E=E, 
                  control.predictor=list(compute=TRUE), 
                  control.compute=list(dic=TRUE,cpo=TRUE))

## the book code is end

and when I run model2, I got Hessian WARNING, like following:

> mod.intII <- inla(formula.intII,family="poisson",data=data,E=E,
+ control.predictor=list(compute=TRUE),
+ control.compute=list(dic=TRUE,cpo=TRUE))

*** WARNING *** Eigenvalue 0 of the Hessian is -27.7435 < 0
*** WARNING *** Set this eigenvalue to 0.897788
*** WARNING *** This have consequence for the accurancy of
*** WARNING *** the approximations; please check!!!
*** WARNING *** R-inla: Use option inla(..., control.inla = list(h = h.value), ...)
*** WARNING *** R-inla: to chose a different `h.value'.

I also check mod.intII$logfile,and found warning info like above. 

and fixed result of model2 and model3 is similar，mean is close to zero, and sd is very big.

> mod.intI$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03615165 0.01521281 0.006234724 0.03614904 0.06604964 0.03614552 4.183286e-12
> mod.intII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.001716422 4055.897 -7963.092 -0.11257 7956.45 0.001848365 9.058292e-16
> mod.intIII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.002882901 2556.503 -5019.437 -0.1489356 5015.27 0.03745706 3.448638e-09
> mod.intIV$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03317507 0.006451598 0.02043958 0.03318634 0.04583411 0.03321004 1.353453e-11
>

in adition, I restarted the R session and run all model again, I found model2 still tip Hessian WARNING.

but fixed result of model3 is different as following:

################################################################################################
# in centos 7

> mod.intI$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03616816 0.01520576 0.006263853 0.03616563 0.06605255 0.03616192 4.216472e-12
> mod.intII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -5.384162e-05 4056.925 -7965.113 -0.1144125 7958.465 0.0001806746 2.716519e-15
> mod.intIII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.01806566 0.1333061 -0.2436593 0.0180619 0.2795721 0.01806566 1.967613e-16
> mod.intIV$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03318039 0.006491015 0.0203669 0.03319124 0.04591975 0.03321439 1.099425e-11
>

#by the way, in this situation I get the same picture result as in the book. so I can't understand....

#cause I think the fixed mean of the four models should be similar....

############################################################################################
# in windows10
> mod.intI$summary.fixed mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) 0.03615165 0.01521281 0.006234724 0.03614904 0.06604964 0.03614552 4.183286e-12 > mod.intII$summary.fixed mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) 0.001716422 4055.897 -7963.092 -0.11257 7956.45 0.001848365 9.058292e-16 > mod.intIII$summary.fixed mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) -0.005845292 7922.847 -15585.55 -0.7158791 15569.51 0.06985049 1.003045e-08 > mod.intIV$summary.fixed mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) 0.03317507 0.006451598 0.02043958 0.03318634 0.04583411 0.03321004 1.353453e-11 >
#############################################################################################
Is anyone have suggestions on how to deal with Hessian WARNING and

how to troubleshoot why different computers produce different results?

Thanks very much!

Best
Dong


Helpdesk

unread,
Sep 29, 2018, 10:28:07 AM9/29/18
to wqdon...@gmail.com, R-inla discussion group, Marta A G Blangiardo, Michela Cameletti
Hi.

I cc to M&M for whom are the authors of this.

If you get this 'Hessian warning' then the results will differ. This
warning has to be fixed.

H
> --
> 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 post to this group, send email to
> r-inla-disc...@googlegroups.com.
> Visit this group at
> https://groups.google.com/group/r-inla-discussion-group.
> For more options, visit https://groups.google.com/d/optout.

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

wqdon...@gmail.com

unread,
Sep 30, 2018, 1:30:16 AM9/30/18
to R-inla discussion group
Hi,H.

I slightly changed the formulas of models 2 and 3 as following:

#--- Type II interaction ---#
formula.intII<- y ~ f(ID.area,model="bym",graph=Georgia.adj,constr=TRUE) +
f(ID.year,model="rw2",constr=TRUE) + f(ID.year1,model="iid",constr=TRUE) +
f(ID.area.int,model="iid", group=ID.year.int,control.group=list(model="rw2"),constr=TRUE)

#--- Type III interaction ---#
formula.intIII<- y ~ f(ID.area,model="bym",graph=Georgia.adj,constr=TRUE) +
f(ID.year,model="rw2",constr=TRUE) +
f(ID.year1,model="iid",constr=TRUE) +
f(ID.year.int,model="iid", group=ID.area.int,
control.group=list(model="besag", graph=Georgia.adj),constr=TRUE)

I added the constr=TRUE option.

The result is model2 do not tip Hessian WARNNING, and the fixed mean of the four models seem to similar, 

and The picture results of models 2 and 3 are no longer the same as in the book, and the pictures of models 1 and 4 without change are the same.

the summary() is following

## fixed result

> mod.intI$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03616816 0.01520576 0.006263853 0.03616563 0.06605255 0.03616192 4.216472e-12
> mod.intII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03114779 0.0108907 0.008937012 0.03261278 0.04512435 0.03323905 0.0001116881
> mod.intIII$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03666893 0.01498052 0.00720342 0.03666576 0.06611844 0.03666131 3.422427e-12
> mod.intIV$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.03318039 0.006491015 0.0203669 0.03319124 0.04591975 0.03321439 1.099425e-11


## summary

> summary(mod.intI)

Call:
c("inla(formula = formula.intI, family = \"poisson\", data = data, ", " E = E, control.compute = list(dic = TRUE, cpo = TRUE), control.predictor = list(compute = TRUE))" )

Time used:
Pre-processing Running inla Post-processing Total
0.7241 7.4754 0.4748 8.6744

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.0362 0.0152 0.0063 0.0362 0.0661 0.0362 0

Random effects:
Name   Model
ID.area BYM model
ID.year RW2 model
ID.year1 IID model
ID.area.year IID model

Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for ID.area (iid component) 28.20 4.715 18.76 28.40 36.89 29.43
Precision for ID.area (spatial component) 629.41 2155.627 42.27 210.05 3776.92 74.73
Precision for ID.year 30489.22 24265.992 5411.06 23992.72 94690.84 13999.68
Precision for ID.year1 32397.74 25974.011 5687.32 25417.98 101159.81 14741.97
Precision for ID.area.year 647.91 231.963 333.43 600.64 1227.24 519.18

Expected number of effective parameters(std dev): 297.12(34.37)
Number of equivalent replicates : 5.887

Deviance Information Criterion (DIC) ...: 11565.65
Effective number of parameters .........: 299.20

Marginal log-Likelihood: -5933.44
CPO and PIT are computed

Posterior marginals for linear predictor and fitted values computed

> summary(mod.intII)

Call:
c("inla(formula = formula.intII, family = \"poisson\", data = data, ", " E = E, control.compute = list(dic = TRUE, cpo = TRUE), control.predictor = list(compute = TRUE))" )

Time used:
Pre-processing Running inla Post-processing Total
8.2420 77.0691 2.6888 87.9999

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.0311 0.0109 0.0089 0.0326 0.0451 0.0332 1e-04

Random effects:
Name   Model
ID.area BYM model
ID.year RW2 model
ID.year1 IID model
ID.area.int IID model

Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for ID.area (iid component) 2195.94 1956.41 276.49 1652.93 7377.32 776.61
Precision for ID.area (spatial component) 2231.14 2126.26 217.27 1622.13 7911.49 619.02
Precision for ID.year 26639.74 20174.40 4376.01 21539.43 79574.48 12042.05
Precision for ID.year1 33977.39 27370.87 6453.85 26482.25 106367.99 15975.79
Precision for ID.area.int 5672.59 3335.60 2078.76 4771.95 14466.76 3564.23

Expected number of effective parameters(std dev): 375.29(16.85)
Number of equivalent replicates : 4.66

Deviance Information Criterion (DIC) ...: 11653.92
Effective number of parameters .........: 376.98

Marginal log-Likelihood: -6323.24
CPO and PIT are computed

Posterior marginals for linear predictor and fitted values computed

> summary(mod.intIII)

Call:
c("inla(formula = formula.intIII, family = \"poisson\", data = data, ", " E = E, control.compute = list(dic = TRUE, cpo = TRUE), control.predictor = list(compute = TRUE))" )

Time used:
Pre-processing Running inla Post-processing Total
3.8572 4777.6689 2.1177 4783.6438

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.0367 0.015 0.0072 0.0367 0.0661 0.0367 0

Random effects:
Name   Model
ID.area BYM model
ID.year RW2 model
ID.year1 IID model
ID.year.int IID model

Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for ID.area (iid component) 29.33 5.309 19.01 29.42 39.53 30.04
Precision for ID.area (spatial component) 409.76 1162.180 32.73 155.13 2357.75 59.05
Precision for ID.year 22394.57 21137.813 2589.28 16373.78 78684.34 7225.31
Precision for ID.year1 18606.70 18367.417 1262.98 13184.21 67168.32 3451.10
Precision for ID.year.int 2157.34 3373.200 255.22 1188.53 9998.78 549.77

Expected number of effective parameters(std dev): 191.76(34.08)
Number of equivalent replicates : 9.121

Deviance Information Criterion (DIC) ...: 11615.01
Effective number of parameters .........: 195.53

Marginal log-Likelihood: -7080.11
CPO and PIT are computed

Posterior marginals for linear predictor and fitted values computed

> summary(mod.intIV)

Call:
c("inla(formula = formula.intIV, family = \"poisson\", data = data, ", " E = E, control.compute = list(dic = TRUE, cpo = TRUE), control.predictor = list(compute = TRUE))" )

Time used:
Pre-processing Running inla Post-processing Total
0.9360 39.7960 0.4168 41.1488

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.0332 0.0065 0.0204 0.0332 0.0459 0.0332 0

Random effects:
Name   Model
ID.area BYM model
ID.year RW2 model
ID.year1 IID model
ID.area.int Besags ICAR model

Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for ID.area (iid component) 2190.24 2041.61 255.03 1611.95 7627.23 716.19
Precision for ID.area (spatial component) 2188.65 2132.51 198.35 1571.58 7887.98 558.49
Precision for ID.year 25394.39 19778.24 4122.66 20276.64 77161.77 11241.19
Precision for ID.year1 32840.10 26694.32 5934.72 25557.82 103378.92 15026.24
Precision for ID.area.int 5772.33 7141.77 997.37 3630.22 23501.62 1960.15

Expected number of effective parameters(std dev): 356.23(15.64)
Number of equivalent replicates : 4.91

Deviance Information Criterion (DIC) ...: 11664.36
Effective number of parameters .........: 358.46

Marginal log-Likelihood: -7354.32
CPO and PIT are computed

Posterior marginals for linear predictor and fitted values computed
####################################################################


by the way, you can get code & dataset in this url :




Best,
Dong






在 2018年9月29日星期六 UTC+8下午10:28:07，help help写道：
Reply all
Reply to author
Forward
0 new messages