Both µ,d_i T_k,(d_T)ik are fixed effects, whereas h_ij and e_ijl are random effects.

The assumptions are

the h_ij are iid N (0, s)

the e_ijkl are iid N (0, s²)

the h_ij and e_ijl are independent

Thus, the random effect h_ij describes again the effect of a particular herd on the response variable.

Comparing two drugs at two different doses: a split plot design

The data structure of Example 1.2.5 can be described by the following statistical model

Y_ijk = µ + d_i + r_j + h_ij + T_k + (_dT)_ik + e_ijk

where, for i = 1, 2; j = 1, ..., 6; k = 1, 2,

Y_ijk is the difference in PCV (=PCV_after - PCV_before)

µ is the overall mean

d_i is the effect of the Berenil (d₁) or Samorin (d₂)

r_j is the random effect of region j

h_ij is the random effect of the herd within region j to which level i of the drug has been assigned

T_k is the effect of the high (T_i) or low (T₂) dose

(d_T)ik is the interaction effect between drug and dose

e_ijk is the random error of the average response of animals from region j to which the ik drug-dose treatment has been applied

Both µ., d_i, T_k and (d_T)_ik are fixed effects, whereas r_j, h_ij and e_ijk are random effects; the assumptions are

the r_j are iid N(0, s )

the h_ij are iid N(0, s)

the e_ijk are iid N(0,s²)

the r_j , h_ij and e_ijk are independent

Thus the random effect r_j describes the effect of a particular region on the response variable. The random effect h_ij describes the effect of herd within a region on the response variable.

Variation in ELISA data

The data structure of Example 1.2.6 can be described by the following statistical model

Y_ijk = µ + a_i + t_j + e_ijk

where, for i =1, ..., 4; j = 1, ...,11; k = 1, 2,

Y_ijk is the OD value

µ is the overall mean

a_i is the effect of animal i

t_j is the effect of time j

e_ijk is the random error of measurement k at time j on animal i

In this case a_i is assumed to be a random effect; provided there is no systematic pattern of changes in OD with time, we can also assume that t_j is a random effect. Thus, for the random effects a_i, t_j and e_ijk we assume

the a_i are iid N(0, s)

the t_i are iid N(0, s)

the e_ijk are iid N(0,s ²)

the a_i, t_i and e_ijk are independent

Since the overall mean is the only fixed effect, this is a random effects model. The interest of the investigator is in the magnitudes of the different sources of variability on the OD value.

Helminth-resistance in sheep

The data structure of Example 1.2.7 can be described by the following statistical model

Y_ij = µ + s_i + e_ij

where, for i = 1, ..., 8; j = 1, ..., n_ii,

Y_ij is the weaning weight

µ is the overall mean

s_i is the effect of the sire i

e_ij is the random error of lamb j of sire i

In this case, there is only one fixed effect, the overall mean µ. The other effects s_i and e_ij are random effects.

The assumptions are

the s_i are iid N(0, s )

the e_ij are iid N(0,s²)

the s_i and e_ij are independent

The model is a random effects model. For this example, interest is in the sources of variability as well as in actual responses of the specific sires. It will be shown later how these two goals can be met together.

Breed differences in susceptibility to trypanosomosis

The data structure of Example 1.2.8 can be described by the following statistical model

Y_ikl = µ + ,ß_i0 +ß_il ^t + a_ik + e_ikl

where, for i = 1, 2; k = 1, ..., 6; l = 1, ...14,

Y_ikl is the PCV

µ is the overall mean

ß_io is the intercept of the regression line for breed i

ß_il is the slope of the regression line for breed i

a_ik is the effect of animal k belonging to breed i

e_ikl is the random error of measurement l on animal k belonging to the breed i

There are thus five fixed effects, the overall mean and the intercepts and the slopes for the two breeds. Both a_ik and e_ikl are random effects. The assumptions for the random effects are

the a_ik are iid N(0,s)

the e_ikl are iid N(0,s²)

the a_ik and e_ikl are independent

2.3.3 The mixed model description in matrix form

The mixed model examples, given in Section 2.3.2, can also be written in matrix notation. The matrix form of a mixed model collects the fixed effects in a vector , and the random effects in a vector u, and finally the random error term, which is also a random effects factor, in the vector e. A formal definition is

Y = Xß+Zu + e

with X the known design matrix for the fixed effects and Z the known design matrix for the random effects.

