:
~—.J...~‘
_ ‘fi'
1
1

‘1
M; ,1
2;?-

., m .
., .

x s

._._
“er

:1.

I‘M.

1'1
.

,.
fl .,
r < ‘
Iii"

\ 'u
"1 v

47 . V
‘ gull-1L if:
I; 1 ‘-

4:
1,:

, Tu“
a «"12;
rt‘x'é-‘f'L.-.‘%".}" I "H

.73

g:
5:29!

.. ‘ r-a, 54 "105‘,
«5.142? v.5, u 5. "w t

 

I’K'r.
,.

 

'  willWilli\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ L

This is to certify that the

dissertation entitled

THE APPLICATION OF GIBBS SAMPLING T0 NESTED VARIANCE
COMPONENTS MODELS WITH HETEROGENOUS WITHIN-GROUP

VARIANCE

presented by

RAFA M. KASIM

has been accepted towards fulfillment
of the requirements for

Ph.D. Measurement & Qunt. Mthd.

degree in

 

 

Dr. Stephen W. Raudenbush

 

Date ”M /S/ /ffl7"

 

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE
JW— 7
i .
rm 0?) i7;

 

 

 

 

8594012272003
W953 2 80110

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
abimmmamn

 

 

 

 

THE APPLICATION OF GIBBS SAMPLING TO
NESTED VARIANCE COMPONENTS MODELS
WITH HETEROGENOUS WITHIN-GROUP VARIANCE

BY

Rafa M. Kasim

A DISSERTATION

Submitted to
Michigan State University
College of Education
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Counseling, Educational
Psychology and Special Education

1994

 

 

 

 

ABSTRACT
THE APPLICATION OF GIBBS SAMPLING TO

NESTED VARIANCE COMPONENTS MODELS
WITH HETEROGENOUS WITHIN-GROUP VARIANCE

BY
Rafa M. Kasim

Research in social sciences such as education,
psychology, and sociology often involves analyzing
hierarchically structured data. The use of Bayesian
procedures ‘with. hierarchical linear' regression :models ‘to
analyze such data helped researchers obtain answers to
questions where standard approaches seem to fail. It also
improved the estimates of the regression coefficients by
obtaining their empirical Bayes estimates.

Empirical Bayes estimates of the regression coefficients
at two levels of hierarchy are usually conditioned on the
maximum likelihood estimates of the variance components, and
often assuming homogeneity of within-group variance. There
are times however, where research interest is focused on
studying the variance components, especially when there exists
a clear evidence of heterogeneity of variance, and there is a
great concern about the effect of the uncertainty in
estimating these variances on the.empirical Bayes estimates of
the regression coefficients. A mixed model with random
intercept, which permits heterogenous within—group variances, {03-}
for j=1,---,k, across the k units is introduced within the

Bayesian approach. Information about the within—group

 

 

 

 

variance and the intercept from these k units represent their
prior distributions in a second level of the analysis. The
marginal posterior distributions of the variance components
and regression coefficients parameters are obtained via Gibbs
sampling. The goal is to obtain Bayes estimates of the
parameters of interest, especially those of the within-group
variances, and study the effect of adjusting for the
uncertainty in estimating the variance components on the
estimation of the regression.coefficientsu IBayes estimates of
the parameters for the specified model, with heterogenous
within-group variance obtained via Gibbs sampling, are
compared to their empirical Bayes estimates obtained via HLM
analysis assuming homogeneity of variance.

The process of Gibbs sampling applied on several
artificial data sets, and on one real data set to obtain the
marginal posterior distributions of the parameters for the
specified model. The real data set represents a nationwide
random sample from the high schools in the U.S. Eighteen
artificial data sets wereigenerated.to represent three models.
The variation between the 18 data sets covers the complexity
of the model used (number of predictors at the two levels),
the degree of heterogeneity of within-group variance, and the
number of groups in level-2. These data sets were generated
to reflect the different models used in analyzing the real
data set.

There appears to be no substantial difference between

Bayes estimates (posterior means) and empirical Bayes

 

 

 

estimates of the regression coefficients. Note that empirical
Bayes estimates are based on the assumption of homogeneity of
‘within-group variance, o§=m£. When it comes to the estimation
of the variance components, HLM estimates of the within-group
variance component, of are found to be positively biased,

especially when there exists clear evidence of heterogeneity

of variance.

 

 

 

ACKNOWLEDGEMENT

I would like to thank my dissertation director Dr.
Stephen Raudenbush for his constant help, support, and advice
through my study and the writing of this dissertation. I am
very fortunate to have him as a teacher, as an advisor, and
best of all as friend. I would also like to express my
appreciation to Dr. Dennis Gilliland for his valuable
suggestion and comments on my work, Having him.as a member in
my dissertation committee helped me learn a lot from him. I
would also like to thank. the rest of the dissertation
committee members Dr. Betsy Becker and Dr. Richard Houang for
their support and encouragement" I would also like to express
my appreciation to Dr. Michael Seltzer for his insightful
comments on my work. My wife Hripsime and my kids Nader,
Neda, and Nabeel deserve special thanks for their
understanding and patience so I can finish my dissertation.
Finally, but most importantly, I would like to express my
deepest gratitude to my parents and.my brothers and sister for

their endless love and support.

 

II

III

TABLE OF CONTENTS

Introduction . . . . . . . . . . . . . .
Purpose of the Study . . . . . . . . . .

Objectives of the Study . . . . . . . .
Bayes Solution . . . . . . . . .

Importance of the Study . . . . . . . .

Organization of the Study . . . . . . .

Statement of the Problem . . . . . .

Problems Associated with the Estimation
of 02 and T . . . . . . . . . . . . . .

Problems Associated with Estimating T . . . .
Problems Associated with an Invalid Assumption
of Homogeneity of Variance . . . . . . . . . .
The Analysis of Many Groups Variances a; . . . . .
Bayesian Approach . . . . . . . . . . . . .
Comments on the Analysis of Many Variances . .

Adjusting for the Uncertainty in

Estimating a? and T . . . . . . . . .
Model Specification . . . . . . . . . .
Statistical Model . . . . . . . . . . .
Assumptions of the Model . . . . . . . .

The Joint Density Function of the Parameters

in the Model . . . . . . . . . . . . . .

vi

10
11

13

17

19

21

21

22

26

28
36

37

39

39

42

49

 

 

ififi-T.‘ mr’_:—— --_——_-lw:'.F— —,-_-—-— v-_

IV

Obtaining Marginal Posterior Distributions
Via Gibbs Sampling . . . . . . . . . . . . . .

Data Augmentation . . . . . . . . . . . . . .
Example . . . . . . . . . . . . . . . .

Gibbs Sampling . . . . . . . . . . . . . . . .
The Application of Gibbs Sampling to The
Model of The Study . . . . . . . . . . . . . .

Obtaining Marginal Posterior Distributions of
Parameters of the Model in 3.1 . . . . . . . .

Finding the Conditional Distributions for the

Parameters of the Model in 3.1 . . . . . . . .
The Conditional Distribution of LG
Given 1, 1:2, {03.}, and Y . . . . . . . .

The Conditional Distribution of 0;
Given A, 0, 0.2, {Uj}, and Y . . .
The Conditional Distribution of A

Given HI}, {0%, and Y'. . . . . . . . .
The Condhtional Distribution of 12
Given U7} and Y’ . . . . . . . . . . .
The Condhtional Distribution of 03
Given 6 and Rig} . . . . . . . . . .
The Conditional Distribution of 6
Given a. and k5} . . . . . . . . . . . .

Empirical Applications of Gibbs Sampling . . .

Artificial Data Analysis . . . . . . . . . . .
Data Specifications . . . . . . . . . . .
Data Creation . . . . . . . . . . . . . .

Initial Estimates for Gibbs Sampling . .
Assessing the Heterogeneity of Variance .
Random Number Generation . . . . . . . .
Specifying the Criteria for Convergence .

Real Data Analysis . . . . . . . . . . .

vii

52

52
58

63

67

67

71

72

73

76
78
8O
82
85
86
87
90
97
99
100
101

102

 

 

 

VI Results . . . . . . . . . . . . . . . . . . . . . 106
Required Number of Iterations for Convergence . . 106
Artificial Data . . . . . . . . . . . . . . . . . 118
High School and Beyond Data . . . . . . . . . . . 119

Comparing Estimates of Variance Components . . . . 120
Estimates of of . . . . . . . . . . . . . . 120
Estimates of 6 . . . . . . . . . . . . . . . 122
Estimates of t2 . . . . . . . . . . . . . . . 123
Estimates of kfi} . . . . . . . . . . . . . . 123

Comparing Estimates of the Regression

Coefficients . . . . . . . . . . . . . . . . . . . 124
Estimates of A . . . . . . . . . . . . . . . 125
Estimates of an} . . . . . . . . . . . . . . 127

VII Discussion . . . . . . . . . . . . . . . . . . . . 182
Bayes and Empirical Bayes Estimation . . . . 182
Suggestions and Recommendations for
Future Research . . . . . . . . . . . . . . . 185

References . . . . . . . . . . . . . . . . . . . . . . 191

viii

 

s.__ -

‘11.;

._ ~—.—-~_:_._~u—»—

1_T .-—......- .v— .H-

 

 

CHAPTER 1

Introduction

Social and educational research is usually conducted in
natural settings such as schools or communities rather than in
a controlled experimental setting. Natural settings often
have a hierarchical structure. In education, for example,
students are nested within classrooms, and classrooms are
nested within schools. This hierarchical structure is
reflected in the data collected from these nested levels.
Observations at one level share common characteristics. In
teacher-effectiveness studies, for example, students within a
classroom share common characteristics of the teacher and
his/her teaching method.

Student learning may be viewed as a result of social
interaction within the classroom and the school system
(Bandura, 1977). Students are assigned to classrooms, and
they are taught in a planned and structured manner. Even
though individual students respond differently to the same
teaching' process, their’ responses 'will have commonality.
Therefore, observations on students learning cannot often be
assumed independent. However, independence is one of the
assumptions of many statistical procedures that researchers

often use in their analyses.

 

 

 

2

Ignoring the hierarchical structure of the data leads to
other statistical problems (Burstein, 1980; Cooley, Bond and
Mao, 1981; Cronbach and Webb, 1975; Knapp, 1977; and
Raudenbush and Bryk, 1988), such as aggregation bias.
Aggregation bias occurs when the relationship between two
variables varies (in magnitude and sometimes even in
direction) for different levels oftdata.analysis (Cooley, Bond
and Mao, 1981; Cronbach and Webb, 1975; Robinson, 1950). For
example, the relationship between socioeconomic status (SES)
and academic achievement is more likely to be higher in school
level data than it is in student level data. This difference
in the relationship is attributed to the effect of aggregation
of the. data in each school. In igeneral, school level
membership is related to SES; schools with higher levels of
SES are likely to have higher quality educational programs
which lead to higher achievement in those schools.

Information loss is another substantive problem that
results from ignoring the hierarchical structure in data
analysis. After the research or the evaluation design is
determined, a type of data analysis is usually chosen to
answer questions related to one level of data. If the
research questions were asked about students, it is likely
that they will be answered with a student-level analysis. If
questions were asked about classrooms, then it is likely that
they will be answered with a classroom-level analysis. When

the analysis is done at the classroom-level (i.e., using the

 

 

 

3
classroom as unit of analysis), it is very possible that some
important relationships between variables that were measured
on the student level will be ignored.

However, there has been recent development of statistical
methods which are appropriate for analyzing a hierarchically
structured data based on multi-stage (often two-stage) linear
models (Aitkin and Longford, 1986; Burstein, 1980; Goldstein,
1986; Lindley and Smith, 1972; Smith, 1973; Rao, 1972;
Rosenberg, 1973; Raudenbush and Bryk, 1986). General linear
mixed models or hierarchical linear models (HLM) are some of
the names used to describe the statistical models developed
for analyzing multilevel data.

The multi-stage methods share a common characteristic of
treating the within-unit regression parameters as outcomes for
the. between—unit ‘modelsu They differ, however, in the
procedure for estimating variance components. Some of these
procedures use the iterative generalized least squares method
suggested by Goldstein (1986) , some use the EM algorithm
method as applied by Raudenbush and Bryk (1986), and other use
the Fisher scoring approach (Longford, 1987). The estimates
produced by all these numerical procedures are maximum
likelihood point estimates. They are consistent, efficient,
asymptotically unbiased, and normally distributed (Harville,
1977).

A two stage HLMann.be presented in the following general

form. Within unit j, for j =i1, ,k units, we have

 

 

 

_ (1.1)
Yj" Xij + 53' ,

where Y]. is an njxl vector of observations, X]. is an njxp
matrix of within-unit predictor variables, Bj is a pxl
parameter vector which captures the relationships between Yj
and X]. within each unit, and ej is an njxl vector of random
errors which is assumed to have a multivariate normal
distribution with mean vector 0 and variance covariance matrix
10$; that is ej~N(o,Ia§). For Xj with rank (Xj)=p the
Ordinary Least Square (OLS) estimate for Bj is given by
§j=(XJ{Xj)-1XJ’-Yj with a sampling variance of
Var(DjIBj) = o§(XJI-Xj)‘1.

By allowing the within—unit regression model coefficient
to vary as a function of between—unit variables, we have a

multi—level model that describes hierarchically structured

data. Therefore, in the second (between-unit) stage, we have

B,- = ij + Uj , (1.2)
where W]. is a pxq matrix of known between—unit predictors and
y is a qxl vector of parameters that capture the effect of
the between-unit predictors W]. on the within-units parameters.
Uj is assumed to be N(O,T) , where T is the residual
variance—covariance matrix for [3]. after accounting for the

effects of Wj. Thus, we can have a more general

 

 

 

5
representation of the nudiilevel model by substituting 1.2

into 1.1 to obtain

Yj = XjoY + 53' . (1.3)

where Ej==AGU3+ej. One might estimate the effects in 1.3 by
using an ordinary least squares (OLS) multiple regression
approach to the estimation of multiple simultaneous effects.
However, one of the basic statistical assumptions for OLS is
the homogeneity of variance for Ej. Because Ej is equal to
X3U34-ej, this assumption can only be valid when either LG is
set to equal zero or if X3=X’for all j. Setting either [5
equal to zero or Xj=X , however, implies holding all the
within-unit regression model parameters as a constant across
all units and ignoring the hierarchical structure of the data.
This leads to the violation of the assumption of independence.

The total variance var(fij|wp of the within-unit
regression estimates is :made of 'their sampling ‘variance
var(fij|Bj) plus their residual variance var(Bj|W3). That is

Vaufijmj) = o§(X;Xj)-l + T. (1.4)

Many applications have been found which can utilize the
statistical model in 1.1 and 1.2. In longitudinal studies
(Bryk and Raudenbush, 1987; Laird and Ware, 1982; Strenio,
Weisberg and Bryk, 1983) the growth model parameters for each

subject become the multivariate outcomes for the between-unit

 

 

 

 

6
models. Subject characteristics (either observed or
controlled) can serve as predictors for these outcomes.

In research synthesis (Raudenbush and Bryk, 1985), a
simple measurement model for the observed effect size is set
equal to the true effect size plus some error in each study.
This measurement model is a special form of the within-unit
model in 1.1. The parameter fig Ibecomes the true effect size
in study j and .Kj becomes a scaler with a value of‘one. The
true effect sizes vary across studies as a function of known
study characteristics plus error in the form of the between—
unit model in 1.2.

A major benefit of fitting models 1.1 and 1.2 to hierar-
chical data is the possibility of getting empirical Bayes
estimates for the within-unit regression model parameters.
Empirical Bayes estimates are more stable and outperform the
classical estimates with regard to expected mean squared
error. Information from other groups can be used to get
improved estimates for the within-unit regression model
parameters. In providing a prediction equation, Braun, Jones,
Rubin, and Thayer (1983) presented a situation where there are
several predictors and few cases in subgroups of the
population ‘which. makes it hard. to obtain least squares
prediction equations. The objective of their study was to
provide separate predictive equations for the Graduate
Management Admission Test (GMAT) for the white and minority

students in each business school. Only 4% of 8500 students

 

 

 

 

 

7

sampled from 59 business schools were minority students. The
number of minority students within each school ranged from 2
to 10; this made it hard to obtain a prediction equation for
minority students for a particular school. With the HLM
analysis, information from other schools was borrowed to
obtain the predictive equation for the minority students
within each school.

Assuming that 03: and T are known, empirical Bayes

estimates of Bj and y are given by

[33- = 13-5,- + (1-Aj) ij" , (1.5)
where

Aj = T(Vj+T)'1 , (1.5)

Vj - o§(X§Xj)‘1 (1.7)
and

y. = {§w§[Var(fij|Wj)]’1Wj)‘1{)J;wj’-[Var(§j|wj)]'lfij} . (1.3)

The first component of B; in 1.5 is the OLS estimate of Bj
from the data within each unit weighted by its reliability
Aj. The second component, 7* is an empirical Bayes estimate
of the between-unit parameters, y from 1.2. More

comprehensive presentations of the HLM and different methods

of estimation of the parameters of the model are given by

 

 

 

 

 

8
Lindley and Smith (1972), Smith (1973), and Raudenbush (1984,
1988).

Although empirical Bayes methods have been around for a
long time, they became increasingly popular and were applied
to many types of problems in the early 1970s. Efron and
Morris (1972, 1973, 1975, 1977) discussed many applications of
empirical Bayes methods. Empirical Bayes estimates of the
parameters in hierarchical linear'models are.obtained.when the
true values for a; and T are known. In practice, however,
these values are rarely known to practitioners and often are
estimated from the data. Common estimates often used for
these two parameters are their maximum likelihood (ML)
estimates. Since estimates of the regression coefficients of
the HLM are functions of ML estimates of a? and T, they are
also considered to be ML estimates (Raudenbush, 1988).

In classical ANOVA procedures 03. and T represent the
within- and between-group variance components, respectively.
In balanced designs these variance components.can.be estimated
by solving a set of simultaneous linear equations obtained
from equating the observed mean squares (which are quadratic
forms of the observations) to their expected values (Searle,
1971). For unbalanced designs, there is more than one set of
quadratic forms to use for estimating these variance
components. Depending on the particular set used,
inconsistent estimates of the variance components may be

obtained (Searle, 1971).

 

 

 

 

9

Introducing HLM to the analysis of educational data can
improve the quality of the data analysis in education. It,
also helps applied statisticians obtain sensible answers to
questions where standard approaches seem to fail. However,
there continue to be obstacles related to HLM applications and
theoretical derivations.

In some applications where the assumption of homogeneity
of variance cannot be verified, drawing inferences on the
variance components becomes part of the research interest,
especially the variances of the within—unit residuals
(Leonard, 1975; Rao, 1970; Raudenbush and Bryk, 1987). The
motivation for making such inferences on the variances can be
attributed to the intrinsic interest in the variances
themselves and to the fact that they are part of the standard
errors of estimates for the regression.coefficientsu In other
applications, where we have small samples within each group
and.a small number of groups, maximum likelihood estimates for
the variance components 03- and T might become unstable.
Substituting these unstable estimates as true values for the
variance components will distort the regression effect

estimates.
Purpose of the Study

The purpose of this study'istto obtain marginal posterior

distributions for all effect parameters and variance

 

 

,__.. .___, _ ____..‘ _--__.~.
.- v—._..-- _- ._._ ,A._

 

 

 

 

estimate for the variance. Further,

10

components 111:3 linear mixed model with random intercepts
without relying on the assumption of homogeneity of the
residual variance 0; across groups. Non-informative priors
will be posed on the hyper-parameters that describe the
distributions for random parameters. The marginal posterior
distribution for each parameter will then be approximated
using the Gibbs sampling method (Tanner, 1993).

When research interest is focused on making inferences on
variances, finding the entire posterior distribution for each
of the variance components is more informative than having
onLy a point estimate of the variance. Having the entire

posterior distribution at her/his disposal will allow the

researcher to have the flexibility in choosing an appropriate

s/he can establish

probability intervals around the estimate.

Objectives of the Study

The idea of approximating the marginal posterior

distributions for all parameters in the model leads to the

question of how feasibly to do so. Therefore, the questions

of this study are:

1-

How do the posterior mean (Gibbs) estimates of the
parameters cufaa mixed model with heterogenous within-
group variance differ from their empirical Bayes

estimates which are based on the homogeneity of within—

 

 

 

 

 

 

 

 

11

group variance assumption? (Empirical Bayes estimates
for the within- and between—unit regression coefficients
were conditiOned on ML estimates of the variance
components with homogeniety of within-group variance).
How do inferences about regression coefficients change
when taking into account the uncertainty in estimating
the variance components?

Can a single estimate be used as typical value for all
the different residual variances kfii?

How precise are the posterior mean (Gibbs) estimates of

the residual variances {0%.} in estimating their single

typical value?

The last two questions focus on the parameters of the prior

distribution of the residual variances {03-}. Inferences can be

made on these two parameters from their approximated marginal

posterior distributions.

Bayes Solution

If reasonable priors are posed for all the unknown

parameters in the model, theoretically, the marginal posterior

distributions of the parameters of interest can be obtained by

integrating out the other nuisance parameters in the model.

Bayes estimates for these parameters can then be derived from

their respective marginal posterior distributions based on a

certain loss function.

 

 

 

 

 

 

12

From a Bayesian VieWpoint, estimating a parameter 6,
given its sufficient estimate a from the sample implies
selecting a decision function, say D, so that D(é) is a
predicted value of 0 when both the computed value 8 and the
posterior density function f(0|8) are known. Often, the
researcher predicts an experimental value of any random
variable to be its expected value; however, the median or the
mode can also be used as predicted values. For many cases, it
is desirable that the choice of the decision function
D should depend upon the loss function say L[6, D(é) ], that we

try to minimize (Hogg and Tanis, 1977) . For example, when the

loss function is given by:

L[e,D<é)] =[e —17(?))]2 (1.9)

then Bayes‘ solution that minimizes this loss function is

 

given by D(8) =E(6|é) . It is also possible to obtain Bayes'

probability confidence intervals for the estimates directly

 

 

from their respective marginal posterior distributions.

A problem associated with the application of Bayesian ;
methods to hierarchical models is the difficulty of the
mathematics involved in the derivation of the posterior
distributions (Lindley and Smith, 1972). The multiple
dimensions of the parameter space in the model make it
difficult to express some of the required equations in closed

form. Therefore, certain marginal posterior distributions

 

 

 

 

13
cannot be obtained analytically through direct evaluation of
an integral, and in many applications, integrating out the
nuisance parameters is a difficult job. Consequantlly,
empirical Bayes estimates are used. to approximate truly
Bayesian estimates.

In some situations we might find that the focus of the
research is centered around making inferences about a? and T
rather than the regression coefficients in the model.
Therefore, obtaining the marginal posterior distributions for
these parameters is more informative than obtaining their
single point estimates. Even when concern focuses on the
regression coefficients it will be useful to insure that
posterior uncertainty regarding these coefficients fully
reflects posterior uncertainty about the variance-covariance

components.

 

Importance of the Study

 

In educational settings, where we have several
classrooms or several schools, interest is mainly focused on
studying the factors that affect students' academic
achievement. Traditionally, student achievement is evaluated
by comparing schools or classrooms on their means on some
measure of academic achievement, usually students‘ test
scores. Often, important decisions such as school funding or

teacher promotion are based on this evaluation. The fairness

 

 

 

 

 

 

 

 

 

 

14
of these decisions depends on whether the test score means in
different schools or classrooms truly represent students'
academic achievement. Many factors that could affect these
means can be controlled statistically through model
specification. However, it is still possible that for a
slightly higher test score mean, one school will get more
funds than others or a specific teacher may get promoted. It
is also possible that the school that got extra funds or the
teacher who got promoted produced higher variability in their
students' academic achievement.

A successful educational program or an effective teacher
should arguably strive not only to increase and facilitate
students' academic achievement, but also to reduce the gap
between students' learning and produce equality in their
achievement (Bloom, 1984). Therefore, the evaluation process
should be based not only on average achievement, but also on
the equity of the achievement. Consequently, any improVement
on the methods of estimating these two criteria certainly
helps the decision maker in his/her evaluation of the
educational system.

Another example is drawn from research syntheses
(Raudenbush and Bryk, 1985; Rubin, 1981) where a simple model
for the observed effect size is set equal to unknown true
effect size plus some error for each study. The objective is
to obtain an efficient estimate of each study effect size as

well as its precision. It is also common in some cases to

 

 

 

 

 

 

15

have a number of groups with a small number of observations
for each group which were obtained under the same conditions.
For example, in longitudinal studies, the growth of subjects
from the same age can.be observed.over several points of time.
One might become interested in studying the growth rate of a
particular subject in comparison to other subjects. The
replicate observations on the subject can then be used for
variance estimation.

In all of above examples, there are two sets of
parameters (means and variances) with several parameters of
the same type in each set. Within each set, these parameters
are related by common circumstances; schools or classrooms in
the first example, studies in the second example, and the
individual subject in the last example. Within the Bayesian
framework, it is logical to assume that these many parameters
that share common circumstances would have a common parametric
prior distribution. This distribution summarizes the informa-
tion about these parameters prior to data collection. This
allows the researcher to use such prior knowledge about these
parameters to get their improved estimates.

With the help of modern computers and developments in
simulation'theory, research.interest.is being directed towards
approximating the posterior distributions of the parameters in
the model being investigated and obtaining their Bayesian
estimates (Tanner and Wong, 1987). In his work on

hierarchical linear models, Seltzer (1988) adopted the data

 

 

 

 

 

l6

augmentation method to obtain Bayesian estimates for variance
components. He also investigated Bayes' estimates for the
within-school regression coefficients when they are assumed to
have a t-distribution (see also Seltzer, 1993). Fotiu (1989)
also applied the method of data augmentation in estimating
the joint posterior distributions for the effect parameters
(within- and between-units) when analyzing many groups. Both
Seltzer and Fotiu approximated the posterior distribution of

2

O .

3 under the assumption that it is equal to (9 across all

groups. This assumption is usually made for convenience,
rather than because it necessarily holds true. By making this
assumption, the researcher can pool all the observations from
all the groups to get a large-sample point estimate of 02.
However, when the marginal posterior distributions of all the
parameters in the model can be approximated, this assumption
becomes no longer essential to the analysis. In fact,
sometimes the main objectives of the statistical analysis is
to estimate the posterior distribution for each 0?.

In conclusion, the objective of this study is to
investigate the use of a mixed linear model with heterogenous
variances {03-} across groups. A linear regression model with
a random intercept will be used to represent the relationship
between the criterion variable and the predictors in each
group. By obtaining the marginal posterior distribution for
all the parameters in the model, the practitioner can obtain

Bayes estimates for the effect parameters as well as variance-

 

 

 

 

 

 

17
covariance components. Based on the theoretical formulation
of the model a number of simulated data sets will be generated
to verify the accuracy of the estimation under different
parameter specifications. These specifications reflect the
degree of heterogeneity of variance, the different numbers of
group sizes, and the complexity of the regression model for
each group. The estimation process will be applied to a
national sample of US high schools where math achievement has
been studied as a function of school and student
characteristics. One important feature of this data set is
the heterogeneity of the residual variance across schools.
Radenbush and Bryk (1987) investigated the heterogeneity of
variance in this data set through the application of

hierarchical linear models.

Organization of the Study

Chapter 2 will present a statement of the problem and
review of literature conducted on the topic. That chapter
will also highlight some of the problems in variance
estimation, and its relation to the HLM. In chapter 3, a
general linear mixed model with a random intercept will be
presented to represent the several applications presented in
this study. In addition, a comprehensive description of the
assumptions associated with the model will be provided within

the Bayesian framework. Chapter 4 will explain the Gibbs

 

 

 

 

 

l8
sampler and its application in achieving the objectives of the
study. An example will be provided to illustrate the
iteration process of Gibbs sampling. Chapter 5 will present
'the derivation of the conditional distributions for the
parameters needed in the iteration process. It will also
provide a: description to the computational steps for the
iteration process presented in chapter 4. Chapter 6 will
presents the results of the small simulation study as well as
the results of an analysis of mathematic achievement data from
US high schools. Discussion of the results for the generated
data and the US high school data, as well as conclusions, will

be presented in chapter 7.

 

 

 

 

 

CHAPTER 2

Statement of the Problem

Variance estimation plays an important role in almost
every kind of quantitative research. A basic requirement of
nearly all forms of analysis is that a measure of precision be
provided for each estimate derived from the data. The most
commonly used measure of precision of an estimator is the

reciprocal of its variance. The sample variance is used to

 

estimate the precision of the mean or other location
estimates. It is also used to check for gross errors
affecting' a single observation. Further, it is used to
provide an estimate of a variance component such as a pooled
within sample variance in several kinds of statistical
procedures.
In HLM analysis, the two variance—covariance components
2

0' .

J and T are used in obtaining empirical Bayes estimates for

the effect parameters in equations 1.1 and 1.2. Compared to
the Bayesian approach, the empirical Bayes approach does not
require the specification of prior distributions for a; and
T. The two-level hierarchical linear model and its
assumptions for the empirical Bayes approach can be restated

as

19

 

 

 

20

Y]. = Xjflj + 6]. , where 6- ~ N(0,Io§-) , (2-1)

for the first level, and

