Nothing
################################
## COPAS-LIKE SELECTION MODEL ##
################################
# Obtains maximum likelihood estimates based on the Copas-like
# selection model in Ning et al. (2017)
# INPUTS:
# y = given vector of observed treatment effects
# s = given vector of within-study standard errors
# init = optional initialization values for (theta, tau, rho, gamma0, gamma1)
# If the user does not provide, these are estimated from the data.
# maxit = maximum number of iterations in the optimization
# OUTPUTS:
# theta.hat = estimate for theta
# tau.hat = estimate for tau
# rho.hat = estimate for rho
# gamma0.hat = estimate for gamma0
# gamma1.hat = estimate for gamma1
# H = estimate of Hessian matrix for (theta.hat, tau.hat, rho.hat, gamma0.hat, gamma1.hat)
# Square root of diagonal entries can be used to estimate standard errors
# conv = convergence. "1"=converged, "0"=failed to converge
CopasLikeSelection <- function(y, s, init = NULL, tol=1e-20, maxit=1000){
# Check that y and s are of the same length
if(length(y) != length(s))
stop("Please ensure that 'y' and 's' are of the same length.")
# Check that all values in s are > 0
if(any(s <= 0))
stop("Please ensure that all observed standard errors 's' are greater than zero.")
# Bind 'y' and 's'
data = data.frame(cbind(y,s))
max.s = max(s)
# If no initial values were specified
if(is.null(init)){
### for initial estimates of (1) theta (2) tau (3) rho (4) gamma0 and (5) gamma1
init = rep(0,5)
# initial estimate of theta
init[1]=mean(data[,1])
# moment estimator for tau
init[2]=sqrt(abs(sd(data[,1])^2-(sd(data[,2])^2+mean(data[,2])^2)))
rho.test=seq(-0.98,0.98,by=0.05)
test = vector(mode = "list", length = length(rho.test))
for(k in 1:length(test)){
test[[k]]=optimal.gammas(data,theta=init[1],
tau=init[2],rho=rho.test[k])
}
# Find the rho that maximizes the log-likelihood
loglik.max.index = 1
loglik.max = as.double(test[[1]][2])
for(k in 2:length(test)){
if(as.double(test[[k]][2]) > loglik.max){
loglik.max = as.double(test[[k]][2])
loglik.max.index = k
}
}
# initial value for rho
init[3] = rho.test[loglik.max.index]
# initial value for gamma0
init[4] = as.double(test[[loglik.max.index]][1])
# initial value for gamma1
init[5] = as.double(test[[loglik.max.index]][2])
}
# If initial values were specified but not of correct length
if(!is.null(init)){
if(length(init) != 5)
stop("Please enter a vector of length 5 for initial values of (theta, tau, rho, gamma0, gamma1) in this order.")
if(init[3]>=1){
init[3] = 0.98
} else if(init[3]<=-1){
init[3] = -0.98
}
}
par.new = init # Initialize values
counter = 0 # start counter
while (counter <= maxit) {
counter = counter + 1
par.old = par.new
# E-step for EM algorithm
myM= E.m(par.old,data)
# M-step for EM algorithm
output = stats::optim(par.old,E.loglik,data=data,m=myM,method="L-BFGS-B",lower=c(-Inf,1e-3,-0.99,-Inf,0),
upper=c(Inf,2,0.99,Inf,Inf), control = list(maxit = 1000),hessian=TRUE)
par.new = output$par
#cat("iter = ", counter, "P = ", par.new, "\n")
if (max(abs(par.new-par.old)) < tol){ ###if convergence achieved
# cat("\nSuccessfully Converged\n")
# cat("iter = ", counter, "par = ", par.new, "\n")
theta.hat = par.new[1]
tau.hat = par.new[2]
rho.hat = par.new[3]
gamma0.hat = par.new[4]
gamma1.hat = par.new[5]
# Extract standard errors for theta and tau
H = solve(output$hessian)
return(list(theta.hat = theta.hat,
tau.hat = tau.hat,
rho.hat = rho.hat,
gamma0.hat = gamma0.hat,
gamma1.hat = gamma1.hat,
H = H,
conv=1))
}
}
# print("Convergence Failed")
theta.hat = par.new[1]
tau.hat = par.new[2]
rho.hat = par.new[3]
gamma0.hat = par.new[4]
gamma1.hat = par.new[5]
H = solve(output$hessian)
return(list(theta.hat = theta.hat,
tau.hat = tau.hat,
rho.hat = rho.hat,
gamma0.hat = gamma0.hat,
gamma1.hat = gamma1.hat,
H = H,
conv=0))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.