SDALGCPMCML: Parameter estimation for SDA-LGCP Using Monte Carlo Maximum...

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

View source: R/dm_functions.R

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
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
)

Arguments

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.

Details

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.

Value

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.

Author(s)

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

Emanuele Giorgi e.giorgi@lancaster.ac.uk

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

References

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.

See Also

Aggregated_poisson_log_MCML, Laplace.sampling, summary.SDALGCP

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
### 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)

olatunjijohnson/SDALGCP documentation built on March 20, 2021, 4:24 a.m.