B]. = ij + Uj , where Uj ~ N(O,T) , and 1 ~ N(G),l"), (2.2)

for the second level. The joint posterior distribution for
the effect parameters Dj and y is conditioned on the values
of a; and T (Lindley and Smith, 1972). The joint density

function f1(fij,le3,o§,T) for this posterior distribution can

be expressed as
f1(i3j,lej,o§-,T) o< f2(Yj|pj,o§) f3((3j|o§, y,T) f4(y), (2.3)

where f2(yj|(3j,o§) is the likelihood function of the data,
f3(Bj|o§-, y,T) represents the conditional prior density
function of (33., and f4(y) is the density function for a
noninformative prior for y. The values for a? and T are
often estimated. by one of several methods of numerical
estimation. One common method is the EM algorithm (Dempster,
Laird and Rubin, 1977) which.produces their maximum likelihood
estimates. These estimates are asymptotically normally

distributed, unbiased and efficient estimates.

 

 

 

21

Problems Associated with the Estimation of 02 and T

In some applications, the sample size within each group
as well as the number of the groups might not be large enough
to allow the specific estimates of a; and T to obtain the
asymptotic properties. Nevertheless, practitioners often use
maximum likelihood estimates of a; and T to represent their
true values to obtain the conditional posterior distributions

for the effect parameters.

Problems Associated with Estimating T

In 1.2 of the HLM, T represents the residual variance-
covariance matrix for the within-group parameters Dj across
all the groups. The shape of the sampling distribution of
this variance—covariance matrix depends on the number of
groups in the study. Consider the model investigated in this
study, for example, where only the intercept is considered
random. The matrix T , then, becomes scalar 12, which is the
variance of the intercept across all the groups. When the
number of groups is relatively small, say kle, it is
possible that the distribution of an estimate of t2 becomes
highly skewed. Using the mode of that distribution (maximum
likelihood estimate) as an estimate for 12 might not be a good
representation of its true value.

In a validation study conducted on eight law schools,

Rubin (1983) pointed out the problem of using maximum

22
likelihood estimates for 12 when there are small number of
schools in the study. He stated that

"The specific problem here is that the likelihood

of 12 is not at all symmetric about the maximum

likelihood estimate, and thus this estimate is not

representative of reasonable values for 12.

Integrating over :2 in such cases is a much more

reasonable way to summarize evidence than to fix 12

at some value" (Rubin, 1983, p. 15).

Notice that the estimated value of 12 is being
substituted in 1.4 to 1.8 to obtain the precision weight for
the empirical Bayes estimates for the regression coefficient
of the HLM model. Consequently, invalid estimation of this
parameter causes the estimates of the regression coefficients

of the HLM model in 1.5 and 1.8 and their reliability

estimates in 1.6 to be distorted.

Problems Associated with an Invalid Assumption of Homogeneity
of Variance

When maximum likelihood estimates are obtained for
variance components, {0?' for ;i= 1, ,k7 are often assumed
homogeneous. Researchers often assume homogeneity of
variances across all groups out of convenience rather than
conviction. This assumption allows the pooling of the
observations from all the groups to get one large sample
maximum likelihood estimate of 02 for all the groups. Berlin
(1984) stated that

"large data sets are rarely homogeneous in their

precision and that usual statistical analysis fails

if it does not take such differences into
consideration.“ (p. 209).

23
When the number of observations in each group is large enough
to allow for consistent maximum likelihood estimation, Mason,
Wong and Entwisle (1984) advised that a separate variance
estimate should be obtained for each group.

However, when the number of observations in each group is
not large enough to allow for consistent maximum likelihood
estimation and the parameter variances are not equal, pooling
the observations from all the groups to obtain a single ML
estimate for 02 might give an inaccurate estimate of the
variance. Assuming homogeneity of within—group variance
Bassiri (1988) showed that the within—group variance estimate
becomes unstable when the groups' sample sizes are relatively
small. This could lead to invalid estimates of the regression
coefficients Bj and their precision estimates (see equations
1.4 and 1.5). Consequently, inferences about or confidence
intervals for the regression effects would be invalid.

Equation 1.4 shows when the true values of o; and T were
being used, part of the variability in the estimated variance
of the regression effect Bj stems from the variability inX
across groups. Another part stems from the variability in
og's when no homogeneity of variance assumption is being made.
By assuming homogeneity of variance,o§= 02, we are ignoring
part of the variability in the variances of the regression
effects. The part being ignored is attributable to the

variability in the residual variances.

 

 

 

24
For illustrative purposes, consider the case where A; is
assumed equal to X across the k groups for j=14 m, k. Thus,
the within—unit regression coefficient estimates become
[’3‘]. = (X’X)‘1X’Yj. When homogeneity of variance holds, the

variance—covariance matrix for these estimates can be derived

from 1.4 as
Var(i3j|Wj) = 02(X/X)‘l +T. (2.4)

Equation 2.4 implies that the within—unit effects are
estimated with the same precision across all the k groups.
However, if there is a considerable variability in the
residual variances across groups, the estimation of this
precision becomes invalid and those estimates of the
regression effects no longer have the same precision.

In a classical analysis, when. researchers study any
treatment effect, they compare the means of the criterion
variable from different treatments. One of the assumptions
they often make is homogeneity of variance. The justification
of this assumption is to get a single pooled variance estimate
for the within treatment variance to represent the error term
in the analysis so that treatment effect can be tested. Bryk
and Raudenbush (1988) showed that heterogeneity of variance in
the traditional analysis of treatment effects can be seen as
evidence of an interaction between treatment and subject—
specific characteristics. They warned researchers against the

problem of obtaining biased estimates for treatment effect

 

 

 

25
when ignoring heterogenous variances in the classical
analysis.
Consider the classical case, often referred to as the
Behrens-Fisher case, of testing the hypothesis pa "ih>= 0 with
unequal population variances, oisoi. Using §;-X; as an

estimate for pa—pb, the sampling distribution of the statistic

z: (Ea—Eb) - (pa-uh)
2 2 (2.5)
0 U
(—a)+ —‘3>
na H):
is 1V(0,1). The given hypothesis cannot be tested without

knowing o: and oi. When the sample sizes ha and nb are both
large enough, we can substitute the unbiased estimates of the
variances, s: and SE, for the corresponding parameters, and
the resulting statistic has a sampling distribution
approximated by N(O,1). In practice, however, when 12a and nb
are both relatively small and noticeably unequal, the
resulting statistic

Z: = (Xa_Xb) _ (pa-p19)
S2 32 (2-6)
(—a)+(—5)
a nb

would have what is known by Behrens—Fisher distribution
(Winer, 1971). Therefore, when there are both heterogeneity
of variance and a noticeable difference in the sample sizes
between treatments, the homogeneity of variances assumption

becomes questionable.

 

26

The Analysis of Many Group Variances a;

A common situation where many variances need to be
estimated arises when there is a heterogeneity of variance in
linear models. Rao (1970) introduced a method called Minimum
Norm.Quadratic Unbiased Estimation (MINQUE) for estimating the
residual variances in linear models when these variances are
found to be heterogeneous. MINQUE estimates are a linear
function of the distinct variances under certain conditions.
These conditions are related to the choice of the matrix that
creates a quadratic form in the outcome variable. If some of
the variances are the same, their corresponding coefficients
in the linear function can be chosen to be the same. The
variances estimated using the MINQUE method can be used to
obtain improved estimates of the coefficients of the linear
model. Further, they can be used in obtaining the estimate of
the precision of the simple least squares estimator of the
regression coefficients or any linear function of these
coefficients.

Cook and Weisberg (1983) modeled the variability in
variances as a function of some explanatory variables. The n
residuals from a linear regression model are assumed to have
a multivariate normal distribution, with mean zero and
variance—covariance matrix 02W, where W is a diagonal matrix
with all diagonal elements wq>0 for i=14. ,n. Estimating the

residual variances implies estimating the elements of W.

 

27
The variability in up, which represents the variability
in the variances, is expressed as a function of a 1 x q row
vector of known explanatory variables zi=(zfi) forj=1, ,q
and q><1 vector of unknown parameters A as in the following

two models

q
%=en{ahflfi)l (L7)
J:

or

q
w.= exp(§:Ajlog(zfi)). (2.8)
j=1

1

Note that the variables zfi can have negative values in the
first model but not in the second one. The predicted values
of hp, from the above models, can then be used to estimate the
residual variances from their variance—covariance matrix 02W.

Cox and Solomon (1986) have also dealt with the issue of
estimation of many variances. They provided examples
illustrating the problem of estimating many variances from
small samples. One of the examples covers the situation where
there is a systematic difference in variance between samples.
Observations within a sample are assumed to come from. a
N(pi,o§) population where 1:14. ,k. The ith variance 0? is
a function either of an explanatory variable zi,
characterizing the ith population, or of pi. They adopted
similar versions of the models in 2.7 and 2.8 to represent the
systematic change in the variances.

They also presented another example for the case where

 

28
the changes in the variances are considered random. They
consider it a complement to the above case (where changes in
the variances are systematic) and called it the
"overdispersion model", (Cox and Solomon, 1986, p. 544). In
this case, different population variances are considered to be
independent unobserved values of a random variable 02. This

variable has an inverse gamma distribution with density
1 -1
a — Va

(ivoof) (032v. exp —O—:° {F [2t°]} , (2.9)

where v, represents the degrees of freedom, and of some

 

constant.

Given the inverse Gamma density in 2.9 and considering
the ith. population :mean pi as ea nuisance parameter, the
estimation of the ith population variance 0% is based on the
marginal likelihood of the ith sample. The underlying idea
of the above model suggested by Cox and Solomon (1986)

coincides very closely with Bayesian thinking.

Bayesian Approach

Lindley (1971) used the Bayesian approach in the
estimation of many different means and variances. In the case
where both means and variances were unknown but have
exchangeable distributions, he used the conjugate priors for
these distributions to derive the joint posterior density for
the means and variances. The joint density of all the

parameters was based on a three—stage model of hierarchy. The

 

 

~§~As‘)—_ 1: .. A—

29

first stage describes the data, and its likelihood function
given the means and the variances. That is, given {ufl and
{0;} the data Kg are independent and normally distributed with
.E(Yfi) = uj and var(rfi) = a? for i = 1, ,nj; j = 1, ,k. In the
second stage, the means in; and the variances is? are assumed
to be independent with two different exchangeable
distributions. The exchangeable distribution for the means
{uj} is assumed to be normal with mean 6 and variance 1. The
exchangeable distribution for the variances hfi} is an inverse
chi—square, where the variable voof/ofi is distributed as a x2
with v° degrees of freedom and of as a typical value for 0?.
The third stage describes the prior knowledge about the
parameters in the second stage. For the parameters in the
normal distribution a conjugate prior distribution for the
mean 6 with vague prior knowledge produces locally uniform
prior. Similar to 0%, the variance T is assumed to have an
inverse chi—square where vrof/t distributed as x2 with vt
degrees of freedom and of as a typical value for I. Since a?
is distributed as an inverse chi—square, a conjugate prior for
of given V0 is a function of x2 with r and A as two constants
describing the distribution. The parameter V0 is assumed to
have a uniform prior on (O,W).

The joint posterior density for the means in} and the
variances {0? was derived by integrating out all other
parameters (which are 6, T,V° and of) from the overall joint

posterior density function. The derived estimates were based

 

 

30

on a simple model where there were several groups with
different means and variances. If a more complex linear model
with several predictors were to be used the estimation process
then becomes more complex. The added coefficients for the
predictors in the linear model make integrating the joint
density function of the parameters analytically, extremely
difficult.

Leonard (1975) provided a procedure for modeling the
variability in the means and the variances from several normal
populations using the Bayesian approach. His model for the
means, when the variances are known, represents a general case
of the one presented by Lindley (1971). When populations
means are known, Leonard expressed the log—transformed
variances as a function of some explanatory variables using a
linear model similar to the one for the means. The log-trans—
formed variances were assumed to have an asymptotic normal
exchangeable distribution similar to the one for the means.
Expressing the variances or any transformation of them as
function of other variables is similar to what Cook and
Weisberg (1983) did as shown in equation 2.7. The major
difference here is that Cook and Weisberg consider Aj in 2.7
to be fixed unknown parameters that need to be estimated,
while Leonard proposed the use of prior knowledge about a
similar parameter for the linear model of the log—transformed
variances.

The resulting estimates of the means and variances are

 

 

 

31

shrinkage estimates. They represent a weighted average of the
standard estimate from the sample plus a weighted average of
all the samples estimates. The weights are related to the
reliability of the standard estimate from the sample. The
shrinkage estimates of the variances are found by taking the
exponential of the shrinkage estimates of the log—transformed
variances.

When both means and variances are unknown, Leonard
proposed an iterative method for estimating the joint
posterior mode vectors of the means and the variances. He
substitutes the shrinkage variance estimates in the expression
for the means to obtain their new estimates. The new
estimates of the means are then used to find new shrinkage
estimates for the variances.

Raudenbush and Bryk (1987) tackled the issue of
estimating variances from many groups as a special case of a
more general problem of estimating the parameters of a two—
stage HLM. They provided two different methods for obtaining
empirical Bayes estimates for the variances. The
conceptualization of their models is similar to those
presented by Lindley (1971) and Leonard (1975). In the first
method, estimates were derived using the exact distribution
(chi—square) of the variance estimate within each group. That
is,

sf ~ oixifl/(ni—l) . (2.10)

 

 

 

32

Assuming that the variances have an exchangeable prior
distribution, the variable vgof/oi has a chi-square
distribution with vo degrees of freedom and 03 represents the
typical value for the variances prior to observing the data.
This distribution is the natural conjugate prior for the one
in 2.10. Thus, the variance 0% is distributed as vqofxfij,
where xfij denotes a ‘variable with an inverse chi-square
distribution.

In the second method, however, estimates were derived
using the asymptotic normal approximation to the sampling
distribution of the logarithmic transformed variance estimate

within each group. For the first stage of the HLM
di=6i+ei, where ei~.N(0,(ni—1)”/2) for i=1, ,k, (2-11)

where cQ=-%[log(sf)-ci], Ci is a bias correction, and

6i=-%log(o§). For the second stage of the HLM

5. = W/y + U1, where U1” N(O,‘;1) for i=1, ,k , (2.12)

.1

 

where W1 is a Mx1 vector of known predictors, y is a Mx1
vector of effect parameters, and [Q is a random error normally
distributed with mean equals zero and variance equal ¥;. The
use of the normal approximation to the sampling distributions
of the variances permits the use of the above formulation and
allows for the hypothesis testing associated with it.
Estimates from both methods (exact distribution and

normal approximation) are shrinkage estimates. They are

 

 

 

 

33

weighted averages of two components. For the exact
distribution, the conditional posterior mean estimate for the
variance is given by
32 = W , (2.13)
ni+ v0— 3
The first component is the usual estimate of the variance
based on the observations within each group weighted
proportionally by its degrees of freedom. The second is the
typical value of of the variances in their prior distribution
weighted by the concentration of these variances va around
that typical value.

For the normal theory, however, the empirical Bayes

estimate for the logarithmic transformed variance is given by

(2.14)

~

51 = Xidi + (1_Xi)Wi/T r
where Ai=(ni-1)/(ni+9-1), and 9 is estimated numerically.

The first component is the ordinary estimate of the
logarithmic transformation of the variance estimate within a
group, weighted proportionally by its degrees of freedom. The
second component is a predicted value based on information
from all groups, weighted proportionally by the concentration
of the parameters estimated by the first component around that
predicted value.

The resulting estimates fronlboth.methods are conditional
estimates. In the method where exact sampling distributions

of the variance estimates are being used, the empirical Bayes

 

 

 

 

34
estimates for the variances are conditioned on knowing the
true 'values for the parameters of and v0. In .practice,
however, the true values of these two parameters are usually
unknown. To obtain the shrinkage estimates of the variances
in 2.13, of and v0 have to be estimated. Based on the EM
algorithm approach by (Dempster, Laird and Rubin, 1977) for
numerical estimation, Kasim (1986) developed a procedure for
estimating of and v“ The complete data for estimating 03 and
v0 are {S§L{o?‘ for j =1” ,k, where only {3;} have been
observed. The parameters of and v, are estimated by
maximizing their likelihood function .Zflo?)v°,ofi (M-step),
assuming that (a? have been observed. The sufficient
statistics for estimating of and v0 are functions of (0?.
Since {0? cannot be observed, the posterior expectations (E-
step) of these sufficient statistics are used to estimate of
and Va. Finding the posterior expectations of the sufficient
statistics, however, depends on knowing the values of of and
v0 as well as the data {sfih This dependency between the
posterior expectations and of and v0 is used in an iterative
process to estimate of and v0. Therefore, given the observed
data {Si} and initial estimates of of and vo, the values of
sufficient statistics are estimated by their posterior
expectations. The estimated sufficient statistics are then
used in maximizing ldofilvucf) to obtain new estimates of of
and v0. Going back and forth between the expectation step and

the maximization step, we get reasonable estimates of 03 and

 

35

v In fact, each step of this iteration process increases

1({09lv°,03), which provides us with better estimates of of and
v0 than those from the step before. The iteration process is
terminated when the absolute change in the values of of and
v0, between any two steps, is sufficiently small.

In the method where the normal approximation to the
sampling distribution of the variance estimates is being used,
the empirical Bayes estimates of the variances are conditioned
on the residual variance V in the second stage model (see
2.15). The true value of this parameter is unknown, and it
must be estimated. Similar to the way that of and v. are
estimated, the maximum likelihood estimate if is obtained
numerically via the EM algorithm.

The work of Lindley (1971), Leonard (1975) and Raudenbush
and Bryk (1987) on estimating variances of many groups share
the characteristic of borrowing information from other groups
to get an improved estimate of the variance for a particular
group. As an application of this idea Singh and Sedransk
(1988) provided an example for obtaining improved estimates of
strata variances when strata sample sizes are small. The
improved variance estimator uses data borrowed from other
strata, initially thought to be similar. The resulting
estimate is a shrinkage variance estimate for the stratum,

which is similar to the one in 2.14.

36

Comments on the Analysis of Many Variances

There are still several theoretical and practical
weakness in estimating and investigating the variability in
many variances. Empirical Bayes estimates of the variances,
which were suggested by' Raudenbush and Bryk (1987), for
example, are not adjusted for the uncertainty in estimating
the conditioning parameters of and v6 in the exact method, and
V in the normal approximation method. This is similar to the
problem (Tanner and Wong, 1987) of obtaining empirical Bayes
estimate for the effect parameters in HLM without adjusting
these estimates for the uncertainty in estimating the variance
components.

The focus on obtaining conditional point estimates for
the variances in the model rather than deriving the entire
posterior distribution for each variance can cause some
problems. In some cases, the misrepresentation of the true
values of the conditioning (given) parameters will lead to
obtaining inaccurate conditional estimates of the variances.
For example, in the case of the exact theory of Raudenbush and
Bryk (1987), the variability of the group variances is
inversely proportional to the conditioning parameter v° plus
some constant. Therefore, if there is a clear evidence of
heterogeneity of variance then the likelihood of V0 is not at
all symmetric around its maximum likelihood estimate.

Consequently, this estimate is an inaccurate representation

 

 

 

37

for the true value of v“ Conditioning the variance estimates
for the different groups on the estimated ‘value of this
parameters will distort these variance estimates. This
problem is similar to the one presented earlier by Rubin
(1983) for estimating 12 when the number of groups is small.

Normal approximation for the distribution of the
estimated variances is another crucial point in the analysis
of many variances. Bartlett and Kendall (1946) reported that
the normal approximation is sufficiently accurate when the
sample size in each group is at least equal to 10. For
smaller samples, the normal approximation to the distribution
of the estimated variances ‘might be less accurate. The
inaccuracy in the normal approximation can lead to
inappropriate representation of the model. A severe loss of
efficiency in estimation can occur when sample sizes are less

than 5 (Cox and Solomon, 1986).

Adjusting for The Uncertainty in Estimating a; and T

Standard methods do not exist for adjusting the
conditional posterior distribution for the effect parameters
to reflect increased variability due to the uncertainty of
estimating a? and T (Tanner and Wong, 1987). The estimates
of the regression coefficients given in (1.5) to (1.8) are
conditioned on knowing the true values of o? and T . In

practice, however, these parameters are usually unknown, and

 

 

 

38

their ML estimates are often used to allow for empirical Bayes
estimation for the regression coefficients. Therefore, part
of the variability in the empirical Bayes estimates for the
regression coefficients should be attributed to the
uncertainty in estimating 03 and T. Having the marginal
posterior distributions for the regression coefficients will
facilitate obtaining their estimates without conditioning on
a? and TH Consequently, the estimates will be correctly
adjusted for the uncertainty of estimating the variance
components a? and T.

In summary, the problem in getting empirical Bayes
estimates for the effect parameters in multilevel models is
centered on getting appropriate estimates for a? and T. When
the sample sizes and the number of groups are large enough,
maximum likelihood estimates of a; and T are a good solution
to the problem of getting those empirical Bayes estimates.
However, when there is heterogeneity of variances and the
number of observations within each group as well as the number
of groups are not large enough to allow asymptotic estimation,
maximum likelihood estimates of 02 and T may not be
appropriate; it is essential that one should move toward being

fully Bayesian in this analysis.

 

 

 

 

 

CHAPTER 3

Model Specification

This chapter provides a description of the statistical
model used in this study and its underlying theory. The basic
assumptions for this model are provided within the Bayesian
framework. The joint posterior distribution for the
parameters in the model is also presented.

A common characteristic often shared by evaluation
studies that involve social interventions, like a new
educational program, drug rehabilitation, or health care
promotion, is that the data collection process is based on
more than one level of hierarchy. The most common case is
where there are measures for program or group characteristics
and other measures for characteristics of subjects who are
nested within the groups or programs. Any model set up for
investigating relationships between measures from this kind of
data needs to acknowledge the hierarchical structure of the

data.

statistical Model

This section presents a general linear model that
captures the hierarchical structure of nested data. The model

and its assumptions can be described in three stages. The

39

 

 

 

 

 

 

 

40
first stage describes data obtained from k groups, where

j=1,-~, k. Given the set of parameters A, U. and 03:, let Yj

J!
be a vector of nj independent observations from a normal
distribution with a mean ZjA + ljUj and a variance Ijog. A

linear model with random intercept and fixed effect parameters

for both group and subject level variables can be specified as

where

2

Uj is a random error for the intercept,

F

 

 

Y.-
71
l le Wj lel .. lep '
z]. = s s 2 s s , and A = «(Q
1 W71 qu anjl anjp pl (3.2)
Lfipi

The part [27]., qu] of the matrix Zj represents q known group

level variables, the part [X1- ] represents p known fixed

7'1 Xijp
effect subject level variables, and A is p+q+1 vector of
parameters that capture the effect of group and subject level
variables on the outcome variable Yj.

In the second stage both the intercept y°+Uj and the
residual variance 0; are assumed to be independent and vary
randomly across the k groups. Given A and t2, the random

part of the intercept U]. is distributed as normal with mean

zero and variance 12. Given 0.2 as a typical value of 03, the

 

 

 

 

 

 

 

41

02

variable 67'; is distributed as x2,%) with % degrees of
freedom.

The last stage represents vague prior knowledge about the
hyper-parameters A, 12, 6, and of for the two distributions
presented in the second stage. The prior for each one of
these four parameters is assumed proportional to some

constant.

The model in 3.1 and its associated assumptions can be

summarized as follows:

1- The observed data

Given A, Hi and c@ the data I} is distributed as

 

2 / 2

2- The parameters of exchangeable distributions
Here a; and 09 are assumed independent,

Ule,1:2 ~ N(0,1:2) , (3.4)
and

0.2 2

" xii) ' (3.5)

60] 9

3- The prior knowledge about the hyper-parameters.

The joint prior distribution p(A,tz,6,o§) can.be written

as

._....._. ____.‘—-—

 

 

 

1». ...__-._IH~.—-_..—-——.-q.~.- —— _... -. —»— “illuililita "'—

42

p(),.2,e,e§) « p(1)p(r2)p<e,of), where
pm «X c,
Lfirz) « Cb
p<6,of) .. p<of e)p<e) , where
_p(of 6) a ($3! and
p(6) cc Ce .

 

(3.6)

 

The hierarchical structure of the models in (3.1) and
(3.2) fits nicely within the Bayesian approach. The next

section provides more discussion about the specification of

 

the parametric forms for the.distributions of the exchangeable
parameters LG and a? and.the hyper-parameters A,‘fi% 6, and.of.
These specifications are required to obtain Bayes estimates

for the parameters of interest.

 

Assumptions of the Model

There are three sets of assumptions associated with the
model in 3.1. The first assumption is about the data. Given
A,Cg and 03, we assume that the data Y} are normal with mean
equal to Z3A-+]jLG and variance equal to 1}o§.

In the second stage of the model, the intercept yo+CG and
the residual variance 0? are allowed to vary randomly across
the groups. The parametric forms of the distributions for 03
and a? are based on the hierarchical structure of the model.
More explicitly, there are two sets of parameters with several

parameters of the same. type: in each. set (i.e., several

 

 

 

 

43
intercepts and several variances) . The parameters within each
set usually describe a common variable that is shared by all
the groups. In the example of evaluating educational program
in different schools this common ‘variable could. be the
effectiveness of an educational program adopted by each school
within a district. For the meta-analysis example, the common
variable is the effect size that varies in value across all
studies used in meta—analysis. Within each set of parameters

(intercepts or residual variances), one can think of the value

 

of each parameter in the set as a single realization of a
random variable that has a common distribution with the rest
of the parameters in the set. When this conceptualization is
translated to a probability structure, each parameter can be

viewed as a random variable having an exchangeable

 

distribution with the rest of the parameters in the set.
Another way of saying this is that these parameters are
independently identically distributed random variables with a
prior distribution yet to be specified.

Under the Bayesian approach, the prior knowledge about
the parameters can be used to specify the parametric form of
their prior distributions. Bretthorst (1988) provided general
guidelines for choosing prior distributions. He stated that

"There are two questions one may consider to help

in this. First, one should ask ‘Are the parameters

logically connected?‘ That is, if we gain

additional information about one of the parameters
does it change the estimates we would make about

the others? If the answer is yes, then the

parameters are not logically independent. It will
be useful to find a representation where the

 

 

44
parameters are independent.
Another useful question is ‘What are the
invariances that the prior probability must obey?‘
That. is, what 'transformation.‘would. convert. the
present problem into one where we have the state of
prior knowledge." (Bretthorst, 1988, p.183).
The first question deals with. issue of independence between
the parameters for which.we seek.prior distributions. For the
model in 3.1 we assumed CG and a? independent. It implies
that knowing the residual variance in each group does not
depend on knowing the random error associated with intercept.
The second question deals with the choice of conjugate prior
distribution for the parameter of interest. Therefore, given
A and 12, we assume that the conjugate prior for the random
part of the intercept.LG is normal with mean equal to zero and
variance equal to wk. Further, given 6 and of,'we assume that
(fi/(Bofi) has a chi-square distribution with -% degrees of
freedom. The hyper-parameter of describes the value that one
expects the residual variance 0? to have, prior to observing
the data. Therefore, 03- is distributed as gxé), where xé)
is an inverse chi-square with %-degrees of freedom. WhileO
is a function of the degrees of freedom, it can also be

thought of as a measure of variability of the variances kfi}.

The probability density function p(o§|0,o§) is found to be

 

2 2 (g5)E 2 —(-%+1) ‘03 (3 7)
p(0j|9.0.) = ———1——- (oj) exp 2 . .
IKE) 260j

The above density function is similar to the one given by

Lindley (1965), with.a minor difference of 6=1¢Hu. As Lindley

 

 

 

 

 

 

 

45

described it,

"For all v°>0 the above density function increases
from zero at 03‘- = O to a maximum ate?- = 03v,/(v,+2)
and then tends to zero as og-ow. The density has
therefore the same general shape for all degrees of
freedom as the chi-square distribution itself has
for degrees of freedom in excess of two." (Lindley,

1965, p. 28).

In his interpretation of the above density as a prior density
for 0%, Lindley said:

"To take it as a prior density for o? is equivalent
to saying that values that are very large or very
small are improbable, the most probable value is
o§==o§vJ%v,+2) and, because the decrease from the
maximum as oj increases2 is less than the
corresponding decrease as 05 diminishes, the values
above the maximum are more probable than those a
similar distance below the maximum." (Lindley, 1965

p. 28).

 

The above interpretation to the prior distribution of og‘makes

 

choosing such a prior more realistic than assuming equal
likelihood for all values of 03- (i.e., a locally uniform
prior).

The mean and the variance of 0% can be easily found by
realizing that the density in 3.7 is an inverse Gamma function

