mcll (version 1.2)

mcll_se: Standard error estimation using MCLL

Description

Standard error estimation using Monte Carlo local likelihood

Usage

mcll_se(data, par, H.prior, alp=0.7,   method="Nelder-Mead",  lower = -Inf, upper = Inf, control=list() )

Arguments

data

posterior samples of model parameters. a matrix or data.frame of size $m \times p$ ($m$: sample size, $p$: dimension of the parameters).

par

MCLL parameter estimates on the original scale.

H.prior

Hessian matrix of the prior evaluated at par.

alp

a real value between 0 and 1. $\alpha$ takes a value between 0 and 1, which is the nearest neighbor bandwidth with the $k$th smallest distance $d$ where $k = \lfloor n \alpha \rfloor$ and $d(x, x_{i}) = | x - x_{i} |$ with the sample size $n$.

method

an optimization method to be used to find the coefficients of the polynomial approximation to the log-posterior at the MCLL estimates par. Options from optim are Nelder-Mead, BFGS, CG, L-BFGS-B, and SANN.

lower, upper

bounds on the variables for the L-BFGS-B method in optim.

control

a list of control parameters. See control options for optim.

Value

mcll_se returns a vector containing standard error estimates for the MCLL parameter estimates par.

Details

Standard error estimation in the Monte Carlo local likelihood method. For details, see Section 3 in Jeon et al. (2012). The posterior samples and paramter values should be on the real line (e.g., variance parameters should be in the log-scale).

References

Jeon, M., Kaufman, C., and Rabe-Hesketh, S. (2014). Monte Carlo local likelihood for approximate MLE of complex models. Under revision.

See Also

mcll_est

Examples

## example 

# data preparation
data(samp) 

# prior function
prior.func <- function(vec.t) {
    sum(dnorm(vec.t, m= c(0,0,0,0, -0.9870405, -0.9870405) ,
                sd=c(100,100,100,100, 1/0.766672, 1/0.766672) , log=TRUE))
}

## parameter estimation
run1 <- system.time(
    result1 <- mcll_est(data=samp, prior.func= prior.func, alp=0.7, 
        method = "BFGS", control= list(maxit=10000) )
)

par <- result1$par # original scale 

## standard error estimation 

# H.prior: analytical solution 
p.var = c(100,100,100,100, 1/0.766672, 1/0.766672)^2
H.prior <- -diag(1/p.var)

# H.prior: numerical solution 
# library(numDeriv)
# H.prior <- hessian(prior.func, par)

# SE estimation (NOT RUN)
#run2 <- system.time(
#    se <- mcll_se(data=samp, par=par, H.prior = H.prior,  alp=0.7, 
#        method= "Nelder-Mead" , control=list(maxit=20000)) 
#)

#se
#       b0        b1        b2        b3      tau0      tau1 
#0.4057844 0.5640063 0.4907643 0.6663096 0.3022842 0.2999727