NIMBLE provides the ability to specify and fit conditional autoregressive (CAR) models, both intrinsic and proper models.

Here we’ll use an intrinsic CAR (ICAR) model for disease mapping. We’ll illustrate with data on hospital admissions due to respiratory disease in 2010 from the 134 Intermediate Geographies (IG) north of the river Clyde in the Greater Glasgow and Clyde health board. The data are available in the CARBayesdata package and can be transformed into the neighborhood information needed by NIMBLE using functions from the spdep package. In particular, we need to provide vectors indicating which regions are neighbors of which other regions (adj), the weights for each pair of neighbors (weights), and the number of neighbors for each region (num).

As usual, the mean of the Poisson counts includes an offset for the expected count. We’ll also include a covariate, the percentage of people defined to be income deprived in each spatial region.

library(nimble, warn.conflicts = FALSE)
## nimble version 0.9.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.
library(CARBayesdata, quietly = TRUE)
library(sp, quietly = TRUE)
library(spdep, quietly = TRUE)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
data(GGHB.IG)
data(respiratorydata)

We handle the spatial analysis here with nb2WB from the package spdep.

respiratorydata_spatial <- merge(x = GGHB.IG, y = respiratorydata, by = "IG", all.x = FALSE)
W.nb <- poly2nb(respiratorydata_spatial, row.names =  rownames(respiratorydata_spatial@data))
## Determine neighborhood/adjacency information needed for neighborhood-based CAR model
nbInfo <- nb2WB(W.nb)

# A vector of indices indicating which regions are neighbors of which.
nbInfo$adj
##   [1]   2   3   5  87  91  93   1   3  91 113   1   2   5  11 113   6   7
##  [18]  10  98 101 105 106   1   3   9  11  19  87  90 104 109   4  10  12
##  [35] 105 108   4  10  20 106 109  11  15  16 113 114   5  19  20 109   4
##  [52]   6   7  12  13  14  20   3   5   8  16  19  22 113   6  10  13  17
##  [69]  26 108 111 112 122 124 126  10  12  14  17  10  13  17  20   8  16
##  [86]  18  24 113 114 116   8  11  15  18  22  12  13  14  20  25  26  15
## [103]  16  21  22  24   5   9  11  20  22  28   7   9  10  14  17  19  23
## [120]  25  28 109  18  22  24  11  16  18  19  21  24  27  28  20  25  15
## [137]  18  21  22  17  20  23  26  12  17  25 126  22  28  19  20  22  27
## [154]  30  33  37  38  29  31  34  38  30  34  35  41  44  39  40  29  36
## [171]  37  39  30  31  35  38  41  42  31  34  41  33  37  39  43  29  33
## [188]  36  38  43  29  30  34  37  42  43  32  33  36  40  43  46  47  32
## [205]  39  46  55  31  34  35  42  44  45  34  38  41  43  45  52  36  37
## [222]  38  39  42  47  52  54  31  41  45  50  51  58  41  42  44  50  52
## [239]  59  39  40  47  48  55  39  43  46  48  54  56  46  47  55  56  53
## [256]  57  61  44  45  58  59  70  44  58  70  42  43  45  54  59  63  49
## [273]  55  61  65  43  47  52  56  63  40  46  48  53  56  62  65  66  72
## [290]  47  48  54  55  62  63  74  49  60  61  44  50  51  70  45  50  52
## [307]  63  70  91 113  57  61  64  67  68  73  49  53  57  60  65  68  71
## [324]  55  56  63  66  74  76  83  52  54  56  59  62  74  91 113  60  67
## [341]  69  75  78  53  55  61  71  72  81  55  62  72  76  60  64  73  75
## [358]  60  61  71  73  77  79  64  78  80  82  50  51  58  59 113  61  65
## [375]  68  79  81  55  65  66  76  81  88  90  60  67  68  75  77  85  56
## [392]  62  63  83  91  93  64  67  73  78  85  62  66  72  83  87  90  68
## [409]  73  79  85  86  64  69  75  80  82  85  68  71  77  81  84  86  69
## [426]  78  82  92  95  65  71  72  79  84  88  69  78  80  85  95  62  74
## [443]  76  87  93  79  81  86  88  89  73  75  77  78  82  86  94  95  98
## [460] 101  77  79  84  85  89  96  98   1   5  76  83  90  93  72  81  84
## [477]  89  90  96  97 100 104  84  86  88  96   5  72  76  87  88 104   1
## [494]   2  59  63  74  93 113  80  95  99 103 117   1  74  83  87  91  85
## [511]  95 101 102 105  80  82  85  92  94  99 102  86  88  89  98 100 106
## [528]  88 100 104 109   4  85  86  96 101 106  92  95 102 103 107  88  96
## [545]  97 106 109   4  85  94  98 102 105  94  95  99 101 105 107  92  99
## [562] 107 117 119   5  88  90  97 109   4   6  94 101 102 107 108   4   7
## [579]  96  98 100 109  99 102 103 105 108 110 119   6  12 105 107 110 111
## [596] 112   5   7   9  20  97 100 104 106 107 108 112 119 122  12 108 112
## [613]  12 108 110 111 122 124   2   3   8  11  15  59  63  70  91 114   8
## [630]  15 113 115 116 114 116  15 114 115  92 103 118 119 117 119 120 121
## [647] 123 103 107 110 117 118 120 122 118 119 121 122 118 120 122 123 126
## [664]  12 110 112 119 120 121 124 126 118 121 126 127  12 112 122 126 127
## [681] 129 130  12  26 121 122 123 124 127 132 123 125 126 129 132 131 125
## [698] 127 130 132 125 129 131 132 128 130 132 133 126 127 129 130 131 133
## [715] 134 131 132 134 132 133
# A vector of weights. In this case, all weights are 1.
head(nbInfo$weights)
## [1] 1 1 1 1 1 1
# A vector of length N. num[n] indicates how many neighbors region n contains.
# This helps map the adj vector to the starting region.
nbInfo$num
##   [1]  6  4  5  7  9  5  5  5  4  7  7 11  4  4  7  5  6  5  6 10  3  8  2
##  [24]  4  4  4  2  4  4  4  5  2  4  6  3  4  5  6  7  4  6  6  8  6  6  5
##  [47]  6  4  3  5  3  6  4  5  9  7  3  4  7  6  7  7  8  5  6  4  4  6  4
##  [70]  5  5  7  6  6  5  6  5  6  6  5  6  5  5  5 10  7  6  9  4  6  7  5
##  [93]  5  5  7  6  4  6  5  5  6  6  5  5  7  6  7  7  8  5  3  6 10  5  2
## [116]  3  4  5  7  4  5  8  4  4  3  8  5  1  4  4  4  7  3  2

