Learn R Programming

mclcar (version 0.2-0)

mcl.dCAR: Monte Carlo likelihood calculation for direct CAR models.

Description

This help page contains functions for calculating the Monte Carlo likelihoods and variances for direct CAR models.

mcl.prep.dCAR generates the Monte Carlo samples from the importance sampling distribution to be used in the Monte Carlo likelihood functions.

mcl.dCAR calculates the Monte Carlo log-likelihood ratio for a given parameter value to the importance sampler value.

mcl.profile.dCAR calculates the Monte Carlo profile log-likelihood ratio of rho.

vmle.dCAR calculates the variance at the Monte Carlo MLE.

Avar.lik.dCAR calculates the asymptotic variance of the Monte Carlo likelihood for a direct CAR model on a given size of torus and given size of Monte Carlo samples.

Usage

mcl.prep.dCAR(psi, n.samples, data)

mcl.dCAR(pars, data, simdata, rho.cons = c(-0.249, 0.249), Evar = FALSE) mcl.profile.dCAR(rho, data, simdata, rho.cons = c(-0.249, 0.249), Evar = FALSE)

vmle.dCAR(MLE, data, simdata) Avar.lik.dCAR(pars, psi, data, n.samples, Log = TRUE)

Arguments

psi

the value of the importance sampler parameters

n.samples

the number of Monte Carlo samples

pars

the parameter value for the Monte Carlo likelihood function to be evaluated at

rho

the value of rho for the Monte Carlo profile likelihood function to be evaluated at

data

a data object same as described in loglik.dCAR

simdata

a list object in the same format as the output of mcl.prep.dCAR

rho.cons

rho domain interval

Evar

if TRUE return the estimated variance of the Monte Carlo likelihood

MLE

a list object with

par

the estimated Monte Carlo MLE

hessian

the hessian at the Monte Carlo MLE

Log

if log = TRUE, then the variance is of the log-likelihood ratio; otherwise it is of the likelihood ratio

Value

mcl.prep.dCAR returns a list object containing the following following elements:

  • simYa matrix of the Monte Carlo samples

  • t10 and t20statistics pre-calculate for the importance distribution

mcl.dCAR and mcl.proifle.dCAR return a numeric value of the Monte Carlo likelihood when Evar = FALSE; if TRUE it returns an array of the Monte Carlo likelihood and the corresponding estimated variance.

vmle.dCAR returns the estimated covariance matrix of the Monte Carlo MLE.

Avar.lik.dCAR returns a numeric value of the asymptotic variance.

Details

The asymptotic and estimated variance is derived for a direct CAR model on a \(n \times n\) torus with `rook` style binary neighbourhood matrix. Details of the derivation is given in Sha (2016) ch. 3.

References

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

See Also

mcl.glm, OptimMCL, rsmMCL

Examples

Run this code
# NOT RUN {
## Simulate some data to work with
set.seed(30)
n.torus <- 10
rho <- 0.15
sigma <- 1.5
beta <- c(1, 2)
XX <- cbind(rep(1, n.torus^2), log(1:n.torus^2))
mydata1 <- CAR.simTorus(n.torus, n.torus, rho, 1/sigma)
y <- XX %*% beta + mydata1$X
mydata1$data.vec <- data.frame(y=y, XX[,-1])

## Choose the importance sampler value to the the MPLE
psi1 <- mple.dCAR(mydata1)

## Prepare the Monte Carlo samples
mcdata1 <- mcl.prep.dCAR(psi = psi1, n.samples = 100, data = mydata1)

## Calculate the Monte Carlo likelihoods
mcl.dCAR(c(rho, sigma, beta), data = mydata1, simdata = mcdata1, Evar =
TRUE)
mcl.profile.dCAR(rho, data = mydata1, simdata = mcdata1, Evar = TRUE)

## Do a direct optimization of the Monte Carlo likelihood function
## to find the MLE but it takes very long time to run and may not converge
# }
# NOT RUN {
opt <- optim(psi1, fn = mcl.dCAR,
             data = mydata1, simdata = mcdata1,
             hessian = TRUE, control = list(fnscale = -0.5))
# }
# NOT RUN {
## Assume the we have obtained the opt with the following values
opt<- list(par = c (0.08013547, 1.70294099, 0.01571957, 2.23203089),
           hessian =  matrix(c(  -190.8791352,  -1.5682773,    0.1863733,   -4.844151,
 -1.5682773,  -16.1605186,    0.4911009,    4.403844,
 0.1863733,   0.4911009,  -41.1496774, -156.686631,
-4.8441507,   4.4038442, -156.6866312, -634.316017), 4, 4))

## Calculate the variance of the MC-MLE
vmle.dCAR(opt, data = mydata1, simdata = mcdata1)

## Calculate the asymptotic variance of the likelihood at MC-MLE
Avar.lik.dCAR(pars = opt$par, psi = psi1, data = mydata1, n.samples =
100)

# }

Run the code above in your browser using DataLab