CopasLikeSelection: Copas-like selection model

In RobustBayesianCopas: Robust Bayesian Copas Selection Model

Description

This function performs maximum likelihood estimation (MLE) of (θ, τ, ρ, γ_0, γ_1) using the EM algorithm of Ning et al. (2017) for the Copas selection model,

y_i | (z_i>0) = θ + τ u_i + s_i ε_i,

z_i = γ_0 + γ_1 / s_i + δ_i,

corr(ε_i, δ_i) = ρ,

where y_i is the reported treatment effect for the ith study, s_i is the reported standard error for the ith study, θ is the population treatment effect of interest, τ > 0 is a heterogeneity parameter, and u_i, ε_i, and δ_i are marginally distributed as N(0,1), and u_i and ε_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: γ_0 controls the overall probability of publication, and γ_1 controls how the chance of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through ρ. If ρ=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 (γ_0, γ_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 (θ, τ, ρ, γ_0, γ_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 θ.
tau.hat	MLE of τ.
rho.hat	MLE of ρ.
gamma0.hat	MLE of γ_0.
gamma1.hat	MLE of γ_1.
H	5 \times 5 Hessian matrix for the estimates of (θ, τ, ρ, γ_0, γ_1). The square root of the diagonal entries of H can be used to estimate the standard errors for (θ, τ, ρ, γ_0, γ_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

RobustBayesianCopas documentation built on Jan. 13, 2021, 12:50 p.m.