Example 2.3.3.1

The matrix representation for Example 1.2.2 is as follows

The example contains the mean, the dose factor which is a fixed effects factor at two levels, the herd factor which is a random effects factor at 3 levels and finally the random error term. The overall mean µ and the two fixed dose effects, t_l and t₂, are contained in the ß vector, whereas the three random herd effects h_l, h₂ and h₃ are contained in the u vector.

Example 2.3.3.2

Example 1.2.4 contains the dose factor and the drug factor which are both fixed effects factor at two levels, the herd factor which is a random effects factor at 4 levels and finally the random error term. The overall mean µ., the two fixed dose effects, T_l and T₂, the two fixed drug effects, d₁ and d₂ and the interaction effects (d_T)₁₁, (d_T)₁₂, (d_T)₂₁ and (d_T)₂₂ are contained in the ß vector, whereas the four random herd effects h₁₁, h₁₂, h₂₁ and h₂₂ are contained in the u vector.

This data structure can again be put in matrix form Y = Xß + Zu + e. The different parts are shown on the two next pages.

The X matrix can be divided in four parts, the first part is the vector Xµ, consisting of ls. The second part is the X_d matrix containing the two vectors corresponding to d₁ and d₂. The third part is the X_T matrix containing the two vectors corresponding to T_l and T₂. The last part is the X_y matrix containing the four vectors corresponding to (d_T)₁₁, (d_T)₁₂, (d_T)₂₁ and (d_T)₂₂.

2.4 Variance component estimation and blup

2.4.1 Defining variance components

It has been demonstrated how different variance components describe at different levels the variability in observations with a mixed model structure. The variance-covariance matrix of Y where

Y = Xß  + Zu + e

is obtained as follows.

The first part of the model, Xß , does not contribute to the variance of Y as it only contains fixed effects. Furthermore, the assumptions are made that each element u₁ in the vector u is a random effect which comes from a normal distribution with mean 0 and variance s, and that the elements are independent from each other and that the covariances among the elements of u are thus zero. The variance-covariance matrix for u is therefore given by a diagonal matrix D(u). All the elements in the vector u are assumed to be independent from the elements in e. We further have that the elements in e are also normally distributed with mean 0 and variance s² and are independent from each other. Given these assumptions, the variance-covariance matrix of Y is then given by

D(Y) = D(Zu) + D(e) = ZD(u)Z^T + s² I_N

Example 2.4.1.1

We revisit Example 1.2.5. The vector u consists of 18 random effects, 6 random effects for the regions (r_l .... r₆), 12 random effects for the herds (h₁₁, h₂₁, h₁₂, h₂₂, ..., h₁₆, h₂₆) whereas the Z matrix is a (24 x 18) matrix.

Note that the matrix product can also be written as the sum of two different matrix products. The first 6 entries of the u vector and the first 6 columns of the Z matrix can be defined as u_r and Z_r respectively and the remaining entries of a and Z as u_h and Z_h, respectively. Thus, u_r and Z_r contain all the information with regard to the regions, whereas u_h and Z_h contain all the information with regard to the herds within regions. Thus

For the vector of observations sorted by region, given by

Y^T = (Y₁₁₁ , Y₁₁₂ , Y₂₁₁, Y ₂₁₂, ..., Y₁₆₁, Y₁₆₂, Y₂₆₁, Y₂₆₂), the matrix product Zu is shown below.

As it is assumed that u_r and u_h are independent from each other, the variant of Zu can be obtained as

D(Zu) = D(Z_ru_r) + D(Z_hu_h) = Z_rD(u_r)Z + Z_hD(u_h)Z

Given that the different entries in u_r are independent from each other, the variance-covariance matrix of u_r. is given by

In a similar way, the variance-covariance matrix of u_h can be obtained as

D(u_h) = sI₁₂

Therefore, the variance-covariance matrix of Y is given by

Both Z_r and Z_h can be written in a simpler form as (Appendix C)

Define Y_j as the vector containing the four observations from region j, Y = (Y_ljl,Y_lj2,Y_2j1,Y_2j2)

For D(Yj) = W, the variance-covariance matrix of Yj, we have

Thus, the variance of an observation corresponds to s²+s+s, the covariance between two observations from the same herd equals s +sand finally, the covariance between two observations from the same region, but from different herds, corresponds to s.

The full variance-covariance matrix is then equal to

