inla_bym.R

#http://www.r-inla.org/examples/volume-1/code-for-bym-example
#
# Installer le package R-INLA
#
install.packages("INLA", repos="http://www.math.ntnu.no/inla/R/stable")

## 
## The downloaded binary packages are in
##  /var/folders/hk/d08tl8356d71jk4yj2dfsw4c0000gn/T//RtmplOfOqS/downloaded_packages

require(INLA)

## Loading required package: INLA

## Warning: package 'INLA' was built under R version 3.1.2

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.1.2

## Loading required package: Matrix
## Loading required package: splines

data(Germany)
g = system.file("demodata/germany.graph", package="INLA")
source(system.file("demodata/Bym-map.R", package="INLA"))
summary(Germany)

##      region            E                 Y                x        
##  Min.   :  1.0   Min.   :  1.211   Min.   :  0.00   Min.   : 3.00  
##  1st Qu.:136.8   1st Qu.:  5.164   1st Qu.:  4.00   1st Qu.:38.00  
##  Median :272.5   Median :  9.130   Median :  9.00   Median :51.00  
##  Mean   :272.5   Mean   : 13.388   Mean   : 13.39   Mean   :49.98  
##  3rd Qu.:408.2   3rd Qu.: 15.690   3rd Qu.: 16.25   3rd Qu.:61.00  
##  Max.   :544.0   Max.   :184.463   Max.   :215.00   Max.   :97.00

## just make a duplicated column
Germany = cbind(Germany,region.struct=Germany$region)

# standard BYM model (without covariates)
# pour connaitre tous les modèles dispos 
# names(inla.models()$latent)
#
formula1 <- Y ~ f(region.struct,model="besag",graph.file=g)+f(region,model="iid")

#
# names(inla.models()$likelihood)  pour toutes les family
#
result1<-inla(formula1,family="poisson",data=Germany,E=E)

## Warning in INLA::f(region.struct, model = "besag", graph.file = g):
## Argument 'graph.file' in 'f()' is obsolete; please use the more general
## argument 'graph' instead.

# with linear covariate
formula2 <-Y ~ f(region.struct,model="besag",graph.file=g) + f(region,model="iid") + x

result2 <-inla(formula2,family="poisson",data=Germany,E=E)



par(mfrow=c(2,1))
Bym.map(result1$summary.random$region.struct$mean)
Bym.map(result2$summary.random$region.struct$mean)