with or = 1/26 and [3 = 26/03. Therefore,

 

 

 

 

 

 

46
2
2 2 _ l = O.
E(°jle’°°) _. [3(a — 1) 1 — 20 ’ and 4 (3.8,
Var(o§|9,o§) = 1 = 200°
32(a - 1)2(a - 2) (1 - 28)2(1 — 48)

The results in 3.8 reveal that 6 has to be less than .25 for
the variance of a? to be defined» As 6 gets smaller, the mean
and the variance are approximated to of and 2603.

In many analyses cfiz is used as a measure of precision
for each estimate derived from the data. From 3.7, the prior
distributionofcr}2 is Gamma with a = 1/26 and B = 26/03. Its

mean and variance are found to be

E(o}2|6,o§) = ab = 0:2 , and

(3.9)
Var(o§2|6,o§) = «[32 = 260:4.

Another measure that describes 03- and has the

characteristics of being invariant to the change in the mean
of its distribution, and depending only on 6, is the
coefficient of variation (Linhart, 1965). From 3.8 and 3.9

this coefficient is found to be 1:46 for a? and J23 for 032.

 

Small values of 6 indicate precise knowledge about a? prior
to observing the data.

The last set of assumptions concerns the prior knowledge
about the parameters A,1:2,6, and of. They are the hyper-
parameters for the exchangeable distributions for the
parameters {Uj} and {03.}.

The two questions presented earlier by Bretthorst (1988)
might serve as general guidelines for choosing prior

distributions. However, when. the model in 3.1 is being

 

 

 

 

 

 

47

applied, prior knowledge about the state of the hyper-
parameters A, 1:2, 6, and of is usually limited or non-existent.
That is, the higher we move in the hierarchical structure of
the parameters in the model, the less information we have
about those parameters prior to observing the data. In other
words, one is usually more informed about the distribution of
the {Uj} and the {0%}, prior to observing the data, than about
the distribution of A,t2,6, and 03. Therefore, the prior
distributions for A, t2, 6, and of should have the
characteristic of being non-informative priors. This implies
that the posterior distributions for these parameters should
rely heavily on the information obtained from the data, with
little or no weight given to the prior knowledge about the
parameters.

For the parameter A, we take its conjugate prior
distribution, (which is normal) with vague prior knowledge
about the parameter itself. Therefore, A is distributed

uniformly over its parameter space R that is p(A) °< CA,

q+p+1'
where C). is constant. The subscript A is being used with the
constant to identify it with the parameter.

Similar to 0;, a natural conjugate prior for 12 would be
an inverse chi-square with v12 degrees of freedom. A vague
prior for “t2 implies v12 = 0, which produces what is often
called Jeffery's prior (i.e., p(tz) °<1/1:2). This prior

distribution is found to cause some difficulties when the

number of groups is relatively small (Seltzer, 1988). The

 

 

 

 

 

 

48

effect of this prior with t2 near zero has a great influence
on the density function of the data such that the marginal
posterior distribution of 1:2 becomes heavily concentrated
around 12 = 0 (Morris, 1983 and Lindley, 1983). As an
alternative, a uniform prior distribution is chosen for‘l:2
(Seltzer 1989). That is p(tz) °< 0,2 for (1:2 > 0), where C12 is
a constant. The prior distributions for both A and 1:2 will
then become insignificant in their contribution to their
posterior distributions.

For 6 and of, we write their joint prior distribution
proportional to the product of the conditional prior
distribution of of given 6 and the prior distribution of 6.

That is

6)pfl6). (3.10)

 

p(of, 6) °< p(o.2

We specify p(of [6) to have conjugate prior to the density

in 3.7, which is of a chi-square form as

 

 

_ 2
PM.2 6,I,w) °< (01.2)“1 exp[ (mo) , (3.11)

26

where r and w are two parameters describing the prior for of
(Lindley, 1971). A vague prior for of implies setting r=0
and (0:0 which will produce what is known as Jeffrey's prior.
However, the fact that the marginal posterior distribution of of
is heavily dominated by the data undermines the effect of the
prior of 0.2 on its marginal posterior distribution. Thus, for

mathematical convenience, a constant prior Co? is specified

 

 

 

 

49

for 0.2 .

In the second set of the assumptions, the hyper-parameter
6 is presented as a function of the degrees of freedom for the
exchangeable distributions for {03-}. It also, represents the
variability of {03.}, and mathematically there is no reason why
it should not be any positive number (Lindley, 1971).
Therefore, we assume that p(6) has a uniform distribution on

the real line. That is p(6) °< Ce: where C6 is a constant.
The Joint Density Function of the Parameters in the Model

Applying the Bayesian approach to the model in 3.1
allowed the parameters in the last two stages of the model to
be treated as unknown parameters. In the first stage Hid and
{03-} for j = 1,---,k are the parameters with exchangeable priors.
The second stage covers the hyper-parameters A, 1:2, 6, and of,
where A is defined in 3.2. They are the parameters of the
prior distributions for id; and (0?. Vague prior
distributions for A” 12,6L and.o§ are taken from their usual
conjugate families.

Treating all the parameters in the model as unknown
parameters, allows their joint distribution along with the
data I“ to be distributed as a product of conditional

.7

distributions as follows:
2 2 2 2 2
jfi1p(yjilr thoj)p(Uj|A"1'-2)p(ojle’ WPW‘ '9'°°)- (3.12)

Substituting the density functions for the assumed

 

 

 

 

 

50
distributions from the previous section for each of the above
conditional distributions produces a joint density function

proportional to

-n

fi ( 2)Tj 1f / -2
.=1 Oj exp -§j=1(Yj—ZjA—ljUj) Oj (Yj—ZjA—ljUj)

 

___k l
x (1:2) 2 exp - U]?
21:2 j=1

 

(3.13)

 

Theoretically, the marginal posterior distribution of any
parameter in the model can be obtained by integrating the
joint density function in 3.13 over the space of the other
parameters.

Being able to obtain the marginal posterior distribution
of any parameter in the model implies being able to obtain
interval and point posterior estimates of that.parameteru For
example, an estimate of A based on its marginal posterior
distribution will not be affected by the uncertainty in
estimating the variance components, avoiding an undesirable
characteristic of empirical Bayes estimates of A. Further,
individual group estimates of {03-} for j =l,-~-,k become useful
when we are investigating groups' heterogeneity of variance.
Also, individual estimates of‘iqfi become useful for checking

the normality assumption for the error term associated with

each group intercept. Moreover, inferences on the hyper-

 

 

 

51

 

parameters 12, of, and 6 can be made using their respective
marginal posterior distributions.

With the help of modern computers and developments in
simulation theories, we are now able to approximate these
marginal posterior distribution numerically, to the desired
degree of accuracy. The methods of Data Augmentation (Tanner
and Wong, 1987) and Gibbs sampling (Tanner, 1993) are being
used for numerical integration to obtain marginal
distributions. Morris (1987) shows how one can use these

methods for hierarchical Bayes models.

 

 

 

CHAPTER 4

 

Obtaining Marginal Posterior Distributions

Via Gibbs Sampling

This chapter presents the method of Gibbs sampling as an

procedure for calculating marginal posterior

iterative
The process of Gibbs sampling is best

distributions.
understood within the context of data augmentation (Tanner and

the basic idea of the data

An example, using
A

Wong, 1987). Therefore,

augmentation procedure will be explained.
the normal distribution, will illustrate its application.

 

simple modification to the idea of data augmentation will

facilitate the understanding of Gibbs sampling procedure.

This procedure will then be used in chapter 5 to approximate

the marginal posterior distributions of the parameters in the

model in 3.1.
Data Augmentation

(1964) made in their

The argument that Box and Tiao
investigation of the importance of the assumptions applied.to
the comparison of variances, using the Bayesian approach, can
facilitate our understanding of data augmentation. Very often

the distribution of observations Y depends on more than one
and Z2, one of which, say

say two parameters Z1
For'non-Bayesianmmethods,

parameter,
2&, is of immediate interest to us.

52

 

 

53

 

this could cause an extremely difficult problem in dealing

However, under the Bayesian

with the other parameter Z2.

approach making inferences about Z1 is simplified by finding

the marginal posterior distribution LMZIIY). This marginal

posterior distribution can be obtained by integrating out the

of parameters Z2 from the joint posterior

other set
«distribution p(Zl,Z2|Y). In the data augmentation process, Z2

latent variable or missing data

can be thought of as
1987).

augmenting the observed data Y (Tanner and Wong,

In the following discussion all distributions are

posterior distributions, therefore, the term "posterior" is

deleted from the names of the distributions for simplicity.

If we write the joint distribution p(Z1,Z2|Y) as a product of

the conditional distribution of Z1 and the marginal

distribution of Z2,

p(Z1,ZzlY) =p<ZliZ2.Y> p<Z2|Y> . (4.1)

the marginal distribution of Z1 can then be written as

p(ZIIY) =fTP(Z1iZ2I Y) p(ZzlY) aZ2 ’ (4-2)
where T is the parameter space for Z2. As Box and Tiao (1964)

pointed out:
the marginal posterior distribution of the
parameter ,p(Z2|Y) acts as a weight function
multiplying the conditional distribution
P(Z1IZ2r Y) of the parameter of interest. It is
frequently helpful in understanding the problem and
the nature of the conclusions which can safely be
drawn to consider not only p(Z1|Y) but also the

components of the integral on the right-hand side
One is thus led to consider the

of Z1 for particular

of equation [4.2].
conditional distribution

 

 

 

 

54

 

values of the nuisance parameter Z2 in relation to
the probability of occurrence of the postulated

values of the nuisance parameter." (p. 153).

These two basic ideas, integrating out the nuisance

parameter Z2 from the joint distribution p (Z1, Z2 I Y) , and using
the marginal distribution of the nuisance parameterp(Z2|Y)
as a weight function in that integral, can be generalized to

the process of data augmentation when a direct solution to the

integral in 4.2 can not be found. Tanner (1993) defines

p(ZZIY) in the integral in 4.2 as the predictive distribution

which can be used with the observed data Y to obtain the

posterior distribution p(leY) . In many cases however,

p(ZZIY) is not known, which makes it impossible to obtain

 

p(leY) . The joint distribution of Z1 and Z2 as in 4.1 and
4.2 can then be used to obtain the marginal distribution

,p(Z2[Y). In other words, the joint distribution of Z1 andZ2
in 4.1 can be written as:

p(Z1.Z2IY) =p<zziZ..Y> p(Z.|Y> (4.3)

The marginal distribution for the second parameter p (Z2\Y) can

be obtained by integrating the joint distribution in 4.3 over

the parameter space R of the first parameter le

p(ZZIY) =pr(Z2!Z1,Y)p(Z1|Y)6Z1. (4.4)

the marginal distribution of the

Just as in 4.2,
in 4.4 can be thought of as a weight

parameter p(leY)
function multiplying the conditional distribution of the

parameter P(Z.2(Z1r Y) . Integrating this weighted conditional

 

 

 

 

 

 

 

55
distribution over all admissible values of Z1 will give us the
marginal distribution of Z2.

Carrying out the integration in equation 4.4 analytically
requires knowing the parametric form of the marginal
distribution of the first parameter p (leY) . Likewise,
carrying out the integration in equation 4.2 analytically
requires knowing the parametric form of the marginal
distribution of the second parameter p(ZZIY) . But these two
marginal distributions are unknown to us and we are interested
in finding them. Therefore, this dependency between the two
marginal distributions for the two parameter becomes the key
point of the iteration process of data augmentation.

In 4.1 and 4.3, the joint distribution of Z1 and Z2 is
being expressed as a function of four other distributions,
two of which have unknown parametric forms. Therefore, it is
hard to sample from them. They are the marginal distributions
p(leY) and p(Z2|Y) . The other two have known parametric
forms, and one can easily sample from them. They are the
conditional distributions p (Z1|Z2, Y) and p (Z2 IZ1, Y) . The
iteration process, therefore, involves repeated sampling of Z1
and Z2 from their conditional distributions p(lezz, Y) and
p(ZZIZI, Y), which accomplishes the numerical integration of
the joint distributions in 4.2 and 4.4, and produces numerical
approximations to the marginal distributions for Z1 and Z2.

To exploit the dependency between the two distributions

in 4.2 and 4.4, we start with an initial estimate of Z1. One

 

 

 

 

 

56
can think of this initial estimate as being sampled from a
poor approximation of its marginal posterior distribution,
p(leY) . A sample of M values of Z2 can then be drawn from
its conditional distribution p(ZZIZ1,Y) , where the initial
estimate of Z1 is used along with data Y. For each of theM
values of Z2, we form a conditional distribution of Z1 given
Z2=Z2m defined as p(ZIIZZm, Y) , for m=1, «-,M. The weighing of
each of these conditional distributions is done empirically by
the sampling process from the marginal posterior distribution
p(ZzlY) in 4.2. That is, a large proportion of M of these
conditional distributions are conditioned on values of Z2
which are near the mode of p(ZZIY) . Rubin (1987) refers to
the process of sampling Z2 from p (ZZIY) , in order to form the
conditional distribution of the first set of parameters
p(ZIIZZ, Y) as a multiple imputation process. The mixture of
the M conditional distributions of Z1 given Z21", m=1,~~-,M
represented by their weighted average, where the weighing
process is done empirically (as explained previously),

produces the marginal posterior distribution p(leY) at the
tth iteration
M
P"<Z.IY> x M‘1 Zp<leZZml Y) - W5)
m=1
To simulate the marginal distribution of Zl given in 4.5,
we simply sample M values of Z1; one from each of the M

conditional distributions given in the mixture. These M

values of Z1 represent an approximation to the marginal

 

 

 

 

57

distribution of the first parameter p1(Z1)Y) at time t=1 of
the iteration process. However, this approximation is a poor
one since the sampled values of Z2 came from a rough guess of
its marginal distribution p(ZzlY) . Representing the M
weighted conditional distributions by their average captures
the idea of integrating the joint distribution of Z1 and Z2 in
4.2 over the parameter space of Z2 to produce the marginal
distribution of Z1, when M is relatively large (say M =
4000). An estimate of the density function for the marginal
posterior distribution of Z1 can be obtained by averaging the
conditional densities in 4.5 (Gelfand and Smith, 1990).

For each of the M values of Z1, which are sampled from
4.5, a conditional distribution of Z2," given Z1=Z1m, can be
formed as p(ZZIZm, Y) , for m=1,-«,M. By adopting the same
logic for approximating p (Z1 | Y) , the marginal distribution for

the second parameter can be approximated by

pt(Z2lY) x M‘1 f: p(Zzlzlm, Y) , (4-6)
m=1

where p1(Z2|Y) is an approximation of the marginal

distribution of Z2 at iteration t=1. We then sample M new

values of Z2 for iteration t=2 from the resulting

distribution in 4.6 to get a new approximation to the marginal
distribution of Z1 in 4.5.

By iterating between 4.5 and 4.6, the two marginal
distributions p(leY) and p(ZzlY) stabilize (Tanner and Wong,

1987) . Sampling the values of Z1 and Z2 from their respective

 

 

 

58
stabilized distributions is virtually identical to sampling
from the true marginal distributions of Z1 and Z2.

Tanner and Wong (1987) pointed out that one can select
any value for M'(the sample size for each Z1 and Z2 in each
iteration) to carry out the iteration process. They stated:

"Even when it is as small as 1, the iteration is

still 'in the right direction' in the sense that

the average of p(ZIIZZm' Y) over the augmented data

patterns generated across iterations will converge

to p(leY) ." (Tanner and Wong, 1987, p. 530).

Moreover, they argued that the value of M can be changed
between iterations. In fact they recommended that one Should
start with a small value for M'in the first few iterations and
increase this value as the number of iterations increases.
They stated:

"In practice, however, it is inefficient to take A!

large during the first few iterations when the

estimated posterior distribution is far from the

true distribution. Rather, it is suggested that A!

initially be small and then increased with

successive iterations." (Tanner and Wong, 1987, p.
539).

Example

Let Y'be a random sample of n independent observations
drawn from a normal distribution with mean u and variance 02,
where both parameters are unknown. The sample mean ysrrigbg
and the sample variance 32 = (n—l)“11:31071-17)2 are jointly

sufficient for (u,cfi). The likelihood for (u,02), therefore,

is given by (Box and Tiao, 1973)

 

 

 

 

59

UHIOZIFISZ) °< (02)“2- exp{-2—_lz-[(n—1)s2 + n(u—37)2]} . (4-7)
a

The likelihood in 4.7 reveals that the distribution of Y
depends on the two parameters (1 and 02. From the Bayesian
viewpoint, making inferences on any one of these two
parameters requires finding its marginal posterior distribu-
tion by integrating out the other parameter from their joint
posterior distribution.

The joint posterior distribution of p. and 02 can be
written proportional to the product of their likelihood

function Up, 0257, $2) and the prior knowledge about u. and 02.

Mathematically stated

P(|J,02|Y) °< 0(u,02|37,sz)p(p,02) (4.8)

The second part of the right-hand side of 4.8 represents the
prior knowledge about the two parameters )1 and 02 . It can be
assumed priori that p. and o2 are independent, and the form of
this prior knowledge can be specified by the usual conjugate
priors of these two parameters. A vague prior knowledge about
)1 produces a prior distribution which is proportional to a
constant C. Similarly, a vague prior knowledge about 02

produces a prior distribution proportional to 0‘2 (Box and

Tiao, 1973).

Substituting the assumed forms for the likelihood

function Up, 02 |37, 52) and the prior distributions for p. ando2

in 4.8 produces the joint posterior density for u and 02.

 

6O

p(il,02lY) .. (.2)-<%+1>.xp{2-_12[(.-n.2 + nut—7m} . (4.9)
0

Based on 4.9, and where 02 is known the conditional
distribution p(IJ-loz, Y) is normal with mean 57 and variance
02/11. When 02 is unknown, however, P(l-|v|Y) represents a t-
distribution with mean )7 and variance Sz/n and (n-1) degrees
of freedom, resulting from the integration of 4.9 with respect
to 02 (Box and Tiao, 1973).

When p. is known, the sample variance is defined to be
5'2 = 1nd};1 (yi - (1)2. The variable n.5'3/o2 is therefore
distributed as chi-square with 11 degrees of freedom.
Consequently, ozlu, Y is distributed as nszxfi) where xii, is an
inverted chi-square variable with 11 degrees of freedom. When
u is unknown however, p(ozlY) represents the probability
density of (n-1)Szx2,2,-1) , where x7534) is an inverted chi—square
variable with (n-l) degrees of freedom, resulting from the
integration of 4.9 with respect to n.

For illustrative purpose, let us assume that the forms
for the distributions of p(ply) and p(ozlY) are unknown, and
need to be approximated. Our objective in this example is to
show how these two distributions can be numerically
approximated, using the process of data augmentation presented
in the previous section.

The joint posterior distribution of p. and 02 is written

as a product of the conditional distribution of p. and the

marginal posterior distribution of 02

61

p(wzm =p(p.|02,Y) p(ozly) . (4.10)

Then the marginal posterior distribution of u can be defined

by
p(plY) =fozp()1|02,Y) p(ozlY)d02. (4-11)

Similarly, the joint.posterior distribution.in 4.10 is written

as a product of two new distributions

p(p.02|Y) =p(02|u,y) P(|J,|Y) (4.12)

Then the marginal posterior distribution of 02 can be defined

by
p(ozlY) =f p(ozlp,Y)p(plY>du. (4-13)
p.

Examining 4.11 and 4.13 reveals that obtaining' the
marginal posterior distribution of one parameter depends on
obtaining the marginal posterior distribution of the other
parameter. Therefore, this dependency between the two
equations can be used to ShOW' how the method of data
augmentation can be used to approximate the two marginal
posterior distributions in the following steps.

1. As a starting point for the iteration process, u can be
estimated from the data Y by )1 = gin/fl. This estimate
can be thought of as p, which has been sampled from the

current approximation of its marginal posterior

distribution p(ulY).

 

62
Using the current estimate ii to represent the true value
of u, find the sample variance as S2 =n’1£1(yi - (1)2.

1-

This variance estimate will be used, in conjunction with
the sampled values of a chi-square random variable, in
the next step to sample 02 from its conditional
distribution p(02 In, Y) .
Sample M values of x201)”, for m=1,~,M from a chi-square
distribution with 11 degrees of freedom. For each
sampled value of )(Zm)m find 63,, as 83,, = nSz/xinm. Applying
steps 2 and 3 resembles the process of sampling M values
of 02 from its conditional posterior distribution
p(ozlil ,y) . The mixture of these M conditioned values of
02 forms an initial approximation to the marginal
posterior distribution p(ozlY) . This approximation,
however, is an inaccurate one since the sampled values of
02 are conditioned on a poor estimate of p.
Given the data Y and the M values of 0,3, which are
sampled in step 3, we sample M new values of pm from
p(plog, Y) , which is N(iog/n) , where 57 is the sample
mean from the data. The mixture of the M sampled values
of the mean pm represents an approximation to its
marginal posterior distribution p(p. | Y) .
These new sampled values of pm from step 4 can then be

used to get M new estimates of 8; as in step 2 for a new

cycle of approximation. As we continue iterating between

steps 2 to 5, the mixture of the M values of 0,3, which

 

 

 

63
are sampled in step 3, and the mixture of the M'values
of p” which are sampled in step 4 become increasingly
accurate in representing the marginal posterior

distributions of u and 02.
Gibbs Sampling

The process of data augmentation was presented for the
case:of approximating the marginal.posterior distributions for
only two parameters, Z1 and Z2 with M >1. When M=1, Tanner
(1993) defines the iteration process as "chained data
augmentation". In the normal data example given above, the
two parameters were represented by p and 02‘with.At>1. There
are cases however, where there are more than two parameters
for which we require numerical approximations to their
marginal posterior distributions. To obtain them, a simple
modification to the logic of data augmentation process can be
used.

When. the iteration. process. of data augmentation is
generalized to more than two parameters with M'set equal to
one, the process is called "multivariate chained data
augmentation” (Tannery 1993) or' Gibbs sampling" (Gelfand,
Hills, Racine-Poon and Smith, 1990).

Consider the previous case, where the distribution of the
observations Y depends not only on Z1 and Z2 but also on a

third parameter say Z3. Under the Bayesian approach, making

 

 

 

 

64
inferences about any of the three parameters is simplified by
finding its marginal posterior distribution, by integrating
out the other sets of parameters from the joint distribution
of Z1, Z2 and Z3.
The joint distribution of Z1, Z2 and Z3 can be written as
a product of the conditional distribution of Z1 and the joint

distribution of Z2 and Z3,
p(Z1,Z2,Z3|Y) =p(Z1|Z2,Z3,Y) p(Z2,Z3|Y) . (4.“)

Then, the marginal posterior distribution of Z1 can be written

as

P(Z1)Y) = foWp(Z1]Z2,Z3,Y) p(Z2,Z3lY) 6Z3 8Z2 ’ (4-15)

where T and W are the parameter spaces for Z2 and Z3,
respectively.

Similar to 4.2, the joint distribution p(Z2,Z3|Y) can be
thought of as a weight function multiplying the conditional
distribution of the parameter of interest p(ZIIZ2,Z3, Y) . In
other words, we consider the conditional distribution of Z1
for particular values of parameters Z2 and Z3 in relation to
the probability of getting those values of Z2 and Z3.

Similarly, the marginal posterior distribution for each
of Z2 and Z3 can be obtained by integrating out the other
parameters from the joint.posterior distribution of Z1, Z2.and

Z3 as follows :

p(ZzlY) = IR pr(Z2lZ1,Z3. Y) P(Z1IZ:))Y) 6Z3 621 ' (4,15)

and

 

 

 

 

 

65

p(ZBIY) = f. pr(z3|zl,z2,y) p(Z1,Z2IY) az2 azl. “.17,

As in data augmentation, the three equations 4.15, 4.16
and 4.17 define the iteration process of Gibbs sampling.
Further, it is assumed.that there are only three distributions
that are easy to sample from. They are the three conditional
distributions of Z1, Z2 and Z3, represented by the first part
of the right-hand side in each of 4.15, 4.16, and 4.17. The
iteration process therefore involves repeated sampling of Z1 ,
Z2, and Z3 from these conditional distributions to accomplish
the. numerical approximations of 'their’2marginal posterior
distributions.

Starting with initial values of ZS”, ZS”, and the data
Y, sample ZFJ from its conditional distribution
p(z1|z2‘°’, Z3)”, Y) . Given the values of Z)” , ZS” , and the data

2(1) from its conditional distribution

Y, sample
P(Z2(ZFJ,Z§”,IO. Finally, given the values of Z9), ZEJ, and
the data Y, sample Z?) from its conditional distribution
P(Z3IZ1(1),Z2(1), Y). The three sampling processes complete one
iteration of Gibbs sampling. After a large number of
iterations say X, sampling' any' of’ the three ‘parameters
resemble sampling that parameter from its marginal posterior
distributitwn To simulate a marginal distribution for Z3, for
example, after an X iterations, sample m values of Z3 from
(Z3 |Z{X), Z2(X),Y) . The mixture of these m sampled values of Z3
represents a numerical approximation to its marginal

distribution. That is, the sampled values Z3901, ...... , Z3“), can be

 

 

 

66

 

viewed as a numerical approximation to the marginal posterior
distribution _p(Z3|Y). Similarly; the :marginal. posterior
distributions for Z1 and Z2 can be approximated.

Gibbs sampling can be generalized to as many parameters
as the investigated model has. The basic idea is to write the
joint distribution of all parameters in forms similar to 4.15
to 4.17. Repeated sampling of the parameters from their
corresponding conditional distributions results in
approximation of their marginal distribution.

A special case of Gibbs sampling is found for certain
applications, where some of the parameters depend only on a
selected number of other parameters. For instance, in the
case of approximating the marginal posterior distributions for
the three sets of parameters Z1, Z2 and Z3, it is possible
that the conditional distribution of say, Z3 depends only on
Z2 and the data 1%. While equations 4.15 and 4.16 for the
marginal distributions for Z1 and Z2 do not change, equation

4.17, however, becomes

p(Z3IY) =pr(Z3|Z2,Y) P<Z2|Y) 5Z2 - (4.18)

A repeated sampling of Z1, Z2 and Z3 from their conditional
distributions in 4.15, 4.16 and 4.18 will approximate their
marginal posterior distributions. Chapter 5 will show how the
special case of Gibbs sampling is being used for approximating

the marginal posterior distributions of the parameters in the

model presented in 3.1.

 

 

«« .__ _,_...-~—

 

 

CHAPTER 5

 

The Application of Gibbs Sampling to The Model of The study

Two sections make up this chapter. The first section
presents the application of Gibbs sampling for obtaining the
marginal posterior distributions of the parameters of the
model in 3.1. The conditional distributions of these
parameters are also derived so they can be used in the
application of Gibbs sampling. The second section presents
the steps taken to empirically test the application of Gibbs
sampling. These steps include setting the model, specifying
and generating artificial data sets, assessing the
heterogeneity of variance, obtaining initial estimates for the
parameters to start the iteration process, assessing
convergence of the iterative program, and describing the real

data used in this study.

Obtaining Marginal Posterior Distributions of Parameters of

the Model in 3.1

The proposed model in 3.1 involves two sets of

parameters: those which have exchangeable prior distributions

{UT and {0? for j=1, ,k (i.e., {My =(Q,LQ,W,CQ and
{02}: 0i,0§,m,oi) and the hyper-parameters 12, A, 6, and 03
where AH%Y°,Y1,m,Yq,Blf”,Bp}. Our objective is to approximate,

67

 

 

68
numerically, the marginal posterior distribution for each one
of these parameters using Gibbs sampling.

In general, the marginal posterior distribution for each
one of these parameters can be found by integrating its joint
posterior distribution with the other parameters over the
parameter spaces of all other parameters. Mathematically
stated, given the joint distribution in 3.12, the marginal
distribution for each of the parameters of the model in 3.1

can be written as follows:

p({Uj}lY) = £!{!§)p({Uj}|t2,A,{0§},Y)p(1:2,A,{0§-}|Y)

x 6{0§}6(A)6(12) , (5.1)

p({o§}|y) = [ff
6A

fl

[Whig-HOE, 6, A, {Uj}, Y)p(o§, 6, A, {Uj}l Y)
W)

J

x a{Uj}a<1)a(e)a(o§) , (5.2)

p(AlY) = ffp(A|{0§),{UJ-}, Y)p({o§-},{Uj}|Y)6{Uj} diofi} ,

(5.3)
(03.1%}
p(rZIY) = fp(12){Uj}, Y)p({Uj}lY) ain} , (5.4)
(g)
p(6 Y) = p(6 0.2,{02}, Y)p(0§,{02-} Y) 6{02-} 6(03) ,
l ff I J J I J (5.5)

of (a?)

 

 

69

mean = Hurling), Y)p(e,{o§}m aieg} m)

5.6
6 (a?) ( )

Clearly, any one of the above marginal posterior distributions
depends only on some of the parameters in the joint posterior
distribution in 3.12. More specifically, the marginal
posterior distribution for {U} in 5.1 for example, does not
depend on 03 and 6. Similarly, the marginal posterior
distribution for t2 in 5.4 depends only on {U}. This is true
for the rest of the parameters.

The first term of each integral in 5.1 to 5.6 represents
the conditional distribution for one of the parameter(s) of
interest, which usually has a known parametric form and from
which we can sample. The second term of the integral can be
thought of as a weight function multiplying the conditional
distribution of the parameter of interest. Integrating the
weighted conditional distribution over all the admissible
values of the parameters in the second term will produce the
marginal posterior distribution of the parameter of interest.

To exploit the dependency between the marginal posterior
distributions in 5.1 to 5.6, two things must be known: the
parametric forms of the conditional distributions in the

integral in 5.1 to 5.5 and initial approximations to the

values of the parameters A”), 02w), hpOW, and hfi‘m}. The
superscript "(0)" used with each of those parameters
represents the initial value of the parameters. More