Thus, the different variance components prevail in the variance-covariance matrix of the vector Y. Actually, the only variable elements in the matrix V are the different variance components. In order to derive an estimate of the matrix V, estimates of the variance components axe needed for insertion into V. In the next subsection it is discussed how estimates of variance components can be obtained.

2.4.2 Estimating variance components

There exist different methods to estimate variance components. Only likelihood based estimators are considered here (Harville, 1977). An excellent monograph on likelihood is written by Edwards (1972). Estimation methods for variance components based on other criteria are well described by Searle et al. (1992).

We assume the following distributional assumptions for the observation vector Y

Y ~ MVN(Xß, V)

The likelihood of Y is then given by

where ? V?  means the determinant of the matrix V (see Appendix C).

The only unknowns in V are the variance components s, ..., s, with d the number of variance components.

Thus the log-likelihood can be rewritten as

l_Y (, V) = l_Y (, s, ..., s)

Maximum likelihood estimators can then be obtained by maximising this expression with regard to ß and s, , s Differentiating the log-likelihood with respect to the different parameters and setting the result equal to 0 leads to the maximum likelihood estimator for each parameter.

Differentiating the log-likelihood with respect to ß gives

Thus a solution for ß can be obtained as a function of the unknown variance components. This solution can then be inserted into the remaining differential equations in order to lead to a solution for the variance components. In some cases, however, the value that maximises the log-likelihood turns out to be negative (Searle et al., 1992). Although sometimes a negative variance component makes sense (Henderson, 1953; Nelder, 1954), e.g., if the variance component corresponds to a correlation (see Example 3.3), in most cases the parameter space for the variance components is restricted to the interval containing only positive numbers and zero. If a variance component becomes negative, it is usually set to zero and it can be shown (Searle et al., 1992) that, within the restricted region, the zero estimate corresponds to the maximum likelihood estimate. SAS does this by default (SAS Institute Inc. Cary, NC, 1996; Littell et al., 1996). In what follows, the term maximum likelihood estimators is used for the maximum likelihood estimators in the restricted parameter space.

The maximum likelihood estimators of the variance components are biased. The estimate will on average be smaller than the population value. The reason for this is that in the estimation procedure the existence of fixed effects is ignored, and the degrees of freedom used in deriving the estimators do not adjust for this.

The most straightforward example is the estimator of the variance of a sample Y^T = (Y₁,...,Y_N) of a univariate normal distribution N(µ,s ²). The unbiased estimator of s² is given by

If, however, the maximum likelihood criterion is used, the estimator of s² is given by

Thus, this estimator underestimates s². For the unbiased estimator, one degree of freedom is subtracted from N to account for the one degree of freedom used in calculating Y to estimate µ.

Patterson and Thompson (1971) developed a method to derive unbiased estimates of variance components based on the maximum likelihood principle; it is called restricted maximum likelihood estimation. The restricted maximum likelihood method is based on the likelihood of a vector whose components are independent linear combinations of the observations. The basic idea is to end up with a random vector that contains all the information on the variance components but no longer contains information on the fixed effects parameters. We now give a technical description of this method.

Since X is the design matrix for the overparameterised model, it is not of full rank. The rank of X is the number of cell means in the model (the total degrees of freedom associated with the fixed effects). We call this rank r, i.e., rank(X)=r.

A set of linear combinations, denoted as KY, is derived, where K is a (N. r, N) matrix of full rank (see Appendix C)

and which obeys the following equation

E(KY) = 0

As E(KY) = KXß, it follows that for each individual linear combination

k_j X = 0

The individual linear combinations kjY axe called error contrasts (Harville, 1974; Harville, 1977).

As Y is assumed to be normally distributed MYN(X,V), it follows that

KY ~ MVN(0, KVK^T)

The restricted maximum likelihood criterion is based on the likelihood of KY. This likelihood does not contain fixed effects and furthermore contains fewer terms (N . r) as compared to the likelihood of the observation vector Y. Maximising the likelihood of KY will yield the restricted maximum likelihood estimators for the variance components.

Example 2.4.2.1

Assume that a sample is taken from a large herd to obtain an estimate of the average body weight of animals in that herd. Apart from the estimate of mean body weight also its variability is of interest. The weight of an animal can be assumed to be normally distributed N(µ,s²). Denote the sample as Y = (Y₁, ..., Y_N)^T; the observations are assumed to be independent. Therefore the variance covariance matrix of Y is

