CopasLikeSelection {RobustBayesianCopas}

Copas-like selection model

Description

This function performs maximum likelihood estimation (MLE) of (\theta, \tau, \rho, \gamma_0, \gamma_1) using the EM algorithm of Ning et al. (2017) for the Copas selection model,

y_i | (z_i>0) = \theta + \tau u_i + s_i \epsilon_i,

z_i = \gamma_0 + \gamma_1 / s_i + \delta_i,

corr(\epsilon_i, \delta_i) = \rho,

where y_i is the reported treatment effect for the ith study, s_i is the reported standard error for the ith study, \theta is the population treatment effect of interest, \tau > 0 is a heterogeneity parameter, and u_i, \epsilon_i, and \delta_i are marginally distributed as N(0,1), and u_i and \epsilon_i are independent.

In the Copas selection model, y_i is published (selected) if and only if the corresponding propensity score z_i (or the propensity to publish) is greater than zero. The propensity score z_i contains two parameters: \gamma_0 controls the overall probability of publication, and \gamma_1 controls how the chance of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through \rho. If \rho=0, then there is no publication bias and the Copas selection model reduces to the standard random effects meta-analysis model.

This is called the "Copas-like selection model" because to find the MLE, the EM algorithm utilizes a latent variable m that is treated as missing data. See Ning et al. (2017) for more details. An alternative funtion for implementing the Copas selection model using a grid search for (\gamma_0, \gamma_1) is available in the R package metasens.

Usage

CopasLikeSelection(y, s, init = NULL, tol=1e-20, maxit=1000)

Arguments

y

An n \times 1 vector of reported treatment effects.

s

An n \times 1 vector of reported within-study standard errors.

init

Optional initialization values for (\theta, \tau, \rho, \gamma_0, \gamma_1). If specified, they must be provided in this exact order. If they are not provided, the program estimates initial values from the data.

tol

Convergence criterion for the Copas-like EM algorithm for finding the MLE. Default is tol=1e-20.

maxit

Maximum number of iterations for the Copas-like EM algorithm for finding the MLE. Default is maxit=1000.

Value

The function returns a list containing the following components:

theta.hat

MLE of \theta.

tau.hat

MLE of \tau.

rho.hat

MLE of \rho.

gamma0.hat

MLE of \gamma_0.

gamma1.hat

MLE of \gamma_1.

H

5 \times 5 Hessian matrix for the estimates of (\theta, \tau, \rho, \gamma_0, \gamma_1). The square root of the diagonal entries of H can be used to estimate the standard errors for (\theta, \tau, \rho, \gamma_0, \gamma_1).

conv

"1" if the optimization algorithm converged, "0" if algorithm did not converge. If conv=0, then using H to estimate the standard errors may not be reliable.

References

Ning, J., Chen, Y., and Piao, J. (2017). "Maximum likelihood estimation and EM algorithm of Copas-like selection model for publication bias correction." Biostatistics, 18(3):495-504.

Examples

####################################
# Example on the Hackshaw data set #
####################################
data(Hackshaw1997)
attach(Hackshaw1997)
# Extract the log OR
y.obs = Hackshaw1997[,2]
# Extract the observed standard error
s.obs = Hackshaw1997[,3]

##################################
# Fit Copas-like selection model #
##################################

# First fit RBC model with normal errors
RBC.mod = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="normal", seed=123, burn=500, nmc=500)

# Fit CLS model with initial values given from RBC model fit.
# Initialization is not necessary but the algorithm will converge faster with initialization.
CLS.mod = CopasLikeSelection(y=y.obs, s=s.obs, init=c(RBC.mod$theta.hat, RBC.mod$tau.hat,
                                                       RBC.mod$rho.hat, RBC.mod$gamma0.hat,
                                                    RBC.mod$gamma1.hat))

# Point estimate for theta 
CLS.theta.hat = CLS.mod$theta.hat  

# Use Hessian to estimate standard error for theta
CLS.Hessian = CLS.mod$H
# Standard error estimate for theta
CLS.theta.se = sqrt(CLS.Hessian[1,1]) # 

# 95 percent confidence interval 
CLS.interval = c(CLS.theta.hat-1.96*CLS.theta.se, CLS.theta.hat+1.96*CLS.theta.se)

# Display results
CLS.theta.hat  
CLS.theta.se  
CLS.interval   

# Other parameters controlling the publication bias
CLS.mod$rho.hat 
CLS.mod$gamma0.hat
CLS.mod$gamma1.hat