discussion will be given later on the derivation. of the

 

 

7O

parametric forms of the conditional distributions of the
parameters and the method of obtaining initial values for the
specified parameters. Let us now assume that we know the
parametric forms for the conditional distribution and we can
sample from them. Given current estimates (start with initial
values) of the parameters A“), 03”), RED”, and {03wU, the
iteration process starts as follows:

1- Given HMO” and the data Y, sample one value of 12”) from
_p(12HUjm},Y) in 5.4. Notice that the superscript "(1)"
used with 12 refers to the first cycle of the iteration
process.

2— Given {0§WU, 03m) and the data Y, sample one value of
6”) from its conditional distribution p(6l03m),hfi(m},Y)
in 5.5.

3- Given the sampled value of 6”) from step 2, {0§WU, and

0.3 (1)

the data Y, sample one value of from its

conditional distribution p(og 6“’,hfi(m),Y) in 5.6.

 

4- Given the sampled value of 12”) from step 1, {ofimd,
A‘m, and the data Y, sample one set of values of
{0910-= H41),U§”,o~,0fi“} for j=1q-~,k from its conditional

distribution p(QGHA(M,IzuH{0§mU,Y) in 5.1.

5— Given the sampled values of Qfilw, the sampled value of
6”), the sampled value of of”) from steps 4, 2 and 3
respectively, the initial value of A‘m, and the data Y,

sample one set of values of{0§”U ={0iu),0§u),m,0iuU for

 

 

 

 

 

71
i=1, ,k from its conditional distribution
p({0§}|0.2(1),6‘1),A‘O’,{Uj(1)}, Y) in 5.2.

6- Given the set of Hg”), the set kfidu}, and the data Y,
sample one value of A”) from its conditional
distribution p(Alin‘“},{o§-‘“}, y) in 5.3.

Going through these six steps of sampling from the given

conditional distributions finishes one cycle of the iteration

process in Gibbs sampling. This cycle produces the first
sampled value of each of the parameters, which can then be
used to produce new sampled values of the same parameters in

the second cycle. In the second cycle, starting with step 1

again and given {Ufu}, generate new value of t2”) where the

superscript "(2)" used with 12 refers to the second cycle of
the iteration process. As we continue iterating between the
above six steps, the mixture of the values of any one of these
six sets of‘ parameters becomes increasingly accurate in

representing its marginal posterior distribution.

Finding the Conditional Distributions for the Parameters

of the Model in 3.1

One of the requisites for obtaining the marginal
posterior distributions in 5.1 to 5.6 through the iteration
process of Gibbs sampling is to be able to sample from the
conditional distributions for the respective parameters. This

implies the need for identifying the parametric forms for

 

 

 

 

 

 

72
these conditional distributions. Within the Bayesian approach
the elements of the set {U} =CQ,Lg,m,Uk are assumed
independent and identically distributed with respect to their
conditional distribution as defined in (5.1).. Identifying the
conditional distribution of any element of this set therefore,
can be generalized to the other elements of the set” The same
logic applies in determining the conditional distribution for

{0;} = of, 03, m, 0%.
The conditional distribution of U]. given A, 12, {02}, and Y

From lines 1 and 2 of the joint density function in 3.13,
the conditional distribution of Ule, 172, {03-}, Y is found to have

a density function proportional to

"j

__ 1 _

 

5 “fUJ? (5.7)
j=1

x (1:2) 2 exp
21:2

Since Uj appears only in the exponential parts of the

expression in 5.7, the density function in 5.7, can be written

proportional to

2

l -2 U7"
exp --§ng (dj-ljUj)/oj (dj-ljUj) + :5. ,

(5.8)

where Yj-ZjA = dj. Ignoring the terms that do not depend on

the above expression can be simplified to

Uj,

 

 

 

 

 

 

73

exp{-%)§(—2(UJ/-1§dj)o;2 + U]? (133,032 + 14”} . (5.9)

l

Completing the square in 5.9 results in

_l _ -2 _ - / -22 -2 -2 5.10
exp( Egajj (HjO'j + T 2) lljdjO'j) (DJ-0:] + T )) . ( )

Substituting for 1§dj==nj(Y3—Z5A) in 5.10 results in

exp(-%$<Uj-njo;2 (111.032 +1—2)‘1(Yj—§]/.A))2(njo;2+1:‘2)). (5.11)

The expression in 5.11 represents a kernel of‘a normal density

function. Therefore, £5|A,t2,kfih Y is distributed as normal

with mean U? and variance Vé, where

. _ —2 -2 _ _ _ —/ __ -2 — —
Uj - 133-0: (njoj +1: 2) l(Yj—zji) . Va, — (njoj +1: 2) 1. (5.12)

The conditional distribution of 0% given A, 6, cf, {@3, and Y

From lines 1 and 3 of the joint density in 3.13, the
conditional distribution of 0§|A,v°,03,&gLIY is found to have

a density function proportional to

 

 

 

 

74

l _
exp(-§ (Y,- - zjk - 1,-ija,2(yj - zji - 13%))

(5.13)

 

1 02 % -(—1-+1) —02

1 (26) (03-) 2" exp —-°—

rim)
_ _ / _ _ _ 2

rewritten as

5.13 can be

 

2 2
2__ j 26 ) exp _ anj + 0.
]TEL) 20? 260;

 

(5.14)

Ignoring the terms that do not involve 0;, results in writing

the expression in 5.14 proportional to

2

n' 2 2
(0%)“(73—4'31-54'1) 1(I1jsje + 0'.)
Ci

exp-—— 26 (5.15)

In general, a continuous random variable X with a density

function f’is said to have an Inverse Gamma distribution with

parameters a and B if

f(X)==——¥L——.X*“”)exp(~3;), Osxsw

F(a)pa KB (5.16)

The expression in 5.15 therefore, represents a kernel of a

Inverse Gamma density with the two parameters

.n. 1 _ 26
a = .7 + I 2' — I (5.17)
a; 2 26 (303 an§6 + 0..2

 

To sample 03. from its conditional distribution, the

following strategy can be adOpted: define 03- = h(0}2) , where h

75
is a one to one function defined as 11(032) = (032)‘1. The

conditional density function of 032 is found to be

 

 

'2
g<o§2|1,e,os,{u,},n «f(a§.|),e,oi,w,},y> ah‘—‘i;) , (5.1.,
6(0j)
6h 72
where (—C::) = (032)”. Therefore,
6(0j)

 

 

31.1- - I? 3
g(0}2|A,6,0§,{Uj}, Y) °< (032) 2 2° leXP{_°72(szg+i)} I (5.19)

In general, a continuous random variable X with a density

function f is said to have a Gamma distribution with

parameters a and [5 if

f(x) = X“—1 exp(——pi{) , O s X 3 9° . (5.20)

__1_
PM) B“
The expected value and the variance of X are given by

(5.21)
E(X) = a6 , Var(X) = «(32.

The expression in 5.21, therefore has the kernel of Gamma

density function. Thus, given A, 6, 0.2, {Uj}, and Y, the parameter 032

is distributed as Gamma with

a -2 = E + _1_ , (30.. = ’2—6—7 (5.22)
“i 2 26 J 11].st + 0,,

To sample 03-, simply sample 0? from 5.21 then use the

_ — _ 2
relation h(0j2) = (of) 1 to get 03-.

' - a
The conditional expectation and the variance of oj are

easily found by applying 5.21,

76

 

E(O}2|1,0,03,{Uj}, Y) =((u)jSJa + (1—(i)j)0..2)'1 , and
Var(072|A e 02 {U} Y) = 2"” 32 + (1 — ) 2'2
J I I 0] j] n. [(1)-7' j Qj 0.] , (5023)
.7
where wj =-——£fi——
nj + 6‘1

The conditional expectation in 5.23 has a form similar to
Stein's (1956) shrinkage estimate. It is a form of a weighted
estimate of the variance Si, and the asymptotic overall mean
of the variances 03. The weights depend, in a natural way,
on the degrees of freedom for obtaining S? and 03. As 6
approaches zero, indicating homogenous variances, so does wj
leaving the expected value of of equal to 0?. When 6 gets
very large however, the value of wj approaches one, leaving
the expected value of 032 equal to the inverse of the sample
variance 392.

The operational definition of reliability is the ratio of
the "true" variance to the observed variance. Therefore, one
——£::—— as a reliability index for 032,
nj + 6'1
expressed in terms of its precision.

can think of 1 - wj =

It represents the

. - 4
proportion of the precision in the OLS estimates of oj that

is parameter precision.

- 2
The conditional distribution of A g1ven {U}, {0fi, and Y

From line 1 of the joint density in 3.13, the conditional

distribution of AIQGL{O?,Y' iS found to have a denSIty

function proportional to

77
exp()5—.l (ir'.—z.A—1.U.)/0’.2 (y —Z )_1 (1)) (5.24)
1 2 .7 .7 .7 J J j j j j '

Let Yj-lj Uj = Dj. The expression in 5.24 can now be written

as

eXp{-%€3[(Dj—zji)’o;2(Dj—ZjA)]} . (5.25)

Multiplying the terms inside the exponential, and ignoring the
terms that do not depend on A, the conditional distribution

of A becomes proportional to

eXp(--%€3 032 (A’Zj-ZjA — zit/25%)) . (5.26)

By completing the square, 5.26 is rewritten as

-1 I
exp [—%[A - (g (Zj-Zj) 0332) )::(z,’.o,) of] $(zj4zj) of
(5.27)

-1
x {A - ($2,421.) of) $2,419,) 052)]
1

When Yj—lj Uj is substituted back for Dj in 5.27, the resulting
expression represents a kernel of a normal density. Thus
A |{Uj}, {03-}, Y is distributed as normal with mean A“ and variance
V,“ where

-1
1* = (gm/.21.) of] (ganja?) — $Z§1jUj<732J ,

.7

-1
V, = (fir/.21.) of) .

1

(5.28)

 

78
The mean in 5.28 is a weighted least square estimate with
Yj-ljUj as an outcome, Z- as predictors, and 032 as the

.7

weighing factor.

The conditional distribution of t2 given UL; and 1’

From line 2 of the joint density in 3.13, the conditional
distribution of tzlflgh Y is found to have a density function

proportional to

1. 9:50,? (5.29)

(12) 2 exp

 

212

Compared to the general form of an Inverse Gamma function in
5.16, the expression in 5.29 has a similar density function.

Therefore, tzlhgh Y is distributed as an Inverse Gamma

variable with

“:2 =

k
——1 fi2=—— -

' t (5.30)
2 iv

The conditional expectation and the variance of 12 are

easily found to be

79

EU]?

E(1:2|{Uj}, Y) = 1:_4 , and

2

262%?)
1

(k-4)2(k-6)

 

(5.31)

Vartc2 I{Uj), Y) =

 

To sample 1:2 from its conditional distribution, adopt the

same strategy used for sampling 0?; that is, let

 

 

 

1:2 = h(’E-2) = (.c-Z)-1 . (5.32)
The density function for t“ can then be found as
- 6h(t’2)
(1: 2 {U.}, r) = rd?- {0.}, Y) ,_ (5.33)
where M = (t‘2)‘2. Therefore,
6(1‘2)
_1_<_2 “17—sz3g
h(t‘2|{Uj}, Y) °< (t’z)2 exp 21 (5'34)

Compared to the general form of a Gamma function in 5.20, the

expression in 5.34 is similar. Therefore, t‘2|{Uj},Y is

distributed as Gamma variable with

 

- 1 , [3.4 = - (5.35)
{‘Uj

The conditional expectation and the variance of 1‘2 are

 

80

E(t‘2|{Uj}, Y) = k—2 and

 

(5.35)
Var(t‘2|{Uj}, y) = 2(k-2)

($412

To sample 1:2, simply sample 1'2 from 5.34, then use the

relation in 5.32 to get 12.

The conditional distribution of of given 6 and k@}

From line 3 of the joint density in 3.13, the conditional

distribution of 03

 

(L{0? is found to have a density function

proportional to

 

 

.2 a 1

20 exp -0. 1_2 fi(0§_)“—26*1’ _ (5.37)
0.,iL))k 26 1 o§ 1

26

Let G and H represent the geometric and the harmonic means of

0§'s as follows:

G: fi(0 971‘, Thus flog-=61c
l

(5.38)

k
H: k , Thus 20
1

2.5

 

ale
I
ml»

Q
eels

Substituting the equalities of 5.38 into 5.37 produces

 

81

J”

(ale-)2

we)" 9"“

 

 

1 -k02
3 ” ° (5.39)
(0 ) exp( 26H) .

The first term of the above expression does not involve 03,
therefore, the above expression can be rewritten proportional
to

(of)?!£6

 

‘kOEJ, (5.40)

exp( 26H

The above expression represents the kernel of a Gamma density

with

k 26E!
“0.2 = 2—6 + l I pa? = T . (5041)

c o c c 2
The cond1t1onal expectation and variance for (n

therefore, are

 

E(of e,{o§}) = H+ 3; , and
(5.42)
Var(0§|6,{0§}) = (H + 3;)3-2-5 .

Notice that the first.part of the above expectation represents

information depicted by the harmonic mean of the groups
variances {03.}. The second part involves 6, which reflects the
variability of these variances. As 6 approaches the zero,

implying homogeneity of variance situation, the second.part in

that expectation diminishes, and the conditional expectation

 

 

82

of 0.2 becomes equal to the harmonic mean of the group

variances.

Similar to 03- and 1:2, we can find the conditional

distribution of 0:2 by letting

 

 

 

 

02 = 9(0‘2) = (0-2)-1 (5.43)
The density function for 0:2 can then be found as
'2 5.44
h(0§2|{0§-},6) = f(0§ (03),.) Ml < >
am?)
-2
where M =(O;Z)‘2. Therefore,
6(0Z2)
k- _
h(o;2l{o§},e) o. (0J2)7’32exp(——k—:-2-) _ (5.45)
26H0.

The expression above has the same general form as the inverse
Gamma function in 5.16. Therefore, 0Z2l{0§},6 is distributed
as an Inverse Gamma variable with «0:2 and (30:2 given in 5.41.

o '2
The mean and variance of 0. are

E (0:2 H03}, 6) , and

_1_
H (5.46)

H

Var(0.f2 H02), 6) -——-l-———
1 His-1)
26

e o 2
The Conditional Distribution of 6 g1ven 03 and {oj}

From line 3 of the joint density in 3.13, the conditional

distribution of 6|03,{0§-} is found to have a density function

proportional to

83

k
(03) 2—0 1 2 1 -(i+1)

__ 0 ex —2§_ fi 02. 23 . (5047)
——P ( ) p 2—6 1 a; 1< 3)
Let G and H be defined as in 5.38. Further, to simplify the
term 1"(2—10) in 5.47 we use the first term of Stirling's
approximation. The general expression for Stirling‘s

approximation for a variable P is given by

l + 1 _ 139 ,(5.48)
12p 288p2 51840p3

 

NP) a (21:) %p (p-é) exp (-p) [1 +

Substituting 1/26 for p in the first term of 5.48 will

produces the following approximation.

k
(F(-§%))k z ((21):)3 (2_16.)2_°3 exp(%6.)) _ (5.49)

Substituting G, H, and 5.49 in 5.47 results in

 

k
i 2" (<33)E - i+ —ko§
(20) G 1((29 1) exp( 26H . (5050)

The expression in 5.50 reduces to

2i
——1 [°°)”(-1-)
1.“; G 6

(26)

uln-

 

—.I<o..2 k
_ 5.51
exp( 26H + 26] ( )

The first term of 5.51 does not involve 6, therefore, the

above expression is rewritten proportional to

(5.52)

 

02
(é) 2 kexp[—k (0° H — 1 - log(03) + 109(G))]

 

 

84
The above expression has kernel of an Inverse Gamma density.

Therefore, 6lof,{o§-} has an inverse Gamma distribution with

 

 

 

 

k 2
a - — — 1, =
9 2 B9 03 2 (5.53)
k H - log(o.) + log(G) — 1
The conditional expectation and variance of 6 are found
to be
k 02
2 2 ° 2
E(6|0.,{oj}) = k-4[H —log(o.) +log(G) -1) , and

02 2 (5.54)
2k2(7° — 109m?) + log(G) — 1)

V 6 §,{2.} =
ar( lo 0]) (k-6)(k-4)2

 

Notice that as k-wo, the above conditional expectation and
variance will be

E(6|o§,{o§}) = H“ . and (5 55)

Var(6la§,{o§-}) = 2—(%:)—2- ,
where u.” = log (G) —log (H) .

Sampling 6 from its conditional distribution is done
using the same strategy used in sampling the groups' variances
from their conditional distributions. That is, define
6 =h(6'1) , where h is a one to one function defined as

h(6‘1) = (6‘1)‘1. The conditional density function

9(6"1|o§,{o§}) is found as follows:

 

 

 

 

 

 

 

 

 

85
g(6'1|03,{a§-}) °< f(6|0‘3,{o§}) lgfigt—L , (5.55)
where -Q§19:2-==(6‘U‘2. Therefore,
6(6'1)
f(6’1|o‘3,{o§}) °< (6‘1)§-2exp[ —6:k((: -1-log(o§) +log(G)]] (5'57)
The above expression is the kernel of Gamma density. Thus,

eflficf,{o? has a Gamma distribution with

k.-

2

 

a6—1=——1, 56-1:

2 (5.58)

 

02
1((1’1: — log(o§) + log(G) - 1)

To sample 6, simply sample 6'1 from 5.51, then use the

relation h(6‘1) = (6‘1)‘l to get 6.
Empirical Applications of Gibbs Sampling

The procedures presented in the previous section for
obtaining the marginal posterior distributions for the
parameters of the model in 3.1 were next tested empirically.
This section presents the processes for this test in two
phases. The first phase involves the use of artificial data
sets with pre—specified values of parameters in a FORTRAN

computer program that implements the controlled statistical

 

sampling technique. One of the objectives at this phase was
to validate the theoretical aspects of the statistical

techniques used. The :resulting' :marginal posterior

 

 

 

 

 

86
distributions of the parameters in the model can be examined
against their true values, which are used in generating the
artificial data set and their parameter. spaces. Another
objective is to check the integrity of the computation of the
computer program so that it can be used with confidence on
real data in the second phase.

{A real data set used by Raudenbush and Bryk (1987) is
analyzed in the second phase of the empirical testing of the
procedure. The data set comes from a sample of 160 U.S. high
schools. It includes 83 catholic schools and a random sample
of 77 public schools drawn from the High School and Beyond
data base. Sample sizes range from 14 to 66 student per

school with an average of 44.5 student per school.
Artificial Data Analysis

The general model presented in 3.1 allows for Q'between-
unit variables W’ and p within-unit variables X)
Theoretically, the statistical procedure presented in the
previous chapter works with any number of between and within-
unit variables as long as k>q and nj>p. However, for
practicality; three versions of that model will be presented.
In the first version, no between and within-unit variables
were specified” This is the simplest case, where the model in
3.1 reduces to a simple ANOVA model. One between and one

within-unit.variable will be included.in the second version of

 

 

 

 

 

87
the model with a random intercept and fixed effects
predictors, i.e., q==g3= 1. The last version covers the case
of two between and three within-unit variables with random
intercept and fixed effects predictors , i.e., q = 2 and
;)= 3 (see Table 5.1).

The above three models range in their complexity to
represent. conditions often found. in. educational research
projects. The first is a simple version of the model, used to
facilitate the understanding of the statistical procedure and
its application, The third version of the model is similar to

the example of the real data used in the study.

Data Specifications

Two primary but not completely mutually exclusive
criteria were used to create the data. The first deals with
the degree of heterogeneity of variance. The second criterion
deals with the number of groups k to be generated. For each
combination of these two criteria along with each of the three
versions of the model a data set was generated.

Based on the first criterion (degree of heterogeneity of
variance) and given that the coefficient of variation (C.V.)

Ikur{o? is equal to where O<6<O.25, two alternatives

_1_
1-46’
were chosen. The first reflects the case of heterogenous
variances based on a large value of 6:0.2 with C.V. = 1.41.

The other alternative reflects the case of not so heteroge-

 

 

 

 

 

88
nous variances (i.e., small 6=0.02) with C.V. = 0.21. This
data set is used to test the applicability of holding 6 to a
very small value (close to zero) to reflect the homogeneity of
variance.

For the second criteria, number of groups (k) generated
three values of k, 100, 40 and 15, were chosen. These values
reflect large, moderate, and small numbers of groups.

The combination of the three versions of the model with
the two choices of the first criterion and the three group
sample sizes in the second criterion creates 18 data sets.
Variations of the values of the other parameters (like A, 12,
and of) in the model were arbitrarily controlled based on how
many between- and within—unit variables exist in the model
(i.e., different versions of the model). Table 5.1 presents
a summary of the three models used with eighteen data sets.
Except for the intercept, regression coefficients in these

three models are assumed to be fixed effects.

 

 

 

89

Table 5.1

Models used in generating the artificial data

 

 

 

 

 

 

 

 

 

 

 

 

 

Model (1) k=1oo
Heterogeneity
case _
(6=0.2) k'40
Yj = Zj); + ljUj + Ej I fOI j = l, ',k k=15
1 . k=100
_ Homogeneity
Zj = : , and A = Y. case
1 (6=0.02) k=40
k=15
. k=100
Model (2) Heterogeneity
case =
_ A f <0=0-2> “0 ......
1 Wj le Yo k=100
. Homogeneity
Zj : : I and A : Yl case k=40
1 W5 X%J p (6—0.02)
k=15
k=100
Model (3) Heterogeneity
f case k=40
Yj = sz + ijUj + ej, or (6=0.2)
(Yo. k=l5
71
l le sz X1j1 X1j2 X1j3 72 k=100
Zj : . I . : 1 and l: B I H mogeneity
1 o
l le sz anjl anjg anj3 p2 (eggsgz) k=40
.ij

 

 

 

 

 

 

 

 

k=15

 

 

 

90

Data Creation

Given the model in 3.1, all eighteen data sets were
generated in two general steps. The first step covered the
generation of data for the group level variables. The second
step covered the generation of data within each group based on
the data generated for the group level variables in the first
step. These two steps reflect the inherent nesting of the two

levels of the model.

step one - Between-group data

The following random vectors with k elements (k = number
of groups) in each vector were generated:
1. Within-group sample sizes nj vector, for j = 1, ..., k.
First, k random numbers were generated from a UNIFORM
distribution between 5 and.60. Each number was rounded to the
nearest integer to represent a sample size nj for each of the
k groups. Variation in sample sizes are considered to include
a wide range from small (nj=5) to large (raj-=60). The

following is an example of sample size vector:

Group Sample size nj
1 n1=40
2 n2=15

 

91
2. Within-group variances {0;} vector. For pre-specified
values of 6 and 0.2 the variances {03-} were generated from their

respective distribution.

f(o§-|6,o§) ~ I“1(a,[3), where
(5.59)
a z i = 2_6
26 ’ a:

3. Based on 3.9, the k elements of the {Uj} vector were
generated from a normal distribution with mean equal zero and
variance equal 12.
4. Now, q vectors of W = [W1, W2, ---, Wq] and p vectors of
2?: [§1,)?2,~~,5{_p] each has k elements were generated from a
multivariate normal distribution with mean vector of zero and
pre-specified variance-covariance matrix 2. The values of
the p vectors of §= [§1,)?z,m,fp] , were then used, in step
two, as means of the p vectors for the within-group predictors
Xj = [Xj1,Xj2, -~,ij] in the design matrix Z in 3.2.
Step two - Within-group data

The generated data in the second step is for the within-

group variables. It depends on the values of the between-

group variables generated in the first step. For each j,
where j = 1, ~--, k groups, an p variables ofXj = [Xj1,Xj2,o--,ij]
within-group predictors and an error term ej were generated.
The number of observations for these variables within each

group is equal to the sample size nj generated in the first

data set.

. -._—._———_ —____._ —.—

 

92

1. Based on the model in 3.1, e1.j ~ N(O,Io§-) . Therefore,
given {03-} for j = 1,2,---,k, which are generated in the first
step, nj values of eij, for i = 1,2,---,nj, are generated for

every one of the k groups. The eij's are normal with mean

zero and variance 0?.

2. For each of the k elements of the p vectors of
I? = [§1,)?2,---,)?p] generated in the first step, 12]. vectors of
observations were generated for each of Xij = [Xv-1, X1721“! Xijp] ,

where 1' = 1, 2, ~--, nj and j = 1, 2, --~, k, from normal distribution
with mean fjp and a pre-specified value for the variance air,»-
3. For pre—specified values of the parameter Vector
A’ = [yu yl, ..., yq, [31, [32, ..., [3p] the values of the outcome variable
Y1-j in a particular data set were then generated by
substituting the generated values of all the between- and
within-group variables in the right hand side of equation 3.1.

In summary, one artificial data set for a model with one

between and one within-group variables, for example, can have

the following form:

1- k = 40 number of groups , nj ~ Uniform(5, 60)
2- 0:0.2, 03:50.
1 26
_ 2, ~ _1 -_- :I I a = __ I = —
3 a] I‘ (a 2.5 , [3 008) y 26 0.2

2 - Var(Uj) = 9.0 , Thus Uj ~ N(0,9.0)

 

__ 4 1 93
2 = COV(W,X) =( ), oi“, = 81 , thus W~N(0,4) ,

81) , and eij~N(0,o§-)

6- y°=5.0, 71:3.0, B=O.5, where p=q=l

7- Yij ~ N( 5+3wj+.5Xij+Uj , 03- )

Tables 5.2, 5.3 and 5.4 presents the true values of the
parameters used to generate the 18 artificial data sets used
in this study. Note that the within-group sample sizes nj is

a random sample from uniform (5, 60) for all the artificial

data sets.

 

94

Table 5.2

The true values of the parameters used to generate the
artificial data sets with k = 100

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

k = 100
e = 0.2 e = 0.02
MOdel (1) Model (1)
03:30 12:6.25 A=yg=6 03:15 12:9 A=Yg=5
Model (2) Model (2)
03:30 12:4.00 03:35 12:6.25
Y°=3 . 00 Yo=12 . 5
A=yl=l.50 A=yl=6,00
31:3-50 81:2.50
4.00 8.00
E: 2: ]
0.21 3.00 0.70 4.00
0i = 400 a; = 400
Model (3) Model (3)
03:30 12:6.25 03:30 12:3.0625
ry,,=8.00' 3158.00)
71:3-50 yl=3.50
A =YZ=2'7S A::72=2.75
01:1.00 01:1.00
82:3.50 82:3.50
03:0.75) _ _p3=0.75)
V4.00 “ 3.00
0.21 3.00 0.11 1.50
E==0.72 0.33 7.00 2==0.32 0.35 2.50
0.62 0.21 0.20 4.00 0.09 0.22 0.20 3.00
_0.30 0.20 0.21 0.11 5.00) _0.42 0.20 0.14 0.11 1.00“
2_ 2_
og=676 02:400 02=484 02:225 ofi—289 qa—loo

 

O 0 ' —'
12 is the variance—covariance matrix between the W s and.X s

 

95

Table 5.3

The true values of the parameters used to generate the

artificial data sets with k = 40

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

k = 40
6 = 0.2 6 = 0.02
Model (1) Model (1)
03:25 12:2.25 A=yy=8 03:15 12:9 X=Yg=5
Model (2) Model (2)
03:40 12:2.25 03:30 12:9
Yo=3.00 Y°=l2
A==yl=1.50 1: 71:4
fi1=3.50 61:2
_4b0 2_6.00 J
0.21 3.00 0.30 4.00
a; = 169 a; = 400
Model (3) Model (3)
03:25 13:2.25 03:15 12:6.25
Fy°=8.oo' 'y.=8.00'
Yl=3-50 Y1=3.50
=2.75 Y =2 75
A=Y2 A: 2
01:1.00 31:1.00
02:3.50 02:3.50
p3=o_75 _ _B3=O.75‘
F6.00 ' ‘ 2.00
0.11 3.00 0.13 2.50
2==0.32 0.35 3.50 2==0.12 0.15 1.50
0.12 0.22 0.20 4.00 0.09 0.05 0.10 3.00
0.30 0.20 0.21 0.11 5.00, _0.02 0.20 0.14 0.11 1.00,
2 2‘_
02:121 02:169 02:196 02:64 oa=121 ofi—lOO

 

2 is the variance-covariance matrix between the W's and X's

 

96

Table 5.4

The true values of the parameters used to generate the

artificial data sets with k = 15

 

 

 

 

 

 

 

 

 

 

 

 

 

 

= 15
6 = 0.2 0 = 0.02
Model (1) Model (1)
03:10 12:1 A=Y°=5 03:10 t2=300625 l=yo=5
Model (2) Model (2)
03:20 12:16 03:10 12:36
Yo=8'00 Y°:8.0
A=yl=3.50 1:71:35
51:1.50 61:1'5
3.00 3.00
2: 2: ]
0.11 2.00 0.11 2.00
2 2 _
ox1 = 49 0x1 — 49
Model (3) Model (3)
03:10 12:2.25 03:30 12:4
'Y.=8 . 00' (756 . 00'
7123-50 Y1=4.50
72: .75 lY2=2.50
A: A: _
pl: .00 81—2.00
82:3.50 82:3.50
p3=0.75) . _BB=1.75) 1
'6.00 ' 2.00
0.11 3.00 0.13 2.50
2==0.32 0.35 3.50 2==0.12 0.15 1.50
0.12 0.22 0.20 4.00 0.09 0.05 0.10 3.00
0.30 0.20 0.21 0.11 5.00) _0.02 0.20 0.14 0.11 1.00,
2 2.. 2..
02:121 og=169 ag=196 ofi=25 qh—49 ofi—ie

 

 

9 ' ' —'
2} is the variance-covariance matrix between the W s and X S

 

 