V =s²I_N

and

Y = MVN(µ1_N, s²I_N)

The log-likelihood of Y is given by

Differentiating this expression first with respect to µ, we obtain

and setting it equal to 0 gives

(Y . µ1_N)1_N = 0

leading to

Differentiating the log-likelihood of Y with regard to s², we obtain

When using the restricted maximum likelihood criterion, a matrix K must be found that satisfies the aforementioned specifics. One such K (of dimension (N . 1, N)) is given by

K = (I_{N . 1} . N^{. 1}J_{N. 1} .N ^{. 1}1_{N . 1})

It can easily be shown that the (N - 1) rows are linearly independent, and that each row satisfies k_j X = 0 and therefore

Differentiating this log-likelihood expression with regard to s² and setting it equal to 0 gives

It can be shown (Searle et al., 1992) that, in general,

K^T (KK^T)^{. 1} K = I_N . X(X^TX)^{. 1}X^T

for any appropriate choice of K.

In our example for body weight, the X matrix consists simply of a vector of 1's of length N, 1_N.

Thus

X^T X = N and (X^TX)^{. 1} = 1/N

and furthermore

XX^T = J_N

Therefore

K^T (KK^T) ^{. 1} K = I_N .

Substituting in the formula, for s² gives

Example 2.4.2.2

We revisit Example 1.2.5. There are three different variance components. There is a first variance component associated with the variability among regions, s. A second variance component, s, is associated with the variability among herds within regions. The third variance component,s ², models the variability of the random error component.

Estimates of these three variance components can be obtained using the maximum likelihood (ML) or the restricted maximum likelihood (REML). The table below shows the ML-estimates and the REML-estimates. It provides an illustration of the underestimation of the variance components by the maximum likelihood method. The ML- and the REML- estimates can be obtained by the SAS program shown in Appendix B.2.2.

Table 2.4.2.1 Variance components estimates for Example 1.2.5

Variance component	ML-estimates	REML-estimates
s	4.292	5.151
s	0.322	0.387
s2	1.747	2.096

2.4.3 Best linear unbiased prediction of random effects

In some experiments, the vector of random effects itself is of interest as well as the variance component based on this vector. Example 1.2.7 on the inheritance of helminth-resistance in sheep can be used to illustrate this. In the initial analysis estimation of the variance components s and s² is of main interest. The larger s is relative to s² the greater is the evidence for genetic inheritance.

Subsequently, we need to compare sires We are thus interested in prediction of the values of the random effects s_i, i = 1, ..., 8.

There are two different sources of information contributing to the random effect s_i, the observed values of the lambs of sire i,Yi1,... Y_ini and the distributional assumption s_i ~ N(0,s ). In order to use both sources of information, a random effect is predicted from its expected value, conditional on the mean of the relevant observations.

The random effects model of Example 1.2.7 was given by

Y_ij = µ + s_i + e_ij

If the mean response value of the offspring from sire i is given by , then the expected value of the random effect, conditional on the relevant observations, is given by

The joint distribution of s_i and . is given by

It can be shown (see Appendix D), given this joint distribution, that

As s² tends to zero or n_i s tends to infinity, this expression is approximately equal to (Y_i, . µ). Otherwise, the expression requires REML-estimates of the variance components and the ML-estimate of µ, in order to obtain an estimate _i. As will be shown in the next section, the ML-estimator for µ can be obtained by

Example 2.4.3.1

The random effects model presented above is fitted to the data of Example 1.2.7. The REML estimates for the variance components are given by s = 3.685 and s² = 3.542. The ML-estimate of the mean µ is given by

The best linear unbiased prediction for the random effect of sire 4 is thus given by

The best linear unbiased predictions for the 8 sires in the data set can be obtained by the SAS program shown in Appendix B.2.3

2.5 Statistical inference for fixed effects

2.5.1 Estimation of fixed effects and their variance

Estimates of the fixed effects can be obtained by maximising the log-likelihood of Y

Differentiating the log-likelihood with respect to and setting it equal to 0

If the variance components contained in V are known and can be substituted into the expression, this estimator is called the generalised least squares estimator.

The variance-covariance matrix of ß is given by

However, the variance components contained in V are usually not known. Therefore, the REML-estimates need to be substituted, giving the estimate V, and the estimate based on the maximum likelihood criterion can then be obtained by the expression

