## Pseudolikelihood

This method can be used to fit both random effects models and covariance pattern models. It was proposed by Wolfinger and O'Connell (1993) and is used within PROC GLIMMIX (see Section 9.3). Pseudo-likelihood maximises the quasi-likelihood by iteratively analysing a linearised pseudo variable (i.e. a transformation of y onto the linear scale) using weighted normal mixed models. The method is referred to as 'pseudo-likelihood' because the likelihood function maximised at each iteration is that of the pseudo variable and not that of the original data. The 'pseudo variable' introduced in Section 3.1.5 based on a first-order Taylor series expansion is again used:

z has variance

Vz = var(Xa + ZP) + B—1var(y — ^)B—1 = ZGZ' + B—1RB—1.

By rewriting the residual matrix R as a product of a correlation matrix on the linear scale, P, and AB-1, Vz can be re-expressed as

In random effects and random coefficients models the residuals are uncorrelated and P = I,so Vz then simplifies to

z|P has the following multivariate normal distribution:

(conditioning the z on P allows ZGZ' to be omitted from the variance matrix formula).

z can now be analysed as a weighted normal mixed model with residual matrix P (defined on the linear scale), and diagonal weight matrix A-1B (inverse product of pre- and post-multipliers of P, A1/2B-1/2B-1/2A1/2). Because z and B are dependent on the estimates of a and P, the normal mixed model needs to be fitted iteratively. Any of the methods described in Section 2.2 can be used to do this. However, again, variance parameters will be biased downwards to some extent if ML or IGLS are used, but this problem is not expected with REML and RIGLS. We will now define the iterative procedure more explicitly.

The iterative procedure The raw data, y, can be taken as initial values for z and the identity matrix, I, can be used for the initial weight matrix. Alternatively, g(y) can be taken as initial values for z, with an adjustment if necessary to prevent infinite values. A weighted normal mixed model is then fitted using a method such as REML. This completes the first iteration and provides initial estimates for a and P. New values of z and B are calculated using a and P, and from these a second weight matrix, A-1B. A second weighted mixed model is then fitted and the process is continued until the parameter estimates converge (i.e. until parameter values change very little between successive iterations). This method is computationally costly because normal mixed models, themselves requiring iterations, are fitted iteratively.

Fixed and random effects variance Once the model has converged, the fixed and random effects and their variances can be calculated using the formula specified for normal mixed models in Sections 2.2.3 and 2.2.4:

a = (X'VZ 1X)-1X'VZ 1y, var(a) = (X'V-1X)-1, and for random effects

P = GZ'V-1(y - Xa), var(P) = GZ'Vz-1ZG - GZ'Vz-1X(X'Vz-1X)-1X'Vz-1ZG.

Just as in normal mixed models, because Vz is estimated and not known there can be a small amount of downward bias in var(a) and var(P) (see Sections 2.4.3 and

The results obtained from generalised estimating equations and pseudo-likelihood are expected to be similar; however, there are often slight differences due to the different computational approaches.