Publication Bias in Research Synthesis: Sensitivity Analysis Using a Priori Weight Functions

Sensitivity Analysis

An S-Plus Implementation

Here, we present computer code for estimating the model using Splus. Each group of commands is preceded by a notation of the assumed location in a folder hierarchy; if the user elects to alter that organization, the required changes will be self-evident.

Code for random-effects models. We assume that files are located in a folder named “Swork” on the desktop of a user named “user.” To estimate the model, invoke Splus, and type:

source(“c:/documents and settings/user/desktop/control.sin”) <enter>

The output will consist of:

1.) A listing of parameter estimates (intercept, slopes, variance component);

2.) Their standard errors; and

3.) Adjusted parameter estimates.

The following text should appear in a file in Swork named “control.sin.”

# Get effect sizes, variances, and predictor matrix

# Note to user: enter effect sizes, conditional variances

# and predictor matrix in the file "data.sin"

source("C:/documents and settings/user/desktop/Swork/data.sin")

# Get information about weights

# Note to user: enter information about weights and

# cutpoints in the file "weights.sin"

source("C:/documents and settings/user/desktop/Swork/weights.sin")

# Set a starting point for the vector of regression

# coefficients using conventional fixed-effects model

source("C:/documents and settings/user/desktop/Swork/functions/beta_setup.sin")

# Calculate bij, and establish Bij() function

source("C:/documents and settings/user/desktop/Swork/functions/B_and_b_setup.sin")

# Establish likelihood functions for unweighted and

# weighted models

source("C:/documents and settings/user/desktop/Swork/functions/like.sin")

# Establish standard error function for unweighted model

source("C:/documents and settings/user/desktop/Swork/functions/se.sin")

# Fit unadjusted model by maximum likelihood

out1 <- nlminb(objective=like1, start=params,lower=c(rep(-Inf,p+1),0.0))

par1 <- out1$parameters

se1 <- se(out1$parameters)

# Fit adjusted model

out2 <- nlminb(objective=like2, start=out1$parameters,lower=c(rep(-Inf,p+1),0.0))

par2 <- out2$parameters

print("Unadjusted parameter estimates:")

print(par1)

print("Standard errors:")

print(se1)

print("Adjusted estimates:")

print(par2)

The following text should appear in a file in Swork named “data.sin.”

# Enter number of predictors (not including intercept)

npred <- 1

# Enter effect sizes.

# Here, effect sizes are for the second example described in the paper

# (teacher expectancy data).

myT <- c(

.03, .12, -.14, 1.18, .26,

-.06, -.02, -.32, .27, .80,

.54, .18, -.02, .23, -.18,

-.06, .30, .07, -.07)

# Enter conditional variances.

# (Here, standard errors of effect sizes were supplied; hence, conditional

# variances were entered as the squares (^2) of those standard errors. If

# data are supplied as variances, eliminate the square.

v <- c(

.125^2, .147^2, .167^2, .373^2, .369^2,

.103^2, .103^2, .220^2, .164^2, .251^2,

.302^2, .223^2, .289^2, .290^2, .159^2,

.167^2, .139^2, .094^2, .174^2)

# Enter matrix of predictors, with first column always 1 for the intercept.

# (Here, the second column contains a dummy indicator of familiarity status:

# 0 = teacher has know child two weeks or less, 1 = teacher has known child

# longer than two weeks.)

X <- matrix(

1, 0,

1, 1,

1, 0,

1, 1,

1, 0,

1, 1,

1, 0,

1, 1

),length(myT),npred+1,byrow=TRUE)

# Do not modify below this line.

# Set n = number of effect sizes

n <- length(myT)

if(length(v) != n) stop("Number of conditional variances must equal number of effects.")

for(i in 1:length(v)) if(v[i] <0) stop("Variances must be non-negative.")

s <- sqrt(v)

if(dim(X)[1] != n) stop("Number of rows in predictor matrix must match number of effects.")

# Set p = number of predictors

p <- dim(X)[2]-1

The following text should appear in a file in Swork named “weights.sin.”

# In this example, cutpoints and weights are for the condition denoted

# “severe two-tailed selection” in the paper.

# Enter cutpoints for p-value intervals