97

Initial Estimates for Gibbs Sampling

We recall from the section titled "Empirical Application
of Gibbs Sampling" on page 87, that initial estimates for A,
of, {Ufi' and. (a? are required to start the iteration process
of Gibbs sampling. Reasonable initial estimates for these
parameters could be their empirical Bayes estimates. However,
to avoid the dependency in comparing empirical Bays estimates
of these parameters with those produced by Gibbs sampling, the
formal estimates will not be used as initial estimates for
Gibbs sampling. Alternatively, least square estimates are
used as initial estimates for these parameters as follows:
1- Using the data Y and the design matrix Zlin 3.2, The
parameter A. is estimated by A“) =(J§1Z§Zj) 124le. The
superscript (O) is used with the parameter A to
represents the initial estimate of the parameter.
2- Given the data Y and the computed value of A”) from the

previous step, each element of the set {U9 is estimated

n- n
m)__ l 2 _ i m)
b U- - —— Y- Z-A ).
Y 3 nj 14 1 id 3
3- Similarly, each element of the set {0? is estimated by
03(0) = nj‘l(yj _ 21.1w) — Uj‘°>)/(yj — zjl<°> — Uj‘”) . The values

of A”) and (#0) were computed from the previous two

steps.
4- Finally, given the set of estimated values {ojmn from
step 3, of is estimated by{fhe harmonic mean H'of {ofiww
as follow 02(0) = k 2092(0)) .

The above estimates can be thought of as sampled values

 

 

 

98

of the parameters from the current approximation to their
marginal posterior distributions. These sampled values are
then used to simulate the distributions for the rest of the
parameters. The simulated distributions however, are poor
ones because they are based on poor estimates of A, of, {Uj},
and {03-}. Going through several iterations should improve the
approximations of all the marginal posterior distributions.

Another parameter estimated from the data is 6. While
the estimate of this parameter was not needed to start the
iteration process of Gibbs sampling, its estimated value,
however, is used for empirical confirmation of its Gibbs
(posterior mean) estimate. Except for 6, empirical Bayes
estimates for all other parameters of the model in 3.1 can be
easily found and compared to their counter part estimates from
Gibbs sampling. Therefore, an estimate of 6 from Gibbs
sampling can be compared to a one that is based on the log
transformation of the estimated group variances (Raudenbush
and Bryk, 1987). That is, when S]? for j=1,~-,k is used to
estimate the group residual variance with vj degrees of

freedom, the transformation dj, where

1 _ 5.60

is used to estimate 6]. =—:—109(O§-). Each dj is approximately
distributed as N[5j, 1/ (2Vj—1)], for j=1, ..., k. Similarly, using
the assumption in 3.5, where (IE/(60%) is distributed as a chi-

square variable with 6‘1 degrees of freedom, 5. is distributed

J

as NHL-g) , where A is the average of all 6]..

 

 

 

')

 

 

 

 

99
A total variance estimate of dj approximately consists

of a sampling variance and a residual parameter variance

 

 

Var(dj|A) = Var(dj|6j) + Var(6j|A) , and
f<d - (1)2 (5°61)
E[1 j J = l + .9

When A and 1/(2VJ71) in 5.67 are approximated by

A

A = k“1fdj and V = k'lfz—t— respectively, 6 can be estimated by
1 j

 

f(dj — A)2 (5.62)
6“” = 2 1 k — V .
— 1

Assessing the Heterogeneity of Variance

The primary interest of this research is to study the
analysis of hierarchically structured data when there is

evidence of heterogeneity of variance. Therefore, a measure

for assessing homogeneity of variance was needed. Based on
the results from the previous section, 6]. = glogkfi) is

approximately distributed as N(A, g) , where A is an average
of 6]. for j=1, -~-, k. Therefore, a natural way of assessing the

homogeneity of group variances is to test the hypothesis

Ho:6j = A , j=1,---,k by the statistic

A 5.
zfvjm] A)2 < 63>
j=1

which has a large sample chi-square distribution withk-l

degrees of freedom (Raudenbush and Bryk, 1987).

 

 

 

 

 

 

100

Random number generation

