RobustBayesianCopas (version 2.0)

CopasLikeSelection: 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 \(i\)th study, \(s_i\) is the reported standard error for the \(i\)th 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

# NOT RUN {
####################################
# 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

# }