1
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(
c(
1, 0,
1, 1,
1, 1,
1, 0,
1, 0,
1, 1,
1, 1,
1, 1,
1, 0,
1, 0,
1, 0,
1, 0,
1, 0,
1, 0,
1, 1,
1, 1,
1, 0,
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)
}