R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Standard error estimation using MCLL

mcll_se

R Documentation

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. α takes a value between 0 and 1, which is the nearest neighbor bandwidth with the kth smallest distance d where k = lfloor n α floor 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`.

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).

Value

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

Author(s)

Minjeong Jeon <jeon.117@osu.edu>

References

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

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