Model

Our goal is to analyse of diease count data, more specifically when disease cases are aggregated over a partition, say \((\mathcal{R}_{1}, \ldots, \mathcal{R}_{n})\), of the area of interest, \(A\), which can be written mathematically as \[\begin{eqnarray} \label{eq:data} \mathcal{D} = \left\{(y_{i}, d_{i}, \mathcal{R}_{i}): i=1,\ldots,n\right\} \end{eqnarray}\] where \(y_{i}\) and \(d_{i}\) are the number of reported cases and a vector of explanatory variables associated with \(i\)-th region \(\mathcal{R}_{i}\), respectively. Hence, we model \(y_{i}\) conditional on the stochastic process \(S(X)\) as poission distribution with mean \(\lambda_i= m_{i} \exp\{d_{i}\beta^* + S_{i}^*\}\). Then we assume that \(S^* \sim MVN(0, \Sigma)\), where \[\Sigma_{ij} = \sigma^2 \int_{\mathcal{R}_{i}} \int_{\mathcal{R}_{j}} w_i(x) w_j(x') \: \rho(\|x-x'\|; \phi) \: dx \: dx'\], where \(w(x)\) is population density weight. There are two classes of models in this package; one is when we approximate \[S_i^* = \int_{\mathcal{R}_{i}} w_i(x) S^*(x) \: dx \] and the other is \[S_i^* = \frac{1}{\mathcal{R}_{i}} \int_{\mathcal{R}_{i}} S^*(x) \: dx. \] ## Inference We used Monte Carlo Maximum Likelihood for inference. The likelihood function for this class of model is usually intractible, hence we approximate the likelihood function as \[\frac{1}{N}~ \sum_{j=1}^N~\frac{f(\eta_{(j)}; \psi)}{f(\eta_{(j)}; \psi_0)}.\], where \(\psi\) is the vector of the parameters.

This vignette walk you through how to analyse spatial and spatio-temporal dataset using package. Two illustrative examples were provided; application to primary biliary cirrhosis in Newcastle-upon-tyne, UK (static spatial case) and Lung cancer mortality in Ohio, USA (spatio-temporal case).

SDALGCP I (population weighted): application to primary biliary cirrhosis in Newcastle-upon-tyne, UK

Here we estimate the parameters of the model

Discretise the value of scale parameter \(\phi\)

phi <- seq(500, 1700, length.out = 20)

estimate the parameter using MCML

my_est <- SDALGCPMCML(data=data, formula=FORM, my_shp=PBCshp, delta=200, phi=phi, method=1, pop_shp=pop_den, 
                      weighted=TRUE, par0=NULL, control.mcmc=NULL)

To print the summary of the parameter estimates as well as the confidence interval, use;

summary(my_est)
#and for confidence interval use
confint(my_est)

We create a function to compute the confidence interval of the scale parameter using the deviance method. It also provides the deviance plot.

phiCI(my_est, coverage = 0.95, plot = TRUE)

Having estimated the parameters of the model, one might be interested in area-level inference or spatially continuous inference.

1. If interested in STRICTLY area-level inference use the code below. This can either give either region-specific covariate-adjusted relative risk or region-specific incidence. This is achieved by simply setting in the function.

Dis_pred <- SDALGCPPred(para_est=my_est,  continuous=FALSE)

From this discrete inference one can map either the region-specific incidence or the covariate adjusted relative risk.

#to map the incidence
plot(Dis_pred, type="incidence", continuous = FALSE)
#and its standard error
plot(Dis_pred, type="SEincidence", continuous = FALSE)
#to map the covariate adjusted relative risk
plot(Dis_pred, type="CovAdjRelRisk", continuous = FALSE)
#and its standard error
plot(Dis_pred, type="SECovAdjRelRisk", continuous = FALSE)
#to map the exceedance probability that the covariate-adjusted relative risk is greter than a particular threshold
plot(Dis_pred, type="CovAdjRelRisk", continuous = FALSE, thresholds=3.0)

2. If interested in spatially continuous prediction of the covariate adjusted relative risk. This is achieved by simply setting in the function.

Con_pred <- SDALGCPPred(para_est=my_est, cellsize = 300, continuous=TRUE)

Then we map the spatially continuous covariate adjusted relative risk.

#to map the covariate adjusted relative risk
plot(Con_pred, type="relrisk")
#and its standard error
plot(Con_pred, type="SErelrisk")
#to map the exceedance probability that the relative risk is greter than a particular threshold
plot(Con_pred, type="relrisk", thresholds=1.5)

SDALGCP II (Unweighted)

As for the unweighted which is typically by taking the simple average of the intensity an LGCP model, the entire code in the weighted can be used by just setting in the line below.

my_est <- SDALGCPMCML(data=data, formula=FORM, my_shp=PBCshp, delta=200, phi=phi, method=1, 
                      weighted=FALSE,  plot=FALSE, par0=NULL, control.mcmc=NULL, messages = TRUE, plot_profile = TRUE)