A RANDOM EFFECTS MODEL

In the simple random effects model treated here, measurements of some
quantity are available for individuals in a number of groups.  Each
group has associated with it a mean for individuals in the group
(these are the "random effects").  Note that this is the population
mean, not the mean for those individuals in the group that were
actually measured.  The measurements on individuals in the group are
equal to the group mean plus random individual variation.  The
distribution of the group means is modeled as being normal with some
unknown mean and unknown variance.  The individual variation within a
group is also modeled as being normal, with an unknown variance that
is the same for all groups.  The number of individuals in each group
for which there are measurements is considered to be externally
determined, and is thus not modeled.

The number of individuals in each group together with the sample means
and sample variances for each group are sufficient statistics for this
model; only these are provided in the data file 'edata' (one line per
group, in that order).  This data was synthetically generated using
the S-Plus program in sgen.s.  There are 18 groups, with the following
numbers of individuals in each group:

       1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4

The overall mean used to generate the data was 25, the variance of the
group means was 12^2, and the variance within a group was 6^2.


Specifying the model.

We can specify this random effects model as follows:

    > dist-spec elog.met \
       "u~Normal(0,10^2) + v1~Normal(0,2^2) + v2~Normal(0,2^2)" \
       "t0~Normal(u,Exp(v1)+Exp(v2)/i) + \
       (1-Delta(i-1)) * ExpGamma2(Log(t1),i-1,Exp(v2))"

There are three model parameters: the overall mean (u), the log of the
variance of the group means (v1), and the log of the variance of the
measurements within a group (v2).  The variances are represented by
their logs so that they will not have to be constrained.

The argument after the name of the log file to use ("elog.met") is the
prior specification.  The prior for the overall mean is normal with
mean zero and standard deviation 10.  The priors for the logs of the
variances are both normal with mean zero and standard deviation 2.
These parameters are independent under the prior, so the prior
specifications are just added together.

Each group will correspond to a "case", in the terminology used
elsewhere, since given values for the parameters the data on each
group is independent of the other groups.  The number of measurements
for the group will be regarded as an "input", since it is not being
modeled, while the sample mean and sample variance for the group will
be regarded as "targets".  (This distinction is currently just a
convention, but might be essential in future versions of the software.)

The likelihood for the data in a group is given by the last argument
(which is split between two lines).  The sample mean for the group
(target 0, written as "t0") has a normal distribution whose mean is
given by the overall mean parameter and whose variance is the sum of
the variance of the group means (Exp(v1)) and the variance of the
measurements within a group (Exp(v1)) divided by the group size (the
input, written as "i0", or "i").  As is well known, the sample
variance is independent of the sample mean, and has a gamma
distribution.  This is specified by the last term in the likelihood
argument.  There is a complication due to the possibility of a group
with just one measurement, for which the sample variance is undefined
(it will be zero in the data file).  This is handled by including a
factor of (1-Delta(i-1)), which is zero if "i" is one, and is one
otherwise.  The remainder of this term, ExpGamma2(Log(t1),i-1,Exp(v2)), 
evaluates to the log of the density for Log(t1), the log of the sample
variance.  This specifies that t1 has the gamma distribution with
shape parameter (i-1)/2 and mean Exp(v2), which corresponds to the sum
of the squared deviations having a chi-squared distribution with i-1
degrees of freedom.  Note that the syntax using "~" cannot be used
for a transformed variable such as Log(t1).  


Specifying the data source.

The source of the data is specified as follows:

    > data-spec elog.met 1 2 / edata .

For this model, each "training case" is a group.  The above command
says there is one "input" for each group (the number of measurements)
and two "targets" (the sample mean and sample variance for the group).
The inputs come from 'edata', and the targets come from there as well,
following on the same line as the input.


Sampling with Metropolis updates.

We can now specify what Markov chain operations to use in each
iteration.  The following command specifies that each iteration should
consist of 25 Metropolis updates using a proposal distribution in
which each of the three parameters are changed by an amount that is
normally-distributed with standard deviation 0.5:

    > mc-spec elog.met repeat 25 metropolis 0.5

We can now sample for 10000 iterations with the following command:

    > dist-mc elog.met 10000

This takes 2.2 seconds on the system used (see Ex-test-system.doc).

A scatterplot of the posterior distribution for v1 and v2 can be
obtained using a command such as

    > dist-plt v1 v2 elog.met 100: | plot-points

The first 100 iterations are discarded, as possibly not being
representative of the posterior distribution.  This plot should show
most of the points in an ellipse centred at v1=4.5, v2=3.8, which
corresponds to a between-group standard deviation of about 9.5 and a
within-group standard deviation of about 6.7.  However, a few points
are far outside this ellipse, having much smaller values of v1, and
values for v2 around 5.  These points correspond to the possibility
that the group means are almost the same, with the variation in the
sample means for the groups being almost entirely due to sampling
variation resulting from the within-group variance.  Given the data
seen, this is unlikely, but it is not completely excluded.


Slice sampling.

We can also use slice sampling for this problem.  One variation of
slice sampling is specified as follows:

    > mc-spec elog.slc repeat 5 slice-1 2 1
    > dist-mc elog.slc 10000

Each iteration consists of 5 repetitions of a single-variable slice
sampling update for each parameter.  These updates are based on
randomly positioning an interval of width 2 around the current point,
and then sampling successively from this interval until a point within
the slice is found.  The last argument of "1" specifies that the
maximum interval size is just one times the width parameter.  Since
there is no possibility of widening the interval, there is no need to
evaluate the energy function at the interval's endpoints.  This saving
may make this approach beneficial in some circumstances.

These 5 slice sampling operations take about the same time as the 25
Metropolis updates of the previous section.  The resulting sampling
efficiency is also about the same for this example.