a <- c(.005,.010,.050,.100,.250,.350,.500,.650,.750,.900,.950,.990,.995,1.00)

# Enter fixed weight function

w <- matrix(

c(1.0,.99,.90,.75,.60,.50,.25,.25,.50,.60,.75,.90,.99,1.0),

ncol=1)

# Do not modify below this point

# Set k = number of intervals

k <- length(a)

if(length(w) != k) stop("Number of weights must match number of intervals")

The remaining code is assumed to be located in a subfolder of “Swork” called “functions.” No file in “functions” should require alteration by the user.

The following code should appear in a file in Swork\functions called “B_and_b_setup.sin.”

# Here we set up Bij and bij

# First, we need predicted values

xb <- X%*%beta

sig <- sqrt(v + vc)

# Now bij = -si probit(aj)

b <- matrix(rep(0,n*(k-1)),nrow=n,ncol=k-1)

for(i in 1:n)

for(j in 1:k-1)

b[i,j] <- -s[i] * qnorm(a[j])

# And now Bij (Equation 5)

Bij <- function(params) {

B <- matrix(rep(0,n*k),nrow=n,ncol=k)

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

for(i in 1:n)

for(j in 1:k) {

if(j==1) B[i,j] <- 1.0 - pnorm( (b[i,1]-xb[i]) / sig[i])

else if(j<k) B[i,j] <- pnorm( (b[i,j-1]-xb[i]) / sig[i]) - pnorm( (b[i,j]-xb[i]) / sig[i])

else B[i,j] <- pnorm( (b[i,k-1]-xb[i]) / sig[i])

}

return(B)

}

The following code should appear in a file in Swork\functions called “beta_setup.sin.”

# Do weighted least squares. (Don’t do this at home; the algorithm is unstable

# for general regression problems, but it’s easy and unlikely to cause problems

# in this context.)

varmat <- diag(v)

tmat <- matrix(myT,nrow=n)

beta <- solve( t(X)%*%solve(varmat)%*%X)%*%t(X)%*%solve(varmat)%*%tmat

# Also, set variance component to balanced method of moments start

rss <- sum((myT-X%*%beta)^2)

vc <- rss/(n-p-1) - mean(v)

if (vc < 0.0) vc <- 0.0

# Set starting value for parameter vector

params <- c(beta,vc)

The following code should appear in a file in Swork\functions called “like.sin.”

# Unadjusted likelihood function

like1 <- function(params) {

vc <- params[p+2]

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

vall <- v+vc

-2*(-1/2 * (sum( (myT-xb)^2/vall ) + sum(log(vall))) )

}

# Adjusted likelihood function (Equation 6)

like2 <- function(params) {

vc <- params[p+2]

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

vall <- v+vc

ll <- -1/2 * (sum( (myT-xb)^2/vall ) + sum(log(vall)))

B <- Bij(params)

Bw <- B%*%w

ll <- ll - sum(log(Bw))

return(-2*ll)

}

The following code should appear in a file in Swork\functions called “se.sin.”

# Second derivatives, unadjusted model

se <- function(params) {

vc <- params[p+2]

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

vall <- v+vc

hessian <- matrix(rep(0,(p+2)*(p+2)),p+2,p+2)

for(i in 1:n) {

for(j in 1:(p+1))

for(m in 1:(p+1))

hessian[j,m] <- hessian[j,m] - X[i,j]*X[i,m]/vall[i]

for(j in 1:(p+1)) {

hessian[j,p+2] <- hessian[j,p+2] - (myT[i]-xb[i])/vall[i]/vall[i]

hessian[p+2,j] <- hessian[p+2,j] - (myT[i]-xb[i])/vall[i]/vall[i]

}

hessian[p+2,p+2] <- hessian[p+2,p+2] + 1/vall[i]/vall[i]/2 -

(myT[i]-xb[i])*(myT[i]-xb[i])/vall[i]/vall[i]/vall[i]

}

temp <- solve(hessian)

retval <- rep(0,(p+2))

for(i in 1:(p+2)) retval[i] <- sqrt(-temp[i,i])

return(retval)

}