Three computer programs, using the FORTRAN language, were
written to implement the generation of the data and the
iteration process of Gibbs sampling. The first two programs
are for creating the artificial data sets in two steps. The
third program applies the iterative steps of Gibbs sampling.
In this program sufficient statistics are calculated from an
input data file (either the artificial data created by the
first two programs or the real data) and used to calculate the
initial estimates of the parameters A, of, HG},anui {0%, the
chi-square statistic for testing the hypothesis of homogeneity
of variance as well as the moment estimate of 6. The initial
estimates of the hyper-parameters are then treated as being
sampled from the current approximation of their marginal
posterior distributions and. they are 'used to start the
iteration process between six subroutines. Each subroutine is
written to sample one value of each of the parameters 6, 03,
12, A and one set of {Uj} and {03-} for j=1,---,k from its
corresponding conditional distribution. The iteration process
is then terminated when the number of iterations required for
convergence is reached.

The conditional distributions of the parameters involved
in the iteration process of Gibbs sampling have three general
parametric forms: normal, gamma and inverse gamma. The

International mathematical and Statistical library (IMSL),

 

 

 

 

 

101
Version 10 is used with the programs to generate random
numbers from these distributions. Sampling a random variable
from an inverse gamma distribution with certain parameters was
established by inverting a random variable that was sampled
from gamma distribution with the same parameters.

Only 600 observations are printed out for each of the
parameters as a final representation of its marginal posterior
distribution. To avoid the dependency between two successive
iterations, the 600 observations are drawn from the last 3600
iteration of the program in a systematic order with a 6

iterations jump between selected observations.

Specifying the Criteria for Convergence

Gelfand, Hills, Racine-Poon and Smith (1990) have listed
several procedures for checking convergence in Gibbs sampling.
One of these procedures, they recommend, is the overlay
plotting of the estimated density of the simulated
distribution at several points of the iteration process. One
can assume convergence when these plots become equivalent.
Another procedure is the use of a Quantile-Quantile plot. Two
equal size samples, each drawn several iterations away from
the other, are ordered and plotted. As the number of
iterations before drawing the two samples increases, the plot
of the ordered samples moves toward a 45° line as an

indication of convergence. The above two procedures were used

 

 

 

 

 

102
in this study for checking the convergence in approximating
the marginal posterior distributions. The number of
iterations required for convergence for model 3 with 6=0.2

and k=15 was used for the subsequent runs.

Real Data Analysis

Description of the variables available in the High School

and Beyond data set is given in table 5.5.

Table 5.5

Description of Variables in High School and Beyond data

 

Variable Description

 

Student Characteristics

MATHACH Mathematic achievement (outcome).
SES A measure of socioeconomic status.
MINORITY An indicator of minority status.
GENDER An indicator for females.

 

School Characteristics

SECTOR School affiliation, Public VS. Catholic.

MEANSES Mean Of SES.

 

 

 

 

 

 

 

 

 

 

103

Previous research (Raudenbush and Bryk, 1987) has
demonstrated the use of hierarchical linear model in studying
the effect of the organizational characteristics of schools on
dispersion in mathematics achievement (MATHACH). They have
used the model for estimating the residual dispersion in
mathematics achievement for each school using information from
the whole sample. The resulting estimates of the residual
dispersions are empirical Bayes estimates, which are
conditioned on ML estimates of the hyper-parameters of their
assumed distribution.

This study goes beyond getting a single estimate of the
residual dispersion for each school by approximating not just
the marginal distribution for the dispersions but also the
marginal distribution for every parameter in the 3.1 model,
including the hyper-parameters. Clearly, this will give the
researcher more flexibility in doing her/his inferences about
any parameter of the model. Table 5.6 presents the different

models used with the High School and Beyond data set.

 

 

 

104

Table 5.6

Models used with the High School and Beyond data set

 

    

BASE + SECTOR + MEANSES + MINORITY

+ GENDER + SES

 

2 MATHACH = BASE + SECTOR + GENDER + SES

 

3 MATHACH = BASE + MEANSES + GENDER + SES

 

4 MATHACH = BASE + SECTOR + SES

 

5 MATHACH = BASE + SECTOR + GENDER

 

6 MATHACH = BASE + MEANSES + SES

 

7 MATHACH = BASE + SECTOR

 

8 MATHACH = BASE + SES

 

9 MATHACH = BASE

 

 

 

 

 

 

105
Similar to the models in table 5.1, all school and
student level variables in any of the models in table 5.6 are
assumed to be fixed effects except the intercept, it is

assumed random.

 

 

 

 

I 4

 

 

 

CHAPTER 6

Results

This chapter provides a discussion of the results of the
empirical application of Gibbs sampling. The discussion
covers the results from the artificial data and the High
School and Beyond data set. .Answers to the questions
concerning the objectives of this study which were presented

in chapter 1 are provided through this discussion.

Required Number of Iterations for Convergence

Several procedures have been suggested for checking the
convergence of the iteration process of Gibbs sampling
(Casella & George, 1992; Gelfand, Hills, Racine-Poon and
Smith, 1990 and Iewis and Orav, 1989). In this study two
graphical criteria were used in deciding about the convergence
of the iteration process of the program. The first criterion
is the overlay plotting of the density function for each
parameter approximated by different samples. These samples
were generated from several runs of the program on the same
data set with different number of iterations in each run" The
second criterion is based on "Quantile-Quantile" (Q-Q) plots.
Observations from two different samples based on the same

model, but different number of iterations are sorted and

106

 

 

 

 

 

107

plotted against each other in.a scatter plota “Using the first
criterion, one can assume convergence when the graphs of the
density, from different samples, for the marginal posterior
distribution of a particular parameter become very close or
identical to each other. Using ithe second criterion,
convergence can be detected when the scatter plot of the
sorted observations from two samples forms a 45 degree line.

Model 3 from table 5.1 with 6=0.2 and k=15 was used as
a reference to decide about the appropriate number of
iterations required for convergence for all data sets. It is
assumed.that this model represents the scenario of the largest
number of iterations required for convergence because it has
the largest number of parameters and the smallest number of
groups, k=15, with.heterogeneous variances.

Using the specified model the program was run 4 times.
The total number of iterations in each run were 4600, 8600,
10600 and 15600 iterations respectively. Four samples, one
from each run of the program, with 600 observations in each
sample were generated. To avoid dependency between
consecutive observations within each sample (Casella & George,
1992) the 600 observations were systematically drawn from the
last 3600 iterations in each run. Thus, the first sample is
made of 600 observations drawn systematically from the last
3600 iterations after 1000 iterations in the first run. The
second sample is made of another 600 observations drawn

systematically" from. the last 3600 iterations after 5000

 

 

 

 

 

 

108
iterations in the second run. The third sample is also made
of 600 observations drawn systematically from the last 3600
iterations but after 7000 iterations in the third run. And
finally the fourth sample is made of another 600 observations
drawn systematically form the last 3600 iterations after 12000
iterations in the last run.

Figures 1 and 2 show the overlay graphs of the densities

for the marginal posterior distributions of the hyper-
parameters in model 3 of table 5.1. Four densities, each
based on one sample, were drawn for each parameter. These
densities are smoothed by the "kernel" method (Silverman,
1986), which superimposes a univariate nonparametric kernel
density estimator. The estimator shows areas where the
observations are most concentrated in the sample.
Figures 3 and. 4 jpresent. two Q-Q jplots for each. hyper-
parameter. One plot is based on the samples generated from
the runs of 8600 and 10600 iterations and the other plot is
based on the samples generated from the runs of 8600 and 15600
iterations.

Density graphs for the exchangeable parameters {U} and
{0? for j=fiL,~,15 of all the 15 groups in the chosen data set
are considered in deciding about the convergence of the
iteration process. However, only 4 groups, selected at
random, were presented.here for illustrative;purposeq IFigures
5 and 6 provide the overlay graphs of the densities of the

marginal posterior distributions and the Q—Q plots for a; for

 

 

 

 

 

 

109
the 4 selected groups. The overlay graphs of the densities
and the Q-Q plots for 09 are provided in figures 7 and 8.

In general the overlay graphs show that densities from
the four samples are almost identical for' most of the
parameters. Densities for 6, of, 12 and 63 based on the 4600
iterations sample show little divergence from the rest of the
samples. Considering the overlay graphs and the Q-Q plots of
all parameters, it is decided that 8600 iterations is
sufficient to achieve convergence. The Q-Q plots show that
observations from the sample based on 8600 iterations form a
45 degree line when they are plotted against the observations
of the samples based on 10600 and 15600 iterations.
Divergence from the 45 degree line occurs only for small
number of observations that fall at the tail of the
distribution for every parameter. These observations can be
thought of as outliers. They have much lower probabilities of

being in the distribution than the majority of the

observations.

 

 

 

 

 

110

Figure 6.1

Overlay graphs of the estimated densities of the marginal

posterior distributions for the hyper-parameters

6, 0.2, 12, and Y.

      

 

 

II—I'II—j

v I r _ 7 I I I I I f I 1 I , i. W ’
0.00 025 0.50 0.75 1.00 1.25 1.50 1.75 2.00 0 2 4 6 8 10 12 14 16 18
2
6 o

O

 

 

 

 

 

 

111

Figure 6.2

Overlay graphs of the estimated densities of the marginal

posterior distributions for the hyper-parameters

“(1, 131, [32, and B3

          

 

 

j

I I I l I I l I fl

I I I l I I I I .
2.8 28 8.0 3.2 8.4 8.6 8.8 4.0 0.95 0.97 0.99 1.01 1.08 1.05

7. fl.

 

 

 

j

. , . . I . - I I I I I I I ' I . I ' I '—|
8.40 8.42 8.44 8.46 8.48 8.50 8.52 8.54 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78

fig B.

 

112

Figure 6.3

 

Q-Q plots of the observation of the marginal posterior

distributions for the hyper-parameters

2
6, 0., 1'2,

 

2.00 ‘ ‘

1.75 L
1.60 —'
1.26 -
1.00 -'

8600 iteration

076-
060-
026-I

 

 

GOO

 

om'oz's 'oéo'av'a '160'12'5 '16:) '17's 200

1 0800 Iteration

...............

 

1e
1e -

14-

8600 iteration

 

2r '-

 

 

' I
0 2

.,..,.,...,.
4 510121416

6
1 0600 Iteration

18

 

8600 Iteration

 

 

 

1O

 

1O

0
I

8600 iteration

‘4
l

/

v

 

 

 

v

I ' I

7 s 6
1 060 Iteration

10

and y.

2.00

1.75 '-

8600 iteration

0m ' l ' I ' I ‘ I ' I ' I ' I
0.00 026 050 0.76 1.00 1.25 1.60 1.76 2.00

8600 iteration

8600 iteration

8600 Iteration

1.60 -
1.25 -. _.. -'

1.00 - 'f .‘ .--

oms - l
050 .‘ I.
025 i :

 

 

 

 

“I 5600 Iteration

 

13
13-.
14‘:
12$
10:

3.‘

a:
4L
2.

 

 

 

D

D

, - .
- n
‘ .
."

u' '

.

.I'
.—
_ .
.—
.
v v

'2' '4 '6 5'10'1'2'1'4'1'515
15600

Iteration

 

 

 

 

5 .1 ' é ' 9'
15600 Iteration

‘Y.

I . l

 

‘ I
.
,.
.4
/ .
" II" -
I

 

 

 

I v r ' I

10

. 7 9
15600 Iteratlo

 

113
Figure 6.4
Q-Q plots of the observation of the marginal posterior
distributions for the hyper-parameters

71' 51! 62! and B3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4‘0 1 . . . M . . . .
as r 1
C .- ‘ 95 ‘ .I
:3 j E r. a
m 4—4
\_ $6 - to so — —
31> 6
— 4:
8 3.4 - O 94 ~ _
LO - ‘3
co 3'" g
5.2 . 32 _ . _
.. 3'0 3'2 6'4 ab 32) 4.0 as 3'2 alt 0'0 (1'9 40
1 0600 Iteration 1 5600 Iteration
B. 131
1.06 I . . - 1.06 . . . .
1°“ ‘ r" 1.0.? - _r‘ _
_i 8 H!
i". *5 ,‘If
1‘01 - r-‘J'r » if) 101 - "Jr _
.J"‘ t J"
099 - [I - O 099 _ r" _
[7‘ Q .‘
_.v| (O - _J
097 ’ . CD 5.97 _ ‘ _
095 . . 1 v 095 . I I i
006 0.97 0.90 101 1.03 1 ( 6 o; 6 097 099 01 1.03 1. )5
1 0600 iteration 1 5600 Iteration
B2 I32
354 . . . . a“ . . . . .
362 — 352 _ ,- .
C :‘ C - .:
:8 860 79 360 - —
c6 ' <0
._ t.
0) a)
:I: 343 - — t 64.9 — —
o O - .
8 846 ~ _ - 8 3,45 _ ." _
CD CD _
544 - 3,44 q _
312 I v I l ' 642 u l I l l
5.1 2 644 3.46 646 550 as 5.5 4 5'42 544 5445 5-45 55°. 552 4
10600 Iteration 1 5600 iteratIon
I33
076 0.73 . . . . . .
077 _ 077 "‘ I
' x
8 are — 8 97° ‘ 1"“
A: a ,.I
g 075 - a 0.75 _ JP. _
:6: t 1‘..-
[A74 — C) 0‘74 ' Ir" ‘
8 8 .
— on — _5' -
8 073 CO :1."
072 - 072 1
071 — 071 i I . . I . .
Q 1 ,3 o: 'I 012 0.7: on or: on 077 Q
1 5600 Iteration

 

 

 

 

 

 

114

Figure 6.5

Overlay graphs of the estimated densities of the marginal

posterior distributions of selected groups' variances {0?

 

"'|'l'l' 'I'T‘] |"""|""""l
4 6 8 1012141618 20 22 3 8 13 18 23 28 33 88

 

l—rl'l'l'l'l'l'l'l‘l'l'l
15 25 85 45 55 65 75 85 95105115125

 

2

2
015

 

8600 iteration 8600 iteration 8600 iteration

8600 iteration

115

Figure 6.6

Q-Q plots of the points of the marginal posterior

distributions of selected groups' variances {03-}

 

 

 

v

‘2'! 1‘2 1'5 2'0 2'4 23
10600 iteration
2
0' 13

 

 

6 6 1'2 135 1‘3 2'1
10600 iteration

2
0'15

1 I l l

26354656657586
10600 Iteration

    

8600 iteration 8600 Iteration 8600 iteration

8600 Iteration

 

 

7 9 11 1G 16 17 19 21
15600 Iteration
2
O- 9

l I I

 

 

 

/.I—.-
I ' 32

 

 

é 1'2 1'5 20 2'4 2}
15600 Iteration

2
0'13

B 0 12 16 19 21
15600 Iteration

2
G15

  

aaaAddUflfi'lafldQfi
15600 Iteration

 

 

 

 

116

Figure 6.7

Overlay graphs of the estimated densities of the marginal
posterior distributions of selected groups' intercept errors

{Uj}

 

 

 

 

 

 

 

 

8600 iteration

8600 iteration 8600 iteration

8600 iteration

 

117

Figure 6 . 8

Q-Q plots of the points of the marginal posterior

distributions of selected groups' intercept errors {Uji

U3

5 #lnggnglA

lAlA

LL

 

 

 

 

 

5- L
1 .
b~ . -I-
81 1". i;
at ..
11 F-
01 f
-13 ..., I
-2 '3.0 ._
«I t
-*4'Tfi-'27—-11'01‘;Iév4 a'o
10600 Iteration
U9
4 A L . L1 '
‘ r
3.1 l; -
. .I' ,_
2..
I t
I- I
J
01 E
1 J"
- .J . -'
1‘ _ i-
-2 I

 

 

 

.31 v 6 . 1, . é ,
10600 iteration

5

IL

 

 

 

 

 

 

 

 

a 441 1.1 .......
21 .P
I1 i
oi t
-11 r
.2: i
‘3‘ I")! i-
-‘i i-
—6.,".t, ,.,.,I,.
-5 -2 -1 0 1 3
10600 Iteration
U15
6 . LI i .......
5-1 r
‘1 ..x' r
31 ' F
2 >-
Ii / t
o1 L
-‘I-: i.
41,. ' L
-a‘r.'..‘ fj,r,vh
-O -2 -1 O 5 a

. 1, vi .
10600 iter

3 4
ation

8600 Iteration

8600 Iteration

8800 Iteration

8600 Iteration

 

 

6—
1
‘-

O
_l

R.-
1-4
i

o-
-‘A
-9;

<
_0—

I
1'"
/.
.
:".

 

.1

Ir!

 

‘jfi'filfi'j—VI—rlfi

 

lfitrl
-4 -8 -2 -1

'6 ' i '2 a '1 T
15600 Iteration
U9

. l A l A l L l A A

Y

 

4

81

2.I

or
-14

-2

II"...
If
.1 f,
2

fir l—I T Y 'fi_

' l

 

 

:3

1

v t‘v’fi T

-1 ' 6 ' i 2 a
15600 Iteration
U

13

 

A l A l A

l A l A

 

5

~21 ./

l'l'l'l'l']

 

«5231' 31'1'6 ' 1 '2
15600 iteration

 

I

U15

 

 

 

 

a A l LJ_Ll A l A I A l A l A l A
L
5- . -
4i .- t
91 t
I I
2-4 _
It i
O- ..
-1. ../ t
-21 _' r
-I I-

‘9 ‘fii fir ' r '
—e -2 -1 O o

1 2 '7: ' Ii ' 5
15600 iteration

 

 

 

 

118

Artificial Data

There are 18 data sets (see table 5.1) created by the
combination of the three models, the two values chosen for 6,
and the three choices of k. The marginal posterior
distributions for 6, of, “:2, A, {03-}, and {Uj} where j=1,~--,k
in each model are calculated for each data set. Posterior
means and standard deviations are calculated from the
approximated marginal posterior distributions. Information
about the exchangeable parameters {03'} and {Uji are presented
for only 10 groups selected in systematic sampling for data
sets with k=1oo or 40. Information about these parameters in
data sets with k=15 groups are presented for only 8 groups .
That is, Iknr.k=100, groupl, groupll, group21, ..., group91
were presented, for k=40, groupl, group5, group9, ..., group37
were presented, and for k=15, groupl, group3, group5, ...,
group15 were presented.

Empirical Bayes estimates for the regression coefficients
(i.e., A) with their standard.error'of‘estimateS'were obtained
for each of the 18 data sets using the HLM program. Maximum
likelihood (ML) estimates for the variance components, t2 and
03 were also obtained from the HLM analysis. (Note that the
HLM estimate of of represents a pooled within-group variance
ML estimate.)

Because of the large cost of analysis, the limited time

and resources available, each of the 18 data sets represents

 

 

 

 

 

 

119
just one sample generated from.a population that is different
from the other data sets' populations. Therefore, there is a
chance that any of the data sets might not be a good
representative to its true population and the associated
parameters. This has the implication that estimates from a
given data set are based on one sample, which limit our
assessment of error. A more complete study would involve
generating multiple data sets from the same population.
Estimates of a given parameter can then be obtained from each
generated sample via Gibbs sampling. Given these different
estimates, mean squared errors for example, could be

calculated for the given parameter.

High School and Beyond data

As pointed earlier, this data is from a probability
sample of 160 U.S. high schools. Table 5.5 describes the
variables in this data set. The outcome variable used in all
of the models is a standardized mathematics achievement score.
Nine models (see Table 5.6) were applied to HSB data.
Variations between these models are based on the number and
type (school vs. student) of variables. Variations in school—
level variables are expected to influence mainly the estimated
value of 1:2, while variations in student-level variables
affect all estimates of variance components. Gibbs sampling

was used with 8600 iterations for all of the models. Similar

 

 

 

 

 

120
to the artificial data, information about the exchangeable
parameters {03-} and {Uj} are presented for only 10 schools
selected in systematic sampling. Empirical Bayes estimates
for the regression coefficients in all models were also

obtained via the HLM program.

Comparing Estimates of Variance Components

Part A of Table 6.1 to Table 6.18 provides estimates for
the variance components of, 6, tzanuito? for each of the 18
artificial data sets. Variance component estimates for the
High School and Beyond data set are provided in part A of

table 6.19 to table 6.27.

Estimates of of

The hyper-parameter 03 represents the typical value that
any of the parameters {03} for j=1, ---, k can take in their prior
distribution. One objective of this study is to find a way of
estimating that typical value without assuming homogeneity of
variance. Taking the fully Baysian approach and letting of
be random facilitated the:derivation.of its:marginal.posterior
distribution as well as the marginal posterior distributions
for each of the groups' residual variances.

The fact that the marginal posterior distribution of of

is not far from being symmetric even when k is small (see

 

 

 

 

121

figure 6.1), makes the difference between the posterior mean
and the posterior mode estimates insignificant. If the
marginal posterior distribution of of is far from being
symmetric, it has to be positively skewed, and the posterior
mean is larger than the posterior mode» ‘With.this in.mind, we
find that HLM estimates of of are larger than Gibbs posterior
mean estimates of of in all data sets including high school
and beyond data set. When they are compared to the actual
values of of (part A table 1 to table 18) HLM estimates of of
are still found to be larger than the actual value of the
parameter in all data sets. Furthermore, the over estimation
of of in HLM analysis is more pronounced for data sets with
extreme heterogeneity of variance (see tables with 0:0.20).
This indicates that the HLM estimate of the within-group
residual variance of is over estimated when there is clear
evidence of heterogeneity of variances.

One obvious implication of over estimating the within-
group residual variance is the effect on empirical Bayes
estimates of the regression coefficients. As we see in
equations 1.5 to 1.7 and under the homogeneity of variance
assumption (i.e. o§=o§ , .j=1,m,k) that larger value of of
causes Aj in equation 1.5 to be smaller therefore giving less

A

weight to B in obtaining B}. Another implication stems from

3'
the homogeneity of variance assumption itself, where some of
the groups with smaller actual values of their residual

variance 0% are assumed to have an inflated residual variance

 

 

 

 

 

 

122

which is equal to of because of the homogeneity of variance

assumption.

Estimates of 6

The hyper-parameter 6 represents a scale parameter of the
prior distribution for the exchangeable parameters {03-} for
j=14-u,k. It is inversely proportional to the concentration
cflfto? around their typical value 03. Finding the posterior
distribution of 6 and its estimate provides an answer to the
question: “How precise are the posterior mean estimates of the
exchangeable parameters.{¢? in estimating the typical value
of?" The marginal posterior distribution of 6 is fOund to
cover its actual value within 95 percent confidence interval
in all of the artificial data sets. In.most of the data sets,
and especially in the High School and Beyond data set, Gibbs
estimates of 8 are found to be of the same magnitude of its
moment estimates, which.areibased.on the log transformation of
the residual variances (see chapter 5 for derivation).
Although.they are close in‘theirtmagnitude, Gibbs estimates of
6 derived from data sets with k=15 are found to be larger
than their moment estimates (see tables with k=15). This
finding is expected because when the number of groups is
relatively small the marginal posterior distribution of 8

becomes more positively skewed which causes the posterior mean

(Gibbs estimate) to get larger.

 

 

 

 

 

123

Estimates of t2

The hyper—parameter t2 represents the variance of the
prior distribution of the exchangeable parameters {U}; In
general, Gibbs estimates of 12 are found to be of the same
magnitude of the HLM estimates for both the artificial data
sets and the High School and Beyond data set. However, in
most of the data sets, Gibbs estimates of t2 are found to be
slightly higher than HLM estimates. In fact, Gibbs estimates
of 12 for artificial data sets with small number of groups
(k=15) are found to be noticeably higher than HLM estimates.
That is because as the number of groups get smaller the
marginal posterior distribution 12 becomes more positively
skewed causing the posterior mean estimate of 1:2 (Gibbs

estimate) to be larger than the posterior mode of HLM.

Estimates of kfii

One advantage of this study is the ability to obtain the
marginal posterior distribution for each of the within—group
residual variance {03-} for j=1, ---,k. In all of the artificial
data sets the actual values of {03:} were found to fall within
the 95 percent confidence intervals which were derived from
their marginal posterior distributions. The posterior means
(Gibbs estimates of {031) , the actual values of {03-}, the

standard deviations, the coefficient of variations (C.V.) , and

 

 

 

 

 

 

124
the sample sizes, nj, for the selected groups in each of the
artificial data sets were presented. When compared to their
actual values, Gibbs' estimates<xfio?'were found to have the
same magnitude in most of the data sets. We also found that
Gibbs estimates of {03-} in groups with small sample size have
higher C.V. than estimates with larger sample sizes. This
indicates that estimates of {03'} which are based on small

samples are less stable than those based on a large sample

sizes.

Comparing Estimates of the Regression Coefficients

Part B of table 6.1 to table 6.18 provide estimates for
the regression coefficients in A’= Ho, 'er“'qu: 131,-",Bp] as
well as estimates for {Uj} for j=1, ~-,k in the selected groups
for the 18 artificial data sets. Regression coefficients and
{Uji estimates for the high school and beyond data set are
provided in part B of table 6.19 to table 6.27.

As mentioned earlier, unlike Bayes estimates, which are
derived from their marginal posterior distributions via Gibbs
sampling, empirical Bayes estimates it are conditioned on
knowing the true values of the variance components of and 12.
In practice, these components are unknown and need to be
estimated. One obvious problem with empirical Bayes estimates
is that they do not account for the uncertainty in estimating

the variance components. Consequently one can ask: "How do

 

 

 

 

 

125
Gibbs estimates of the regression coefficient differ from
their empirical Bayes estimates?" And "How do inferences
about regression coefficients change when taking into account
the uncertainty about the estimation of variance components,

especially when there is a heterogeneity of variance?"

Estimates of A

The vector A has two types of regression coefficients.
The first includes {7“ army, for the intercept and between-
group variable effects. Variables like SECTOR and MEANSES in
the High School and Beyond data are examples of this type.
The second includes {BLV~,flQ, for the within-group variable
effects. Variables like MINORITY, GENDER and SES in the High
School and Beyond data set are examples of this type. Except
for the intercept, between— and within-group regression
coefficients are assumed fixed for both the artificial data
sets and the High School and Beyond data set.

In general, Gibbs estimates (posterior means) and HLM
estimates (posterior Modes) of the regression coefficients are
very close to the actual values of the coefficients for the
artificial data sets. There are few cases where the two
estimates are noticeably different from the actual value of
the coefficients, but when compared to each other, Gibbs and
HLM estimates are found to be very close in their values. In

part B of table 6.8, for example, the actual value for ya is

 

 

 

 

 

126

set to be 3.00, while: Gibbs and. HLM: estimates of that
coefficient are 2.598 and 2.609 respectively. Similarly, part
B of table 6.17 shows that the actual value for Y1 is equal to
3.50, while Gibbs and HLM estimates of that coefficient are
4.603 and 4.579 respectively; Large deviations in the
estimates of the regression coefficients from their actual
values are very few and seem to have a random pattern that
could be attributed to the selected random samples. Similar
findings are obtained when using high school and beyond data
set. Gibbs and HLM estimates for the regression coefficients
are found to be of the same magnitude in all the models for
High Schools and Beyond data set.

A standard deviation computed from the marginal posterior
distribution produced by Gibbs sampling for each regression
coefficient is compared to the standard error of the estimate
from HLM. In general we find that both the posterior standard
deviations of Gibbs sampling (standard error of Gibbs
estimates) and the standard error of the estimates from HLM
are very close in their values and in some cases they are
equal. Since Gibbs estimates of the regression coefficients
are close in their 'values and. their standard error of
estimates to their counterpart estimates from.HLM, inferences
(t-test statistics) about. these. coefficients 'using' Gibbs
estimates are not drastically different from the HLM
inferences. This is found to be interesting especially for

the cases were there is a:noticeable heterogeneity of variance

 

 

 

 

 

127

and of is being over estimated in HLM analysis.

Estimates of {Uj}

Based on their marginal posterior distributions, Gibbs
estimates (posterior means) of {Uj}, where j=l,u-,k and the
standard deviations for the selected groups in each of the
artificial data sets are presented with their actual values.
. Because Gibbs estimates of {Uji are derived from one sample in
each data set, qualities of the estimates are difficult to
evaluate. However, in most cases we find that these estimates
are close in their magnitude and their sign to the actual

values of {Uj} .

 

 

 

128
Table 6.1-A
Generated data: Model 1, 6 = 0.20 and k = 100

Variance components estimates

Hyper-parameters

 

 

6 of 1'2
GIBBS““ 0.242 GIBBS““ 29.700“ 7.079
MOMENT EST. 0.317 HLM 48.512 7.550
ACTUAL 0.200 ACTUAL * 30.000 6.250
GIBBS s.0.* 0.042 GIBBS S.D. 2.313 1.251
CHI-SQUARE 809.349 D.F. = 99 P-value = 0.000

Within-groups residual variances

 

 

 

01 011 0:1 0:1 0:1
GIBBs*“* 12.144 136.379 29.342 16.480 ' 35.496
ACTUAL * 7.018 147.361 23.720 19.950 73.466
GIBB§*S.D. 5.239 26.185 5.600 3.714 20.162
C.V. 0.431 0.192 0.191 0.225 0.568
nj 12 55 59 43 9

2 2

0.3.1 0:1 071 081 0:1
GIBBs*“* 30.851 35.306 30.168 16.077 534.123
ACTUAL 26.693 27.382 30.924 11.768 508.802
GIBB§*S.D.* 6.804 16.158 6.200 3.518 202.966
C.V. 0.221 0.458 0.206 0.219 0.380
n. 36 9 49 42 14

.7

 

Standard deviation of the marginal posterior distribution
Pooled within—group variance estimate

Coefficient of variation

Posterior Mean

****

 

 

Generated data:

129

Table 6.1-B

Model 1, 8 = 0.20 and k = 100

Regression coefficient and random error CG estimates

Regression coefficient

 

 

 

 

 

YO
GIBBSH’" 5 . 978
HLM 5.905
ACTUAL 5_ooo
GIBBS s.D.* 0.293
HLM s.E.** 0,309
Random error CG
U1 U11 U21 U31 U41
GIBBs**** —0.750 1.399 0.627 —3.o44 -3.218
ACTUAL -0.531 2.102 0.289 —2.077 —3.807
GIBBS s.D.* 0.944 1.298 0.718 0.673 1.637
5%1 [Q1 £41 [Q1 £51
GIBBS**** —1.550 0.685 0.405 1.188 -0.793
ACTUAL -2.935 2.678 -0.795 1.525 -2.614
GIBBS 8.0.* 0.917 1.616 0.779 0.647 2.459

 

Standard deviation of the marginal posterior'distribution
H Standard error of the HLM estimate of the parameter

Posterior Mean

 

 

130

Table 6.2-A

 

 

 

 

 

 

 

 

 

Generated data: Model 1, 0 = 0.20 and k = 40
Variance components estimates
Hyper-parameters

6 t2
GIBBS**** 0.178 GIBBs**** 24 . 553" 3 . 248
MOMENT EST 0.181 HLM 34.201 3.524
ACTUAL 0.200 ACTUAL * 25.000 2.250
GIBBS S.D.* 0.052 GIBBS S.D. 2.657 1.124
CHI-SQUARE 229.232 D.F. = 39 P-value = 0.000
Within—groups residual variances
01 a: a: 0:3 017
GIBBs**** 23.474 111.119 36.078 32.839 31.858
ACTUAL 37.935 172.278 39.534 29.618 24.233
GIBBS S.D.* 8.440 37.033 8.530 8.794 11.685
C.v.*** 0.360 0.333 0.236 0.268 0.367
nj 12 17 34 29 16
2
0:1 0:5 0:9 0:3 037
* . 57
GIBBS *** 19.849 66.245 14.515 16.994 18 9
ACTUAL 15.690 53.687 20.660 12.568 16.566
GIBBS S.D.* 3.647 17.411 3.296 7.929 3.351
C.V.*** 0.184 0.263 0.227 0.467 0.177
nj 59 27 42 10 58
* Standard deviation of the marginal posterior distrlbution
" Pooled within—group varlance est1mate

Coefficient of variation
Posterior Mean

 

131

Table 6.2-B

Generated data: Model 1, 6 = 0.20 and k = 40

Regression coefficient and random error CG estimates

Regression coefficient

 

 

 

 

 

YB
GIBBs**** 7 . 7 80
HLM 7.800
ACTUAL 8.000
GIBBS S.D.* 0,349
HLM S.E.** 0.347
Random error CG
U1 U5 U9 U13 U17
GIBBs**** -1.804 -l.785 —2.445 3.204 2.490
ACTUAL -2.259 -2.884 -2.816 3.721 1.602
GIBBS S.D.* 1.143 1.620 0.897 0.948 1.168
U21 U25 U29 U33 U37
GIBBs**** —1.490 —0.469 1.913 —o.903 3.3%
ACTUAL -2.116 -0.850 1.459 —o.961 0.680
GIBBS S.D.* 0.660 1.182 0.646 1.03 .

 

Standard deviation of the marginal posterior distribution
H Standard error of the HLM estimate of the parameter

Posterior Mean

 

132

 

 

 

 

 

 

Generated data: Model 1, 6 = 0,20 and k = 15
Variance components estimates
Hyper—parameters

0 o. :2
GIBBS**** 0 . 309 GIBBs**** 12 . 755 1. 439
MOMENT EST 0.229 HLM 20.544“ 1.036
ACTUAL * 0.200 ACTUAL 10.000 1.000
GIBBS S.D. 0.181 GIBBS S.D.* 2.997 0.882
CHI-SQUARE 90.567 D.F. = 14 P-value = 0.000
Within-groups residual variances
01 03 a: 03 09
GIBBS**** 8.226 16.472 6.023 26.187 31.382
ACTUAL 10.489 14.213 4.842 17.763 35.425
GIBBS S.D.* 3.744 3.477 2.494 6.481 7.892
C.V. ** 0.455 0.211 0.414 0.247 0.251
I1]. 12 42 17 37 34
Oil 0:3 0.5.5
GIBBS**** 36.008 36.058 17.321
ACTUAL 37.924 35.642 18.583
GIBBS S.D.* 6.959 9.575 4.712
C.v.*** 0.193 0.266 0.272
nj 55 29 29
:* Standard deviation of the marginal posterior distribution
Pooled within—group variance est1mate

****

Coefficient of variation
Posterior Mean

 

133
Table 6.3-3
Generated data: Model 1, 6 = 0.20 and k = 15

Regression coefficient and random error CG estimates

Regression coefficient

 

Ya
GIBBS*“* 5.129
HLM 5.051
ACTUAL 5.000
GIBBS s.g.* 0.362
HLM s.E. 0.346

 

Random error [5

 

 

 

 

U 1 U3 U5 U7 U9
GIBBS*“* 1.068 -0.310 —0.854 -0.178 0.125
ACTUAL 2.184 —0.629 -0.526 —0.661 1.001
GIBBS S.D.* 0.722 0.583 0.615 0.699 0.737
U1 1 U13 U15
GIBBS*“* 0.801 —0.281 0.436
ACTUAL 1.491 -1.015 0.205
GIBBS S.D.* 0.713 0.799 0.689
* Standard deviation of the marginal posterior distribution

M Standard error of the HLM estimate of the parameter
Posterior Mean

 

 

 

 

 

 

 

 

 

 

134
Table 6.4-A
Generated data: Model 1, 6 = 0.02 and k = 10°
Variance components estimates
Hyper—parameters
6 0. 12
GIBBS**** 0 . 037 GIBBS**** 15 . 147 9 . 917
MOMENT EST 0.040 HLM 16.162" 9.837
ACTUAL * 0.020 ACTUAL 15.000 9.000
GIBBS S.D. 0.010 GIBBS S.D.* 0.599 1.531
CHI—SQUARE 210.674 D.F. = 99 P—value 0.000
Within-groups residual variances

Oi 011 0:1 031 0:1
GIBBS**** 13.155 14.107 17.298 12.285 14.344
ACTUAL 11.426 13.504 14.683 13.404 18.347
GIBBS S.D.* 3.241 2.280 2.721 2.270 3.791
C.v.*** 0.246 0.162 0.157 0.185 0.264

nj 12 55 59 43 9

2 2

05231 0:1 071 081 091
GIBBS**** 17.973 16.501 18.226 19.889 17.447
ACTUAL 17.665 14.297 20.474 18.957 15.502
GIBBS S.D.* 3.528 4.244 3.030 3.598 4.237
C.v.*** 0.196 0.257 0.166 0.181 0.243
nj 36 9 49 42 14
:* Standard deviation of the‘marginalzposterior'distribution

Pooled within—group varlance estlmate

****

Coefficient of variation
Posterior Mean

 

135
Table 6. 4-B

Generated data: Model 1, 6 = 0.02 and k = 100

Regression coefficient and random error CG estimates

Regression coefficient

 

Y.
GIBBS*“* 4.986
HLM 4.977
ACTUAL 5.000
GIBBS S.D.* 0.281
HLM s.E.“ 0.325

 

Random error CG

 

 

 

 

 

 

U1 U11 U21 U3 1 U41
GIBBS*“* -0.852 2.395 0.589 —3.340 -4.346
ACTUAL -0.637 2.522 0.347 -2.493 —4.568
GIBBS S.D.* 1.053 0.584 0.604 0.582 1.140
U51 Us 1 U7 1 U81 U9 1
. . —3.302
GIBBS*** -2.380 1.828 0.001 1 :12 _3 137
ACTUAL —3.522 3.213 -0.954 3.705 1.013
GIBBS S.D.* 0.745 1.240 0.658 . .

 

 

 

ginal posterior distribution

' ‘ the mar
Standard.dev1at1on.of meter

M Standard error of the HLM estimate of the para
Posterior Mean

 

 

 

 

Generated data:

136
Table 6.5-A

Model 1, 8 = 0.02 and k = 40

Variance components estimates

Hyper-parameters

 

 

 

 

 

 

 

 

0 of 12
GIBBS““ 0.026 GIBBS““ 15.961“ 13.485
MOMENT EST 0.020 HLM 16.777 12.800
ACTUAL 0.020 ACTUAL * 15.000 9.000
GIBBS S.D.* 0.014 GIBBS S.D. 0.934 3.455
CHI-SQUARE 62.049 D.F. = 39 P-value = 0.009
Within—groups residual variances
Ci 05 a: 013 017
GIBBS*“* 14.304 14.116 16.287 15.650 19.158
ACTUAL 13.883 11.068 16.891 12.937 20.906
GIBBS S.D.* 3.056 2.850 2.730 2.931 3.976
C.V.“* 0.214 0.202 0.168 0.187 0.208
nj 12 17 34 29 16
2 2
0:1 0:5 0:9 c‘33 037
17.387
GIBBS*“* 18.359 21.648 13.340 15.458
ACTUAL 16.387 22.196 16.551 11.647 16.846
GIBBS S D * 2 676 4.415 2.269 3.496 2-558
. . . 47
C.V.“* 0.146 0.204 0.170 0.226 0.1
nj 59 27 42 10 58
* Standard deviation of the marginal posterior d15tr1but1on
M Pooled within—group variance est1mate.

Coefficient of variation
Posterior Mean

 

137
Table 6. 5-B
Generated data: Model 1, 6 = 0.02 and k = 40

Regression coefficient and random error CG estimates

Regression coefficient

 

 

 

 

 

 

 

Y.
GIBBS"" 4.890
HLM 4.808
ACTUAL 5 , o 0 o
GIBBS S.D.* 0.589
HLM s.E.“ 0.581
Random error LG
U1 U5 U9 U13 U17
GIBBS*"* -4.450 -5.957 -5.805 7.551 4.774
- 18 -5.768 —5.631 7.441 3.205
ACTUAL * 4.5 0 873 1 178
GIBBS S.D. 1.218 1.045 0.876 . .
U21 U25 U29 U33 U37
** -2.143 —3.455
_GIBBS* * —3.773 -l.536 3.329 _ 3 _3 939
—1 699 2.919 1.92 .
ACTUAL * -4'232 .009 0 801 1.355 0.752
GIBBS S.D. 0.788 1. .

 

e marginal posterior‘distribution

Standard deviation of th parameter

H Standard error of the HLM estimate of the
*“* Posterior Mean

 

Generated data:

138

Model 1,

Variance components estimates

Hyper-parameters

6 = 0.02 and k = 15

 

 

 

 

 

 

e of 1:2
GIBBS””* 0.050 GIBBS 11.288“ 5.771
MOMENT EST 0.039 HLM 12.243 4.747
ACTUAL 0.020 ACTUAL * 10.000 3.063
GIBBS S.D.* 0.051 GIBBS S.D. 1.253 2.778
CHI-SQUARE 24.469 D.F. 14 P-value 0.032
Within—groups residual variances

of 03 O5 03 a:
GIBBS*“* 9.864 12.575 9.934 14.300 9.744
ACTUAL 9.080 11.087 8.585 10.694 8.563
GIBBS S.D.* 2.724 2.228 2.588 2.800 2.007
C.V.“* 0.276 0.177 0.261 0.196 0.206
nj 12 42 17 37 34

Oil Oi}; 02:5
GIBBS*“* 14.705 9.318 10.911
ACTUAL 17.018 6.224 10.676
GIBBS S.D.* 2.530 2.120 2.215
C.V.“* 0.172 0.228 0.203
nj 55 29 29
* ' ' inal osterior‘distribution

Standard dev1at10n of the marg .p

“ Pooled within—group var1ance est1mate

Coefficient of variation
Posterior Mean

 

 

 

139
Table 6 . 6-B

Generated data: Model 1, 0 = 0.02 and k = 15

Regression coefficient and random error CG estimates

Regression coefficient

 

Yo
GIBBS*“* 5.098
HLM 5.092
ACTUAL 5.000
GIBBS S.D.* 0.653
HLM s.E.“ 0.590

 

Random error CG

 

 

 

 

 

U1 U3 U5 U7 U9
GIBBS*”* 2.915 —0.840 —1.489 -0.815 1.228
ACTUAL 3.822 —1.101 —0.920 -1.157 1.751
GIBBS S.D.* 1.013 0.813 0.981 0.878 0.815
U1 1 U1 3 U1 5
GIBBS*“* 2.304 —1.564 0.660
ACTUAL 2.611 -1.776 0.25:
GIBBS S.D.* 0.825 0.830 0.
* Standard deviation Of the marginal posterior'distribut1on
M Standard error of the HLM est1mate of the parame er

Posterior Mean

 

140
Table 6.7-A

Generated data: Model 2, 6 = 0.20 and k = 100

Variance components estimates

Hyper-parameters

 

 

6 of :2
GIBBS*“* 0.216 GIBBS“”* 33.081“ 4.715
MOMENT EST 0.247 HLM 56.951 4.116
ACTUAL 0.200 ACTUAL * 30.000 4.000
GIBBS S.D.* 0.037 GIBBS S.D. 2.603 0.922
CHI-SQUARE 863.006 D.F. = 99 P-value = 0.000

Within-groups residual variances

 

 

 

 

 

 

Oi 011 0:1 031 0:1
GIBBS““ 16.329 41.471 24.863 25.736 21.571
ACTUAL 12.933 32.057 24.293 32.472 17.323
GIBBS S.D.* 7.083 8.120 4.858 5.365 11.159
C-V.“* 0.434 0.196 0.195 0.208 0.517
nj 12 55 59 43 9
2
051 061 071 0:1 091
GIBBS“" 70.313 23.057 414.811 43.554 25.085
ACTUAL 59.509 18.808 341.339 37.248 16.019
* 53 12.659 82.279 9.882 10.227
GIBBS S.D. 16.2 408
C.V.“* 0.231 0.549 0.198 0.227 0.
nj 36 9 49 42 14
* Standard deviation of the‘marginallposterior'distribution
“ Pooled within—group variance est1mate

Coefficient of variation
Posterior Mean

 

Generated data:

141

Model 2,

6 = 0.20 and k = 100

Regression coefficients and random error cg estimates

Regression coefficients

 

 

 

 

 

 

Y. Y1 Bl
GIBBS*“* 2.880 1.623 3.493
HLM 2.927 1.629 3.493
ACTUAL 3.000 1.500 3.500
GIBBS S.D.* 0.254 0.112 0.004
HLM s.E.“ 0.254 0.115 0.005
Random error CG
CG £51 £51 [51 C51
GIBBS*“* -2.281 0.692 3.165 -0.925 3.678
ACTUAL —2.371 0.245 3.314 —0.418 3.625
GIBBS S.D.* 1.118 0.885 0.667 0.782 1.295
[£1 [51 Uh £51 £51
— - . 0.607
GIBBS*“* 1.074 2.349 0.475 1 842 1 560
90 -1.741 —2.601 .
ACTUAL —0.691 2.2 1 109
GIBBS S.D.* 1.179 1.277 1.808 0.984 .

 

* Standard deviation of t .
H Standard error of the HLM est1mate of t

Posterior Mean

 

he marginal posterior distribution

he parameter

 

Generated data:

142

Model 2, 0 =

0.20 and k = 40

Variance components estimates

Hyper-parameters

2

 

 

 

 

 

 

 

 

6 O. t
GIBBS*“* 0.254 GIBBS“”* 39.420“ 1.336
MOMENT EST 0.269 HLM 74.133 1.041
ACTUAL 0.200 ACTUAL * 40.000 2.250
GIBBS S.D.* 0.066 GIBBS S.D. 4.913 0.636
CHI-SQUARE 319.225 D.F. = 39 P—value = 0.000
Within—groups residual variances

oi a: oi oi. 0%
GIBBS*“* 27.456 47.465 21.786 38.281 28.705
ACTUAL 35.656 34.510 14.777 24.560 27.952
GIBBS S.D.* 11.696 16.689 5.734 10.126 10.606
C.V.“* 0.426 0.352 0.263 0.265 0.369
nj 12 17 34 29 16

2 2

0:1 0:5 0:9 033 037
GIBBS““ 75.117 31.557 31.286 51.772 80.495
ACTUAL 76.792 52.517 39.124 85.585 90.053
GIBBS S.D.* 13.517 8.694 7.126 25.840 14.678
C.V.“* 0.180 0.275 0.228 0.499 0.182
nj 59 27 42 10 58
* Standard deviation Of the marginal posterior dlstr1bution
H Pooled within—group var1ance est1mate

Coefficient of variation
Posterior Mean

 

143
Table 6.8-B
e = 0.20 and k = 40

Generated data: Model 2,

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

 

 

Y. 71 pl
GIBBS*“* 2.598 1.412 3.489
HLM 2.609 1.575 3.478
ACTUAL * 3.000 1.500 3.500
GIBBS 8.3. 0.260 0.154 0.009
HLM S.E. 0.306 0.167 0.010
Random error CG
U1 U5 U9 U13 U17
GIBBS*“* —0.792 —0.311 1.713 0.274 0.077
ACTUAL -1.729 —0.517 2.001 0.608 0.416
GIBBS S.D.* 0.946 0.933 0.778 0.834 0.849
U21 U25 U29 U33 U37
GIBBS*“* 1.292 0.502 —0.286 -O.808 0.101
ACTUAL 2.023 1.551 —1.673 —1.383 —1.162
GIBBS S.D.* 0.818 0.827 0.695 1.084 0.841

 

Standard deviation 0 .
“ Standard error of the HLM est1mate o

Posterior Mean

f the marginal posterior distribution

f the parameter

 

 

 

144

Table 6.9-A

 

 

 

 

 

 

 

Generated data: Model 2, 6 = 0.20 and k = 15
Variance components estimates
Hyper—parameters

6 of 12
GIBBs“”* 0.206 GIBBS“" 20.978“ 26.121
MOMENT EST 0.136 HLM 32.276 21.553
ACTUAL * 0.200 ACTUAL * 20.000 16.000
GIBBS S.D. 0.112 GIBBS S.D. 3.699 12.460
CHI-SQUARE 94.728 D.F. = 14 P-Value = 0.000
Within-groups residual variances
01 03 a: 07 09
GIBBS*“* 19.828 14.153 32.483 18.957 18.681
ACTUAL 34.028 9.497 25.433 14.237 14.756
GIBBS S.D.* 8.519 3.090 11.107 4.483 4.261
C.V.“* 0.430 0.218 0.342 0.237 0.228
nj 12 42 17 37 34
2
011 013 015
GIBBS*“* 90.083 38.498 15.047
ACTUAL 75.116 27.346 13.960
GIBBS S.D.* 17.875 9.963 4.210
C.V.”* 0.198 0.259 0.280
n 55 29 29

J

:* Standard dev . . t
Pooled within—group var1ance estlma e
“* Coefficient of variatlon

Posterior Mean

****

iation of the marginal posterior distribution

 

145
Table 6. 9-B

Generated data: Model 2, 6 = 0.20 and k = 15

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

Yo Y1 pl
GIBBS**** 8.360 3.471 1.489
HIM 8.469 3.456 1.496
ACTUAL 8.000 3.500 1.500
GIBBS S.D.* 1.353 0.565 0.032
HLM s.E.“ 1.314 0.498 0.037
Random error CG
U1 U3 U5 U7 U9
GIBBS**** 1.199 —2.392 4.201 —4.396 3.303
ACTUAL 2.455 -2.504 5.519 —4.095 3.339
GIBBS S.D.* 2.744 3.532 2.052 1.553 1.535
U11 U13 U15
GIBBS**** -0.561 2.127 4.4%;
ACTUAL -0.809 2.790 0.47187
GIBBS S.D.* 1.952 1.866 1.

 

Standard deviation o _
M Standard error of the HLM est1mat

Posterior Mean

e of the parameter

f the marginal posterior distribution

 

146

Table 6.10-A

Generated data: Model 2, 6 = 0.02 and k = 100

Variance components estimates

Hyper—parameters

 

 

 

 

 

 

 

6 o. t2
GIBBS**** 0 . 017 GIBBS**** 35.539 7.145
MOMENT EST 0.011 HLM 36.606** 7.018
ACTUAL * 0.020 ACTUAL * 35.000 6.250
GIBBS S.D. 0.007 GIBBS S.D. 1.238 1.265
CHI-SQUARE 147.231 D.F. = 99 P-value = 0.001

Within-groups residual variances

Ci 0:1 5:1 031 521
GIBBS**** 35.437 43.223 35.010 32.267 35.255
ACTUAL 64.585 39.928 34.088 34.930 29.343
GIBBS S.D.* 5.708 6.309 4.788 4.625 6.840
C.v.*** 0.161 0.146 0.137 0.143 0.194
12]. 12 55 59 43 9

2

0:1 0:1 031 0:1 091
GIBBS**** 38.987 34.889 36.209 35.478 40.695
ACTUAL 36.968 43.337 26.444 30.017 49.005
GIBBS S.D.* 5.997 6.333 5.013 4.751 7.505
C.v.*** 0.154 0.182 0.138 0.134 0.184
n, 36 9 49 42 14
.7
* Standard deviation Of the‘marginal;posterior'distribution
M Pooled within-group variance est1mate

Coefficient of variation
Posterior Mean

****

 

 

 

 

Generated data:

Model 2,

147

Table 6.10-B

6 = 0.02 and k = 100

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

y. 7. Bl
GIBBS**** 12.402 6.114 2.494
HLM 12.423 6.103 2.492
ACTUAL * 12.500 6.000 2.500
GIBBS 5.3. 0.303 0.092 0.005
HIM s.E. 0.296 0.094 0.005

Random error CG

U1 U11 U21 U31 U41
GIBBS**** —3.010 0.762 3.940 —1.147 4.319
ACTUAL —2.963 0.306 4.143 -0.523 4.532
GIBBS S.D.* 1.538 0.893 0.824 0.871 1.665

U51 U61 U71 U81 U91
GIBBS**** 0.891 2.876 -1.608 —2.707 0.419
ACTUAL —0.863 2.863 —2.176 —3.251 1.949
GIBBS S.D.* 1.088 1.577 0.826 0.938 1.484

 

Standard deviation of the marginal posterior distribution
M Standard error of the HLM estimate of the parameter

Posterior Mean

 

 

 

 

 

 

 

 

148

Table 6.11-A

 

 

 

 

 

 

 

Generated data: Model 2, 6 = 0.02 and k = 40
Variance components estimates
Hyper—parameters

6 12
GIBBS**** 0.035 GIBBs**** 29.917” 7.263
MOMENT EST 0.032 HLM 32.254 6.721
ACTUAL * 0.020 ACTUAL 30.000 9.000
GIBBS S.D. 0.018 GIBBS S.D.* 1.948 2.043
CHI-SQUARE 77.033 D.F. = 39 P-value = 0.000
Within-groups residual variances
oi oi oi oi. oi7
GIBBs**** 26.545 33.325 27.540 41.700 29.654
ACTUAL 25.764 27.759 20.545 35.204 32.972
GIBBS S.D.* 6.314 7.933 5.057 8.774 6.750
C.v.*** 0.238 0.238 0.184 0.210 0.228
nj 12 17 34 29 16
2
0:1 035 0:9 033 037
32.862
GIBBS**** 26.641 23.280 26.810 27.529
ACTUAL 24.737 23.238 31.347 28.753 36.224
GIBBS S.D.* 4.310 5.255 4.580 7.098 5.119
C.v.*** 0.162 0.226 0.171 0.258 0.156
27 42 10 58
nj 59
* Standard deviation of the marginal posterior distrIbution
H Pooled within-group var1ance est1mate

Coefficient of variation
Posterior Mean

 

 

 

 

149

Table 6.11-B

Generated data: Model 2,

6 = 0.02 and k = 40

Regression coefficients and random error Cg estimates

Regression coefficients

 

 

 

 

 

Y. 71 [31
GIBBS**** 11.437 3.925 1.993
HIM 11.437 3.945 1.992
ACTUAL * 12.000 4.000 2.000
GIBBS 8'2: 0.446 0.205 0.006
HLM s.E. 0.453 0.190 0.006
Random error (G
U1 U5 U9 U13 U17
GIBBS**** —2.712 —0.781 4.310 1.104 0.671
ACTUAL -3.458 —1.034 4.003 1.215 0.832
GIBBS S.D.* 1.580 1.268 1.149 1.192 1.311
U21 U25 U29 U33 U37
GIBBS**** 4.498 2.814 —1.766 -2.587 3.32;)
ACTUAL 4.046 3.102 —3.346 -2.766 .
GIBBS S.D.* 0.884 1.012 0.952 1.756 0.996

 

 

Standard deviation of the marginal posterior distribution
“ Standard error of the HLM estimate of the parameter

Posterior Mean

 

150

Table 6.12-A

 

 

 

 

 

 

 

Generated data: Model 2, = 0.02 and k = 15
Variance components estimates
Hyper-parameters

0 of 12
GIBBS*“* 0.022 GIBBS*“* 12.254“ 57.949
MOMENT EST 0.007 HLM 12.661 50.305
ACTUAL * 0.020 ACTUAL * 10.000 36.000
GIBBS S.D. 0.023 GIBSS S.D. 1.190 26.081
CHI—SQUARE 13.851 D.F. = 14 P-value = 0.266

Within—groups residual variances
of a; 0% 03 a:
GIBBS*“* 11.264 14.066 12.354 12.040 13.844
ACTUAL 10.470 12.036 8.624 8.755 13.937
GIBBS S.D.* 2.423 2.263 2.194 1.866 2.199
c.v.*“‘ 0.215 0.161 0.178 0.155 0.159
nj 12 42 17 37 34
2

0:1 032.3 015
GIBBS*“* 13.817 13.743 11.593
ACTUAL 12.099 10.924 10.002
GIBBS S.D.* 1.990 2.324 1.927
C.V.“* 0.144 0.169 0.166
n. 55 29 29

J

* Standard deviation of the marginal posterior dlstribution

Pooled within—group variance estimate
Coefficient of variatlon
Posterior Mean

 

 

 

 

 

151
Table 6 . 12-B

Generated data: Model 2, 6 = 0.02 and k = 15

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

Ya 71 B1
GIBBS**** 8.485 3.661 1.493
HIM 8.792 3.471 1.493
ACTUAL 8.000 3.500 1.500
GIBBS S.D.* 2.162 0.710 0.025
HIM S.E.** 1.961 0.740 0.024
Random error (9
U1 U3 U5 U7 U9
GIBBS**** 1.899 —5.080 7.603 —6.932 4.577
ACTUAL 3.683 —3.755 8.279 —6.142 5.008
GIBBS S.D.* 3.494 4.444 2.636 2.207 2.178
U11 U13 U15
GIBBS**** -1.349 3.609 «.ng
ACTUAL -1.213 4.185 0.528
GIBBS S.D.* 2.482 2.314 2.

 

Standard deviation of the marginal posterior distribution
Standard error of the HLM estimate of the parameter

Posterior Mean

 

 

 

 

 

152
Table 6.13-A

Generated data: Model 3, 6 = 0.2 and k = 100

Variance components estimates

Hyper-parameters

 

 

3 of 12
GIBBS””* 0.294 GIBBS“” 28.457“ 6.244
MOMENT EST 0.364 HLM 81.360 5.433
ACTUAL 0.200 ACTUAL * 30.000 6.250
GIBBS S.D.* 0.050 GIBBS S.D. 2.542 1.097
CHI-SQUARE 1368.057 D.F. = 99 P—value = 0.000

Within—groups residual variances

 

 

 

 

2
oi oil 0:1 0'31 041
GIBBS*“* 16.973 23.737 11.577 26.024 20.127
ACTUAL 22.335 22.998 16.267 32.926 12.840
GIBBS S.D.* 7.233 4.809 2.224 5.754 9.068
C.v.*“' 0.426 0.203 0.192 0.221 0.451
nj 12 55 59 43 9
2
051 061 071 081 091
GIBBs*"* 16.712 23.608 35.107 46.319 30.096
ACTUAL 18.524 32.182 26.896 30.733 23.624
GIBBS S.D.* 4.018 11.181 7.283 10.596 11.934
C V “* 0.240 0.474 0.207 0.229 0.397
. . 9 49 42 14
nj 36
t ' I ' . I n
* Standard deviation of the marginal posterior dlStrlbuth
M Pooled within-group var1ance est1mate.

Coefficient of variation
Posterior Mean

 

Generated data:

153

Table 6.13-B

Model 3, 6 = 0.20 and k = 100

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

 

Y5 Y1 Y2 pl [32 B3
GIBBS*“* 7.694 3.769 2.856 1.001 3.491 0.748
HLM 7 896 3.743 2.782 1.003 3.494 0.742
ACTUAL 8.000 3.500 2.750 1.000 3.500 0.750
GIBBS S.D.* 0.307 0.131 0.180 0.003 0.004 0.004
HLM s.E.“ 0.299 0.135 0.186 0.003 0.006 0.005

Random error CG

U1 U11 U21 U31 U41
GIBBS““ -2.435 2.122 —2.429 1.271 2.331
ACTUAL -1.797 0.936 —3.079 0.558 3.604
GIBBS S.D.* 1.189 0.792 0.647 0.829 1.399

(%1 5&1 C51 [51 £51
GIBBS““ —1.346 1.175 3.343 3.266 0.095
ACTUAL —2.548 0.713 1.216 3.795 —$.§g§
GIBBS S.D.* 0.916 1.383 0.883 1.037 .
* Standard deviation of the marginal posterior distribution
H Standard error of the HLM est1mate of the parame er

Posterior Mean

 

154
Table 6.14-A

Generated data: Model 3, 6 = 0.20 and k = 40

Variance components estimates

Hyper—parameters

 

 

0 a. t
GIBBS“" 0.313 GIBBS““ 25.840 1.869
MOMENT EST 0.322 HLM 47.353“ 2.340
ACTUAL * 0.200 ACTUAL 25.000 2.250
GIBBS S.D. 0.093 GIBBS S.D.* 3.601 0.940
CHI—SQUARE 405.029 D.F. = 39 P-value = 0.000

Within-groups residual variances

 

 

 

 

01 05 a: 013 017
GIBBS*“* 14.152 35.146 52.942 16.728 28.825
ACTUAL 16.970 33.372 67.623 17.316 21.081
GIBBS S.D.* 6.476 13.091 13.311 4.566 10.611
C.V.”* 0.458 0.372 0.251 0.273 0.368
H, 12 17 34 29 16

2 2

0:1 035 039 033 037
GIBBS*“* 26.171 39.628 190.809 11.235 15.297
ACTUAL 40.242 64.957 184.955 12.052 10.423
GIBBS S.D.* 4.754 10.809 43.059 5.450 2.906
C.v.*“ 0.182 0.273 0.226 0.485 0.190
nj 59 27 42 10 58
* Standard deviation of the marginal posterior distribution
H Pooled within—group variance est1mate
* I I
** Coefficient of variation

****

Posterior Mean

 

155
Table 6.14-B

Generated data: Model 3, 6 = 0.20 and k = 40

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

 

 

7° Y1 Y2 [31 [32 p3
GIBBS**** 7.977 3.932 2.780 1.002 3.490 0.746
HLM 8.112 3.952 2.831 1.000 3.476 0.750
ACTUAL * 8.000 3.500 2.750 1.000 3.500 0.750
GIBBS S.D. 0.285 0.124 0.136 0.009 0.009 0.007
HLM S.E.H 0.333 0.152 0.166 0.012 0.011 0.009
Random error CG

U1 U5 U9 U13 U17
GIBBS**** -1.024 0.357 -1.175 1.889 0.171
ACTUAL -2.831 0.945 -1.059 2.246 0.568
GIBBS S.D.* 1.019 0.950 0.993 0.826 0.973

£51 £55 £59 £53 £57
GIBBs**** 1.289 0.533 0.453 0.529 -0.333
ACTUAL 1.755 2.368 -0.490 1.533 -8.370
GIBBS S.D.* 0.687 0.934 1.079 1.05 .
* Standard deviation of the marginal posterior'distribution
M Standard error of the HLM est1mate of the parame er

Posterior Mean

 

 

 

 

Generated data:

156
Table 6.15-A

Model 3, 6 = 0.20 and k = 15

Variance components estimates

Hyper-parameters

 

 

 

 

 

 

0 03 12
GIBBs**** 0 . 435 GIBBS**** 8 . 719 2 . 064
MOMENT EST 0.319 HLM 16.590“ 2.061
ACTUAL * 0.200 ACTUAL 10.000 2.250
GIBBS S.D. 0.226 GIBBS S.D.* 2.309 1.209
CHI-SQUARE 106.600 D.F. = 14 P-value = 0.000

Within-groups residual variances

0i cg a: 05 03
GIBBS“”* 3.796 10.678 3.858 37.938 12.871
ACTUAL 3.559 6.980 2.339 31.373 16.116
GIBBS S.D.* 1.834 2.302 1.603 9.166 3.378
C.V.“” 0.483 0.216 0.416 0.242 0.262
nj 12 42 17 37 34

CI]. 0:3 015
GIBBS“”* 19.420 9.881 46.783
ACTUAL 19.165 11.071 49.120
GIBBS S.D.* 3.763 2.869 13.431
C.V.“* 0.194 0.290 0.287
nj 55 29 29
* Standard deviation of the marginal posterior distribution
** Pooled within—grOUp var1ance est1mate

‘“* Coefficient of variation

***‘k

Posterior Mean

 

 

Generated data:

157

Table 6 . 15-B

Model 3, 8 = 0.20 and k =

15

Regression coefficients and random error CG estimates

Regression coefficients

 

 

 

 

 

Y” Y1 72 [31 92 p3
GIBBS**** 7.774 3.512 2.947 1.004 3.488 0.743
HLM 7.566 3.525 2.975 1.015 3.484 0.744
ACTUAL * 8.000 3.500 2.750 1.000 3.500 0.750
GIBBS 5'2: 0.490 0.115 0.218 0.013 0.011 0.009
HIM s.E. 0.503 0.130 0.246 0.015 0.012 0.011
Random error CG

U1 U3 U5 U7 U9
GIBBS**** 0.875 0.486 0.557 -0.320 1.079
ACTUAL 0.149 0.915 0.478 -O.556 2.500
GIBBS S.D.* 1.095 1.140 0.714 0.911 0.755

U11 U13 U15
GIBBS**** -1.920 —0.911 1.137
ACTUAL -2.992 -1.568 1.961
GIBBS S.D.* 0.813 1.006 1.040

 

Standard deviation Of the marginal posterior distribution
H Standard error of the HLM estimate of the parameter

Posterior Mean

 

Generated data:

158

Model 3,

6 = 0.02 and k = 100

Table 6.16-A

Variance components estimates

Hyper-parameters

 

 

 

 

 

 

9 12
GIBBS**** 0 . 029 GIBBs**** 30 . 205 3 . 361
MOMENT EST 0.054 HIM 32.231“ 3.290
ACTUAL * 0.020 ACTUAL * 30.000 3.063
GIBBS S.D. 0.010 GIBBS S.D. 1.184 0.683
CHI-SQUARE 206.971 D.F. = 99 P-value = 0.000

Within—groups residual variances

Ci 0:1 031 031 021
GIBBS**** 27.289 28.398 22.138 25.187 30.355
ACTUAL 27.242 26.377 26.066 25.761 29.086
GIBBS S.D.* 5.976 4.280 3.510 4.057 6.672
C.V.*** 0.219 0.151 0.159 0.161 0.220
nj 12 55 59 43 9

2

051 061 031 081 091
GIBBS**** 31.380 28.421 32.066 37.254 33.112
ACTUAL 38.593 28.148 25.195 28.179 32.216
GIBBS S.D.* 5.334 6.337 4.800 6.000 6.805
C.V.*** 0.170 0.223 0.150 0.161 0.206
11]. 36 9 49 42 14
* Standard deviation of the marginal posterior distribution
‘” Pooled within group-var1ance est1mate
‘"* Coefficient of variation

****

Posterior Mean

 

 

 

 

159
Table 6 . 16-B

Generated data: Model 3, 0 = 0.02 and k = 100

Regression coefficients and random error CG estimates

Regression coefficients

 

Ya Y1 Y2 B1 92 p3
GIBBS*“* 7.863 3.724 2.855 1.002 3.490 0.744
HLM 7.864 3.722 2.860 1.002 3.492 0.742
ACTUAL 8.000 3.500 2.750 1.000 3.500 0.750
GIBBS S.D.* 0.225 0.116 0.200 0.005 0.005 0.009
HLM s.E.“ 0.219 0.113 0.193 0.005 0.005 0.009

 

Random error CG

 

 

 

U1 U1 1 U2 1 U3 1 U4 1
GIBBS““ —1.423 1.450 —1.360 0.834 1.150
ACTUAL -1.258 0.656 —2.155 0.390 2.523
GIBBS S.D.* 1.169 0.713 0.665 0.805 1.252

U5 1 U6 1 U7 1 U8 1 U9 1
GIBBs*”* -1.080 0.618 2.432 2.268 0.134
ACTUAL —1.783 0.499 0.851 2.656 -0.328
GIBBS S.D.* 0.937 1.353 0.805 0.875 1.238

 

 

Standard deviation of the marginal posterior distribution
H Standard error of the HLM estimate of the parameter

Posterior Mean

 

 

 

160

Table 6.17-A

 

 

 

 

 

 

Generated data: Model 3, 9 = 0,02 and k = 40
Variance components estimates
Hyper—parameters

e of t2
GIBBS“"' 0.031 GIBBS““ 15.104 5.922
MOMENT EST 0.020 HLM 16.084“ 5.389
ACTUAL * 0.020 ACTUAL 15.000 6.250
GIBBS S.D. 0.017 GIBBS S.D.* 0.974 1.599
CHI-SQUARE 66.218 D.F. = 39 P-value = 0.003
Within-groups residual variances
oi as 09 013 017
GIBBS*“* 13.732 15.412 14.665 14.005 16.535
ACTUAL 13.445 14.186 17.284 13.484 13.841
GIBBS S.D.* 3.229 3.482 2.548 2.665 3.459
C.V.*“ 0.235 0.226 0.174 0.190 0.209
nj 12 17 34 29 16
2
0:1 035 039 0:3 037
GIBBS*“* 13.157 12.930 18.131 13.852 18.624
ACTUAL 18.566 14.631 19.086 19.518 14.904
GIBBS S.D.* 2.113 2.569 3.093 3.288 2.892
C.V.“* 0.161 0.199 0.171 0.237 0.155
nj 59 27 42 10 58
* Standard deviation of the'marginaI;posterior'distribution
H Pooled within—group variance est1mate

Coefficient of variation
Posterior Mean

 

 

 

161
Table 6 . 17-B

Generated data: Model 3, 6 = 0.02 and k = 40

Regression coefficients and random error LG estimates

Regression coefficients

 

7.. Y1 Yz pl pz [33
GIBBS*“* 8.297 4.603 2.837 0.996 3.484 0.749
HLM 8.319 4.579 2.831 0.999 3.483 0.748
ACTUAL 8.000 3.500 2.750 1.000 3.500 0.750
GIBBS S.D.* 0.422 0.326 0.218 0.014 0.010 0.011
HLM s.E.“ 0.397 0.289 0.218 0.014 0.010 0.011

 

Random error LG

 

 

 

 

 

 

U1 U5 U9 U13 U17
GIBBS*“* -2.482 1.192 —1.502 3.093 0.660
ACTUAL -4.718 1.575 -l.765 3.744 0.946
GIBBS S.D.* 1.399 1.032 1.020 0.782 1.098

U2 1 U2 5 U2 9 U3 3 U3 7
GIBBS*”* 1.594 1.856 —1.293 0.369 —0.901
ACTUAL 2.925 3.947 -0.817 2.566 —0.634
GIBBS S.D.* 0.741 0.845 0.826 1.594 1.028

 

 

f the marginal posterior distribution

' 'On 0
Standard dev1at1 he parameter

" Standard error Of the HLM estimate of t
Posterior Mean

 

 

 

 

162
Table 6.18-A

Generated data: Model 3, 6 = 0.02 and k = 15

Variance components estimates

Hyper-parameters

 

 

9 03 t2
GIBBS“”' 0.067 GIBBS“" 32.417 3.363
MOMENT EST 0.060 HLM 35.615“ 3.104
ACTUAL * 0.020 ACTUAL 30.000 4.000
GIBBS S.D. 0.051 GIBBS S.D.* 4.231 2.253
CHI-SQUARE 29.843 D.F. = 14 P-value = 0.006

Within-groups residual variances

 

 

 

 

of a: 0% 03 a;
GIBBS*“* 25.497 55.196 30.492 32.398 28.137
ACTUAL 20.225 42.615 22.852 25.591 30.396
GIBBS S.D.* 8.123 11.692 8.370 6.365 5.662
C.V.“* 0.319 0.212 0.274 0.196 0.201
nj 12 42 17 37 34
Oil 013 015

GIBBS*“* 33.894 37.273 30.061
ACTUAL 33.129 44.566 27.008
GIBBS S.D.* 5.721 8.680 6.977
C.V.”* 0.169 0.233 0.232
nj 55 29 29
:* Standard deviation of the marginal posterior‘distribution

Pooled within-group var1ance est1mate

Coefficient of variation
Posterior Mean

****

 

163
Table 6.18-B

Generated data: Model 3, 6 = 0.02 and k = 15

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

Y“ Y1 Y2 [31 [32 p3
GIBBS*“* 5.487 4.505 2.910 2.040 3.451 1.712
5.504 4.507 2.933 2.035 3.452 1.705
ACTUAL . 6.000 4.500 2.500 2.000 3.500 1.750
GIBBS S'R' 0.637 0.303 0.370 0.054 0.034 0.065
HLM s.E. 0.622 0.281 0.346 0.052 0.033 0.067
Random error [E

cg Lg cg I; I;
GIBBS*“* 1.023 0.680 0.705 -0.583 1.266
ACTUAL 0.199 1.220 0.638 —0.742 3.333
GIBBS S.D.* 1.506 1.566 1.221 1.124 1.045

£51 Uh Uh
GIBBS*“* —2.152 —1.205 1.834
ACTUAL -3.990 -2.091 2.615
GIBBS S.D.* 1.122 1.290 1.167

 

of the marginal posterior distribution

Standard deviation
the HLM estimate of the parameter

Standard error of
Posterior Mean

 

164
Table 6.19-A

HSB data: Model 1: MATHACH = BASE + SECTOR + MEANSES
MINORITY + GENDER + SES

Variance components estimates

Hyper—parameters

 

 

0 O. t
GIBBS“* 0.007 GIBBS“* 35.502 1.824
MOMENT EST* 0.010 HLM 35.908“ 1.790
GIBBS S.D. 0.004 GIBBS S.D.* 0.751 0.319
CHI-SQUARE 233.039 D.F. = 159 P—Value = 0.000

Within-groups residual variances

2 2 2 2 2
01 017 033 049 065

 

35.339 33.639 36.259

 

 

GIBBS“* 39.826 35.557

GIBBS S.D.* 4.557 3.902 3.749 3.637 3.609
nj 47 29 47 53 58
2 2 2 2 2

081 097 0113 0129 0145

GIBBS“* 38.854 36.921 35.065 34.116 35.909

GIBBS S.D.* 3.986 3.927 3.857 3.748 3.998
49 36

nj 66 35 44
_________________________________________________________________
n of the marginal posterior distribution

:* Standard.deviatio _
Pooled within group—variance est1mate

Posterior mean

 

165

Table 6.19-B

HSB data: Model 1: MATHACH = BASE + SECTOR + MEANSES
MINORITY + GENDER + SES

Regression coefficients and random error lg estimates

Regression coefficients

BASE SECTOR MEANSES MINORITY GENDER SES

 

 

 

 

 

Y. Y1 72 I31 [32 Ba
GIBBS*“ 13.310 1.727 2.152 —2 826 —1.245 1.892
HLM * 13.316 1.718 2.106 —2.838 -1.244 1.917
GIBBS S'R' 0.212 0.287 0.357 0.211 0.161 0.115
HLM S.E. 0.202 0.278 0.356 0.201 0.159 0.108

Random error CG

U.1 U17 U33 U49 U6 5
GIBBS“* -0.604 1.747 0.161 —0.425 —1.038
GIBBS S.D.* 0.732 0.857 0.769 0.687 0.755

U81 U97 U113 U129 U145
GIBBS“* —0.605 1.148 -1.742 0.230 1.229
GIBBS S.D.* 0.700 0.826 0.776 0.703 0.833

_//

Standard deviatio

H Standard error of the HLM estimate

Posterior mean

n of the marginal posterior distribution
of the parameter

 

166
Table 6.20-A

HSB data: Model 2: MATHACH = BASE + SECTOR + GENDER + SES

Variance components estimates

Hyper—parameters

2

 

 

0 O. t
GIBBS”* 0.008 GIBBS“* 36.402“ 3.550
MOMENT EST 0.008 HLM * 36.818 3.396
GIBBS S.D.* 0.004 GIBBS S.D. 0.699 0.515
CHI-SQUARE 219.569 D.F. = 159 P—Value = 0.000

Within—groups residual variances

 

 

 

2
01 017 533 0:9 065
GIBBS“* 41.255 36.181 37.185 34.213 36.626
GIBBS S.D.* 4.990 4.247 3.946 3.496 3.931
nj 47 29 47 53 58
2 2 2 02 02
081 097 0113 129 145

 

 

GIBBS“* 39.494 37.160 35.774 34.270 37.458
GIBBS S.D.* 4.030 4.298 3.777 3.416 4.289
n 66 35 44 49 36

j
_________________________________________________________________
Standard deviation of the marginal posterior distribut1on
H Pooled within group-variance est1mate
Posterior mean

 

 

167
Table 6.20-B

HSB data: Model 2: MATHACH = BASE + SECTOR + GENDER + SES

Regression coefficients and random error CG estimates

Regression coefficients

 

 

BASE SECTOR GENDER SES

7. 71 B. (32
GIBBS*” 12.323 2.116 —1.207 2.309
HLM 12.347 2.100 -1.203 2.341
GIBBS S.D.* 0.229 0.336 0.161 0.112
HLM s.E.“ 0.237 0.330 0.164 0.105

 

Random error CG

 

 

 

 

Q. £57 £53 (Q9 [Q5
GIBBS*“ -0.678 3.644 0.643 0.842 —3.512
GIBBS S.D.* 0.871 0.981 0.816 0.778 0.742

U81 U97 U113 U129 U145
GIBBs“* 0.468 -0.151 -4.173 1.087 0.3::
GIBBS S.D.* 0.717 0.917 0.901 0.802 O.
* Standard deviation of the marginal posterior distribution
“ Standard error of the HLM estimate of the parameter

Posterior mean

 

168
Table 6.21-A

HSB data: Model 3: MATHACH = BASE + MEANSES + GENDER + SES

Variance components estimates

Hyper—parameters

 

 

6 03 12
GIBBS“* 0.008 GIBBS“* 36.303“ 2.580
MOMENT EST 0.011 HLM * 36.798 2.496
GIBBS S.D.* 0.004 GIBBS S.D. 0.712 0.395
CHI-SQUARE 238.342 D.F. = 159 P-Value = 0.000

Within—groups residual variances

 

 

 

2
Ci 017 033 0:9 065
GIBBS”* 41.358 35.680 36.938 33.793 36.397
GIBBS S.D.* 5.063 4.028 3.858 3.765 3.979
nj 47 29 47 53 58
2 2
081 0:7 0113 0129 0145
GIBBS“* 40.024 37.082 35.367 34.199 32.33:
GIBBS S.D.* 4.498 4.074 3.812 3.530 . 6
n- 66 35 44 49 3

J

_________________________________________________________.__._

Standard deviation of the marginal posterior'distribution
Pooled within group-variance est1mate
Posterior mean

 

169
Table 6.21-B
HSB data: Model 3: MATHACH = BASE + MEANSES + GENDER + SES

Regression coefficients and random error CG estimates

Regression coefficients

 

BASE MEANSES GENDER SES

Yo Y1 Bl p2

GIBBS"* 13.278 3.668 —1.185 2.123
HLM * 13.280 3.610 —1.183 2.153
GIBBS 3.3. 0.170 0.358 0.167 0.111
HLM S.E. 0.168 0.368 0.163 0.108

 

Random error (G

 

 

 

 

U'1 U17 U33 [J49 [JG 5
GIBBS*“ -0.294 1.026 -0.021 —0.114 -1.339
GIBBS S.D.* 0.869 0.888 0.761 0.717 0.715

U81 U97 U113 U129 U145
GIBBS“* 0.236 0.459 -2.260 1.054 1.001
GIBBS S.D.* 0.737 0.896 0.783 0.759 0.898
* Standard deviation of the marginal posterior'distribution
M Standard error of the HLM estimate of the parameter
***

Posterior mean

 

170
Table 6.22-A

HSB data: Model 4: MATHACH = BASE + SECTOR + SES

Variance components estimates

Hyper-parameters

 

 

6 of 12
GIBBS”* 0.008 GIBBS”* 36.564 3.799
MOMENT EST* 0.008 HLM 37.037" 3.685
GIBBS S.D. 0.004 GIBBS S.D.* 0.724 0.564
CHI-SQUARE 220.466 D.F. = 159 P—value = 0.000

Within—groups residual variances

 

 

Ci 0:7 033 0:9 025
GIBBS“* 41.847 36.353 37.309 34.250 36.847
GIBBS S.D.* 5.183 4.103 4.086 3.503 3.844
nj 47 29 47 53 58

2 2 2 2 2

081 097 0113 0129 0145

 

37.522 35.941 34.406 37.989
3.929 3.513 4.355
36

GIBBS“* 39.946

GIBBS S.D.* 4.159 4.573
n- 66 35 44 49

Standard deviation of the marginal posterior'distribution
H Pooled within group-variance est1mate
Posterior mean

 

171
Table 6.22-B

HSB data: Model 4: MATHACH = BASE + SECTOR + SES

Regression coefficients and random error LG estimates

Regression coefficients

 

BASE SECTOR SES

Yo Y1 pl

GIBBS“* 11.703 2.103 2.350
HLM * 11.719 2.101 2.375
GIBBS S'R‘ 0.222 0.332 0.108
HLM s.E. 0.228 0.341 0.105

 

Random error (G

 

 

 

U1 U17 [J33 U49 U65
GIBBS”* -0.801 3.555 1.095 0.877 -2.967
GIBBS S.D.* 0.884 0.970 0.818 0.735 0.776

U81 U97 U113 U129 U145
GIBBS“* 0.561 -0.311 —4.171 1.593 0.884
GIBBS S.D.* 0.783 0.919 0.873 0.812 0.934

 

* Standard deviation of the marginal posterior distribution
H ' f the arameter
Standard error Of the HLM est1mate O p

Posterior mean

 

172
Table 6.23-A

HSB data: Model 5: MATHACH = BASE + SECTOR + GENDER

Variance components estimates

Hyper-parameters

 

 

0 of 12

GIBBS‘“ 0 016 GIBBS“* 37 998
. . 6.424
MOMENT EST* 0.019 HLM 38.854“ 6.242
GIBBS S.D. 0.005 GIBBS S.D.* 0.851 0.861
CHI-SQUARE 291.573 D.F. = 159 P-value = 0.000

Within-groups residual variances

 

 

2 2 2 2
01 017 033 049 0:5
GIBBS*“ * 46.613 36.517 38.510 32.513 35.701
GIBBS S.D. 6.900 5.451 5.469 4.703 4.778
nj 47 29 47 53 58
2 2
081 037 0113 0129 014s

 

GIBBS*” 45.759 40.558 35.976 34.322 41.787
GIBBS S.D.* 5.768 6.166 4.827 4.446 5.936
n 66 35 44 49 36

j

 

Standard deviation of the marginal posterior distribution
M Pooled within group-variance estimate
Posterior mean

 

HSB data:

Model 5:

173

Table 6.23-B

MATHACH = BASE + SECTOR + GENDER

Regression coefficients and random error LG estimates

Regression coefficients

 

 

 

 

 

 

BASE SECTOR GENDER

Yo Y1 Bl
GIBBS*** 12.080 2.807 -1.337
HLM * 12.113 2.791 -1.371
GIBBS S.E. 0.312 0.448 0.171
HLM S.E. 0.298 0.426 0.171

Random error (G

Ul Ul7 [J33 U49 U65
GIBBS*** -l.388 5.665 1.459 1.940 —4.969
GIBBS S.D.* 0.952 1.068 0.900 0.808 0.835

U81 U97 U113 U129 U145
GIBBS*** 1.023 —0.900 -5.504 1.439 0.6:):
GIBBS S.D.* 0.879 1.023 0.948 0.828 1.0
* Standard deviation of the marginal posterior'distributlon
" Standard error of the HLM est1mate of the parameter

** a
* PosterIOr mean

 

174

Table 6.24-A

 

 

 

HSB da : - —
ta Model 6. MATHACH _ BASE + MEANSES + SES
Variance components estimates
Hyper-parameters

0 03 :2
GIBBS*“ 0.009 GIBBs“* 36.533 2.753
MOMENT EST* 0.012 HLM 37.019“ 2.692
GIBBS S.D. 0.005 GIBBS S.D.* 0.761 0.400
CHI-SQUARE 240.750 D.F. = 159 P—value = 0.000

Within—groups residual variances

2 2 2 2 2
01 017 033 O49 O65

 

GIBBS“*
GIBBS S.D.*

 

42.197 35.813 37.292 34.027 36.388
4.923 4.028 4.079 3.891 3.901
47 29 47 53 58

2 2
0:1 0:7 0113 0129 0145

 

GIBBS“*

40.194 37.526 35.623 34.302 38.145

 

GIBBS S.D.* 4.301 4.543 4.098 3.796 4.424

Hj 66 35 44 49 36

:* Standard deviation of the marginal posterior'distribution
Pooled within group—variance estimate

Posterior mean

 

175
Table 6 . 24-B

HSB data: Model 6: MATHACH = BASE + MEANSES + SES

Regression coefficients and random error LG estimates

Regression coefficients

 

BASE MEANSES SES

Yo Y1 pl

GIBBS“* 12.651 3.740 2.160
HLM * 12.661 3.675 2.191
GIBBS S.E. 0.149 0.403 0.114
HLM S.E. 0.149 0.378 0.109

 

Random error LG

 

 

 

U1 U17 U33 U49 U65
GIBBS*“ -0.315 0.957 0.362 —0.159 -0.817
GIBBS S.D.* 0.895 0.872 0.833 0.734 0.726

U81 U97 U113 U129 U145
GIBBS“* 0.200 0.409 —2.202 1.508 1.028
GIBBS S.D.* 0.740 0.870 0.786 0.745 0.885

 

* Standard deviation of the marginal posterior distribution
M Standard error of the HLM estimate of the parameter

Posterior mean

 

176
Table 6.25-A

HSB data: Model 7: MATHACH = BASE + SECTOR

Variance components estimates

Hyper-parameters

 

 

0 of 2c 2
*** ***
GIBBS 0.017 GIBBS 38.243 6.871
MOMENT EST* 0.020 HLM 39.151" 6.680
GIBBS S.D. 0.005 GIBBS S.D.* 0.925 0.930
CHI-SQUARE 299.157 D.F. = 159 P-value = 0.000

Within-groups residual variances
2 2 2 2 2
01 G17 O33 O49 065

38.249 32.763 36.508

GIBBS“* 47.627 36.521
GIBBS S.D.* 6.670 5.717 5.185 4.435 5.064
n 47 29 47 53 58

i_/,/

2 2 2 2 2
O81 097 0113 0129 0145

__//

34.949 43.323

 

GIBBS“* 46.460 40.872 36.236

GIBBS S.D.* 5.865 5.905 5.293 4.869 7.007

nj 66 35 44 49 36
distribution

ion of the marginal posterior

Standard.deviat .
p-variance est1mate

H Pooled within grou
Posterior mean

_—“-“I

177
Table 6.25-B

HSB data: Model 7: MATHACH = BASE + SECTOR

Regression coefficients and random error LG estimates

Regression coefficients

 

BASE SECTOR

Yo YI

GIBBS“* 11.398 2.790
HLM * 11.393 2.805
GIBBS S.E. 0.294 0.442
HLM S.E. 0.293 0.439

 

Random error LG

 

 

 

U1 U17 U3 3 U49 U65
GIBBS"* —1.455 5.739 2.113 2.037 —4.376
GIBBS S.D.* 1.006 1.105 0.934 0.840 0.830

U81 U97 U113 U129 U145
GIBBS”* 1.115 —1.163 -5.459 2.055 0.568
GIBBS S.D.* 0.877 1.043 0.907 0.847 1.076

 

** Standard deviation of the marginal posterior‘distribution
Standard error of the HLM estimate of the parameter

Posterior mean

 

178
Table 6.26-A

HSB data: Model 8: MATHACH = BASE + SES

Variance components estimates

Hyper-parameters

 

 

9 o. r
GIBBS“* 0.007 GIBBs"* 36.621 4.944
MOMENT EST* 0.007 HLM 37.034“ 4.768
GIBBS S.D. 0.003 GIBBS S.D.* 0.739 0.713
CHI-SQUARE 214.629 D.F. = 159 P-value = 0.000
Within-groups residual variances
01 “i7 033 0:9 0:5

__________________________________________________________________

GIBBS*“ 41.826 36.284 37.445 34.474 36.657
GIBBS S.D.* 4.728 4.176 3.898 3.621 3.474
n. 47 29 47 53 58

J

__________________________________________________________

2 2 2 2 2
O81 097 0113 0129 0145

__//

GIBBS"* 39.659 37.365 36.274 34.900 37.725
GIBBS S.D.* 3.791 4.028 3.803 3.673 3.847
n- 66 35 44 49 36

.7

 

of the marginal posterior distribution

* Standard deviation '
est1mate

** I l l
Pooled w1th1n group—var1ance
Posterior mean

—““i

*fi

179
Table 6 . 2 6-B

HSB data: Model 8: MATHACH = BASE + SES

Regression coefficients and random error [9 estimates

Regression coefficients

 

BASE SES

Yo pl

GIBBS”* 12.640 2.358
HLM * 12.657 2.390
GIBBS S.E. 0.189 0.112
HLM S.E. 0.188 0.106

 

Random error CG

 

 

 

U1 U17 U33 U49 U55
GIBBS“* —1.600 3.137 2.160 1.972 -2.071
GIBBS S.D.* 0.904 1.007 0.810 0.807 0.755

U81 U97 U113 U129 U145
_________,______________________________________________
GIBBS“* 1.586 -1.110 —3.379 2.693 0.163
GIBBS S.D.* 0.757 1.010 0.902 0.809 0.931

 

on of the marginal posterior distribution

d deviati
Standar e of the parameter

“ Standard errorthe HLM estimat
Posterior mean

___TF1____442222444 44———————————————————------‘---II

180
Table 6.27-A

HSB data: Model 9: MATHACH = BASE

Variance components estimates

Hyper—parameters

 

 

e 002 t2
GIBBS“* 0.016 GIBBS“* 38.189 8.943
MOMENT EST* 0.020 HLM * 39.148“ 8.615
GIBBS S.D. 0.005 GIBBS S.D. 0.868 1.146
CHI-SQUARE 299.157 D.F. = 159 P-Value = 0.000

Within-groups residual variances

 

 

 

 

01 017 033 0:9 065
GIBBS“* 47.229 36.496 38.592 32.718 36.002
GIBBS S.D.* 7.184 5.733 5.258 4.268 4.959
nj 47 29 47 53 58

2

0:1 037 0113 0129 0145
GIBBS“* 46.327 41.056 36.334 34.998 42.652
GIBBS S.D.* 5.972 6.164 5.137 5.077 6.185
n. 66 35 44 49 36

.7

 

Standard deviation of the marginal posterior distribution
Pooled within group—variance estimate
Posterior mean

 

 

 

181
Table 6.27-B

HSB data: Model 9: MATHACH = BASE

Regression coefficient and random error CG estimates

Regression coefficient

 

 

 

 

 

BASE
Y.
GIBBS“* 12.638
HLM 12.637
GIBBS S.E.* 0.249
HLM S.E. 0.244
Random error CG
U1 U17 U33 U49 U65
GIBBS“* —2.625 4.817 3.499 3.509 -3.054
GIBBS S.D.* 0.967 1.062 0.842 0.769 0.810
U81 U97 U113 U129 U145
i GIBBS"* 2.568 —2.244 —4.114 3.496 —0.507
GIBBS S.D.* 0.797 1.115 0.943 0.875 1.033

 

Standard deviation of the marginal posterior distribution
Standard error of the HLM estimate of the parameter
Posterior mean

 

 

CHAPTER 7
Discussion

Discussion of the results is presented in this chapter.
It is focused on comparing the application of Bayesian
approach via Gibbs sampling in multi—level analysis to the
empirical Bayes approach via HLM analysis when there is
heterogeneity of variance in the first level. Suggestions and
recommendations for future research are also presented in this

chapter.
Bayes and Empirical Bayes Estimation

The main purpose of this empirical study was to apply the
fully Bayesian approach to the analysis of multi—level data
for the cases where the homogeneity of variance assumption can
not be granted and when interest of the research is focused on
making inferences on some or all of the groups variances.
Available empirical Bayes methods for analyzing multi—level
data often assume homogeneity of variance and concentrate on
obtaining empirical Bayes estimates for the regression
coefficients. The question here is: "How do Bayes estimates

of variance components and regression coefficients behave when

computed via Gibbs sampling?"

 

 

 

183

In general, the results presented in the previous chapter
suggest that there are no substantial differences between the
two approaches in the estimation and inferences about the
regression coefficients. This finding suggests that. HLM
estimates of the regression coefficients, particularly for
group effects are robust to the violation of homogeneity of
variance. However, when it comes to the estimation of the
variance components, HLM estimates of of are found to be
positively biased especially when there exists clear evidence
of heterogeneity of variance. A moderate heterogeneity of
variance with large number of groups seems to have little
effect on the HLM estimate of 03. This was demonstrated in
all models used on the High School and Beyond data set.

Since empirical Bayes estimates of the regression
coefficients are conditioned on estimating the variance
components, it was anticipated that the regression coefficient
estimates will be affected by the uncertainty in estimating
those variance components, especially in the cases when the
HLM estimate of of is quite different from its true value.
However, when they are compared to the Bayes estimates and the
actual values of the regression coefficients they were found
to be within the same range of values with about the same
random error of estimate. What appears to happen is that
larger HLM estimates of of help compensate for the

heterogeneity of variance, 0; when it comes to the estimation

 

 

184
of regression coefficients, and their standard errors of
estimate.

One advantage of applying the Bayesian approach to the
analysis of multi-level data is the ability to obtain the full
posterior distribution for each parameter involved in the
model under study. A wide rage of statistics describing a
particular parameter and its distribution such as mean, mode,
percentile points and. variance can be obtained from ‘the
posterior distribution of the parameter.

Furthermore the ability to Obtain the marginal posterior
distribution for each of the residual variances and the
parameters 0 and of of their prior distribution becomes quite
important when research interest is focused on the residual
variances themselves and their heterogeneity. To assess this

heterogeneity, a coefficient of variation (C.V.) for 'the

 

residual variance can be easily found as 1320 (lindley,
1965). Inferences about a particular group residual variance

also can be made using its marginal posterior distribution.
The application of the Bayesian approach to the analysis
of multi-level data involved the use of Gibbs sampling. One
disadvantage of this procedure is the time and cost involved
in the process. This problem is reflected in this study by
having only one sample for each of the 18 conditions used in
the study. Limiting the analysis to one sample prevented us

from carefully studying the characteristics of the parameter

 

 

 

185
estimates using their mean squared errors and statistical
tests against their true values.

The idea of getting the marginal distribution of a
variable from its conditional distribution on another variable
is very appealing. Gibbs sampling and Data augmentation are
two approaches which adopt that idea in many classical and
Bayesian. analyses. However, developments :h1 their
applications are limited to certain research problems.
General algorithms and software that utilize these approaches
are not widely available to practitioners for use in their

practices.

Suggestions and Recommendations for Future Research

To investigate the differences between the fully Bayesian
approach and the empirical Bayes approach more closely, it is
recommended to use the empirical Bayes model that allows for
heterogeneity of variance, Mfi}, in a simulation study where
many estimates of the same parameters in a given model are
derived from several samples. Based on those many estimates
of the parameters, a probability statement and mean squared
errors can be used to compare the two approaches in parameter
estimation.

Simulating several samples for a given model might seem
straightforward in studying the hyper-parameters 6, 03 , :2 and

K, where their values can be pre—specified for the data

 

 

 

 

186
generation. However, it becomes more complex in studying the
parameters {02} and {Uj} where j=1,-~,k. That is because the
true values for both sets Of parameters, {0?} and {Uj}
represent two random samples (see "Data creation" section in
chapter 5) from their prior distributions with pre-specified
values of the hyper-parameters. This means that we have no
direct control on the actual true values of {0%} and {Uj}. The
only way that one can alter the true values of {03'} and {Uj} is
by either respecifying the values of the hyper-parameters of
their prior distributions or respecifying totally different
parametric forms for their priors. More specifically, the
true values for {0;} represent a'random sample drawn from an
inverse gamma distribution with a=—2% and (i=4? (see equation
5.59) . To alter these true values we need either to respecify
the values of 03 and 6 in the prior distribution of {03} in
5.59, or chose different parametric form for their prior.
Similarly, the true values for {Uj} represent a random
sample drawn from a normal distribution with mean equal to
zero and variance equal to 1:2. To alter these true values we
either respecify the value of t2 in the prior distribution of{Uj}
or chose different parametric form of the prior such as t—
distribution (see Seltzer, 1993).
Realizations of the outcome variable Y and the predictors
in Z are based on the generated values of {02} and {Uj} in the

above steps and the pre-specified values Of A. Thus,

randomness in Y and Z, took place in two steps: in the

 

 

 

 

 

 

187
generation of the true values of {0? and {U}, and in the
generation of Y and Z which are based on the generated values
of kfi}, Hg} in the first step and 1. Therefore, the posterior
estimates of the hyper—parameters, 6, of, t2 and 1 based on
a given sample of Y and Z can only be as "good" (close to the
true values of 0, of, t2 and A) as those estimates of 0, 03,12
and 1 produced by the parameters,{o? and MGL Another way of
explaining this is that if we define 6*, 03*, A‘ and 12* as
estimates of 0, of, A and 12, which are based only on the
true values of {0? and {03; also we define 0’, 0?, 1’ andty
as estimates for the same parameters which are based on the
parameters {03,-} and {Uj} and the data Y and Z, then the
estimates 0’, 0?, 1’ and I” cannot be as good as 0*, 03*,1*
and 12* in estimating the true values of 0, of, A and 12.
That is because Of the extra randomness in Y and Z added
through the generation of (a? and MEL

When several (say G) samples of {0? and {0; are being
generated for the same model, the chance of misrepresenting
the true values of the hyper—parameters by one bad sample will
be greatly reduced by the presence of the other samples offifi}
and.{Ufi. A mean squared error of estimates and a distribution
of the estimates derived from these samples can be used to
make inferences about the hyper-parameters.

If only one set of realization of Y and Z is generated
for each of the G samples of the true values Of{o$ and Ugh

then the problem of having only one set of estimates of (a?

 

 

 

 

188
and.tqfi will be the same as in the hyper-parameters. That is
because for each one of these G samples of {03'} and {Uj} we only
have one set of realizations for the variables Y and Z.
Therefore, we can only Obtain one set of estimates of {0;} and
{QB for each of the G samples of Y and Z.
To overcome this problem we simply generate several (say
.1) samples of the variables Y'and Z, for each of the G sample
of the true values of {0;} and {Uj}. This allows us to obtain
a distribution of J estimates of {03:} and {Uj}, which can be
used to make inferences about the two sets of parameterskfi}
and.{qfi. Therefore, the total number of samples that need to
be generated to estimate all the parameters in the model is
equal to GU'samples.
One of the debatable issues in Bayesian analysis is the
choice of the prior distributions (Deely and Lindley, 1981).

In this study conjugate priors for the second stage

(exchangeable) parameters were normal (0 , 'cz) for the {Uj},
and inverse gamma 1"1(-é%,3‘g) for {03.}. Priors for the third
00

stage parameters (hyper-parameters) were chosen to be
proportional to constants to reflect noninformative priors.
It would be extremely useful to know how inferences about all
parameters of the model change when different priors are being
used in both stages.

In this study, densities of the produced marginal
posterior distributions were approximated using the kernel

method. These estimates were used graphically in determining

 

 

 

 

189

the stopping point for the iteration process of Gibbs
sampling. Another method of estimating a density of a
marginal posterior distribution for a particular parameter is
by using equation 4.5 on the mixture of the conditional
densities for a particular parameter. Zeger and Karim (1991)
recommended the use of this method for getting better
estimates for the tail of the distribution. Gelfand, Smith
and Lee (1992) argued that using the mixture of the
conditional density produces a more accurate representation of
the density function than the kernel estimate. It would be
'Very useful to know if overlay graphs of densities produced.by
(4.5) are more efficient.in determining the stopping point for
Gibbs sampling; This is important because of the large amount
of computer resources required by Gibbs sampling.

This study attempted to model conditions often found in
many education research projects where only intercepts from
regression models of many groups were allowed to be random in
multi-level analysis. Interactions between individual
characteristics and group characteristics were not considered
in this study. As a result of this layout, the variance-
covariance matrix of the regression coefficients T reduced to
a scaler 12 which represents the variance of the intercept.
It is possible that as more elements added to this variance-
covariance matrix (i.e., allowing some of the within-group
effects to be random), the Bayesian estimates of these effects

might then become different from their empirical Bayes

 

 

 

190
estimates. Thus, it would be extremely informative if the
analysis in this study were expanded to allow some of the
within-group regression coefficients to be random with the
possibility of them being related to the group level
‘variables. This will change the scalar:2 to a full variance—
covariance matrix, T. Investigating such models provides some
insight on what conditions the Gibbs sampling is superior and
what conditions other statistical procedures are adequate.
Based on the result Of this study, investigators who are
only interested in the regression coefficients when analyzing
data similar to the one in this study are advised to use
empirical Bayes procedures, with assumption.of homogeneity of
variance, and for the cases with large number of groups, k.
Note in table 6.9—B that when k=15 and 0:0.20, the HLM
standard errors Of the 7's are too small. However, when
k=100 and 8:0.20, (table 6.13—B), the HLM standard errors of
the y's are nearly perfect. For these cases (large k)
empirical Bayes estimates of the regression coefficients and
their associated tests are not that different from those
obtained by fully Bayesian procedures. Also, computer
programs and application software that utilize empirical Bayes
procedure are widely available. When there is a clear
evidence of heterogeneity of variance and.research interest is
focused on studying group variances, applying the fully

Bayesian procedure might become more informative than

empirical Bayes procedure.

 

 

 

 

 

References

 

Aitkin, M., and Longford, N. (1986). Statistical Modeling
Issues in School Effectiveness Studies. Journal of the

Royal Statistical SocietyI Series A, 149, (1), 1-43.

Bandura A. (1977). Social learning theory. New Jersey,
Englewood Cliffs: PRENTICE HALL.

Bartlett, M. 8., and Kendall, D. G. (1946). The Statistical
Analysis of variance-Heterogeneity and the Logarithmic

Transformation. SuppZJournal of Royal Statistical
Societv. 8, 128-138.

Bassiri, IL. (1988). large and Small Sample Properties of
Maximum Likelihood Estimates for the Hierarchical Linear
Model. Ddoctoral Dissertation, IMichigan State
University.

Berlin, K. (1984). Statistical Analysis Of Mixed Linear
Models with Inhomogeneous Variances. COMPSTAT thsica-
Verlag, Vienna for IASCl International Association for
Statistical Computing. 209—213.

Bretthorst, G. L. (1988) . Bayesian Spectrum Analvsis and
Parameter Estimation. Lecture Notes in Statistics (48).
New York, Heidelberg: Springer—Verlag.

Burstein, L. (1980). The Analysis of Multilevel Data in
Educational Research and Evaluation. In D.C. Berliner
(Ed.), Reviey;of Research.in.EducatiLuL 8. ‘Washinqton.
DC: American Educationa Research Association.

Bloom, B. S. (1984). The z-Sigma Problem: The Search for
Methods of Group Instruction as Effective as One-to-One.
Educational Researcher, 1 , (6), 4-16.

Box, G. E. P., and Tiao, G. C. (1964). A Bayesian Approach to
the Importance of Assumptions Applied to the Comparison
of Variances. Biometrika. 51, (1) and (2), 153-167.

Box, G. E. P., and Tiao, G. C. (1973). Bavesian inference in
statistical analysis. Reading, Massachusetts: Addison-
Wesley.

Braun, H. I., Jones, D. H., Rubin, D. B., and Thayer, D. T.
(1983) . Empirical Bayes Estimation in the General Linear
Model with Data of Deficient Rank. Psychometrika. 48,
(2), 171-181.

191

 

 

192

   

Bryk, A. 8., and Raudenbush, S. W. (1987). Application of
Hierarchical Linear Models to Assessing Change.

Psychological Bulletin 101, (1), 147-158.

Bryk, A. S., and Raudenbush, S. W. (1988). Heterogeneity of
Variance in Experimental Studies: A Challenge to
Conventional Interpretations. Psychological Bulletin,
104, (3), 396-404.

 

Casella, G., and George, E. I. (1992). Explaining the Gibbs
Sampler. The American Statigtician. 46, (3), 167-174.

Cook, R. D., and.‘Weisberg, S. (1983). Diagnostics for
Heteroscedasticity in Regression. Biometrika. 70, (1),

 

1-10.

Cooley, W. W., Bond, L., and Mao, B. (1981). Analyzing
Multilevel Data. In Berk, R. A. (Ed.), Educational
evaluation methodology. Baltimore: Johns Hopkins

 

University press.

Cox, D. R., and Solomon, P..J. (1986). Analysis of Vriability
with large Numbers of Small Samples. Biometrika. 73,
(3), 543-554.

Cronbach, In {I., and Webb, N. (1975). Between and Within
Class Effects in Reported Aptitude-by-Treatment
Interaction: Reanalysis of a Study by G. L. Anderson.
Journal of Educational Psychology, 6, 717-724.

Deely, J. J., and Lindley, D. V. (1981). Bayes Empirical
Bayes. Journal of the American Satistical AssociationL
16, (376), 833—841.

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).
IMaximum. Likelihood from. Incomplete Data ‘Via the IBM
Algorithm. Royal Statistical SocieQLJournal. Series. B.
35, 1—22.

Efron, B., and Morris, C. (1972). Empirical Bayes on Vector
Observations. Biometrika, 59.

Efron, B., and.Morris, C. (1973). Stein's Estimation.Rule and
Its Competitors-An Empirical Bayes Approach” Journal of
the American Statistical Association. 68, 117—130.

 

Efron, B., and Morris, C. (1975). Data analysis Using Stein's
estimator’ and Its Generalizations. .Journal of ‘the
American Statistical Association. 72, 311-319.

Efron, B., and Morris, C. (1977). Comments. Journal Of the
American Statistical AssociationL 72, 91—93.

¥

 

193

Fotiu, R. P. (1989). A Comparison of the EM and Data
Augmentation Algorithms on Simulated Small Sample
Hierarchical Data from Research on Education. Doctoral
Dissertation, Michigan State University.

Gelfand, A. E., Hills S. E., Racine-Poon,.A., and.Smith, A. F.
M. (1990). Illustration of Bayesian Inference in Normal
Data Models Using Gibbs Sampling. Journal of the
American Statistical Association, 8 , (412), 972-985.

Gelfand, A. E., and Smith, A. F. M. (1990). Sampling-Based
Approaches to Calculating Marginal Densities. Journal of
the American Statistical Association, 85, (410) , 398—409.

Gelfand, A. E., Smith, A, F. M., and Lee, T. (1992). Bayesian
Analysis of Constrained Parameter and Truncated Data
Problems Using Gibbs Sampling. Journal of the American
Statistical Association. 87, (418), 523-532.

Goldstein, H. (1986). Multilevel Mixed Linear Model Analysis
Using Iterative Generalized Least Squares. Biometrika-
ll, (1). 43—56.

Harville, D. A. (1977). Maximum Likelihood Approachs to
Variance Component Estimation and to Related Problems.
Journal of'the.American.Statistical.Association 72, 320-

 

 

 

338.

Hogg, R. V., and Tanis, E. A. (1977). Probability and
statistical inference. New York: Macmillan Publishing
Co.

Kasim, R. (1986, April). Examining Variances in Hierarchical
Models. Paper presented at the annual meeting of the
American Educational Research Association, San Fransisco
CA.

Knapp, T. R. (1977). The 'Unit of .Analysis Problem. in
Applications of Simple Correctional Research. Journal of
Educational Statistics. 2, (3), 171-186.

Laird, N. M., and Ware, J. H. (1982). Random—Effect Models
for Longitudinal Data. Biometrich 38, 963-974.

Leonard, T. (1975). A Bayesian Approach to the Linear model
with unequal variances. Technometrics.,l7, (1), 95-102.

Lewis, P. A. W., and Orav, E. J. (1989). Simulation
methodology for mtatisticians, operations analysts, and
engineers. Volume 1. Pacific Grove, California:
Wadsworth.& Brooks/Cole, Statistics/Probability Series.

 

 

 

 

 

194

Lindley, D. V. (1965) . Introduction to probability and
statistics from a Bayesian vieypoint, part 2. Cambridge:
Cambridge University Press.

Lindley, D. V. (1971). Foundations of statistical inference
(ed. by V. P. Godmabe and D. A. Sprott). The Estimation of
Many Parameters. Toronto: Holt, Rinehart, and Winston, 435-
455.

Lindley, D. V., and Smith, A. F. M. (1972). Bayesian
Estimates for the Linear Model. Journal of the Royal

Statistical Society, Series B. 34, 1-41.

Linhart, H. (1965) . Approximate Confidence Limits for the
Coefficient of Variation of Gamma Distributions.
Biometrics. 21, 733-738.

 

 

Longford, N. T. (1987) . A Fast Scoring Algorithm for Maximum
Likelihood Estimation in Unbalanced Mixed Models with
Nested Effects. Biometrika. 74, 817-827.

Mason, W. M., Wong, G. Y., and Entwistle, B. (1984) .
Contextual Analysis through the Multilevel Linear Model.
In S. Leinhardt (Ed.) . Sociological Methodology, 72—103.
San Francisco: Jossey-Bass.

 

Morris, C. N. (1987). Comment on Tanner and Wong (1987).
Journal of the American Statistical Association. 82, 542-
543.

Rao, C. R. (1970). Estimation of Heteroscedastic‘Variances in
Linear Models. Journal of the American Statistical

 

Association. 56, 161-172.

Rao, C. R. (1972). Estimation of Variance and Covariance
Components in Linear Models. Journal of the American
Statistical Association- 67, 112-115.

Raudenbush, S. W. (1984). Application of a Hierarchical
Linear Model in Educational Research . Doctoral
dissertation. Harvard University, Graduate School of

Education.
Raudenbush, S. W. (1988). Educational Applications of
Hierarchical Linear Models: A. Review. Journal of

 

Educational Statistics. 13, (2), 85-116.

Raudenbush, S. W., and Bryk, A. S. (1985). Empirical Bayes
Meta—Analysis. Journal of Educational StatisticsL 10,
(2), 75-98.

 

 

 

 

195

Raudenbush, S. W;, and Bryk, A. S. (1986). A Hierarchical
Model for Studying School Effects. Sociology of
Education, 5 , 1-17.

Raudenbush, S. W., and Bryk; A. S. (1987). Examining
Correlates of Diversity. Journal of Educational
Statistics; 12, (3), 241-269.

Raudenbush, S. W}, and Bryk, A. S. (1988). Methodological
Advances in Analyzing the Effects of Schools and
Classrooms on Student Learning. In E.Z. Rothkopf (Ed.),
Review of Research in Education 423-475. Washington, DC:
American Educational Research Association.

Robinson, W. S. (1950). Ecological Correlations and the
Behavior of Individuals. American Sociological Review,
15, 351-357.

Rosenberg, B. (1973). Linear Regression with Randomly
Dispersed Parameters. Biomstrika. 60, 61-75.

Rubin, D. B. (1981). Estimation in Parallel Randomized
Experiments. Journal of Educational Statistics, 6, (4),
377-400.

Rubin, D. B. (1983). Some Applications of Bayesian Statistics
to educational data. The Statistician._32, 55-68.

Rubin, D. B. (1987). Multiple Imputation for Nonrespsnse in
surveys. New York: Wiley and Sons.

Searle, S. R. (1971). Topics in Variance Component
Estimation. Bigmetrics. 27, 1-76.

Seltzer, M. H. (1988). The Use of Data Augmentation in
Fitting Hierarchical Models to Educational Data.
Doctoral Dissertation. The University of Chicago.

Seltzer, M. H. (1993). Sensitivity Analysis for Fix Effect in
Hierarchical Model: A Gibbs Sampling Approach. Journal
of Educational Statisticsy, 6, (4), 377-400.

Singh, B., and Sedransk, J. (1988). Variance Estimation in
Stratified Sampling When Strata Sample Sizes are Small.
Sankhya: The Indian Journal of Statistics. 50- Series B.
Pt. 3, 382-393.

Smith, A. F. M; (1973). A General Bayesian Linear Model.
Journal of the Royal Statistical Society. Series B. 35,
67-75.

Strenio, J. F., Weisberg, H. 1., and Bryk, A. S. (1983).

 

 

 

196

Empirical Bayes Estimation of Individual Growth—curve
Parameters and Their Relationship to Covariates.
Biometrics, 3 , 71—86.

Tanner M. A. (1993). Tools for statistical inference.
Observed data and data augmentation methods. Lecture
Notes in Statistics 67. New York, Heidelberg: Springer-
Verlag.

Tanner, M. A., and Wong, W. H. (1987). The Calculation of
Posterior Distributions by Data Augmentation. Journal of
the American Statistical Association, 82, (398), 528—540.

Winer, B. J. (1971). Statistical principles in experimental
design. Second Edition. NeW‘ York: McGraw—Hill Book

Company.

Zeger, S. L., and Karim, M.R. (1991). Generalized Linear
Moddels with Random Effects; A Gibbs Samplin Approach.

Journal of the American Statistical association, 8 ,
(413), 79-86.

NICHIGQN STQTE

1|)le 1| IHIWINWIWIHM

 

2.... :, .éfi
gw—q- .‘ 1L. —.

MP- 2
‘3’.