**
****
**Schools: ranking schoolexamination

results using multivariate

hierarcical models

**
**Goldstein et al. (1993) present an analysis of examination results from inner London schools. They use hierarchical or multilevel models to study the between-school variation, and calculate school-level residuals in an attempt to differentiate between `good' and `bad' schools. Here we analyse a subset of this data and show how to calculate a rank ordering of schools and obtain credible intervals on each rank.

Standardized mean examination scores (Y) were available for 1978 pupils from 38 different schools. The median number of pupils per school was 48, with a range of 1--198. Pupil-level covariates included gender plus a standardized London Reading Test (LRT) score and a verbal reasoning (VR) test category (1, 2 or 3, where 1 represents the highest ability group) measured when each child was aged 11. Each school was classified by gender intake (all girls, all boys or mixed) and denomination (Church of England, Roman Catholic, State school or other); these were used as categorical school-level covariates.

We consider the following model, which essentially corresponds to Goldstein et al.'s model 1.

Y

m

+ b

+ b

log t

where i refers to pupil and j indexes school. We wish to specify a regression model for the variance components, and here we model the logarithm of t

The fixed effects b

which represents our prior guess at the order of magnitude of

The BUGS code is given below:

model

{

for(p in 1 : N) {

Y[p] ~ dnorm(mu[p], tau[p])

mu[p] <- alpha[school[p], 1] + alpha[school[p], 2] * LRT[p]

+ alpha[school[p], 3] * VR[p, 1] + beta[1] * LRT2[p]

+ beta[2] * VR[p, 2] + beta[3] * Gender[p]

+ beta[4] * School.gender[p, 1] + beta[5] * School.gender[p, 2]

+ beta[6] * School.denom[p, 1] + beta[7] * School.denom[p, 2]

+ beta[8] * School.denom[p, 3]

log(tau[p]) <- theta + phi * LRT[p]

sigma2[p] <- 1 / tau[p]

LRT2[p] <- LRT[p] * LRT[p]

}

min.var <- exp(-(theta + phi * (-34.6193))) # lowest LRT score = -34.6193

max.var <- exp(-(theta + phi * (37.3807))) # highest LRT score = 37.3807

# Priors for fixed effects:

for (k in 1 : 8) {

beta[k] ~ dnorm(0.0, 0.0001)

}

theta ~ dnorm(0.0, 0.0001)

phi ~ dnorm(0.0, 0.0001)

# Priors for random coefficients:

for (j in 1 : M) {

alpha[j, 1 : 3] ~ dmnorm(gamma[1:3 ], T[1:3 ,1:3 ]);

alpha1[j] <- alpha[j,1]

}

# Hyper-priors:

gamma[1 : 3] ~ dmnorm(mn[1:3 ], prec[1:3 ,1:3 ]);

T[1 : 3, 1 : 3 ] ~ dwish(R[1:3 ,1:3 ], 3)

}

Note that school is a 1978 x 3 matrix taking value 1 for all pupils in school 1, 2 for all pupils in school 2 and so on. For computational convenience, Y , mu and tau are indexed over a single dimension p = 1,...,1978 rather than as pupil i within school j as used in equations above. The appropriate school-level coefficients for pupil p are then selected using the school indicator in row p of the data array --- for example alpha[school[p],1] .

Results

A 1000 update burn in followed by a further 10000 updates gave the parameter estimates

The school-specific intercept a