Code for fixed-effects models. To estimate the fixed-effects model, the organization of files in the hierarchy of folders is identical to the organization for the random-effects model. Effect sizes, variances, and weights are entered in the same fashion. The following file replacements must be made.

Replace control.sin with the following:

# Get effect sizes, variances, and predictor matrix

# Note to user: enter effect sizes, conditional variances

# and predictor matrix in the file "data.sin"

source("C:/documents and settings/user/desktop/Swork/data.sin")

# Get information about weights

# Note to user: enter information about weights and

# cutpoints in the file "weights.sin"

source("C:/documents and settings/user/desktop/Swork/weights.sin")

# Set a starting point for the vector of regression

# coefficients using conventional fixed-effects model

source("C:/documents and settings/user/desktop/Swork/functions/beta_setup.sin")

# Calculate bij, and establish Bij() function

source("C:/documents and settings/user/desktop/Swork/functions/B_and_b_setup.sin")

# Establish likelihood functions for unweighted and

# weighted models

source("C:/documents and settings/user/desktop/Swork/functions/like.sin")

# Establish standard error function for unweighted model

source("C:/documents and settings/user/desktop/Swork/functions/se.sin")

# Fit unweighted model by maximum likelihood

out1 <- nlminb(objective=like1, start=params)

par1 <- out1$parameters

se1 <- se(out1$parameters)

# Fit weighted model

out2 <- nlminb(objective=like2, start=out1$parameters)

par2 <- out2$parameters

print("Unadjusted parameter estimates:")

print(par1)

print("Standard errors:")

print(se1)

print("Adjusted estimates:")

print(par2)

Replace functions\B_and_b_setup.sin with the following:

# Here we set up Bij and bij

# First, we need predicted values

xb <- X%*%beta

sig <- sqrt(v)

# Now bij = -si probit(aj)

b <- matrix(rep(0,n*(k-1)),nrow=n,ncol=k-1)

for(i in 1:n)

for(j in 1:k-1)

b[i,j] <- -s[i] * qnorm(a[j])

# And now Bij (Equation 5)

Bij <- function(params) {

B <- matrix(rep(0,n*k),nrow=n,ncol=k)

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

for(i in 1:n)

for(j in 1:k) {

if(j==1) B[i,j] <- 1.0 - pnorm( (b[i,1]-xb[i]) / sig[i])

else if(j<k) B[i,j] <- pnorm( (b[i,j-1]-xb[i]) / sig[i]) - pnorm( (b[i,j]-xb[i]) / sig[i])

else B[i,j] <- pnorm( (b[i,k-1]-xb[i]) / sig[i])

}

return(B)

}

Replace functions\beta_setup with the following:

# Do weighted least squares (Don’t do this at home; the algorithm is unstable

# for general regression problems, but it’s easy and unlikely to cause problems

# in this context.)

varmat <- diag(v)

tmat <- matrix(myT,nrow=n)

beta <- solve( t(X)%*%solve(varmat)%*%X)%*%t(X)%*%solve(varmat)%*%tmat

# Set starting value for parameter vector

params <- beta

Replace functions\like.sin with the following:

# Unadjusted likelihood function

like1 <- function(params) {

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

-2*(-1/2 * (sum( (myT-xb)^2/v ) + sum(log(v))) )

}

# adjusted likelihood function (Equation 6)

like2 <- function(params) {

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

ll <- -1/2 * (sum( (myT-xb)^2/v ) + sum(log(v)))

B <- Bij(params)

Bw <- B%*%w

ll <- ll - sum(log(Bw))

return(-2*ll)

}

Finally, replace functions\se.sin with the following:

# Second derivatives, unadjusted model

se <- function(params) {

beta <- matrix(params[1:(p+1)],ncol=1)

xb <- X%*%beta

hessian <- matrix(rep(0,(p+1)*(p+1)),p+1,p+1)

for(i in 1:n) {

for(j in 1:(p+1))

for(m in 1:(p+1))

hessian[j,m] <- hessian[j,m] - X[i,j]*X[i,m]/v[i]

}

temp <- solve(hessian)

retval <- rep(0,(p+1))

for(i in 1:(p+1)) retval[i] <- sqrt(-temp[i,i])

return(retval)

}