= (X^TV^{. 1}X)^. X^TV^{. 1}Y

Furthermore, D() is estimated by (X^{T. 1}X) ^. X^{T . 1}X (X^{T. 1}X) . This estimator is biased downwards (i.e. is on the average smaller than the population value (X^TV^{. 1}X) ^. X^TV ^{. 1}X (X^TV^{. 1}X) . ) since the uncertainty introduced with the REML-estimates of the variance components in V is not taken into account in the approximation for (X^TV^{. 1}X) (Kackar and Harville,l984).

Example 2.5.1.1

We revisit Example 1.2.5. The mixed model is given by

Y_ijk = µ + d_i + r_j + h_ij + T_k + (d_T)_ik +e_ijk

Since the model is overparameterised, the solution for ß is not unique; there exist an infinite number of solutions (see Section 2.2). This can also be seen by the fact that the matrix X^TV ^{. 1}X is not of full rank and that we thus have to use generalised inverses. One solution is to set some of the parameter levels equal to 0, thereby reducing the number of parameter levels. In this example, the following restrictions are used

d² = 0; T₂ = 0; (d_T )₁₂ = 0; (d_T)₂₁ = 0;(d)₂₂₌ 0

Therefore, a maximum likelihood estimate of the parameter vector

ß^T = (µ, d₁, d₂, T₁, T₂, (dT)₁₁, (dT)₁₂, (dT)₂₁, (dT)₂₂)

can be shown, following substitution of  for V, to be

^T = (17.13, 7.15, 0, 4.35, 0, . 3.18, 0, 0, 0)

and the variance-covariance matrix

These results can be obtained by the SAS program shown in Appendix B.2.4.

2.5.2 Estimation of linear combinations of fixed effects and their variance

Interest is not in the individual parameters of an overparameterised model (as they are not estimable), but in estimable functions of those parameters c^Tß so that the solution is independent of the generalised inverse. A simple example of such an estimable function is the mean of one level of a parameter (e.g. a treatment) averaged over other parameter levels. Both the estimate and the variance of the estimate of a linear combination of the parameters can easily be obtained if the vector of parameter estimates and its estimated variance-covariance matrix are available.

Assume that c^Tß is an estimable function. It is easy to show that the maximum likelihood estimator of c^Tß is given by c^Tß; the estimated variance (see Appendix E) is

Example 2.5.2.1

As already mentioned, the cell means, being linear combinations of the parameters, are relevant and estimable parameters. In Example 1.2.5 the expected value for the first cell mean is given by the following linear combination

It represents the mean response of animals treated with a high dose of Berenil.

With c = ( 1 0 1 0 1 0 0 0 1 ) we obtain

This cell mean represents the mean response of animals treated with the low dose of Samorin.

Estimates for the cell means are then obtained by substituting the parameter

estimates obtained in Example 2.5.1.1

The variance of a cell mean ₁₁  can also be obtained.

For instance, the estimated variance of ₁₁ is given by