Now we have the three pieces of information we need. We’re ready to use the dcar_normal distribution in a nimbleModel.

nregions <- nrow(respiratorydata_spatial)

code <- nimbleCode({
    # priors
    beta ~ dnorm(0, sd = 100)
    sigma ~ dunif(0, 100)   # prior for variance components based on Gelman (2006)
    tau <- 1 / sigma^2
    # latent process
    s[1:N] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N], tau, zero_mean = 0)
    # likelihood
    for(i in 1:N) {
        log(lambda[i]) <- log(expected[i]) + beta*x[i] + s[i]
        y[i] ~ dpois(lambda[i])
    }
})

x <- respiratorydata_spatial$incomedep
x <- x - mean(x)  # center for improved MCMC performance

set.seed(1)

constants <- list(N = nregions, L = length(nbInfo$adj), 
               adj = nbInfo$adj, weights = nbInfo$weights, num = nbInfo$num,
                        x = x, expected = respiratorydata_spatial$expected)
data <- list(y = respiratorydata_spatial$observed)
inits <- list(beta = 0, sigma = 1, s = rnorm(nregions))
      
model <- nimbleModel(code, constants = constants, data = data, inits = inits)
## defining model...
## Adding adj, weights, num as data for building model.
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.
cModel <- compileNimble(model)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
conf <- configureMCMC(model, monitors = c('beta', 'sigma', 's'))
conf$printSamplers()
## [1] RW sampler: beta
## [2] RW sampler: sigma
## [3] CAR_normal sampler: s[1:134]
MCMC <- buildMCMC(conf)
cMCMC <- compileNimble(MCMC, project = cModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
samples <- runMCMC(cMCMC, niter = 10000, nburnin = 1000)
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Here are density plots for some posterior samples of interest.

plot(density(samples[,"sigma"]))

plot(density(samples[,"s[1]"]))

plot(density(samples[,"beta"]))