SDALGCPMCML: Parameter estimation for SDA-LGCP Using Monte Carlo Maximum...
In olatunjijohnson/SDALGCP: Spatially Discrete Approximation to Log-Gaussian Cox Processes for Aggregated Disease Count Data

Description Usage Arguments Details Value Author(s) References See Also Examples

This function provides the maximum likelihood estimation of the parameter given a set of values of scale parameter of the Gaussian process, phi.

SDALGCPMCML(
  formula,
  data,
  my_shp,
  delta,
  phi = NULL,
  method = 1,
  pop_shp = NULL,
  weighted = FALSE,
  par0 = NULL,
  control.mcmc = NULL,
  plot = FALSE,
  plot_profile = TRUE,
  rho = NULL,
  giveup = NULL,
  messages = FALSE
)

`formula`	an object of class `formula` (or one that can be coerced to that class): a symbolic description of the model to be fitted.
`data`	data frame containing the variables in the model.
`my_shp`	A SpatialPolygons or SpatialPolygonsDataFrame object containing the polygons (i.e each regions).
`delta`	distance between points
`phi`	the discretised values of the scale parameter phi. if not supplied, it uses the default, which is 20 phis' which ranges from size of the smallest region to the one-tenth of the size of the entire domain.
`method`	To specify which method to use to sample the points, the options are 1 for Simple Sequential Inhibition (SSI) process, 2 for Uniform sampling and 3 for regular grid. 1 is the default
`pop_shp`	Optional, The raster of population density map for population weighted approach
`weighted`	To specify if you want to use the population density, default to FALSE, i.e population density is not used.
`par0`	the initial parameter of the fixed effects beta, the variance sigmasq and the scale parameter phi, specified as c(beta, sigma2, phi). Default; beta, the estimates from the glm; sigma2, variance of the residual; phi, the median of the supplied phi.
`control.mcmc`	list from PrevMap package to define the burnin, thining, the number of iteration and the turning parameters see `controlmcmcSDA`.
`plot`	To display the plot of the points inside the polygon, default to TRUE
`plot_profile`	logical; if TRUE the profile-likelihood is plotted. default is FALSE
`rho`	Optional, the packing density, default set to 0.55
`giveup`	Optional, number of rejected proposals after which the algorithm should terminate, default set to 1000
`messages`	logical; if messages=TRUE, it prints the results objective function and the parameters at every phi iteration. Default is FALSE.

This function performs parameter estimation for a SDALGCP Model Monte Carlo Maximum likelihood. The Monte Carlo maximum likelihood method uses conditional simulation from the distribution of the random effect T(x) = d(x)'β+S(x) given the data y, in order to approximate the high-dimensional intractable integral given by the likelihood function. The resulting approximation of the likelihood is then maximized by a numerical optimization algorithm which uses analytic expression for computation of the gradient vector and Hessian matrix. The functions used for numerical optimization are nlminb. The first stage of estimation is generating locations inside the polygon, followed by precomputing the correlation matrices, then optimising the likelihood.

An object of class "SDALGCP". The function summary.SDALGCP is used to print a summary of the fitted model. The object is a list with the following components:

D: matrix of covariates.

y: the count, response observations.

m: offset

beta_opt: estimates of the fixed effects of the model.

sigma2_opt: estimates of the variance of the Gaussian process.

phi_opt: estimates of the scale parameter phi of the Gaussian process.

cov: covariance matrix of the MCML estimates.

Sigma_mat_opt: covariance matrix of the Gaussian process that corresponds to the optimal value

llike_val_opt: maximum value of the log-likelihood.

mu: mean of the linear predictor

all_para: the entire estimates for the different values of phi.

all_cov: the entire covariance matrix of the estimates for the different values of phi.

par0: the initial parameter of the fixed effects beta and the variance sigmasq used in the estimation

control.mcmc: the burnin, thining, the number of iteration and the turning parameters used see controlmcmcSDA.

call: the matched call.

Olatunji O. Johnson o.johnson@lancaster.ac.uk

Emanuele Giorgi e.giorgi@lancaster.ac.uk

Peter J. Diggle p.diggle@lancaster.ac.uk

Giorgi, E., & Diggle, P. J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software, 78(8), 1-29. doi:10.18637/jss.v078.i08

Christensen, O. F. (2004). Monte Carlo maximum likelihood in model-based geostatistics. Journal of Computational and Graphical Statistics 13, 702-718.

Aggregated_poisson_log_MCML, Laplace.sampling, summary.SDALGCP

### Prepare the input of the model
data(PBCshp)
data <- as.data.frame(PBCshp@data)  #get the data
### Write the formula of the model
FORM <- X ~ propmale + Income + Employment + Education + Barriers + Crime +
Environment +  offset(log(pop))
### set the discretised phi
phi <- seq(500, 1700, length.out = 10)
#### get the initial parameter
model <- glm(formula=FORM, family="poisson", data=data)
beta.start <-coef(model)
sigma2.start <- mean(model$residuals^2)
phi.start <- median(phi)
par0 <- c(beta.start, sigma2.start, phi.start)
# setup the control arguments for the MCMC
n <- 545
h <- 1.65/(n^(1/6))
control.mcmc <- controlmcmcSDA(n.sim = 10000, burnin = 2000,
                 thin= 8, h=h, c1.h = 0.01, c2.h = 1e-04)
###Run the model

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