Vâr(c = c()c₁=1.27

This shows that the estimates for the cell means model can easily be obtained from the analysis of the overparameterised model. The SAS program that derives the cell mean estimates with their estimated variances is shown in Appendix B.2.5.

2.5.3 Confidence intervals

We often like to determine an upper and lower confidence limit for an estimable parameter, say c^Tß.. Confidence intervals provide an appropriate probabilistic tool for this. It leads to the determination of a lower and an upper confidence limit for the parameter function under consideration. The distributional properties of a specific linear combination, needed to derive confidence intervals, follows easily from the fact that the estimators are linear combinations of the observations Y. Since Y is multivariate normal,. the estimator c^T also follows a normal distribution.

c^T ~ N(c^T, c^Tß, cT (X^T V^{. 1} X) ^. c)

In general, however, the variance components in V must be estimated. So c^T  is assumed to follow a t-distribution with v degrees of freedom, where v corresponds to the degrees of freedom available to estimate the variance c^T (X^T V^{. 1}X) ^_ c.

In some situations, the calculated variance corresponds to the expected mean square of one of the random effects factors. In this case, the degrees of freedom are the degrees of freedom of the mean square itself, which can easily be determined (see Example 2.5.3.1).

In general in a mixed model, however, the variance of a linear combination is a complicated expression of different expected mean squares, and the degrees of freedom can only be approximated.

The Satterthwaite procedure (Satterthwaite, 1946), based on quadratic forms (mean squares typically are quadratic forms), provides a way to obtain approximate degrees of freedom and is given by

In this expression, the variance term in the denominator must be estimated. An approximation for this variance term has been proposed by Giesbrecht and Burns (1985) and McLean and Sanders (1988). The derivation of both the expression for the degrees of freedom and the approximation for the variance is discussed in detail in Appendix F.

If c^T (X^TV^{. 1}X)^. c can be written as a linear combination of mean squares, i.e. a₁MS₁ +...+ a¹MS_l with a_l, ..., a_l constants and MS_j the mean square due to the j^th source of variation in the analysis of variance, then the Satterthwaite procedure simplifies to give

with df_j the degrees of freedom for MS_j.

The (1 - a)100% confidence interval for c^Tß is then given by

Example 2.5.3.1 Confidence intervals

In the analysis of the data of Example 1.2.5, different linear combinations might be of interest. Two linear combinations cß ~  and cß have already been shown in Section 2.5.2.

If the parameter vector is given by

ß^T = ( µ d₁ d₂ T₁ T₂ (d_T)₁₁ (d_T)₁₂ (d_T)₂₁ (d_T)₂₂ )

some other examples of linear combinations are

1. The mean response for animals treated with Samorin, µ2,, given by

cß = ( 1 0 1 0.5 0.5 0 0 0.5 0.5) ß

2. The mean difference between animals treated with Samorin and Berenil, µ₂ . µ_1., given by

cß = (0 -1 1 0 0 -0.5 -0.5 0.5 0.5) ß

3. The mean response for animals treated with a high dose of Samorin, µ₂₁, given by

cß = ( 1 0 1 1 0 0 0 1 0) ß

4. The mean difference between animals treated with a high versus a low dose of Samorin, µ₂₁ . µ₂₂, given by

cß= (0 0 0 1 . 1 0 0 1 . 1 ) ß

It can easily be shown that an estimate for cß is given by cß= ₂.., with _2.. the mean of the 12 observations on Samorin. The corresponding variance is

Similarly, variances for the other linear combinations can be obtained and are given by

By substituting the REML-estimates of the variance components ( = 5.15; = 0.387; ² = 2.096) in these expressions, the following estimates for the variances of the linear combinations are obtained

Vâr(c) = 1.0975

Vâr(c) = 0.478

Vâr(c) = 1.272

Vâr(c) = 0.698

In order to derive confidence intervals for these linear combinations, based on the t-distribution, the degrees of freedom associated with the variance term need to be calculated.

An appropriate tool to derive these degrees of freedom is the table of the different mean squares with their corresponding degrees of freedom and expected values, which is given below for Example 1.2.5.

Table 2.5.3.1. Table of mean squares and their degrees of freedom an expected values foe Example 1.2.5

where f (factor) is a term in the expected mean squares due to an individual fixed effect. The expected mean squares for this example can be obtained by the program given in Appendix B.2.6.

The derivation of the degrees of freedom is only straightforward for cß and cß. They have 5 and 10 degrees of freedom as their variances are only dependent on the expected mean square for the herd (5 df) and error (10 df), respectively.

The variances of the other linear combinations can be obtained by combining the expected mean squares of the different factors in the following way.

The variance of c= is (2s + 2s + s²)/12. The numerator can be obtained from 0.5E(MS_region)+0.5E(MS_herd). Using the simplified Satterthwaite formula, with a_l = a₂ = 0.5, we obtain

For c. the variance is (s + s + s²)/16; the numerator can be obtained from

0.25E(MS_region)+0.25E(MS_herd)+0.5E(MS_error). Therefore

Thus, 95% confidence intervals are given by

From these confidence intervals, it can be concluded that there is a large and significant difference between animals treated with Samorin and Berenil (cß.). There is, however, no difference between the high and low dose of Samorin (cß.). These estimates and their standard errors can be obtained by the SAS program given in Appendix B.2.7.

2.5.4 The appropriate inference space

From the example in the previous section it is clear that the variance of c^Tß. with c^Tß. estimable, is a linear combination of the variance components present in the mixed model. In some situations it is meaningful to narrow the scope of the study and to think about one or more of the random effects as fixed effects. As a consequence, the corresponding variance components will not appear in the expression of the variance for the estimator c^Tß.. If we think of one or more of the random effects as fixed effects we say that we consider a narrow inference space. If we assume all effects to be random we talk about the broad inference space. Based on examples we show that the confidence limits for c^Tß. are quite different for broad and narrow inference spaces. If we consider all factors as fixed effects the only variance component is s², the error term variance. It typically follows that the narrower the inference space the smaller will the confidence interval be.

Example 2.5.4.1 Confidence intervals in the fixed effects and mixed model

It has been shown for the data of Example 1.2.5 that the variance of the mean response of animals treated with Samorin, µ₂, is given by

The variance of ₂ ..., the estimate of µ₂, in the fixed effects model, i.e., the model with region and herd as fixed effects, is

This estimate is much smaller than that calculated previously for cß. when r_j and h_2j are assumed random, and therefore the width of the confidence interval for µ₂. is much narrower.

The estimated standard errors for the linear combination µ₂, can be obtained by the SAS program given in Appendix B.2.8 for both the fixed effects model and the mixed model.

2.6 Hypothesis testing

Testing general linear hypotheses for fixed effects by an appropriate F-test is much more complex in mixed models than in fixed effects models. As has already been shown, there is only one source of random variation in fixed effects models. The F-statistic consists of a ratio of two mean squares, the denominator being the error mean square. For that reason the derivation of the degrees of freedom of the denominator of the F-test is straightforward, as shown before. In the mixed model, however, there are different sources of random variation, and the denominator of the F-test consists in most cases of a combination of these different sources of random variation. The derivation of the appropriate degrees of freedom for the denominator in the F-test is therefore not straightforward, and can, in most practical cases, only be approximated by the Satterthwaite procedure as described in Appendix F.

General linear hypotheses for the fixed effects take the form

H₀ : C^Tß. = 0 versus H_a : C^Tß.? 0

with ß a (b x 1) vector and C a (b x q) matrix of v_i independent and estimable functions. Thus (rank(C)=v_l).

To estimate C^Tß we use C^Tß and since

D(C^Tß) = C^T D(ß)C = C^T (X^T V ^{. 1}X) ^.C

we use C^T (X^T V ^{. 1}X) ^. C as the estimated variance-covariance matrix.

It can be shown that the statistic

approximately follows a F-distribution under the null hypothesis with v₁ as numerator degrees of freedom.

As can be seen from this expression, F is no longer a ratio of two mean squares as in the fixed effects model. The role of the denominator is taken over by (X^TV^{. 1}X). The denominator degrees of freedom, v_2, thus correspond to the degrees of freedom associated with (X^TV^{. 1}X) and can be approximated by the Satterthwaite procedure in a similar way as shown for the t-test. Further details on these approximate F-tests for general linear hypotheses are in Fai and Cornelius (1996).

The P-value for testing the null hypothesis is then given by the probability that Fv_1, v₂ is larger than the corresponding value of the F-distribution, i.e.

Prob (Fv_1,v₂ = F)

Example 2.6.1 Hypothesis tests

The null hypothesis and alternative hypothesis of no difference between Samorin and Berenil in Example 1.2.5 are given by

H₀ : µ₂. = µ₁. and H_a : µ₂ ? µ₁.

Thus in this simple case, where C is a (q x 1) vector, the null hypothesis C^Tß = 0

can be written as

follows a F-distribution with 1 and 5 degrees of freedom (with 5 being the degrees of freedom associated with region as calculated before for cß). The estimate for C^Tß equals . 5.55 and an estimate for its variance is given by 0.478.

 Therefore, the P-value for the null hypothesis is given by

In testing a more complex hypothesis in which C is a matrix rather than a vector, assume that, for Example 1.2.5, interest is in the following hypothesis

H₀ : µ₁₁ = µ₁₂ = µ₂₁ = µ₂₂ versus H_a : not all µ_ij equal

We can rewrite this hypothesis by building up three simultaneous independent estimable functions and testing them together

The first estimable function tests whether there is a difference between drugs (=µ₁. . µ₂.), the second between doses (=µ_.1 . µ_.2) and the third whether there is an interaction between dose and drug (=(µ₁₁. µ₁₂) . (µ₂₁ . µ₂₂)).

It can be shown that the value of the F-statistic is 31.2 and that the denominator degrees of freedom equal 7.14. The rank of C and thus the numerator degrees of freedom equal 3. Therefore, the P-value is Prob(F_3,7 = 31.2) = 0.0002. The null hypothesis is clearly rejected. The corresponding SAS program is given in Appendix B.2.9.