I have used R-INLA in scientific work to be published; how should i cite it?It would be great if the report describing the methodology could be cited, which isH. Rue, S. Martino, and N. Chopin. Approximate Bayesian inference for latent Gaussian models using inte- grated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B, 71(2):319{392, 2009. and mention also the web-cite for where the R-INLA package is located, www.r-inla.org, I want to use a prior for my hyperparameter that seems not to be implemented, what should I do?If you want to use a prior for the hyperparameter that is not yet implemented there are two choices. If you think that your prior should be on the list and that other might use it to, please send us the expressions and we will add it. Alternatively, you can define your own prior using prior = "expression: ...."; see this help-page.Does INLA support the use of different link-functions?Yes, the type of link function is given in the 'control.data' statement using `link="<link>"', and the type of link-functions implemented are listed on the documentation for each likelihood. The default option `link="default"' is a macro for the second entry on the list. Here is an example using the binomial.n = 100 If you miss other link-functions, or want to enable more of them for likelihoods where they are not yet enabled; please let us know.z = runif(n) eta = 1 + 0.1*z N = 20 ## logit p = exp(eta)/(1+exp(eta)) y = rbinom(n, size = N, prob = p) r = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials = rep(N, n), control.data = list(link = "logit"), control.predictor = list(compute=TRUE)) ## probit p = pnorm(eta) y = rbinom(n, size = N, prob = p) rr = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials = rep(N, n), control.data = list(link = "probit"), control.predictor = list(compute=TRUE)) ## cloglog p = exp(-exp(-eta)) y = rbinom(n, size = N, prob = p) rrr = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials = rep(N, n), control.data = list(link = "cloglog"), control.predictor = list(compute=TRUE)) How can I do predictions using INLA?In INLA there is no function ``predict'' as for glm/lm in R. Predictions must to done as a part of the model fitting itself. As prediction is the same as fitting a model with some missing data, we can simply set y[i] = NA for those ``locations'' we want to predict. Here is a simple examplen = 100 n.pred = 10 y = arima.sim(n=n, model=list(ar=0.9)) N = n + n.pred yy = c(y, rep(NA, n.pred)) i = 1:N formula = yy ~ f(i, model="ar1") r = inla(formula, data = data.frame(i,yy), control.data = list(initial = 10, fixed=TRUE)) ## no observational noise which gives predictions > r$summary.random$i[(n+1):N, c("mean", "sd") ] mean sd 101 -2.6021569048 1.397605790 102 -2.2801830807 1.639556612 103 -2.0008859575 1.807935005 104 -1.7635334748 1.932743951 105 -1.5569410367 2.025845322 106 -1.3790793264 2.097800368 107 -1.2261449461 2.153969553 108 -1.0913806143 2.197121102 109 -0.9767910767 2.232332803 110 -0.8749748925 2.259614345 Quantiles such like r$summary.fitted.values and r$marginals.fitted.values, if computed, use the identity link if y[i] = NA. If you need to transform marginals, you can use the function inla.marginal.transform(), as in the following example (see demo(Tokyo)) ## Load the data data(Tokyo) summary(Tokyo) Tokyo$y[300:366] <- NA ## Define the model formula = y ~ f(time, model="rw2", cyclic=TRUE, param=c(1,0.0001)) - 1 ## The call to inla result = inla(formula, family="binomial", Ntrials=n, data=Tokyo, control.predictor=list(compute=T)) ## need to recompute the fitted values for those with data[i] = NA, ## as the identity link is used. n = 366 fitted.values.mean = numeric(n) for(i in 1:366) { if (is.na(Tokyo$y[i])) { if (FALSE) { ## either like this, which is slower... marg = inla.marginal.transform( function(x) exp(x)/(1+exp(x)), result$marginals.fitted.values[[i]] ) fitted.values.mean[i] = inla.expectation( function(x) x, marg) } else { ## or like this, which is much faster... fitted.values.mean[i] = inla.expectation( function(x) exp(x)/(1 +exp(x)), result$marginals.fitted.values[[i]]) } } else { fitted.values.mean[i] = result$summary.fitted.values[i,"mean"] } } plot(fitted.values.mean) When I run it on a Mac, I get an error "dyld: unknown required load command 0x80000022" ?The binary for the inla-program is built using version 10.6 of MacOSX, and you are likely using version 10.5.x of MacOSX. An upgrade of MacOSX to 'Snow Leopard'/10.6 will resolve this issue. Otherwise, you will have to compile the inla-program yourself. There are scripts for doing this; see inla.googlecode.comSome of the models needs a graph, how do I specify it?Some of the models in INLA needs the user to specify a graph in terms of a graph.file. The graph.file is a filename for the graph and the format of the graph.file is as follows. The first entry is the number of nodes in the graph, `n'. The nodes in the graph are labelled 1, 2, ..., n. The next entries, specify the number of neighbours and the neighbours for each node, like `node number.of.neighbours neighbour.1 neighbour.2 ...'. An example will clarify the format.41 2 3 4 2 03 1 14 1 1This defines a graph with four nodes, where node 1 has 2 neighbours 3 and 4. Node 2 as 0 neigbours, node 3 has 1 neighbour 1, and node 4 has 1 neighbour 1, like this  Note that we need to specify the implied symmetry as well. In this example 4 is a neighbour of 1, then we need to specify that 1 is a neighbour of 4 as well. Some extensions:
# This is a simple graph# The first node 1 2 3 4 2
0 3 1 1 4 1 1which specify the same graph as above. You will find more examples of graph-specifications in the examples. Hints: You can read a such graph into R using inla.read.graph("graph.dat"), and if you think you have errors in your specification, you can try inla.debug.graph("graph.dat") which help you spot the error. If you have a correctly specified graph, but its numbering starts from 0 instead of 1, you can covert it using g = inla.read.graph("graph.dat"); inla.write.graph(g, "new.graph.dat"); How does inla() deal with "NA"s in the data.frame?For a formula likeformula = y ~ x + f(k, model= <some model>) then NA's in either y, x or k are treated differently. NA's in the response yIf y[i] = NA, this means that y[i] is not observed, hence gives no contribution to the likelihood.NA's in fixed effect xIf x[i] = NA this means that x[i] is not part of the linear predictor for y[i]. For fixed effects, this is equivalent to x[i]=0, hence internally we make this change: x[is.na(x)] = 0NA's in random effect kIf k[i] = NA, this means that the random effect does not contribute to the linear predictor for y[i].NA's in a factor xHere, we assume that 'x' is a factor; is.factor(x) == TRUE. NA's in a factor is not allowed, and must be avoided. If an entry in the data.frame or list(), is a factor, then this will be detected and and error produced (unless NA is a level itself, of course). ## INLA detect the NA in the factor
> x = as.factor(c(1,2,NA)) > y = c(1,2,3) > r = inla(y~x, data = data.frame(x,y)) Error in function (x, ...) : Error in the dataframe. INLA does not support factor(s) 'x' having NA's. Please rewrite the factor(s) into fixed effects. Unfortunately, the internal check fail to detect the following case, hence this formulation should be avoided!!! ## Here, INLA does not detect the NA in the factor. Please avoid a such formulation! > x = c(1,2,NA) > y = c(1,2,3) > r = inla(y~ as.factor(x), data = data.frame(x,y)) Warning message: In `[<-.factor`(`*tmp*`, thisvar, value = 0) : invalid factor level, NAs generated Can INLA deal with missing covariates?No, INLA has no way to `impute' or integrate-out missing covariates. You have to adjust your model to account for missing covariates. Sometimes, you can formulate a joint model for the data and the covariates, but this is case-specific.How can I compute cross-validation or predictive measures of fit?inla() provides two types of ``leave-one-out'' predictive measures of fit. It is the CPO value, which is Prob(y_i | y_{-i}), and the PIT value Prob(y_i^new \le y_i | y_{-i}). To enable the computation of these quantities, you will need to add the argumentcontrol.compute=list(cpo=TRUE) to `inla()'. Since the CPO/PIT quantities are computed keeping the integration-points fixed, you may want to use an improved version of the grid integration, for example you may additionally add the argument control.inla=list(int.strategy = "grid", diff.logdens = 4) In some examples, you may also want to increase the accuracy of the tails of the marginals, as the tail behaviour is sometimes important for the CPO/PIT calculations (see the variable `result$cpo$failure' below), like If the resulting object is `result', then you will find the predictive quantities as `result$cpo$cpo' and `result$cpo$pit'. control.inla = list(strategy = "laplace", npoints = 21) where add more evaluation points (npoints=21) instead of the default npoints=9. Another issue, are the implicit assumptions made in the inla-program. There are internal checks in the inla-program about these assumptions, and these check will appear as `result$cpo$failure'. In short, if `result$cpo$failure[i] > 0' then some assumption is violated, the higher the value (maximum 1) the more seriously. If `result$cpo$failure[i] == 0' then the assumptions should be ok. It might be a good idea to recompute those, if any, CPO/PIT values for which `result$cpo$failure>0'. However, this must be done `manually' by removing `y[i]' from the dataset, fit the model and try to predict `y[i]'. To provide a more efficient implementation of this, we have provided a special version of inla() inla.cpo( .... ) which works like `inla()', but will recompute (in an efficient way) the cpo/pit for which failure > 0, afterwards. If you want to do this yourself, see the source for inla.cpo() for details. WARNING: inla.cpo() is an experimental feature and will likely change in the future. If you have opinions on this issue, please let us know. Even though inla() is fast, is it possible to make it run faster?The R-INLA library contains precompiled binaries of the inla-program. To make the program run for all various processors, its is compiled using generic compiler options and 32bit, for Windows, MacOSX and Linux. There is certainly possible to compile the inla-program particulary for your computer/cpu. On my laptop, I get a speedup of about 25% easily. There are 'build' files available here, which require you to install mercurial which is availble in the package-system for all major Linux distros and for MacOSX prebuilt binaries are availble. You may want to play around with ATLAS/GOTO as well, and ATLAS is bundled in the build files. For MacOSX you'll need a newer compiler than the one bundled in XCode, for example using the ones here.To use a native compiled inla-program, download the appropriate version and install it as <path>/inla, for example /usr/local/bin/inla. The easiest way to override the builtin inla-program in R-INLA, is to do library(INLA) inla.setOption("inla.call", "/usr/local/bin/inla") I have access to a remote Linux server, is it possible to run the computations remotely and running R locally?Yes! This option allow INLA to use a remote server to do the computations. In order to use this feature, you need to do some setup which is different from (various) Linux distributions, Mac and Windows. In short the
inla.setOption("inla.call", "remote") or per inla-call inla(...., inla.call="remote") HOWTO setup ssh-keys: For the unexperienced user, this is somewhat tricky; sorry about that. Any Google search on "How to setup ssh-keys" willl reveal that more than one found this a bit tricky. Well, on some Linux distributions there are now tools which do this automatically (in Ubuntu at least), but a manual setup is not hard as ssh-agent is used by default when you log into your Desktop in Linux (Make sure you have installed `ssh' packages; do a search and install all of them if you're not sure what is the least requriement.) For MacOSX, the easiness depends on your version of the OS, but for MacOSX version prior to 10.6 this involves installing the program SSHKeychain (but make sure to consult this page to get some hints how to activate the ssh-agent at login for older MacOSX's , or alternatively, you may also use the same approach as for Windows using a fixed socket). For the newer MaxOSX version 10.6, this feature is builtin. Here is how to setup everything manually. If you have already generated your key or use a separate key for various machines, please make appropriate changes. I assume the remote Linux/MacOSX host, which shall do the computations, is called `mc.math.ntnu.no' and the remote username is `hrue'. You have to change these accordingly. LINUX/MacOSX: First time only: (all the following commands you'll do in the terminal window)
First time only:
/cygdrive/c/PROGRA~1/R/R-29~1.1/library/INLA/bin/remote/ssh-copy-id hrue@inla.math.ntnu.no Every time you start your desktop:
replacing `hrue' with your username. Then run the command `ssh-add' and give your password. You R-session should now be able to connect to the remote server.
I want to compare the results from inla() with those obtained with MCMC; is there any way to do this using inla()?Even though inla() use the integrated nested Laplace approximations to compute the marginals, it has also a builtin MCMC-sampler. The sampler is the ``one-block MCMC sampler'' described in chapter 4 in the book by Rue and Held (2005). There is only one type of sampler, and if the model is either to large or to far away from a conditional Gaussian one, the accept-rate can be to low to get some reliable results out. In this case you may want to try to reduce the scaling of the proposal (option -S <scale>) but this does not always help and you have to find other ways to sample from your model.In the demo-directory of INLA, there is an example how to compare the inla-results for the Tokyo-example using the builtin sampler, and you can run it using demo("Tokyo-compare") and the R-code is found as the output of the following command system.file("demo/Tokyo-compare.R", package="INLA")
When running the example, you'll enter a menu displaying Chose directory to view These are the (set of) nodes that is available, and entering for example `1' displays the density estimates from MCMC and inla(), and `-1' display the density estimate from inla() together with the histogram of the MCMC trace. [1] hyperparameter-random.effect00000001-parameter [2] hyperparameter-random.effect00000001-parameter-user-scale [3] random.effect00000001 (<0: histogram, >0: density: q:quit) : I want to fit survival models using INLA, how do I do this? Easiest way is, first to check the given examples. OR follow these steps:
I have a complicated ``state-space model'', can this be fitted using INLA?Good question! There is no `state-space' model in the list of latent models, it does not mean that you cannot define your own latent model in a state-space form. It will be more complicated though, at least until some of these model are properly implemented and included in INLA. It would be sufficient to give an example only, then you'll figure out (if you're able to do it)...## This is an example how to implement the following state-space model ## without using the builtin model for {x_t} ## ## Observational equation: ## y_t = x_t + eps_t, t=1..n ## ## State equation: ## x_t = 2*x_{t-1} - x_{t-2} + v_t, t=3...n ## ## The 'trick' is to rewrite this as ## y_t = x_t + eps_t, t=1..n ## 0 = x_t - 2*x_{t-1} + x_{t-2} + v_t, t=3...n ## and then we can implement this using an augmented model two ## different likelihoods. The first n datapoints are Gaussian with ## unknown precision, whereas the last `m' datapoints, which are faked ## to be 0, is observed with a high and fixed precision. Note that the ## model for x_t, we just use model=iid with fixed low precision. n = 100 m = n-2 y = sin((1:n)*0.2) + rnorm(n, sd=0.1) formula = Y ~ f(i, model="iid", initial=-10, fixed=T) + f(j, w, copy="i") + f(k, copy="i") + f(l, model ="iid") -1 Y = matrix(NA, n+m, 2) Y[1:n, 1] = y Y[1:m + n, 2] = 0 i = c(1:n, 3:n) # x_t j = c(rep(NA,n), 3:n -1) # x_t-1 w = c(rep(NA,n), rep(-2,m)) # weights for j k = c(rep(NA,n), 3:n -2) # x_t-2 l = c(rep(NA,n), 1:m) # v_t r = inla(formula, data = data.frame(i,j,w,k,l,Y), family = rep("gaussian", 2), control.data = list(list(), list(initial=10, fixed=T))) I have multinomial data, but I cannot find the multinomial likelihood; isn't it available?Multinomial data implies that one datapoint depends on the latent field through a non-linear combination of some of its elements, and this situation is not yet supported in INLA. The classical way around it, can however, be used. It is possible to rewrite multinomial data into several Poisson data, using a data-indexed intercept. A classical reference is found here, and a Bayesian discussion is found in section 9.7 of the WinBugs manual. Here is the Alligator data example of multinomial analysis taken form the WinBugs manual vol I.I have some linear combinations of the nodes in the latent field that I want to compute the posterior marginal of, is that possible?Yes! These are called 'linear combinations'. There are handy functions, 'inla.make.lincomb()' and 'inla.make.lincombs()', to define one or many such linear combinations. Single linear combinations made by using 'inla.make.lincomb()' can easily be joined into many. Its use is easiest explained using a rather long example...Control variables for linear combinations, are defined in ?control.lincomb. Here is the example, that explains these features. ## A simple model n = 100 a = rnorm(n) b = rnorm(n) idx = 1:n y = rnorm(n) + a + b formula = y ~ 1 + a + b + f(idx, model="iid") ## assume we want to compute the posterior for ## ## 2 * beta_a + 3 * beta_b + idx[1] - idx[2] ## ## which we specify as follows (and giving it a unique name) lc1 = inla.make.lincomb(a=2, b=3, idx = c(1,-1,rep(NA,n-2))) names(lc1) = "lc1" ## strictly speaking, it is sufficient to use `idx = c(1,-1)', as the ## remaining idx's are not used in any case. r = inla(formula, data = data.frame(a,b,y), ## ## add the linear combinations here lincomb = lc1, ## ## force noise variance to be essiatially zero ## control.data = list(initial=10, fixed=TRUE), ## (this is the default option) control.inla = list(lincomb.derived.only=TRUE)) ## to verify the result, we can compare the mean but the variance and ## marginal cannot be computed from the simpler marginals alone. lc1.1 = 2 * r$summary.fixed["a", "mean"] + 3 * r$summary.fixed["b", "mean"] + r$summary.random$idx$mean[1] - r$summary.random$idx$mean[2] lc1.2= r$summary.lincomb$"lincomb.lc1.derived"$"mean" print(c(lc1.1 = lc1.1, lc1.2 = lc1.2)) ## the marginals are available as `r$marginals.lincomb$...'. ## we can derive the lincomb above more accurately at a higher cost, ## by adding a new node to the graph, using derived.only=FALSE. Thus ## the lincomb can be derived using the simplified or full laplace ## approximation. For the Gaussian case here, they are all the same... r = inla(formula, data = data.frame(a,b,y), lincomb = lc1, control.data = list(initial=10, fixed=TRUE), control.inla = list(lincomb.derived.only=FALSE)) lc1.3= r$summary.lincomb$"lincomb.lc1"$"mean" print(c(lc1.1 = lc1.1, lc1.2 = lc1.2, lc1.3 = lc1.3)) ## let us compute the posterior for ## ## linear.predictor[2] + intercept ## ## which we specify as follows (and giving it a unique name). ## ## NOTE: The linear predictor has a predefined name = "Predictor", and ## the intercept is named "(Intercept)" (as lm() does, don't blame ## me!). lc2 = inla.make.lincomb(Predictor = c(NA,1, rep(NA,n-2)), "(Intercept)" = 1) names(lc2) = "lc2" ## same here, `Predictor = c(NA,1)' is sufficient. We need to spesify ## them all the way until the last !is.na(). ## to use more than one linear combination, we can spesify more than ## one using the c() operator all.lc = c(lc1,lc2) r = inla(formula, data = data.frame(a,b,y), lincomb = all.lc, control.predictor = list(compute=TRUE), control.data = list(initial=10, fixed=TRUE)) lc2.1 = r$summary.linear.predictor$mean[2] + r$summary.fixed["(Intercept)", "mean"] lc2.2 = r$summary.lincomb$"lincomb.lc2.derived"$"mean" print(c(lc2.1, lc2.1)) ## There is an another function which is handy for specifying many ## linear combinations at once, that is inla.make.lincombs() [note the ## plural 's']. Here each 'column' define one linear combination ## ## let wa and wb be vectors, and we want to compute the marginals for ## ## beta_a * wa[i] + beta_b * wb[i], for i=1..m. this is done ## conveniently as follows m = 10 wa = runif(m) wb = runif(m) lc.many = inla.make.lincombs(a = wa, b=wb) ## we can give them names as well, but there are also default names, like print(names(lc.many)) ## if we want to join in the previous ones, as well, we can add them ## using c() r = inla(formula, data = data.frame(a,b,y), lincomb = c(lc.many, all.lc), control.data = list(initial=10, fixed=TRUE)) for (i in 1:length(r$summary.lincomb)) print(r$summary.lincomb[[i]]$mean) ## if 'wa' and 'wb' is defined in a matrix, like A = cbind(a=wa, b=wb) ## then we can define the same lincombs using inla.uncbind(), lc.many = inla.make.lincombs( inla.uncbind(A) ) ## we can join this trick with "(Intercept) = wi" using (do not ask ## why it is like this: read the code instead. If someone wants to fix ## it, please do!) wi = runif(m) lc.most = inla.make.lincombs( c( inla.uncbind(A), list( "(Intercept)" = wi)) ) ## terms like 'idx' above, can be added as idx = IDX into ## inla.make.lincombs(), where IDX is a matrix. Again, each column of ## the arguments define one linear combination. There is a further option available for the derived linear combinations, that is the option to compute also the posterior correlation matrix between all the linear combinations. To activate this option, use control.inla = list(lincomb.derived.correlation.matrix = TRUE) and you will find the resulting posterior correlation matrix as result$misc$lincomb.derived.correlation.matrix Here is a small example where we compute the correlation matrix for the predicted values of a hidden AR(1) model with an intercept. n = 100 Here are some of the details: The default option, is the compute the marginal for the lincomb, a^Tx, say, where 'a' are the weights, and 'x' is the latent field, as follows. Conditionally on the hyperparameters, we assume the marginal is Gaussian with mean a^TE(x), where E(x) is the mean for x computed with either one of the approximations (gaussian, simplified laplace, laplace), and the variance is a^TQa, where Q is the precision matrix for the Gaussian approximation. Corrections are needed of'course, it there are linear constraints. So, unconditionally on the hyperparameters, the marginal for a^Tx is a mixture of Gaussians. Note that doing it this way, we achive that the lincomb's does not affect the graph of the latent field, and the marginals for the lincomb's can be computed in parallel. All this means, that you can have thousands of (sparse) linear combinations...nPred = 10 phi = 0.9 x = arima.sim(n, model = list(ar=phi)) * sqrt(1-phi^2) y = 1 + x + rnorm(n, sd=0.1) time = 1:(n + nPred) Y = c(y, rep(NA, nPred)) formula = Y ~ 1 + f(time, model="ar1") ## make linear combinations which are the nPred linear predictors B = matrix(NA, nPred, n+nPred) for(i in 1:nPred) { B[i, n+i] = 1 } lcs = inla.make.lincombs(Predictor = B) r = inla(formula, data = data.frame(Y, time), control.predictor = list(compute=TRUE), lincomb = lcs, control.inla = list(lincomb.derived.correlation.matrix=TRUE)) print(r$misc$lincomb.derived.correlation.matrix) I tried one of the examples on this site, and I get this weird error....The documentation is kept up to date, at least we try to, with the TEST version of the software. The version you'll get installing the default package will likely be an older version (but is likely more stable). To install the TEST version, which is updated often several times a week, doinla.upgrade(testing=TRUE) You may revert back to the older version, doing inla.upgrade(). I want to compute/approximate the joint density of some components of the latent field, how can I do that?The easiest approach is to define these components as linear.combinations, and then use the control.inla = list(lincomb.derived.correlation.matrix = TRUE) option, see the related FAQ question about linear combinations. In this case you get the posterior correlation matrix, which you can use as a Gaussian copula.INLA seems to work great for near all cases, but are there cases where INLA is known to have problems?The methodology needs the full conditional density for the latent field to be "near" Gaussian. This is usually achived by either replications or smoothing/"borrowing strength". A simple example which do not have this, is the following:n = 100 u = rnorm(n) eta = 1 + u p = exp(eta)/(1+exp(eta)) y = rbinom(n, size=1, prob = p) idx = 1:n result = inla( y ~ 1 + f(idx, model="iid"), data =data.frame(y,idx), family = "binomial", Ntrials = 1) For each binary observation there is an iid "random effect" `u', and there is no smoothing/``borrowing strength'' (apart from the weak intercept). If you plot the loglikelihood for eta for y=1, say, then its an increasing function for increasing eta, so the likelihood itself would like eta = infinity. With an unknown precision for `u', we run into problems; INLA has a tendency to estimate a to high precision for 'u'. However, it must be noted that the model is almost singular and you'll have a strong prior sensitivity in the (exact) results as well. There is a similar discussion in here as well for the Salamander data example. I can plot the results using 'plot(result)' but how can I save the plots as pdf or postscript files one by one?There is an option to plot (spesific for iNLA),plot(result, single = TRUE) which will not group the plots. Using this option, you can now just 'print' the plot on a spesific window, using dev.print(pdf, file="plot.pdf") for example. #Load libraries library(spdep) library(INLA) #Load data from a shapefile included in the spdep package nc.sids <- readShapePoly(system.file("etc/shapes/sids.shp", package="spdep")[1]) #Create adjacency matrix nc.nb <- poly2nb(nc.sids) #Compute expted number of cases nc.sids$EXP<-nc.sids$BIR74*sum(nc.sids$SID74)/sum(nc.sids$BIR74) #Compute proportion of non-white births nc.sids$NWPROP<-nc.sids$NWBIR74/nc.sids$BIR74 #Convert the adjacency matrix into a file in the INLA format nb2INLA("nc.adj", nc.nb) #Create areas IDs to match the values in nc.adj nc.sids$ID<-1:100 m1<-inla(SID74~NWPROP+f(nc.sids$ID, model="besag", graph.file="nc.adj"), family="poisson", E=nc.sids$EXP, data=as.data.frame(nc.sids), control.predictor=list(compute=TRUE)) #Get realtive risk estimates nc.sids$RR<-m1$summary.fitted.values[,1] #Display relative risk estimates spplot(nc.sids, "RR") |