rsmMCL: Response surface method for maximising the Monte Carlo...
In mclcar: Estimating Conditional Auto-Regressive (CAR) Models using Monte Carlo Likelihood Methods

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

The function maximise the Monte Carlo likelihood by using the response surface methodology.

1	rsmMCL(data, psi0, family, exact0 = NULL, control = list(), mc.control = list())

`data, psi0, family, mc.control`	are the same as those in `OptimMCL`
`control`	see more in Details
`exact0`	the exact log likelihood ratio values of a given range or parameters to the importance sampler parameter value. This is used for comparing the fitted response surface of the Monte Carlo likelihood to the exact if the exact evaluation is required.

The response surface methodology finds the maximum of a function by using first and second order polynomial approximations to the target function. The target function is usually unknown or expensive to evaluate. The function is approximated by a first order model in regions that are far away from the maximum and by a second order model around the maximum region. Then the stationary point of the second order model is used as an approximation to the maximum point. By changing the design region, the number and location of the evaluation points, the approximated surface can achieve different accuracy and variance properties. Details of the methodology can be found in Box, G.E.P. and Draper, N.R.(2007) and details of the implementation of the method in this package can be found Sha (2016).

The procedure implemented in this function locates the maximum region of the target function by iteratively updating the design centre and design region. Given a starting point as the design centre, the design region is given based on the finite variance region of the importance sampling approximation as described in Sha (2016), chap. 4. Then with the design region, design points are generated by the central composite design (Box, G.E.P. and Draper, N.R. (2007)). The Monte Carlo likelihood and the corresponding estimated variance at these design points are evaluated and used to fit the first or second order polynomials.

If the stationary point is not found or found to be outside the current design region then, either a steepest ascent analysis (for a first order approximation) or a ridge analysis (for a second order approximation) is done to find the next design centre. The updating step is also limited by considering the estimated variance of the Monte Carlo likelihood.

The procedure is defined to be converged once the stationary point is found and within the design region or the difference between the design centres in two consecutive iterations is smaller than a given tolerance. The Monte Carlo sample size is then increased by a multiple of m once the procedure converges. With the new Monte Carlo sample size a few more iterations are done until the procedure achieve convergence again. The result of the second procedure is given as the final output.

The control argument is a list containing the following elements:

mc.samples,

an array of numbers specifying the Monte Carlo sample size(for direct CAR models only), the increasing multiple and the maximum Monte Carlo sample size to use.

n.iter,

maximum number of iterations; default value is 20

time.max,

maximum total time for the computation

K,

an array of two numbers for generating the design region; the first number K[1] decide the initial size of the design region and once the system find stationary point within the design region, K[1] becomes K[1]/K[2] to get a smaller design region for better approximation

mc.var,

if true, the estimated covariance matrix for the MC-MLE is returned in every iteration; default is FALSE

psi.lim,

a list that contains the range for rho and sigma in the CAR covariance matrix; default values are r.lim = c(-0.2499, 0.2499) for rho and s.lim = c(0.1, 8) for sigma

exacts

(for direct CAR only) a list containing the setting for plotting the exact values for comparison with the fitted rsm surface; the default values are eval = TRUE (for plotting the exact values, default for glm is FALSE), rho = c(-0.25, 0.25), sigma = c(0.5, 2) define the plotting region, length = 100 define number of grid point in each coordinate in producing the plotting grid

rsm.fit

a list of numbers that controls the fitting of the response surfaces and the elements are

n01, number of central points in a central design in generating the design points. The central points are used to estimate the variance of the evaluation and the default is 4.
n02, number of central points in a composite design and the default is 2.
Rsq, R-squared in the fitted first order model in judging the goodness of fit of the first order model; the default value is 0.9.
lof, p-value of the lack of fit test for the first order model; the default value is 0.05.
st.diff, control the range of the approximated function; the approximated surface has less variance for values closes to zero, so this should not be too large and the default value is set to be 50.
mcl.diff, maximum value of the approximated function evaluated at the MC-MLE; again this value should be close to zero and the default is 10.

plotrsm,

if true the approximated response surface with steepest or ridge analysis is plotted; defualt is TRUE

trace.all,

if true, each iteration is stored and returned in the final output, default is TRUE

verbose,

if true, an summary message is printed after each iteration, default is TRUE

When trace.all is TRUE, the function returns a list containing the following objects:

psi0,: a vector of the parameter values used as the initial importance sampler parameter values
Psi,: a list of importance sampler parameter values from all iterations
FVbox,: a list containing the finite variance regions for all iterations
DAs,: a list containing the design region for all iterations
DAbox,: a list containing the design region box for all iterations that can be used for plotting
Ttime,: a list containing the computational time used in each iteration
Expts,: a list of design points in each iteration
Expt.Val1,: a list of all values evaluated at the central design points in each iteration
Expt.Val2,: a list of all values evaluated at the composite design points in each iteration
RSM.1,: a list of all the fitted first-order response model
RSM.2,: a list of all the fitted second-order response model
ESAs,: a list of all explorations along the steepest ascent path
Sim.data,: a list object generated by mcl.prep.dCAR or mcl.prep.glm from all iterations;
Exacts,: a list of exact values in each iteration for the direct CAR models
data,: the data object supplied to the function;
N.iter,: the total number of iterations;
total.time,: the total time elapsed;
convergence,: a logical value indicating whether the procedure converged or not;
mcsamples,: an array of two entries for the initial Monte Carlo sample size and the increased Monte Carlo sample size after the first convergence.

When trace.all is FALSE, the function returns the same list with each object in the list containing the result of the final iteration only.

Zhe Sha zhesha1006@gmail.com

Box, G. E. P. and Draper, N. R. (2007) Response Surfaces, Mixtures, and Ridge Analyses, Wiley, New York

Sha, Z. 2016 Estimating conditional auto-regression models, DPhil Thesis, Oxford.

summary.rsmMCL, plot.rsmMCL, OptimMCL

#### Poisson glm with CAR latent variables
## Simulate some data
set.seed(33)
n.torus <- 10
nb <- 30
rho <- 0.2
sigma <- 1.5
beta <- c(1, 1)
pars.true <- c(rho, sigma, beta)
X0 <- cbind(rep(1, n.torus^2), sample(log(1:n.torus^2)/5))
mydata3 <- CAR.simGLM(method = "poisson", n = c(n.torus, n.torus),
                      pars = pars.true, Xs =as.matrix(X0))

## Fit by rsm
rsm.mcmle2 <- rsmMCL(data = mydata3, psi0 = c(0, 1, 2, 2), family = "poisson",
                     control = list(n.iter = 2, trace.all = TRUE),
                     mc.control = list(N.Zy = 1e3, Scale = 1.65/(n.torus^(2/6)),
                                       thin = 5, burns = 5e2,
                                       method = "mala", scale.fixed = TRUE))
summary(rsm.mcmle2, family = "poisson",  mc.covar=TRUE)
plot(rsm.mcmle2, family = "poisson")