AN ORTHODOX BLUP APPROACH TO GENERALIZED LINEAR MIXED MODELS by Renjun Ma B.Sc, Wuhan University of Water Transportation Engineering, China, 1982 M.Sc, Wuhan University, China, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES Department of Statistics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April, 1999 ©Renjun Ma, 1999 In presenting degree at this the thesis in partial fulfilment of University of British Columbia, I agree freely available for copying of department publication this or of reference thesis by this for his thesis and study. scholarly or for her Department of The University of British Vancouver, Canada DE-6 (2/88) Columbia purposes gain shall requirements that agree may representatives. financial permission. I further the It not be is that the permission granted allowed an advanced Library shall by understood be for for the that without make it extensive head of my copying or my written Abstract We introduce a new class of generalized linear mixed models assuming Tweedie exponential dispersion model distributions for both the response and the random effects. This class of models accommodates a wide range of discrete, continuous and mixed data. By letting the random effects enter as weights as well as means in the conditional distributions, the variance matrix may be expressed as a sum of variance components. We consider an orthodox BLUP approach to parameter estimation and random effects prediction for this new class of models based on a predictor of the random effects that is truly best linear and unbiased, in contrast to the conventional BLUP which is the conditional mode. We obtain an optimal estimating equation based on the orthodox BLUP, which is solved by a modified Newton algorithm. This approach facilitates analysis of residuals and allows justification of asymptotic results under realistic conditions through standard estimating equation theory. An important feature of this approach is that the principal results depend only on the first and second moment assumptions of unobserved random effects. The common fitting algorithm based on orthodox BLUP enables us to study this new class of models as a single class, rather than as a collection of unrelated different models. This approach is illustrated with the analyses of seed germination data, epilepsy data and cake baking data. By means of asymptotic justifications, simulations and worked examples, we conclude that the orthodox BLUP approach is of practical value for analysis of clustered non-normal data. ii Contents Table of Contents iii List of Tables vii List of Figures viii Acknowledgements x 1 Introduction 1 2 Preliminaries 7 2.1 Generalized linear models with random effects 7 2.1.1 Definitions 7 2.1.2 Tweedie exponential dispersion models 9 2.2 2.3 Prediction of random effects 11 2.2.1 Likelihoods 11 2.2.2 Penalized quasi-likelihood and h-likelihood approaches 2.2.3 Comparison for different predictors Estimating functions . . . . 12 13 . 15 2.3.1 Unbiased estimating functions 15 2.3.2 Optimal estimating functions 19 iii 2.3.3 2.4 3 4 5 Nuisance parameter case 19 Data examples 20 2.4.1 Epilepsy data . . 21 2.4.2 Seed germination data 23 2.4.3 Cake baking data 23 Tweedie mixed models 26 3.1 Models 26 3.2 Covariance structure •. 31 3.2.1 Derivation 31 3.2.2 Matrix expressions 33 Orthodox B L U P predictors of random effects 37 4.1 Orthodox BLUP 38 4.2 Random effects predictors 39 4.3 Mean squared distance 42 4.4 Consistency 45 Parameter estimation 49 5.1 Estimation of regression parameters 49 5.1.1 Estimated score function 49 5.1.2 Standard errors of regression parameter estimators 52 5.1.3 Newton scoring algorithm 52 5.1.4 Optimality 54 5.1.5 Derivation of matrix form of ib((3) 57 5.2 Estimation of dispersion parameters 61 5.2.1 61 Adjusted Pearson estimators iv 6 7 Asymptotic properties 63 5.2.3 Heterogeneity 65 Residual analysis and computational procedure 67 6.1 Residual analysis 67 6.2 Computational procedure 69 Conventional, semiparametric and binomial mixed models 71 7.1 Conventional Tweedie mixed models 72 7.1.1 Covariance and variance function 73 7.1.2 Random effects predictors 75 7.1.3 Parameter estimation 77 7.1.4 Adjusted Pearson estimators 78 7.2 8 5.2.2 Random effects beyond Tweedie family 78 7.2.1 Nonparametric random effects 79 7.2.2 Log-normal random effects 81 7.2.3 Binomial mixed models 83 Data analysis 86 8.1 86 8.2 Count data 8.1.1 Overview of the previous analyses 8.1.2 Poisson-Tweedie models 89 8.1.3 Model checking 91 8.1.4 Comparison of different approaches 94 8.1.5 Computational aspects 99 Binomial data 8.2.1 . . . 87 99 Paired Poisson-Tweedie models v 100 8.3 9 8.2.2 Model checking 102 8.2.3 Comparison of different approaches . . . 102 Continuous data . 105 Simulation 112 9.1 Results over 100 simulations 112 9.1.1 Summary statistics 114 9.1.2 Confidence and prediction intervals 116 9.1.3 Normality of parameter estimates and random effects 117 9.2 Residual analysis 120 10 Discussion 124 10.1 Conclusion 124 10.2 Further study 127 10.2.1 More than two levels of random effects 127 10.2.2 Crossed designs 128 10.2.3 Survival data analysis 129 Appendix 132 A Data sets 132 Bibliography 136 vi List of Tables 2.1 Summary of Tweedie exponential dispersion models 10 8.1 HGLM parameter estimates for the epilepsy data 89 8.2 Parameter estimates for the epilepsy data 91 8.3 Parameter estimates for the seed germination data based on paired Poisson mixed models 100 8.4 Parameter estimates for the seed data based on binomial models. . . 8.5 Parameter estimates for cake baking data 107 9.1 Averages of parameter estimates over 100 simulations 114 9.2 Simulated and estimated standard errors 115 9.3 Mean squared errors over 100 simulations 115 9.4 Coverage counts of 95% confidence and prediction intervals 117 A.l Seed germination data (a)seeds: O. Aegyptiaca 75 and 73, (b)root extracts: bean and cucumber 102 132 A.2 Epilepsy data: successive two-week seizure counts for 59 epileptics, (a) Trt: treatment 0=placebo, l=progabide); (b) Base: eight-week baseline seizure counts; (c)Age (in years) A.3 Cake baking data: breaking angles (degrees) vii 133 135 List of Figures 2.1 Overdispersion diagnostic plot for epilepsy data 22 2.2 Scatter plot of proportions of seeds germinated versus plates 24 2.3 Boxplots of cake baking data by batches 25 8.1 Profile plot of log transformed epilepsy seizure counts 8.2 Normal plots of level 1, 2 and 3 residuals for epilepsy data 92 8.3 Scatter plots of level 2 and 3 residuals for epilepsy data 93 8.4 Heterogeneity plots for epilepsy data: (a) scatter plot of within subject . 87 sample variances versus sample means; (b) boxplots of within subject responses for all patients 97 8.5 Heterogeneity plot for epilepsy data 98 8.6 Normal plots of residuals for seed germination data 103 8.7 Scatter plots of residuals for seed germination data 104 8.8 Normal plot of ANOVA model residuals for cake baking data 108 8.9 Boxplots of cake baking data by temperature 109 8.10 Normal plots of residuals for cake baking data 110 8.11 Scatter plots of residuals for cake baking data Ill viii 9.1 Normal plots of parameters over 100 simulations: x and y axes are quantiles of standard normal distribution and ordered parameter estimates over 100 simulations, respectively. . ., 118 9.2 Plots for random effects predictors ui and u 9.3 Comparison of normal plots of level 1, 2 and 3 residuals for epilepsy n over 100 simulations. . data and simulated data 9.4 121 Comparison of histograms of level 1, 2 and 3 residuals for epilepsy data and simulated data 9.5 119 122 Comparison of scatter plots of level 2 and 3 residuals for epilepsy data and simulated data 123 ix Acknowledgements Foremost I would like to thank my supervisor, Dr. Bent J0rgensen, for his excellent guidance, for his inspiration, for his constant encouragement and for his great patience. I am very grateful to other members of my supervisory committee, Dr. John Petkau and Dr. Martin Puterman for their valuable comments and assistance. I am particularly indebted to Dr. John Petkau for his careful reading of the manuscript; his helpful discussion and constructive suggestions lead to the substantial improvement of the presentation of this thesis. I wish to thank my program advisor, Dr. Jean Meloche, for helping me to adapt to the new environment. I would like to thank Dr. Nancy Heckman for her invaluable advice and help, for her constant encouragement and support throughout my program. Many thanks go to Dr. Harry Joe for his useful advice and help throughout my program. Many thanks go to Dr. Li Sun, Hongbin Zhang and Jeevanantham Rajeswaran for their valuable help with Latex and computer programming. I am very grateful to Christine Graham for her constant help and patience. I express my gratitude to Department of Theoretical Statistics, Aarhus University for providing me help and access to research facility during my academic visit. I also acknowledge the support of the University of British Columbia through a University Graduate Fellowship. Chapter 1 Introduction In recent years much effort has been devoted to extending regression methodology to clustered non-normal data. A wide range of applications of generalized linear mixed models to clustered data has been investigated by many researchers (Breslow and Clayton 1993; Lee and Nelder 1996). Traditional generalized linear models (McCullagh and Nelder 1989) have unified regression analysis for a variety of independent responses, but have difficulties in providing valid inference for clustered, and hence correlated, data. Clustered data usually exhibit significant heterogeneity and overdispersion as well (cf. Engel 1987; McCullagh and Nelder 1989, p. 124-126, 198-200). Such cluster effects are often modelled by incorporating random effects into generalized linear models. The application of generalized linear mixed models to clustered data has attracted much attention since the seventies (Crowder 1978; Laird 1978). However, the early development of generalized linear mixed models focussed mainly on relatively simple cases such as random intercept models (Hinde 1982; Williams 1982; Breslow 1984; Stiratelli, Laird and Ware 1984; Anderson and Aitkin 1985; Brillinger and Preisler 1 1986). The modelling of more complicated problems has largely been hampered by the intractability of high-dimensional integrals involved in evaluating the likelihood. To avoid these numerical problems, numerous alternative approaches to the analysis of clustered data have been proposed recently. The most popular approaches are generalized estimating equations (Liang and Zeger 1986; Zeger, Liang and Albert 1988), Bayesian methods using computational techniques (Zeger and Karim 1991), penalized quasi-likelihood (Breslow and Clayton 1993) and maximum hierarchical likelihood (Lee and Nelder 1996) methods. We briefly review these approaches and discuss their advantages and disadvantages. Generalized estimating equations (GEE) approach focuses on the marginal relationship between covariates and clustered responses, but only estimates the covariances as nuisance parameters by adopting a so-called "working correlation". This approach enjoys robustness against mis-specification of covariances, but sometimes suffers from a lack of efficiency due to such incorrect specification (Lipsitz et al. 1994; Fitzmaurice 1995). It is of limited use when the correlation structure is of primary interest. This approach mainly deals with marginal models instead of random effects models. On the other hand, in areas such as genetics, it is often the unobserved genetic effects which are of primary interest (Harville 1977; Clayton 1991; Karim and Zeger 1992). Robinson (1991) presented a variety of applications where the random effects themselves are of interest; therefore random effects models are more relevant in such situations. As a general approach to complicated generalized linear mixed models, Zeger and Karim (1991) proposed to cast these models in a Bayesian framework and approximate maximum likelihood estimates using flat or diffuse priors. However, such 2 approximations are often impossible because the posterior may not exist for such priors (Natarajan and McCulloch 1995; Hobert and Casella 1996). This problem of the posterior may not be detected when using computational techniques such as the Gibbs sampler; therefore misleading estimates may result. On the other hand, recent non-Bayesian approaches to generalized linear mixed models have focussed on the explicit or implicit modification of the E-step in the EM algorithm (Dempster, Laird and Rubin 1977) due to the difficulty evaluating the conditional expectation of the random effects given the data. The most widely adopted technique is the generalization of the linear mixed model equations of Henderson et al. (1975), but with various approximations (Gilmour, Anderson and Rae 1984; Harville and Mee 1984; Schall 1991; Breslow and Clayton 1993; Wolfinger 1993; McGilchrist 1994; Lee and Nelder 1996). It can be shown that these approaches essentially modify the EM algorithm with the conditional expectation of the random effects given the data being replaced by the predictors based on the mode of the corresponding conditional distributions. These modal predictors are often referred to as the 'BLUP' (Schall 1991; McGilchrist 1994), although in general it is neither linear nor unbiased for non-normal distributions. In contrast to the modal predictor approaches, Bayesian approaches enjoy great flexibility in modelling the data with normal or non-normal random effects distributions, but are computationally intensive. The computational time required is sufficiently long as to possibly discourage fitting several different models (Karim and Zeger 1992). This drawback may pose serious problems in practice because it is generally difficult to justify a particular distribution for the unobserved random effects. Data will often point with almost equal emphasis at several different models and it is 3 important that we recognize these models and their possible different interpretations. On the other hand, the modal predictor approaches enable us to explore several different models with reasonable computing time. The penalized quasi-likelihood (PQL) approach deals with models with approximate multivariate normal random effects. The maximum hierarchical likelihood (MHL) approach widens the choice of random effects distributions to include conjugate distributions (George et al. 1987), but fails to model the dependence structure of the random effects. The introduction of flexible, yet tractable distribution classes to model both the distributional shape and dependence structure of random effects is certainly needed to facilitate appropriate inference. The implementation of the Bayesian approach is relatively straightforward in contrast to the requirement for manipulation of large matrices for modal predictor approaches when there are a large number of random effects; however the assessment of convergence of computational techniques such as the Gibbs sampler remains an area of debate (Glifford 1993; Smith and Roberts 1993). On the other hand, non-Bayesian approaches provide a natural framework for model checking, but are forced to rely on asymptotics. In fact, the existing approaches to generalized linear mixed models mainly concentrate on the model fitting part, but to a large extent ignore another important ingredient of the modelling process, the model checking component (Lee and Nelder 1996). The justifications of these modal predictor approaches were generally intuitive until Breslow and Clayton (1993) presented some ad hoc justifications for their penalized quasi-likelihood approach and Lee and Nelder (1996) provided rigorous justification for their maximum hierarchical likelihood approach. The asymptotic results for quite 4 general settings are obtained under certain conditions by Nelder and Lee (1996), but with respect to large cluster sizes with fixed number of clusters. This raised concerns about more practical situations where the number of clusters is large, but with relatively small cluster sizes (Clayton 1996; Engel and Keen 1996). In addition, the estimating equations based on the modal preditors are generally biased (Breslow and Lin 1995; Lin and Breslow 1996a). This thesis considers generalized linear mixed models based on the class of Tweedie exponential dispersion model distributions (J0rgensen 1987a) for both the response and the random effects. This gives a very flexible class of models which includes various combinations of Poisson, normal, gamma, inverse Gaussian, compound Poisson and extreme stable and positive stable distributions. In the context of clustered data, the hierarchical structure is very clear so that the modeled covariance structure should clearly reflect those hierarchies. The dominant tradition in accounting for dependence between or within clusters is to explicitly incorporate random effects into a monotonic transform of the conditional mean ignoring the dispersion components. This approach generally does not lead to variance components decomposition structure for covariance matrix of the response. By incorporating correlated or uncorrelated random effects into both mean and dispersion components, the covariance matrix of our model possesses an interpretable variance components decomposition. The novelty of our approach lies in the introduction of a new unbiased estimating equation, which is based on a modification of the EM algorithm where conditional expectations are approximated by an orthodox BLUP, in the context of generalized linear mixed models. An orthodox BLUP is defined as the best linear unbiased predictor in the literal sense, (cf. Brockwell and Davis 1991 p. 64). The unbiased estimating 5 equations introduced via the orthodox BLUP lead to consistent estimators for both regression and dispersion parameters under practical conditions where the number of clusters is large. The estimating equation for the regression parameters is optimal in the sense of Godambe (1976). While the parametric nature of our models facilitates residual analysis, our estimating procedure also allows a semi-parametric interpretation of the models. Our approach does not require manipulation of large matrices; therefore is computationally simpler than the modal predictor approaches. This approach is applicable to a wide range of clustered discrete, continuous and mixed data. The organization of this thesis is as follows. In Chapter 2, besides a brief introduction of Tweedie exponential dispersion model distributions and some estimating function results, we compare the orthodox BLUP and modal predictors and present . a few data examples to motivate our study. In Chapter 3, we propose a class of nested random effect models and derive their moment structures. The prediction of random effects and the consistency properties of the orthodox BLUP are discussed in Chapters 4. In Chapter 5, we discuss estimation for both regression and dispersion parameters as well as asymptotic properties of these parameter estimators. An outline of the residual analysis and computational procedure is presented in Chapter 6. We address in detail a so-called conventional model and its relationship with other models in Chapter 7. Illustrative examples and simulations are presented in Chapters 8 and 9, respectively. We present a discussion in Chapter 10. 6 Chapter 2 Preliminaries 2.1 2.1.1 Generalized linear models with random effects Definitions We study generalized linear models with random effects based on the class of Tweedie exponential dispersion model distributions. Here wefirstdefine exponential dispersion models. A random variable Y is said to follow reproductive exponential dispersion model ED(n, a ) if its probability density functions can be written in the form: 2 P(V, </>, V) = <f>) exp{0[y rj - K(T,)]}, (2.1) where \x = E[Y] and a = 1/0. Let «'(•), the first derivative of K(-), be denoted by 2 r(-). V(/JL) = T'{T~ 1 (fj,)} is called the variance function. Further, d(y; /*) = 2[y{T-\y) - r " ^ ) } - / . { r " ^ ) } + / . { r " ^ ) } ] 1 1 1 (2.2) is known as the (unit) deviance function. The distribution of Z = <\>Y is called an additive exponential dispersion model, 7 denoted by Z ~ ED*(r],(f)). The justification of terminology 'reproductive' and 'additive' can be found in J0rgensen (1997, p l l ) . Now we can define a generalized linear model with random effects. (Yii,Vim,Y i,Y m responses. m n j n ) be an n = YnLi n^-dimensional vector of observed T Let (3 = (Pi,3 ) be a p-dimensional vector of fixed effects and T P U = (Ui,U ) Let Y = be a m-dimensional vector of random effects. T m We suppose that, given U = u, Y u , Y i n ,Y i i,Y m m n j are conditionally in- n dependent and Y \U ij Let fx = u (/z^, = u^ED(^,a ), (2.3) 2 . . . , / / i , ...,Aimi. - > / C n ) n i T a n d l e m We denote ( g ^ ) , g ( t f J , g ( ^ ) \ g ( ^ J ) n t #(•) denote the link function. by #(/i ). T u n Suppose further that g(^) (2.4) = Xf3 + Zu, where X and Z are two known matrices. Then E [ F | U = u] = ^ , (2.5) ij The expectation and covariance matrix of U are denoted by the following equations: E[U] = 77 and Var[U] = D ( ) , 7 (2.6) where D ( ) is a q x q covariance matrix depending on an unknown vector of variance 7 components 7 . 8 The expectations of U or transformed U are often assumed to be known. The transformed U is also often called random effects. 2.1.2 Tweedie exponential dispersion models Many exponential dispersion models have variance functions that are asymptotically of the Tweedie form, leading to a general convergence theorem with the Tweedie models as limiting distributions (J0rgensen et al. 1994). For this reason, Tweedie models occupy a central position among exponential dispersion models. Now we define Tweedie exponential dispersion models. We call (2.1) a reproductive Tweedie exponential dispersion model, denoted by Tw (/i, a ), if E(Y) = p and Var(Y) = 2 p a pP. 2 Here p = 0, 2, 3 and 1 < p < 2 correspond to well known normal, gamma, inverse Gaussian and compound Poisson distributions respectively. The case p = 1 with a = 1 corresponds to Poisson distributions. In fact, using usual notations, we 2 have Twi(/i, 1) = Poisson(^x), Tw2(//, cr )= Gamma(/i, ff ), 2 2 and Tw (/i,a ) = Inverse Gaussian(//, a ). 2 2 3 Furthermore Tweedie exponential dispersion models posses the following scale transformation property: cTw {p, a ) = Tw (cp, 2 p p c ' a ). 2 p 2 (2.7) A complete list of Tweedie exponential dispersion models is given in Table 2.1 where S, f2 and 0 denote the support, mean parameter domain and canonical parameter domain of the Tweedie exponential dispersion model, respectively (J0rgensen 9 Table 2.1: Summary of Tweedie exponential dispersion models. e Distributions P S 0 Extreme stable p<0 R R Normal p= 0 R R [Do not exist] 0<p<l — R + Poisson p=l N R + Compound Poisson K p < 2 Ro R+ R_ Gamma p= 2 R+ R + R_ Positive stable 2 <p < 3 R + R + —Ro Inverse Gaussian p=3 R + R + —Ro Positive stable p Extreme stable p = oo >3 0 Ro + R R+ R+ R R Ro R —Ro R_ Notation: —Ro = (—oo, 0] 1987a). To ease the derivation of the estimation function for regression parameter for our random effects models later, we rewrite the density of the Tweedie exponential dispersion models Tw (/i, a ) as follows: 2 p c,(2/;^ )ex {^(^-^)} if 2 P f (y; A*, P P 7 U,2, o-) = { c (2/;^ )exp{-^ (j + log^)} iip = 2, Ci(y)exp{ylog/u - p,} if P = 1, 2 2 2 where the explicit expressions for c (y;a ) are given by J0rgensen (1987). The fact 2 p that c (y; a ) do not depend on p is crucial to the derivation of our unbiased estimating 2 p function in Chapter 5. The exact expressions for c (y; a ) are immaterial in the 2 p derivation of the unbiased estimating function, thus omitted. For more details about 10 Tweedie exponential dispersion model, see J0rgensen (1987, 1996). 2.2 Prediction of random effects The prediction of random effects plays an important role in random effects models. It is especially useful in the identification of outliers. The distinction between Bayesian and non-Bayesian approaches to generalized linear models with random effects is clear. Bayesian approaches use posterior means or modes as point estimators for both parameters and random effects. The regression coefficients and the random effect variance are assumed to be random vectors and treated symmetrically with the observed responses and unobserved random effects. On the other hand, non-Bayesian approaches such as penalized quasi-likelihood approach of Breslow and Clayton (1993) and the maximum h-likelihood estimation approach of Lee and Nelder (1996) predict random effects using the mode of the conditional distribution of the random effects given the data. We concentrate on making comparisons of these two approaches with our approach. 2.2.1 Likelihoods Before we make comparisons, we need to introduce the concept of the partially observed joint (log) likelihood as follows. Let Y be the response and U be the unobserved random effects. We assume that both the conditional distribution of Y given U and the distribution of U are parametric. The partially observed joint log likelihood with only the response being 11 observed is defined as £(/3, a, 7 ; Y, U) - a; Y|U) + £ ( , U), T where £((3, a ; Y|U) is the logarithm of the conditional density of Y given U with 2 (3 and cc being the regression and dispersion parameters respectively, and £ ( 7 ; U) is that for U with parameter 7 . Thus the marginal likelihood is |exp[£(/3,a,7;Y,U)]dU. In Breslow and Clayton's paper, the partially observed joint likelihood is approximated by the penalized quasi-likelihood defined as q(J3,a, , Y , U ) = r ^n^Ud^Y^p^) - ^u D- ( )u, T 1 7 (2.8) where dY,-(-; •) denotes the deviance function defined in Section 2.1.1. 2.2.2 Penalized quasi-likelihood and h-likelihood approaches Breslow and Clayton consider the model in (2.3), but assume U follows, at least approximately, a multivariate normal distribution with mean 0. They use solutions (/3, u), for fixed 7 , from the following equations flg(y,/3,u,7) d(3 _ n dg(y»ft,u,7) _ du n as the estimators for ((3, u). That is, they obtain estimators through maximizing an approximation of the partially observed joint likelihood. They adopted the Fisher's scoring algorithm to obtain the solutions. The multivariate normality assumption is convenient to incorporate correlation structure among random effects, but quite 12 restrictive to model the distributional shape. On the other hand, Lee and Nelder consider the model in (2.3), that is u~ED(^,o ), Y \\J = 2 ij but with </(/0 = X/3 + v, where v is a monotonic transformed variable of u, denoted by v = v(u), and the distribution of random effects is assumed appropriately. They call the partially observed joint log likelihood of (Y,V) the h-likelihood, denoted by £(Y, V; (3, a, 7 ) . The solution (/3, V), for fixed 7 , from the following equations: M( ,/3,v,7) dp y 0*(y,frv,7) _ n ' av _ n u ( { Z 9 Q - , y ) is used as the estimator for (/3, V). They also adopted the Fisher's scoring algorithm to obtain the solutions. Lee and Nelder's method is feasible when the distribution of the random effect is conjugate to that of the observed response. 2.2.3 Comparison for different predictors To better compare our approach with other current approaches, we make a comparison of conditional expectation, conditional mode (the mode of conditional distribution of U given Y) and the orthodox BLUP. When both Y and U are normally distributed, the conditional mode of U given Y equals the corresponding conditional expectation which is also orthodox BLUP. 13 Hence Lee and Nelder approximated the conditional expectation by the maximum hlikelihood estimates (MHLEs), that is, the solutions of equation (2.9). In fact, these solutions are equivalent to the modes of the conditional distribution of U given Y since £({3, « , 7 ; Y , U) = ^ ( / 3 , a, 7 ; U|Y) + £ ((3, a, 7 ; Y ) . 2 Breslow and Clayton also adopted this approach in their penalized quasi-likelihood method, but they differentiated penalized quasi-likelihood which is an approximation of the partially observed joint likelihood. That is, they used an approximation of the conditional mode to be the predictor of random effects. Clearly the mode is neither necessarily a good approximation to the corresponding mean when the conditional distribution is not approximately symmetric about its mode such as in the case of normal, nor is it easy to evaluate its departure from the conditional mean. Lee and Nelder actually consider the partially observed joint likelihood of the response and a monotonic transform of original random effects. They claimed that their MHLE predictors of random effects are invariant with respect to monotone transformation of the random effects U . However the location of the mode depends on the dominating measure, therefore estimates on the original scale of random effects are preferred. Thus we suggest approximating the conditional expectation by orthodox BLUP on the original scale of random effects where appropriate. The orthodox BLUP is defined as a linear unbiased predictor of U given Y which minimizes the mean square distance between the random effects U and their predictor within the class of linear functions of Y . We call it 'orthodox' BLUP to distinguish 14 it from the modal predictor. Comparing with the conditional expectation which is the best unbiased predictor for U among all functions of Y , the orthodox BLUP is truly the best unbiased predictor for U among all linear functions of Y . We will demonstrate that the orthodox BLUP is often a good predictor for U though the conditional expectation is generally a non-linear function of Y . Unlike the conditional mode, the mean square distance between random effects and the corresponding orthodox BLUP can usually be evaluated easily. 2.3 Estimating functions Our orthodox BLUP approach is based on estimating functions. Some results on estimating functions will be used repeatedly later. We briefly summarize them here. 2.3.1 Unbiased estimating functions Suppose Y i , . . . , Y m are independent random vectors. Let us consider the estima- tion for the parameter 0 = . . . ,9 ) T q based on estimating functions of the form V>(0; ) = E £ i ^ i ( 0 ; Y i ) = £ £ i V \ ( 0 ) , where tp^ i = 1,..., m are unbiased estiY mating functions, that is E i/, (0;Y) g i = 0, where tb(0;Y) is of the same dimension as the dimension of 0. We also assume that ib(9;Y) is regular, (J0rgensen and Labouriau 1994; McLeish and Small, 1988). Let the density function of Y and the range of Y be denoted by 15 p(Y; 0) and X, respectively. We define a regular estimating function as follows: Definition 2.1 An estimating function ib(0;Y) is said to be regular if the following conditions are satisfied for all 0 in the parameter space: 1. The support ofY does not depend on 0; 2. E ih(0;Y) 0 = 0; 3. The partial derivative ^^' ^ a exists for almost every Y in X, j = 1,..., q; Y 4- The order of integration and differentiation may be interchanged as follows: AI j = m Y) (Y; = j P x J - me; Y W Y ; *)} dY, l,...,q. 5. S(0) = Eg d^^' ^ is a nonsingular q x q matrix; Y 6. V(0) = EQ[IP(0; Y)xb(9; Y ) ] is a q x q positive-definite matrix. T S(0) and V(0) are called the sensitivity and variability matrices, respectively. Note that we do not assume that ib^O; Y)s are regular. Actually the corresponding sensitivity and variability matrices for Y) are often singular in practice, due to the presence of categorical covariates. However, we assume that 1/^(0; Y) satisfies all above conditions except the non-singularity requirement of its sensitivity and variability matrices. 16 In the estimating function approach one considers estimators which can be expressed as solutions of the following estimating equation: m tf(0;Y) = 5>i(0) = O. i-l By standard asymptotic theory for estimating functions, we may show that under certain regularity conditions, the sequence of roots, 9^ \ associated with the estimating function ^ ^ ( 6 ) , is consistent for 0 and asymptotically normal. Specifically we have (Artes and J0rgensen 1998): Lemma 2.1 xAr7(0 (m) -ff)-2+N (0, J (0)) as m -> oo, (2.10) _1 where 3(6) can be expressed in terms of the sensitivity cluster i, defined by S*(0) = E { ^ r } and variability matrices for , and V;(0) = E {^(0)V\ (0)} respec- 1 T 0 0 tively: jm-^SiWJ- ^ESiWJjm-^V^)! (2-11) The asymptotic covariance matrix of 0 is therefore given by the inverse of the Godambe information matrix defined by 3(0) = S(0)V(0)- S(0) , 1 where S(0) = ET=i Si(0) and V(0) = Z? =1 T V<(0). The estimating function is a generalization of the score function U(6): 17 (2.12) where £(Y; 0) = logp(Y; 0) is the log likelihood. The sensitivity and variability matrices for the score function have following relationship: S(0) = - V ( 0 ) , which does not hold for regular estimating functions in general. The estimator and Godambe information matrix for the score function are the maximum likelihood estimator, denoted by MLE, and the Fisher information matrix, denoted by 1(0), respectively. Among all regular estimating functions, the score function is optimal in the sense that the estimator associated with the score function attains the minimum asymptotic variance among estimators associated with all regular estimating functions. To state this result precisely, let the Godambe informa- , tion matrix for any given regular estimating function ip be denoted by J^(0). Then (J0rgensen and Labouriau 1994) jj(0)-r\0), is nonnegative-definite for all 0 in the parameter space. However, in the context of generalized linear mixed models, the full score function for both the response and the random effects is not available since the random effects are unobserved. Thus we consider estimating functions other than the score function. We can also consider the optimality property, but within a more restricted class of estimating functions. In the next section, we state an optimality result within a certain linear class. 18 2.3.2 Optimal estimating functions Crowder (1986, 1987) proved, under some regularity conditions, the following Lemma 2.2 m *l> (0) = 2Si(0)V- (0)ib (0) i=i (2.13) 1 opt is an optimal estimating function Z l within the class of all linear estimating functions of the form m EQiWiW. (- ) 2 14 i=i where Qi(0) is a constant weights matrix of appropriate dimensions. Qi(0) is a matrix function of parameter 0, but does not involve More specifically, Y. The solution 0 from the estimating equation xb (0) = 0 is then asymptotically opt normal with the asymptotic mean 0 and asymptotic variance given by the inverse of m •VP) m m = £ - S ; ( 0 ) = £ V , ( 0 ) = £ Ji(0). i=l i=l i=l That is, the Godambe information for the optimal estimating equation is the sum of the Godambe informations for each i. 2.3.3 Nuisance parameter case Suppose that 0 = (0(i), 0( )) where 0(i) and and 0( ) are taken as the parameI T 2 2 ter of interest and the nuisance parameter, respectively. We partition ib (0) into , i {\bf\0), tp[ \0)), where xb^f and ipf 2 1 are of the same dimensions as those of 0( ) 1 X and 0(2), respectively. Let ib^ = Y^Li ipf" k = 1, 2. Then the asymptotic covariance 1 matrix for 0 (1) is given by (2.15) if E d e ^ { 0 ) = 0. (2) Lemma 2.3 Asymptotic Var(0 ) = S^(0)V (0)S^(0) , T {l) n 19 (2.15) where S {0) = kl EQ^O^- andV (0) (0)V> (0) | k,l = 1,2. (/) = EQ kl t Knudsen (1998) obtained this result in his unpublished thesis. Note that since Si (0) = 0, his proof is straightforward, and is based on the following equation: 2 Sn(0) 0 V 1 V (6>) n ^ S (0) S (0) ) 21 \ V (0) 22 sri WVn(0)sr (°" 21 1 1 1 (0) T V (0) 12 -T Sn(0) 0 V (0) j \ S (0) S (0) 22 21 22 * V Clearly the upper left block of the right hand side, Sr (0)Vii(0)Sn (0) , is the 1 1 T 1 asymptotic covariance of 0(i) since the right hand side is the asymptotic covariance of 0. The asymptotic covariance of 0(i) will, in general, be affected by both the variance of the nuisance parameter estimator and the variability for the estimating function for the nuisance parameter. However, under nuisance parameter insensitivity, that i> Efl \ A s U = 0, this result tells us that the asymptotic covariance matrix of 0"(2) the estimator for the parameter of interest will be affected by the variability and sensitivity of the estimating function for the parameter of interest, but not by the remaining parts of the variability and sensitivity matrices for the estimating function. 2.4 Data examples In this section, we describe some clustered data examples where over-dispersion and heterogeneity may exist. We will analyze these data sets in Chapter 8. 20 2.4.1 Epilepsy data Thall and Vail (1990) presented longitudinal data (see Table A.2) from a clinical trial of 59 epileptics who were randomized to a new anti-epileptic drug progabide (Trt=l) or a placebo (Trt=0). Baseline data available at the start of the trial included the number of epileptic seizures during the 8-week period prior to randomization and the patient's age. A multivariate response variable consisted of the counts of seizures during the 2-week periods before each of four visits to the clinic. Count data are traditionally modelled by the Poisson distribution; however the Poisson assumption that E(Y) = Var(F) is generally inconsistent with the empirical behaviour for clustered data. A set of count data is called over-dispersed or underdispersed if Var(y) > E(F) or Var(F) < E(F), respectively. Over-dispersion is very common in biological data. An over-dispersion diagnostic plot for count data relative to Poisson regression models (Lambert and Roeder, 1995) is displayed for the epilepsy data in Figure 2.1. The graph is convex if over-dispersion exists with respect to the corresponding generalized linear model. The clear convexity of the over-dispersion diagnostic graph for the epilepsy data with respect to the Poisson generalized linear model may indicate the existence of over-dispersion. Subject (random) effects are often incorporated into Poisson models to account for heterogeneity, over-dispersion relative to Poisson models, and dependence among the repeated measurements within the same subject. It is natural to introduce a second level of random effects at familial level to account for familial aggregation of epilepsy (Paik, Tsai and Ottman 1994); however we do not have such familial grouping information for this data set. A second level of random effects is often introduced to 21 O 10 20 30 40 50 60 mean parameter Figure 2.1: Overdispersion diagnostic plot for epilepsy data. 22 account for uncontrollable covariates at each visit (Lee and Nelder 1996). We will thus consider generalized linear models with two levels of random effects in next section. Such two levels of random effects arises from many practical situations (Hedeker and Gibbons 1994; Goldstein 1995). 2.4.2 Seed germination data Crowder (1978) presented data (see Table A.l) on the proportion of seeds that germinated on each of 21 plates arranged according to a 2x2 factorial layout by seed variety and type of root extract. In particular, he presented the total number of seeds on each plate and the number of seeds germinated on that plate (Table A.l). He noted that there is between-plate heterogeneity of proportions. Figure 2.2 shows that these proportions vary from 0 to 0.83 among plates. Furthermore, the variability among clustered binary responses also often exceeds what would be expected due to binomial variation alone. It is natural to account for such heterogeneity and overdispersion by means of random effects models. 2.4.3 Cake baking data Cochran and Cox (1957) presented data (see Table A.3) from an experiment in baking chocolate cakes. Three recipes were tested and each recipe was replicated 15 times, giving total of 45 batches of cake mix. Each batch was divided into six cakes and these were baked at six different temperatures. After baking, a breaking angle was measured as the response of the experiment. Firth and Harris. (1991) found the residual plots based on the analysis of variance revealed heterogeneity between batches. This can also be seen from the boxplot of the responses by batches in Figure 2.3. 23 -8 10 plate numbers Figure 2.2: Scatter plot of proportions of seeds germinated versus plates. 24 8H batches Figure 2.3: Boxplots of cake baking data by batches. 25 Chapter 3 Tweedie mixed models In this chapter, we introduce a class of models with nested random effects based on the class of Tweedie exponential dispersion model distributions. We call these models Tweedie mixed models. When the conditional responses follow a specific distribution, say, the Poisson or gamma, we call it a Poisson mixed model or gamma mixed model, respectively. For a model with Poisson responses and gamma random effects distributions, we call it Poisson-gamma model. We will also discuss the covariance structure of the model. 3.1 Models In this section, we consider three-level hierarchical models where each model is composed of m independent clusters indexed by i. Within each cluster i, there are J; correlated sub-clusters indexed by are Then within each sub-cluster there correlated observations. Let the vector of observations be denoted by Y = (Yi ,...,Y n n n i l ,...,Y m J r n l ,...,Y j m ) . T m n m J m Then Y ijk represents the kth observation in sub-cluster (i, j). One such hierarchy would be a multi-center longitudinal clinical trial involving repeated measurements within patients and patients nested within 26 study centers. Another example is teaching method evaluation involving high scool students within classes and classes nested within schools. (Goldstein 1995; Gray et al. 1995). Yijk denotes the A;th measurement for patient j taken in the ith medical center for the former example, whereas the represents the test score for student k in class j of the ith. school. Denote the cluster, sub-cluster and observation covariates by z^Zjj and z ^ , respectively. We assume that there exist cluster and sub-cluster specific random effects. The random effects for the cluster i and sub-cluster (i,j) are denoted by Ui and Uij, respectively. Thus the vector of the random effects can be written as U = (Ui, ...,U ,Un,U ) = (U*, £/**), where U* and T m mJm and (Un, ...,U j ) , stand for (U U) T u m for short. We assume further that, given the random effects, T m m the responses are independent and follow certain Tweedie distributions. Specifically: Al) Given U = u, Y m , Y n n u ,Y i n ,Y i j n ..,Y m J i,F m m J m „ m J m are con- ditionally independent, and the conditional distribution of Y ^ , given U = u, depends on Uij only which is Yi \U = u ~ £w {}i j u j,p u j ~ ) r jk = 2 p i k i l lp i UijTwp ifiijh, — ) , V (3.1) ij) U where / i ; ^ = exp(z^ ./9( )). For case p = 1, namely the Poisson distribution, p = 1. 2 A 3 Furthermore we assume that, given the cluster level random effects, the sub-cluster level random effects are independent and follow Tweedie distributions as follows: A2) Given U* = u*, Un, ...,U j m m are conditionally independent, and the condi27 tional distribution of [/»_,-, given U* = u*, depends on Ui only which is t/ij|U„, = where p tj U* ~ TWq(fJ, jU ,UJ Ui = UiTw (^lij, ^- j . 2 i i 1 ) Q (3.2) q = exp(zJ/3 ). (2) A3) Finally we assume that the cluster level random effects are independent with Ui ~ Tw (^,a ), (3.3) 2 r where ^ = exp(z /3 ). T t (1) 5fc = (7> ij> ijfe)- L e t x z z z T h e n t h e f u l 1 covariate matrix is X = ( x , . . . , x m m J m „ m J m ) The Tweedie distributions Tw , Tw and Tw are called the level 1, level 2 and level r 9 p 3 distributions, respectively. To avoid non-positive random effects, r > 2 and q > 2 are required. Now we discuss some implications of the model assumptions. First, note that (3.1) interprets Uij as weights. Note also that the dispersion component takes a special form p U\~ 2 v for the assumed level 3 Tweedie distribution in assumption Al). This special form will enable us to obtain unbiased estimating equation for the regression parameters later. It also leads to variance component structure of the covariance matrix which will be shown in the next section. This special dispersion component form also gives the linear structure of the conditional variance of the level 3 distributions. Thus we have the following double linearity: 28 Double linearity E[Y \\J] ijk = m^Uij, Var[^ |U] = ifc pp U . 2 p ijk l3 The first statement is clear from the assumption, whereas the second statement is immediately verified as follows: V a r o l i i ] = Ut- i4 U% 2 = p p P jk 2 ^ . That is, both the conditional mean and the conditional variance of the level 3 distribution is linear in the second level random effects. In fact, for any given p, the double linearity is equivalent to the assumption Al). Similarly, for any given q, A2) is equivalent to the following double linearity of the level 2 distributions. Varolii*] = wV'A E [ c y U » ] = fMjUi, We discuss some decompositions of the models. We consider the structure of the conditional mean. In the literature, the modeling of within cluster dependence is usually focused on incorporating random effects into a monotonic transform of the conditional mean assuming constant dispersion components. For multiplicative models, it is usually defined as E[Y \V,V] = ijk p UV . Hjk i ij This implies that Var[Y \\J]=uj V(Li U V ) 2 ijk ijk i ij where V(-) is the variance function for the conditional distribution of Yij , given k (U, V). In addition, Ui and V}j are often taken as independent, but the covariances 29 of the responses in general do not then exhibit a variance component structure (Morton 1987; Thall and Vail 1990; Firth and Harris 1991; Lee and Nelder 1996). In our model assumptions, the conditional expectation of Y given U is linear in the random effects. To compare our model assumptions with others, we give the following multiplicative decomposition for the conditional expectation of Y given U . Multiplicative decomposition EfY^-fclU] = fajkUij = fajkUi-jji- = PijkUiVij, where Vij = ^ji- and Cov[C/j, V~ij] = 0. The latter can be shown as follows: E , | , u . ] - E { ™ } therefore Cov[^,vy = E (ujfy = E (Uij) - - HiHij E(0i)E = 0. Additive decomposition Besides the above multiplicative decomposition, the response of this model also possesses the following additive decomposition with three uncorrelated, but dependent components as follows: Yijk = {Yijk — P-ijkUij) + {p-ijkUij 30 — Hij^ijkUi) + HijllijkUi. (3.4) This additive decomposition will facilitate residual analysis. The two sides of the equation are clearly equal. The three components on the right-hand side can be easily shown to be uncorrelated through the covariance structure described in the next section. Tweedie mixed models with one level of random effects A Tweedie mixed model with one level of random effects is a special case of the Tweedie mixed model with two levels of random effects by setting us = 0 and J; = 1 2 for all i. 3.2 Covariance structure The derivations of both our parameter estimators and random effects predictors are based on the moment structure of the model. Thus we investigate the moment structure of the model here. 3.2.1 Derivation We begin our investigation with the moments of the random effects. The intra- dependence within clusters are clearly reflected by the covariance structure described below. The derivations of the following moment expressions are straightforward using the conditioning technique. Using Kronecker notation 8^) being 1 if s = i, 0 otherwise, the covariance structure can be expressed in the following way: E[Ui] = ^ E[iy.y»] = and Cov[U , s iHjUi Ui] = 8(s, i)o- f4. and Var[Ly[/,] = 31 2 co ^. 2 Cav[U , Uij] = S(s, a i)<T /Zi/iij, 2 E[Uij] = miMj, (3.5) Cov[<7 , Uij] = S(s,i) \o ^iiijHn + 6{t,j)u fMifiij] , 2 (3-6) 2 st Cov[U ,Y ] s = ijk 5(s,i)a pJ n fj, , 2 i Cov[U , Y ] = S(s, i) [a ^pu ijk + 6(t, J)UJ miAj] 2 st ij p, 2 ijk (3.7) ijk E [ ^ | U ] = iMjkUij. Var[F |U] = p tf U , 2 ljfc jk l3 E[Y ] = mpijHijk = exp(xJ /3), ijk (3.8) fc therefore the link function is log. Cov[Y ,Y ] M = 5(s,i) ijk + \o tfiPijiiitPijkPiti 2 SitJ^pirfjPijkPiji + 5(l,k)p »ip tf ]}. (3.9) 2 l3 jk Since all the derivations of these moment structures are similar, we show, as an example, the derivations of covariance of the responses. Given U , Ym,..., • • •, Y i j n i j ,Y m J m l ,..., Y j m st( ijfc ^ Cov[y ,^- |U] sU fc ,..., Y x, nnu are conditionally independent, so we have mnmJjn Cov[r ,r |u] = o if ( s , u ) Y (M, *0, thus = S(s,i)8(t,j)S{l,k)Vai(Y \U) = 8{s,i)S(t,j)5(l,k)p f/ U . t3k 2 ijk i3 Therefore Cav[Y ,Y ] M ijk = E { C o v [ F « , y | U ] } + Cov[E(r «|U),E(F |U)] s = ufe s S(s,i)S(tJ)6(l,k)E\yax{Y \U)] ijk 32 iifc i3 +Cov[p U , fJ.ijkUij] stl st 5(s,i)5(t,j)5(l,k)p ^ E(U ) = 2 ijk ij +tiaiilHjkCov[U t,Uij]. (3.10) a The proof is completed by plugging (3.5) and (3.6) into the last equation. Now we return to (3.4) and verify that the three components on the right-hand side are uncorrelated. The verifications are very similar so we only show that the first two components are uncorrelated. It follows from (3.6) and (3.7) that Cov[Y , Uij] = fjLij Var(Uij); ijk k hence Cov[Y ijk - fj,ij Uij, Pij Uij] k k = = Hij Cov[Y , U^} - fj, Cov[Uij, Uij] 2 k ijk ijk 4-*Var(^-) - p i V&r(Uij) = 0. 2 jk Similarly we have Cov[Yij —fajkUij,Hijfj-ijkUi] = 0. k So the result follows immediately. Finally, note that Corr[W ,]S clearly depends on the mean parameters. 3.2.2 Matrix expressions To facilitate the derivations of the quantities of interest in the rest of the thesis, it is desirable to express the moment structure in matrix form. After introducing 33 some matrix notation, we will state the covariance matrices between random effects and the responses in matrix form. The matrix form will be useful in the derivation of the orthodox BLUP predictors of the random effects. We will then concentrate on the derivation of the inverse of the covariance matrix of the response since this inverse plays a key role in all stages of the derivation of the orthodox BLUP approach. Let us first introduce some matrix notation. Let Y$ denote the responses corresponding to the ith cluster, that is, (Y ,..., Yn ,..., in Mi* = (Mil) • • • > t^ijjy (^it^ijIHjli Y ,..., Y nil Uil Vij* = (Mtjl> • • • J ^ijmjY") Mi** = (Mil*> • • • • i fJ'ifJ'ijH'ijriij) ) . Let T i J i n i J i • • ' MZ/i*) d a n ij* U • Then E(Y ) = ^ i = i K i*,---,^i*) T T In addition, for any a = (ai,..., a ), a is defined as (a[,..., a^). Then we have r n Cov[U ,Y ] i = i Covp7y, Yi] = ^tt-^uj (3.11) a- i4- vJ, 2 1 + ^ ( O , . . . , i/?,, . . . , 0 ). 7 T (3.12) Now we derive an explicit expression for the inverse of Var(Y). Note first that I Var-^Yi) 0 \ Var (Y) _1 V 0 Var-^YJ j since different clusters are independent. Thus it is enough to derive the inverse of Var(Yj). The inversion of Var(Yj) can be further simplified by noting that Var(Y;) is a patterned matrix which reflects the hierarchy. To find a simple expression for 34 = Var(Yi), we define / •A-ij v P P'iP'ij Then Ui 0 X Var(Y ) = i 0 Therefore Var(Yj) can be expressed as a positive-definite matrix plus a matrix product between a column vector and a row vector. Thus we can invert Var(Yj) using the following well-known formula (cf. Rao 1973) A- ab A~ 1 + b Aa ' l (A + a b ^ - = A " 1 1 T l T (3.13) where a and b are vectors and A is invertible. With ( Aii V a = <T //£ 2 Ui and b = i>i, V a r f Y , ) - - (A, + 0 \ 0 we obtain OV-W)" = A."' - where 35 f f ^ f f i ^ . (3.14) A AT - 0 1 1 0 Now where each block A^- is again of the form (3.13). Applying the same matrix inversion formula to the A^- gives < ± A" where ^ 0 \ 1 (fillip = l/(p 2 + to ^ 2 1 £*=i Mi?/)- 2 _ M ' (3.15) Thus we have ^ Ajj ^jj* [ A ^ I J * Mij Z^fc=l p2 H-ijk ij W 1-p Mij* 1-P — ijH>ij* , w Plugging the last two equalities into (3.14) gives the explicit formula for the inverse of the covariance matrix of the responses. We will write out this explicit expression with simpler notation later. This explicit expression will be very useful in our discussion of random effects prediction and parameter estimation in the following chapters. 36 Chapter 4 Orthodox B L U P predictors of random effects In this chapter, we study the orthodox BLUP of random effects. Robinson (1991) presented a variety of applications where the prediction of random effects is of interest. These applications arise from areas such as sports, genetic study, quality and management, time series, image analysis, geostatistics, actuarial science and smallarea estimation. The prediction of random effects is also found to be very useful in the identification of outliers. Fellner (1986) discussed this topic with some examples. He found that looking at the predictors of random effects is often more sensitive for detecting outliers than merely looking at residuals. We derive the explicit expressions for the orthodox predictors of random effects. Let us begin with the general expressions for orthodox BLUP. 37 4.1 Orthodox BLUP Let U and Y be random vectors with finite second moments. We define the orthodox BLUP of U given Y by: U = E(U) + Cov(U, Y ) V a r ( Y ) (Y - E(Y)). (4.1) _1 The mean squared distance between U and U can be evaluated through the following equation (Harvey 1981; J0rgensen et al. 1996b): Var(U-U) = E[(U-U)(U-U) ] T = Var(U)-Cov(U,Y)Var- (Y)Cov(Y,U). (4.2) 1 This mean squared distance can be further decomposed into the sum of the mean squared distance between orthodox BLUP of U given Y and the conditional expectation of U given Y and mean squared distance between the conditional expectation of U given Y and U as follows: E[(U-U)(U-U) ] T = E[(E(U-E(U|Y))(U-E(U|Y)) ] T +E[(E(U|Y) - U)(E(U|Y) - U) ]. T (4.3) In addition, we have the following two desirable orthogonality properties concerning the orthodox BLUP: Cov[U - U , U] = 0 and Cov[U - U , Y] = 0. (4.4) That is, the residuals between the random effects and their orthodox BLUP predictors are orthogonal to both the response and the predictors. These orthogonality properties will be repeatedly used later. 38 The first two moments of the orthodox BLUP are as follows: E(U) = E(U), Var(U) = Cov(U, Y)Var (Y)Cov(Y,U). 1 (4.5) The first equation follows immediately from (4.1) since E[Y — E(Y)] = 0. The second statement can be easily verified by using (4.4) as follows: Var(U) = Var(U - U + U) = Var(U - U) + Var(U) + 2Cov(U - U , U) = Var(U - U) + Var(U) + 0. (4.6) The verification is completed by noting (4.2). Comparing (4.2) with (4.5), we have Var(U) < Var(U). Thus the variance of the orthodox BLUP predictor of random effects is generally smaller than that of the random effects. Thus the orthodox BLUP predictor is also referred to as the shrinkage predictor of random effects. 4.2 R a n d o m effects p r e d i c t o r s Explicit expressions for the orthodox BLUP predictors of random effects U given Y can be derived from (4.1). Since different clusters are independent, we derive the random effects predictors from the following two formulae: Ui = E(Ui) + Cov(U Y )Var(Y )" {Y - E(Y^), (4.7)' Un = E(Uij) + Cov(% Y )Var(Y )" (Y* - E(Y,)). (4.8) 1 u i i t and 1 i 39 i Let us first derive the explicit expression for the orthodox BLUP predictor Ui. It is more convenient to derive the expression using a matrix form based on (3.11), (3.12) and (3.14) as follows: u- - u I . V - W A - I 1+ ^r 2 A ^(Ar^) ) T ( Y v ) (4.9) a^- ujA^Ui 2 Noting that we have Ji riij j=l fc=l Hence where u>ij = l/(p + co /4j Y^=i Pnjk) 2 2 1 a s defined in (3.15). Thus the orthodox BLUP predictor for each cluster i is a linear function of responses within the cluster. A reasonable predictor of Ui should be nonnegative since Ui is. The orthodox BLUP predictor Ui is clearly nonnegative when the responses are. In fact, we mainly concentrate on models with nonnegative responses in this 40 thesis. Similarly we can derive an explicit expression for l \ . We introduce a notation which denotes ( 0 , . . . , T 0 ) . We have t/J^,..., T Cov[tf , y Yt] T = o 2 ^ ^ 1 (4.12) + uj ^ J . 2 e 3 Thus Uij = = (a pr ^vl HjUi „ +^ - ^ V a r ^ ) " x 2 + ^ - n) 1 rftf^elVaxOTiy^Yi-Vi) f F x „ 2 „ H a p\- Al Ui(Ai Ui) \ 2 T | i - i 0 2 l l T / v l + a^-^A- ^)^ 1 = mjUi + ^Vr ( 1 A r %) (Y 1 T Plugging A ^ e y = ( 0 , . . . , T Uij = i - -UilJi). 0 ) T Wij(fj,l~ ) ,..., p /Jy/7i + T into (4.13), we obtain T Pi'kiYijk uj nlJ Wij 2 (4.13) l - HijHijhUi) fc=i = + u fj, [ WijJ^i2 [ Y j , p w jfMi Ui 2 2 i j q 1 1 i (4.14) p i k i k fc=i where the second expression is obtained from the identity Uij 2 pwij i = l - u 2 o-l filj Wij 2^ fc=l 41 2-p The first expression for Uij shows that U^ is PijUi adjusted by a linear function of the responses within sub-cluster (i,j). The second expression shows that Uij is also nonnegative when the responses are. 4.3 Mean squared distance To evaluate the distance between random effects predictors and the corresponding random effects, we study the diagonal elements of Var(U — U). These mean squared distances form the basis of our discussion of consistency of random effects predictors in next section. Using the notation c(i) = E(Ui — Ui) and c(ij) = E(f/y — U^) , it follows from 2 2 (4.6) that Vai(Ui) = Var(Ui-Ui + Ui) = Vai(Ui-Ui) = c(i) + V a r © ) + 0. + Vai(Ui) + 2Cov[Ui-Ui,Ui] Hence c(i) = Viur(Ui)-Vax(Ui) = a ^ - V a r (ft). 2 (4.15) Similarly, we have c(ij) = Var(iy - Var(l\) = <r t*il*ij + u Pi(J, ij 2 2 42 9 - Var(<7y). (4.16) We derive explicit expressions for c(i) and c(ij) via Var(U). The explicit expression for Var(CTj) can be derived from (4.5) as follows: Var(^) = Cov(C/ Y )Var(Y )- Cov(Y ,f/ ) 1 i) i i 2 r—1 T I i - l l -"-t u , 1+ i "i) 1+ J2..t-1,. a^ujAr^i a^-'ujA- ^ 1 Thus l + aMiV r V A " ^ 1 = ^ 1 Hence we have the following proposition: Proposition 4.1 1+ <7 /i£ 2 Efcil 1 iil ijl 'i k w x i 3 Now we derive the explicit expression of c(ij) based on (4.16). As the derivation of a concise expression for using the matrix form in (4.5) becomes much more V&r(Uij) complicated, instead we derive V&v(Uij) directly from (4.4). To ease the derivation, we first rewrite Uj as c(i) Ui = Mi EE J i n u iiViilhrk iik 1=1 k=l w Y Also note that 43 + constant. a 2 Ji r-1 H n 2-p ijVijVijk w j=ifc=i _u 2 . , r A*i -> = 7777 - 1c w Furthermore, it follows from (4.4) that Cov(Uij, Uij) = Cov(C/y, C/y). Now we have Var(l\) = CoviU&Uij) = Cov(p WijHijUi + u \i\~ Wij 2 2 VijkYijk, Uij) x fe=i = (PwijiUjCoviUi, Uij) + <J^ Wij £ /4/Cov(y , C/y). X (4.17) ijfc fc=i But ( = c Ji (j\ nu zZ zZ iiVii^i\k iik w ^ + const, C7y Y i=l fc=l — ££^^/^fc Cov(y^,c/y) W l=ik=l p r^l. (=1 fc=l C R 2 A i" i 1 zZ £ wulMi&b + ^lAj^Wij l-\ k-l Ihjk k=l = c(z)^ |^--l + l-pV | j J = cP-^Hij - p c(i)nijWij, (4.18) 2 and (CT ^A4 + ^MijWtfj u Aj ^ij zZ ^ C o v ( y ; , i y = 2 p fc 2-P a zZ Hk 2 jA: lw k=l k=l = { a 2 ^ l + u ^ j){l-p w j). 2 q (4.19) 2 i i Plugging (4.18) and (4.19) into (4.17), after some simplification we have Var(ftj) = CT A</4 2 + uPmnlj - p Wij Therefore we have 44 2 (w #j/4j + p c(z)wy/Li ) . 2 2 2 Proposition 4.2 c(ij) = o w + ffctywijfij) 2 ij . This expression does not involve a and pij explicitly, but implicitly via c(i) and 2 k w^. In the expressions for both c(i) and c(ij), the p,ij s appear only in the form of k Z)fc=i A i / ) i 7 P c that is the average of p ~ 2 within sub-clusters. p k 4.4 Consistency In this section, we draw consistency results about random effects predictors based on the mean squared distances derived in last section. According to Chebyshev's .—. p p inequality, we have c(i) —» 0 implies Ui —> Ui, and c(ij) —¥ 0 implies Uij —> Uij. To draw our consistency results based on this inequality, we need to find upper bounds for c(i) and c(ij). Since oversimplified upper bounds do not result in desirable consis- tency results, we derive the upper bounds through a sequence of inequalities. After some algebra, we can show that c(i) is bounded above as follows: Lemma 4.1 fr /4 2 c{i)< p + u minj(nij) 2 (f?+ maxj(p,jj 2 u min {n )max (tf7 )min {p]j )) 2 1 j )minj, (p ) k ij jjk + a minj(riij)Ji 2 ijk p j mui k ^ j(l ij) ''' j,k(f 'ij ). m n I m We start our investigation with w^. 1 2 Efc=l Mijfc 2 1 ^ ,2-PV p + w min -(ri 0min /(/j?7 )min /b(^7) . 29 , <- . . 29 _ : _ J 1 iJ 45 j i J k Proof P +U p1j n ( 4 - 2 0 ) Also we have nij 1 .,2-p sr^ij E ijf ijk w 2-p _ Z^fc=l P fc=l -I- W /Xy H'jjk 2^k=l Pijk 1 E fc„ij P=,,2-p l ^yfc \-L0 ufl~ 2 2 n l ^ PlJ I > min (n, )£in J J j f c ( ^) + a ; 2 m a X ^^ 1 ) mi"j(^)minjfc(^ fc ) p + u mm (n )ma,x ([i 7 )mm k(^ r 2 2 q j ij j P ' 1 i j hence 2 r-iv*^ 2- P ^ gV» ^min^^ )mm (n )mm {p ) p j P +w min (n )max (^ i=ifc=i 2 J li i ij jk ijk )mm (n ) p jk ijk Now we have c(i) = o ti ( P + ^ m%(^j) ax (pf~ )min ( 2 < 2 2 m 1 J j;fc p + w min (n )max (pf7 )min - (p 7 ) + ^min^n^) J^^mm^/Jminj,k{lAjk 2 2 1 j ij i 2 J fc J p fe • The derivation of an upper bound for c(ij) is straightforward. Note that Proposition 4.1 immediately implies that c(i) < a t f • 2 In addition, we have i .. p w ?! = 13 o + uPulj E ^ i 2 1 46 < i lAik ) Thus c(ij) = pw ^ 2 {uPliifij + 2 ij I 2 Q , p c{i)w p ) 2 2 ij i] 2 r 2 \ It follows from (4.20) that c(ij) is bounded above as follows: Lemma 4.2 c(ij) < The following 'small dispersion asymptotics' (J0rgensen 1987b) is an immediate consequence of Lemma 4.1 and 4.2: Proposition 4.3 1. Ui as a -> 0 or (UJ , p ) -> (0,0); 2. ftj - A 2 2 2 as p -> 0 or (a , a; ) -)• (0,0). 2 2 2 For large sample asymptotics, clearly we have Proposition 4.4 p Ui as Ji a// Pi, Uij and p,ij k —>• oo and C/j,- p C/jj as mirijinij) — > • oo, are contained in a compact set not containing 47 zero. This result distinguishes the consistency conditions for Ui from those for U^. The former is implied by a large number of sub-clusters within a certain cluster, while the latter is implied by large sizes of the sub-clusters. This result matches our intuition. The above results can be expressed in a more delicate form derived from Chebyshev's inequality as follows: Proposition 4.5 under the same conditions. It follows from (4.3) that c(i) ->• 0 and c(ij) -> 0 also imply Ui and Uij - A E{Uij\Y), respectively. Thus Ui -A E(Ui\Y) and U {j - A -A E(Ui\Y) E(J7„|y) are clearly implied by the conditions stated above. Similarly, we have Proposition 4.6 A more delicate discussion of the consistency of the random effects predictors would involve the relationships among the quantities further results will be discussed in Section 7.1.1. 48 fa, fj,ijk,p,q,r,riij. Some Chapter 5 Parameter estimation 5.1 Estimation of regression parameters We begin our discussion on the estimation for the regression parameters w i t h the case of known dispersion parameters. T h e inclusion of unknown dispersion parameters w i l l be discussed i n next section. 5.1.1 Estimated score function A s the score function is optimal among a l l regular estimating functions, we begin our investigation w i t h the partially observed score function m Ji / X " — £ £ ij~~2~(Uij 30(2)!) i=l j=l 1 9 Z m 3/3(3) £ Ji defined below: ritj — UiPij), 1-p £ £ ijk Z i=l j=l k=l 49 K O^ijk P' UijPijk)- The marginal score function for the response Y is obtained by taking expectation of the above partially observed score function with respect to the conditional distribution of the unobserved random effects U given the responses Y. Since the partially observed score function is linear in the unobserved random effects U, the marginal score function for the response Y is obtained by replacing unobserved random effects U in the partially observed score function by the expectations of the unobserved random effects U given the responses Y, denoted by E(U|Y). To find the maximum likelihood estimate, in principle, we can then solve the marginal score equation obtained by setting the marginal score function to zero.. However, a closed form expression for E(U|Y) is generally difficult or impossible to obtain. Note that E(U|Y) is the best unbiased predictor of unobserved random effects U given the response Y and we showed that the orthodox BLUP predictor, the best linear unbiased predictor of the unobserved random effects U, generally converges to E(U|Y) in probability. Thus we approximate E(U|Y) by the orthodox BLUP predictor of the unobserved random effects in the marginal score function as follows: m m (5.1) IkpiP) </> (/3) = £ £ z (2) (u {(3)-um^M) ; l3 i=i j=i m Ji = iZ^\P) (5.2) ™ij </> (/3) = £ £ £ (Y (3) Z i = l jz=l k = l y f P - UyfflfMikiPJ) = z24 \P)Z ljk c i=i 2 50 (5-3) The three functions in (5.1),(5.2) and (5.3) are clearly linear functions of the responses as the orthodox BLUP predictors of the random effects are. We call ib((3) as defined below the estimated score function (based on the orthodox BLUP): = (V> W,V> W,V> W) ( ( (3 m m m i=l i=l i=l T m i=i m = £V>,(/3), i=l where ipi{0) corresponds to the estimated score function for the zth independent cluster. Clearly t/>j(/3) is an unbiased estimating function as the orthodox BLUP predictors of the random effects are unbiased. If we treat the unobserved random effects as 'missing data', the EM algorithm can then be used to obtain the maximum likelihood estimate. This algorithm iterates between an E-step, which involves evaluating the conditional expectations of unobserved random effects given the response in the marginal score function using current parameter values, and an M-step, which involves obtaining updated parameter estimates by solving the marginal score equation. The orthodox BLUP approach is equivalent to approximating the E-step by replacing the conditional expectations by orthodox BLUP predictors and approximating the M-step by solving m E^(/3)=0. 1=1 The roots of this equation are then used as regression parameter estimates. 51 5.1.2 Standard errors of regression parameter estimators Letting (3^ denote the sequence of roots oiY%L\ rfiifi) = 0, it follows from Lemma 2.1 that $ ^ is consistent for (3 and asymptotically normal as m —> oo. Specifically, the asymptotic variance is given by the inverse of the Godambe information matrix: J(/3) = S(/3)V(/3)- S(/3) , 1 (5.4) T where the sensitivity matrix S(j3) and the variability matrix V(/3) are given by: m m V(/3) = £ V M = zZ /3 {*M1>J(P)} E 1=1 • 1=1 From now on, we simply denote ft ^ by /3. An analogue of Wald's test is available for testing the hypothesis H : (3^ = 0, 0 where (3^ is a sub-vector of (3. The test statistic is: ^ = 3 {j (3)}~ 3 , r 11 ( 1) 1 (1) where J (0) is the corresponding block of the asymptotic covariance matrix of /3. n Asymptotically, this statistic follows a x (A;)-distribution, where k is the size of the 2 sub-vector (3^. 5.1.3 Newton scoring algorithm The orthodox BLUP approach is actually a linearized E M algorithm. Instead of using this approximate E M algorithm, we adopt a more efficient algorithm, the Newton 52 introduced by J0rgensen et al. (1996a) to solve the estimating scoring algorithm, equation ib((3) = 0. The Newton scoring algorithm is defined as the Newton algorithm applied to the equation il)((3) = 0, but with the derivative of t/?(/3) replaced by its expectation S((3). The resulting algorithm gives the following updated value for (3, (3* =(3-S- {(3)ib{P). (5.5) 1 Clearly the sensitivity and variability matrices are crucial in this estimation procedure. USij denotes and denotes E(V> (/3)V> (/3) ), then the sensitivity (i) 0) T and variability matrices can be expressed as follows: Sn s 12 S13 S21 S22 S23 \ S31 S32 S33 v 1 2 v 1 3 V i V 2 2 V 2 3 V31 v 3 2 V 3 3 i m = V n 2 \ / Note that S(/3) is apparently of asymmetric form. For example, S i = 3 whereas S31 = y J^>'• ^r a , Actually S{(3) can be shown to be symmetric. However a block-by-block direct derivation of these matrix blocks from (5.1), (5.2) and (5.3) would be complicated by the nonlinearity of f3 in U (J0rgensen et al. 1996b, 1996c). Instead, we introduce a concise matrix expression for ib(f3) in the next section. This 53 expression not only facilitates derivation of the sensitivity and variability matrices, but also helps to clarify the relationships between the sensitivity and variability matrices. 5.1.4 O pt imality After introducing a matrix expression for we will study the relationships be- tween the sensitivity and variability matrices and the optimality of the estimated score function. The concept of optimality here is in the sense that the estimated score function attains the minimum asymptotic covariance for /3 among a certain linear class discussed below. Recall that the full covariate matrix is X with x7 = (z , zJ- zj ). Let Xj denote T fc : fc the sub-matrix of X corresponding to the ith cluster. We state the following results: Theorem 5.1 For Tweedie mixed models, the estimated score function, sensitivity and variability matrices can be expressed as follows: 1. = X. diag{E{Y)) r 2. V(/3) = Var{xb{f3)) Var~ (Y) l (Y - = X diag(E(Y)) E{Y)), VarCY)' T 1 diag{E(Y)) X, 3. J ( / 3 ) = - S ( / 3 ) = V ( / 3 ) . The first statement gives a global matrix expression for the estimated score func- tion. This expression can be viewed as a multivariate version of the quasi-score function first proposed by Wedderburn (1974) . The second statement gives the expression for the variability matrix. The third statement for the sensitivity matrix is 54 similar to that for the Fisher information matrix. Both the derivation and computing efforts are greatly eased by the relationship given in the second and the third statements. Clearly now the asymptotic covariance matrix of (3 is given by the inverse of the variability matrix alone. The key statement in Theorem 5.1 is the first statement; the second and third statements are immediate consequences. We will leave the derivation of ib((3) to a separate section since it is long and technical. Now we show the second statement. V((3) = Var(V(/3)) = Var(X diag(E(Y))Var- (Y)(Y-E(Y))) = X diag(E(Y))Var- (Y)Var(Y)Var- (Y)diag(E(Y))X = X diag(E(Y))Var (Y)diag(E(Y))X. T 1 T 1 T 1 _1 The third statement can also be easily shown based on thefirststatement. Noting that <9E(Y) d/3 T = diag(E(Y))X and E ^ (Y - E(Y)) = 0, we have E/3 j f d (x diag (E(Y)) VarfY)- ) T 1 | (x diag (E(Y)) Var(Y)- ) T + E / 3 1 ((Y - E(Y)) ® E) 9 ( Y ~ ^ ( Y ) ) } f3(x diag(E(Y))Var(Y)- ) . . ^ (E (Y - E(Y)) ® E) T 1 0 + (x diag(E(Y))Var(Y)- ) T 1 55 (-^^) = 0-X diag(E(Y))Var(Y)- diag(E(Y))X T 1 = -V(/3), ' where <S> denotes the Kronecker product and E is an identity matrix whose order equals the size of /3. Now we can derive the explicit expression for S(/3). Note that S((3) = -V09) = -X diag(E(Y))Var(Y)- diag(E(Y))X = - ( X , . . . , Xjjdiag (E(Y)) Var(Y)- diag (E(Y)) = - £ X d i a g ( E ( Y ) ) V a r ( Y ) - d i a g ( E ( Yi)) Xj. T 1 T 1 T (5.6) 1 I i i i=i Rewrite Var (Y,) as X An / o \ c(i) Var-^YO Vo 1-p Aij, l-p\T\ I / 1 where A y is given in (3.15). Plugging this expression for Var {Y ) 1 into (5.6), we 1 i have j mj zZ S) (zZ zZ ijtojVij~k ijk) (zZ ZZ m (P) S = Ji nij c( t w j=i i=l m ^2^2 + zZ i=l x fc=i Ji nj P {j V m^Jk)CzZ 2 fc=i j=l 56 T n { zZ VijtfPizZ ijto3t£jkXijk) w j-i fc=i fc=l »ij~kXijk) T - £ £ £ -^rffkXijkxJjk- (5-7) i=ij=ifc=iP In next section, we show that tl>M = XTdiag (E(Yj)) VarCY,)- (Y, - E(Y;)). 1 This immediately implies that the results of Theorem 5.1 hold for each cluster. Applying Lemma 2.2 to the estimated score functions for each cluster, we may directly show the optimality of the estimated score function among the class of all linear estimating functions of the form Y%Li Qi(/3)V>j(/3) if the variability matrices Vj(/3)s are nonsingular; however, these variability matrices are often singular in practice. On the other hand, J2iLi Qi(/3)V'i(/3) is a subclass of the class of all linear estimating functions of the form YliLi^iiP) (Yj — E(Y;)) since is a linear function of Yj — E(Yj). Applying Lemma 2.2 to the latter linear class, we have m </>(/3) = = E i i=i X r d i a E Y 1 l -t p{ E g( ( 0)Var(Y )- (Y -E(F )) d { Y l e f Y l ) ) } l V a r ( Y - i ) - 1 ( Y » " E ( Y ' ) } is the optimal estimating function among the linear class. 5.1.5 Derivation of matrix form of ip((3) We will derive the matrix expression for ib((3) in terms of ipi((3). Furthermore we will deal with t/>^(/3), %bf\(3) and ij)f\p) separately. The basic technique of handling V>|^(/3) and ibf\(3) is to convert (5.1) and (5.2) into explicit linear functions of Yj - E(Yj) directly. We will then deal with il>f\(3) by converting (5.3) into explicit linear function of Yj — E(Y*) indirectly due to the technical convenience. 57 The derivation of ib\ \(3) is quite straightforward. First rewrite Ui as l Ui = tM + a V r E ( Y ) V a r ( Y ) " (Y* - E(Y,)). 1 T 1 i i Thus we have # (/3) *£l'Ui-tii) = } = z ^ ( a V r E ( Y ) V a r ( Y ) - (Y, - E(Y,)) = z E(Y ) Var(Y )" (Y, - E(Y )) = (Zi, )diag 1 T i 1 l l T i 1 i l S (E(Yj)) V a r ^ ) " (Y, 1 Zi E(Y0) • Similarly, we may deal with ibf^ (/3) by noting that iHjUi + o ; V r = l e 5 V a r ( i) ( i Y _ 1 Y "0- Thus we obtain Uij - KjUi = ^nf-'eljVMCY^CYi - E(Y0). • Plugging this formula in ib^\/3), we have ib?\(3) = j=i = E v z 3=1 = ^Zij^&j-U^j) ^ ( ^ ^ ^ V a r f Y O - ^ Y , - E(Y,))) U X;^-e5Var(Y0- (Y -E(Y0) 1 i = (E z ^ V a r f Y O - ^ Y i - E(Y,)) J'=I = = . . . ,z i/ lJi T l )Var(Y )- (Y - E(Y,)) 1 Ji (za,..., z a , . . . , z ,..., Ui i i z )diag (E(Y;)) Var(Y )" (Y, - E(Y;)). 1 i7i 58 i A direct derivation of ib\ ((3) would be much more complicated than those of tp^\(3) and xp[ \(3). Instead, we introduce a new matrix notation to ease the deriva2 tion. Let Uj*A*i** — (UnPm,..., UiiHn ,..., nn Uu^ij^,..., Uu^u^j.) This new vector is the same size as the longer, vector p with components iJfif UijHijhS- That is, each component pij of p^ . is matched by a random effect com- k ponent Uij within the same sub-cluster. We may call it a nested product between = (Un,..., U ) and p^ = (p ,Pij .)- Ui in This nested product notation is inu especially convenient in the generalization of Tweedie mixed models with two levels of nested random effects models to higher levels. Similarly we have Ui*Pi** = (Uiipni,..., C/jiMiinjD . . . , Uij^ij^,..., UuifJ-iJimj.) '• 1 With this notation, we can express ibf\(3) simply in terms of Yj—-XJ^p^. Before doing so, we note: Ui^i** = E(Yj) + C o v ( U ^ „ , Yi)Var(Yi)- (Y< - E(Y<)), 1 Now we can rewrite the covariance of the responses in terms of that between Ut*/A I]M , and Yi as follows: P Pin Var(Yi) = H1 U Cov(U A* **,Y ) + i# < i P ViJ n . iJ n j V i 59 iJ l i i J / ' = P— 1 2 M i l l P Cov(U ,/* „,Y ) + i i diag(E(YO) i P^PHJXJ. J Hence Yj — U i * / / ^ has the following simple expression: Yj - U i , / i = Y i - E(Yi) - C o v ( U ^ „ Yi)Var(Yi)- ( Y - E(Y,)) 1 i M = {Var(Y ) - C o v ^ ^ , YO} Var(Y )- (Y< - E(Yj)) 1 i / 2 l P—1 ' P Mill diag (E(Yi)) Var(Y )" (Yi - E(Y,)) 1 l p ^ „ r - i ; Now we have W (0) = — — EE ufc " J " - ^ijMijfc) i=ifc=i P i-p . / - P /'Mill ^iJinjjj s / I ill 5 -•• > • i^niJ,A *» z ( Z i l l , , . . Z 2 P = :pjr v Z 2 • , Zi p J i n i J { v — ' )diag (E(Yj)) Var(Yi)" (Y, - E(Yj)). 1 These equations for ipf\l3), ibf\p) and ibf\/3) give ^! (/3) 2) V ^ (/3) 3) Zi diag(E(Yi)) VartY,)" (Yi - E(Y,)) 1 Zil ZiJi Z i l l ZiJvnv 60 = (x ,..., x = x T d i a g (E(YO) VarCY,)" (Y, - E(Y,)). m )diag (E(YO) VarCY,)" (F, - E(Y,)) 1 UiniJi (5.8) 1 Now Theorem 5.1 follows immediately. Noting that Var (Y) is block diagonal, _1 we have m </>(/3) = £</>i(/3) t=i m = iE =i z X T d i a g ( ( 0 ) Var(Yj)- (Yi - E( F/)) E Y 1 / Y!-E(Y!) ^ (X7,...,X^)diag(E(Y))Var(Y) - l ^ Y = 5.2 m - E(Y ) ) m X'diag(E(Y))Var- (Y)(Y-E(Y)). 1 Estimation of dispersion parameters We now discuss the situation when the dispersion parameters are unknown. We estimate unknown dispersion parameters using the adjusted Pearson estimator (J0rgensen et al. 1996b); that is, the Pearson estimator adjusted by bias correction. 5.2.1 Adjusted Pearson estimators Recall (4.16), or equivalently c(i) a Mi , + E(Ui - Hif (5.9) M i We may thus estimate a by the following adjusted Pearson estimator: 2 2 _ 1 " (Uj ™ i=l //j) ,,r Mi 61 2 1 " c(i) ™ 2-~i T=i m ,,r Mi • (5.10) To obtain an unbiased estimator for to , we consider 2 E(Uij - mftif Var(r7 ) + tfjVaxiUi) = ij -2piijCov(Ui, Uij) = VaxiUij) + tfjVaxiUi) -2fXijCov(U Uij), (5.11) h where the last equality follows from (4.4). Rewrite (4.15) and (4.16) as Var(£/i) = o p\ - c(i), (5.12) 2 and VaxiUij) = o i ix\j + 2 ~ c(ij). r i i (5.13) Plugging (5.12), (5.13) and (4.18) into (5.11), after some simplification we obtain w //i/4 = E(Uij - HijUi) + c{i)ii j + c(ij) - 2p c(i)w iJ j. 2 2 2 2 (5.14) 2 iJ Thus we may estimate u as follows: 2 Lb 2 1 = m 1 (Uij l^ijUi) i=l *^» j=zl VipHj 1 A 1 ^c(i)tfj + | c(ij)-2p c(i) jtfj 2 Wi The dispersion parameter p can be estimated similarly. Noting that (4.4) implies 2 CovCYijk, U^) = Cov(Y , U^), we have ijk E(Yij - UijUijk) = 2 k = VsLi(Yi ) + p j V&v(Uij) - 2pij Cov(Y , l \ ) 2 jk k k ijk Va,v(Y ) + p, j Va,v(Uij) - 2p j Cov(Y j , U^) 2 ijk k 62 Ji k i k (Vax(Yi ) - PijkCov(Y , Uij)) jk ijk +H% (Var(fty) - Var(E^)) k P PiPijP 2 (5.16) + c(ij)n . P 2 ijk jk Therefore p can be estimated as follows: 2 1 ^2 1 1 i=i m m i=i (^ijfc UijP*ij } k j = i i j fe=l PiPijP%k n •'t j=i n MiMij fc=l » j . (5 , 7) Similar to the REML estimator, we may also consider a small sample degree of freedom correction. For example, we may replace m by m minus the number of regression parameters estimated; however the latter correction is suggested on an intuitive basis. 5.2.2 Asymptotic properties Let £ denote (o , u , p ) . Let 2 2 2 T Pi 0( ) (0 |) 2 1 = jr f ^ ~^'^^ Pi i + C ^ ~ P 2 2 c ^ ^l W i | and ^(3)^ 0 = £ ~ E l j=l ij fc=l I n ^^ PiPijPijk Uij ijk + ^fi^Jk PiPij \ - p I 2 Let m W 0 -i = E i = i -p ( ^ O , 0 , <t>? ( f t # ( / 3 > 0 ) ] m 63 } m -I 0- = E i = i m If we replace (a ,u ,p ) in (5.10), (5.15) and (5.17) by (a , cu , p ) , clearly these 2 2 2 2 2 2 T three equations can be rewritten as an unbiased estimating equation m -I i=l m where </>;(/3, £) is unbiased estimating function. Together with m V(/3,0 = £ ^ ( / 3 ) = 0, i=l we obtain a set of unbiased estimating equations. Applying Lemma 2.1 again, we obtain that the sequence of roots, ( / 3 , £ ) , is consistent for ( / 3 , £ ) and asympT T T T T T totically normal with mean (/3 ,,£ ) , as m —t oo. T T T Note that f 3 (x diag(E(Y)) Var-^Y) (Y - E(Y))) | T = = a(x^ia (E(Y))Var^(Y ) = 0. g ) E ( Y ^ Taking /3 and £ as parameter of interest and nuisance parameter respectively and applying Lemma 2.15 to (tb((3, £ ) , <p((3, £ ) ) = 0, we conclude that the asymptotic T T T variance for 0 is still given by S ^ O ^ V ^ S " ^ ) = - S " ^ ) . 1 1 In the literature of generalized linear mixed models, the asymptotic variance of /3 for the unknown dispersion parameters case is generally taken as that obtained for the known dispersion parameters case, but with unknown dispersion parameters being naively replaced by their corresponding estimators. The usage of this naive plug-in 64 estimator generally needs to be more cautious since it ignores the additional variability stemming from the need to estimate the dispersion parameters. Lee and Nelder (1996) proved that this information loss is asymptotically negligible for h-likelihood. Breslow and Clayton (1993) and McGilchrist (1994) did some simulations to investigate this problem for their methods. Now we showed that the asymptotic variance of /3 for our orthodox BLUP approach is exactly the same as the naive plug-in estimator. This result coincides with the noted observation from simulation studies for penalized quasi-likelihood method (Breslow and Lin 1995). That is, inference about regression parameters is mainly affected by the bias, rather than the variance, of the dispersion parameter estimators. 5.2.3 Heterogeneity In practice, heterogeneity is often substantial across clusters. To account for the heterogeneity in the marginal distributions of the responses, we may allow different dispersion parameters for different clusters. One of the choices would be to assume the similar models as those in Section 3.1, but with dispersion parameters (a , u> , p ) 2 being replaced by cluster specific dispersion parameters (a ,u ,p ). 2 2 2 2 2 All the previous derivations of orthodox BLUP random effects predictors and estimated score function remain valid if we replace (a ,u ,p ) by (a , to , pf). The dispersion parameters can 2 2 2 2 2 now be estimated through the following equations: + l_ Ji c(i)/4- + c(ij) - 2p c{i)w p PihAj 3= 1 65 2 2 ij j (5.18) and 1 ~2 Pi ~f J 1 Ji / , (v. ij n / y i j=i ij k=l n fj .. \2 ijh ijk) p \ iP iji ijk u X 1 ^ 1 ^ +yE — E Ji j—i n<ij _ J l c(ij)ii - X 2 p k 3 / • (5-19) Mi Mi? fc=1 That is, the appropriate quantities are now averaged only within each cluster to account for the heterogeneity. Therefore consistency of dispersion parameter estimators (a) , p ) for large number of clusters is not available for this situation. 2 2 In addition, dispersion parameters (a ,uj ,p ) often appear in the literature when 2 2 2 the designs are balanced with respect to the number of sub-clusters, that is, J\ = • • • = J = J- Such dispersion parameters can be estimated as follows: m u (Uij-^Ui) 1_ ™ 2 m 3 fr[ 2 p i p l j and ^2 3 1 j=l i j m + ^> \ ijk I fc=l n ijH'ijk) u P'i^ijPijk IgigcW)^_ m i=l ^ t j fc=l ( 5 2 1 ) MiMij The index j here often represents time in longitudinal data settings where the dispersion parameter estimators are obtained by averaging the appropriate quantities across time. These two estimators are consistent for large number of clusters. 66 Chapter 6 Residual analysis and computational procedure In this chapter, we briefly discuss residual analysis and computational issues. 6.1 Residual analysis Residual analysis for both the responses and the random effects is an important ingredient of our approach. Note that the marginal residuals of the responses can be decomposed into three uncorrelated residual components as follows: Yijk — PiPijPijk = (Yijk — PijkUij) + (Uij — Hij)Hij Ui k + (Ui — HijHijHijk- (6.1) These three components actually correspond to the residuals for the three levels of distributional assumptions given in the model. Thus we may check those distributional assumptions via estimated residuals Yij — HijJJij, Uij — [HjUi k a n d Ui — All three estimated residuals have zero mean and can be standardized to have unit variance. Actually the variances of the three estimated residuals are already available from (5.16), (5.14) and (5.9). We now define the appropriate types of standardized 67 (estimated) residuals for the three levels of distributions. Level 1: r Ui - Hi (i) (6.2) = Level 2: Uij \l^ PiP ij 2 q Pij Ui ~ c{i)H j - c(ij) + 2 (6.3) 2p c{i)wijH j 2 2 Level 3: ijk ^p HiHijH 2 The level 2 and 3 residuals, HijkUij v and ik + (6.4) c(ij)ix\ k , are residuals for conditional distributions, so we will call them conditional residuals. The level 1 residual is the residual for the first level random effects distribution. We can also check the marginal distributional assumptions of the response and the second level random effects through the following marginal residuals: Ui r (3) = The basic idea in residual analysis is then to use plots of standardized residuals in much the same way as in standard generalized linear models. Plots of standardized residuals against each covariate are useful for detecting nonlinearity of the data relative to the model. To check the log link assumption, we plot logYyfc against the 68 log fitted values logPijkUijk- Ideally this plot should show a horizontal linear rela- tionship; curvature or other unusual shape would indicate inadequacy of the log link assumption. As in generalized linear models, we may check the form of the variance functions and hence the distributional forms by plotting the standardized residuals against the corresponding fitted values. In order to check the variance function of the distribution for the random effects, we thus plot the standardized residuals against log fitted values. 6.2 Computational procedure Initial values for regression parameter estimates are taken as the regression parameter estimates obtained from standard generalized linear models techniques assuming independent responses. We take initial random effects predictions for Ui and Uij as the average of the responses within cluster i divided by the average of all responses and the average of the responses within sub-cluster divided by the average of all responses, respectively. That is of = m ^i=\ j 2^j=l { 2"k=\ n i j ijk I and 1 x^ ij n *3 U v. _i_ sr^m J_ y>Ji m ^»=1 Ji £"j=l riij \-^ ij ^k=l _1_ n -y The initial dispersion parameter estimates are calculated from Pearson estimators via (5.10), (5.15) and (5.17), but omitting the bias-correction terms. More specifically, 69 we have i™{Ui- (o) = — a in? E and U) ( o ) - ~ L j L " b i = l i TTTF. j=l J > riF-ij and P O) ( = - m E T E — i-l Ji j=l ij n E k=l ViVijV' ijk The algorithm then iterates between updating the regression parameter estimates via the Newton scoring algorithm, updating random effect predictors via the orthodox BLUP and updating dispersion parameter estimates via the adjusted Pearson estimators. 70 Chapter 7 Conventional, semiparametric and binomial mixed models In the literature of multiplicative random effects models, the marginal means of the random effects are usually taken to be 1 (Morton 1987; Thall and Vail 1990; Firth and Harris 1991; Lee and Nelder 1996). Imposing these conventional constraints on the Tweedie mixed models is equivalent to taking fa = I and faj = 1 for all i and j. We thus call a Tweedie mixed model with such contraints a conventional Tweedie mixed model. After discussing the conventional Tweedie mixed models in the following section, we will discuss mixed models beyond the Tweedie family which are related to the conventional Tweedie mixed models. The first model is a semiparametric mixed model with only first and second monents assumptions about the random effects; the second one is a mixed model assuming log-normally distributed random effects; whereas the last model involves binomial distributions for the conditional responses. 71 7.1 Conventional Tweedie mixed models To be more specific, we give the assumptions for the conventional Tweedie mixed models: BI) Given U = u, Y 1 U , Y U n , Y n i j Y U i j n i ,Y j m J m l ,Y m J j n r i m J m are con- ditionally independent, and the conditional distribution of Y ^ , given U = u, depends on only as Y \U ijk = u ~ Tw (ij,ij Uij,p Uij ~ ) 2 p 1 p k P 2 = where n ijk UijTw (p,ij ,—), p k Uij = exp(xT. /3). fe B2) Given U* = u», Un, ...,U j m m are conditionally independent, and the condi- tional distribution of Uij, given U* = u*, depends on Ui only, which is £/y-|U, = u» ~ Tw (ui, cu 2 Ui - ) l g q to UiTw,(l,—). Hi 2 = B3) Ui,...,U m iid Tw (l,(7 ). 2 r The derived random effects predictors, estimated score function and adjusted Pearson estimators for the more general models in Section 3.1 remain valid with the constraints p; = /iy = 1. Since the detailed expressions for the random effects predictors, estimated score function and adjusted Pearson estimators will be useful to discuss related models beyond the Tweedie family, we will present these expressions in the following sections. 72 7.1.1 Covariance and variance function The first and second moments of the second level random effects possess the following simple structure E[U ] = 1, l3 and ' a +u 2 if (s,t) = 2 cov[[/ ,cy = ^ a if s = i and t ^ j 2 si 0 otherwise. The covariance between the random effects and responses is simply given by Cov[U , Y ] = 5(s,i) [a + 5(t,j)tu } /j, . 2 st 2 ijk ijk The first and second moments of the response for the conventional Tweedie mixed models are then given by E[Vij/t] — pij , k and P f4 2 2 ijk if (s, t) = o VijkViti if s = i and t ^ j 0 otherwise. 2 i]k Cov[Y Y if (s, *, I) = (i, j, k) (CT + uj )ii piji 2 sth + (o + 0J )p? 2 jk andl^k ijk 2 The covariance of the response displays an interpretable variance components structure which clearly shows the variance contribution from each level of random effects. 73 The marginal variance of Y i3k for conventional Tweedie mixed models has a simple form as follows: Var(Yy ) = p ^ + (a + a, )/, .,. 2 2 fc 2 (7.1) 2 k The case p = 1 and p = 2 which correspond to Poisson-Tweedie models and gammaTweedie models are especially interesting. The marginal variances of these two models are 1 Var(Y^) = p + (a + u) )p , 2 ijk 2 2 jk and Var(y ) = (p + 2 ijfc a 2 + to )tf , 2 jk respectively. Therefore the marginal variance functions of these two models coincide with those of the negative binomial and gamma distributions, respectively. It is known that the marginal distribution of Y of the conventional Poisson-gamma models with ijk only one level of random effects (to — 0) follows a negative binomial distribution. 2 This negative binomial distribution is frequently adopted to account for overdispersion relative to Poisson model (Venables and Ripley 1994). It is unclear, in general, if the marginal distributions of the responses of the Poisson-Tweedie and gammaTweedie models follow the negative binomial and gamma distributions, respectively. It is known that the distribution of a random variable is uniquely characterized by its variance function within the exponential dispersion models (J0rgensen 1997); however, the marginal distribution of Y ijk may not follow an exponential dispersion model (J0rgensen 1987). Hence whether the marginal distributions of the responses of these two models follow the negative binomial and gamma distributions, respectively, remains an open question. 74 Finally, we give an expression for the correlation between the responses: Uorr(rj,- , Y ) f c stl = . .— ylftiik + (° 2 + rfHkJPAu — . + (° 2 + <* )tit 2 Clearly the correlation matrix still depends on both the regression and the dispersion parameters for the conventional Tweedie mixed models. 7.1.2 Random effects predictors The random effects predictors for conventional Tweedie mixed models also have simple expressions as follows: r f _ 1 + ° 1 where w = l/(p + u 2 {j 2 2 s £ i ffii , 7 - + ° 2-3=1 Lk=l i3Vijk W 1%$), and Uij = p WijUi + UJ Wij 2 Vijk i3k- 2 Y (7-3) k=l The conventional Poisson-gamma model coincides with the conjugate Poissongamma model studied by Lee and Nelder (1996) when there is only one level of the random effects. For this model both the orthodox BLUP predictors and MHLEs for the random effects are exactly the conditional expectations of the random effects given the responses; therefore both orthodox BLUP and MHLE estimation procedure lead to maximum likelihood estimates for (3 when dispersion parameters are known. In general, the orthodox BLUP predictors and MHLEs are not identical. Lee and Nelder (1996) presented a class of conjugate hierarchical generalized linear models with one level of random effects, and provided an explicit expression for MHLEs of the random 75 effects. As with any multiplicative conjugate hierarchical generalized linear models within this class, a direct calculation shows that the MHLEs for the random effects are linear functions of the response. However, the orthodox BLUP predictors of the random effects are the best linear unbiased predictor for Tweedie mixed models. The estimates for Cov(/3, U ) and Var(U—U) are useful in making inferences about realized or sample values of U (Harville 1976). Lee and Nelder (1996) mentioned that they have not found consistent estimates for these quantities except in few special cases. They further conjectured that consistent estimates may not exist under their regularity condition: the number of clusters is fixed when the number of observations increases. We have not obtained estimates for Cov(/3, U ) for the orthodox BLUP approach yet, but exact expressions for Var(U — U ) are available below: Cov[Ui - U U - U ] = 5( )c(i) h s s Covfft - Ui, U - U ] = st at Gov[Uij - U^, Ust - U ] = 5( ,i) st s Sii S( ,i)P c(i)w 2 s it p c(i)w ]} 2 . 2 [p Wij[5^j)Uj •+ 2 it These are clearly consistent estimates for Var(U — U ) when the parameter estimators are plugged in. Besides the general small dispersion and large sample asymptotics established in Section 4.4, we give another small dispersion asymptotics result: c ^ a ( p + ^^in^nijQmiiijjfe^ ^)) 2 < ~ p 2 2 +a; min (ny)min (tfr ) p 2 j ^ ^ 2 :/fc fc and 76 + a mmj{nij)mmj {p ~ )Ji 2 2 k p k P> c(« < . 2. . + CT ) 2 2 ,2 ^ , _ \ " J _ ,>- (7-5) P + ^ min (n )min (^) j ij R V jfc Obviously, when p ^ 2, we also have Uij - A Uij as m i n ^ / x / ) oo. 2 For Poisson-Tweedie (p = 1) and compound Poisson-Tweedie (1 < p < 2) models, m\Uj (iJ? ~ ) —> oo is equivalent to all /^fc -> oo. For positive stable Tweedie p k i3 k models (p > 2), including the inverse Gaussian, min (p r ) 2 jfc p oo is equivalent to all fc Vijk ->• 0. Clearly when p = 2, the same large sample asymptotics conclusion as in Section 4.4 holds without any restrictions on /XjjfcS. Actually, the large sample asymptotics conclusion upon the second level random effects now holds for any p under much relaxed conditions such as mmj(nij)mmj (fi ~ ) 2 k —>• oo as min,-(njj) —> oo. A sufficient p k condition is that m i n ^ / i ^ ) > clog(minj(ny))/minj(nij) for a positive constant c. 2 That is, the only restriction is that p,ij should not tend to zero too fast for the Poisk son and compound Poisson cases, whereas min (p 7 ) should not tend to infinity too 2 jfc J p fc fast for positive stable case. Actually, for conventional Tweedie mixed models, only the behavior of Uij affects the orthodox BLUP estimation procedure. This will become apparent in Section 7.2. 7.1.3 Parameter estimation The estimated score function now reduces to m W ) = Ji nij E jE= \ E *m^t{Yij 1=1 l-p k k=l p 77 - U i m k \ (7.6) whereas the matrix form of the estimated score function is still given by ip(P) = X diag (E(Y)) Var (Y) (Y - E(Y)). T _1 In addition, the expression for sensitivity matrix S(/3) also has the following simplified version: m S(/3) Ji nij J £ c ( z ) ( £ zZ ijPiik^ijk)(zZ = i=\ m j=lk=l Ji . . 2 zZ i3hhik*ijk) w T j=lk=l ritj "y " i=lj=l fc=l ^h^ijk i = l j = l jfc=l P 7.1.4 Uij t w fc=l ^ijk^ijk Adjusted Pearson estimators The adjusted Pearson estimators have simpler forms: m 1 * = -i " ™ 1 ^ 1 * = - E i=l (-) 7 E {(ft* - ft) + c(tj) + c(i) - 2p c(i)w \ j=l 2 1 ? = ltjt^"±{ 7.2 t=i 2 ij T j=i • 2 m .2 m 1) + - £ ' ( 0 , 2 J i j=i 'Hj =l k " ~S" ' (V I k Mijfc lU i)2 , + cdfinl-A 8 (7.9) . (7.10) J Random effects beyond Tweedie family For the general Tweedie mixed models defined in Section 3.1, the estimation procedure of the orthodox BLUP approach depends on the specific parametric forms of the random effects distributions among Tweedie family of distributions through 78 Pi, Uij, q and r. Since we have //; = pij = 1 for conventional Tweedie mixed models, this dependence of the estimation procedure on the specific parametric random effects distributional forms forms is then via q and r only. However, for conventional Tweedie mixed models, q and r actually do not enter the estimation procedure in either the estimated score function, the sensitivity matrix, the random effects predictors or the adjusted Pearson estimators. That is, the estimation procedure of the orthodox BLUP approach to conventional Tweedie mixed models is robust against misspecification of the random effects distributions within the Tweedie family. This robustness indicates that the orthodox BLUP approach to conventional Tweedie mixed models depends on the random effects only via the first and second moments of the second level random effects. Therefore we can extend the orthodox BLUP approach to deal with conventional models with only the first and second moment assumptions about the second level of random effects. These models are usually referred as semiparametric random effects models with nonparametric random effects. 7.2.1 Nonparametric random effects To discuss this semiparametric model in more detail, we assume that Cl) Given U = u, F m , ...;Y ,Yyi,Yij ,Y j x,Y j nnn nij m m are con- m mnmJm ditionally independent, and the conditional distribution of Y^, given U = u, depends on only as Yi \U jk = u ~ Twp(pi Uij,p Uij ~ ) 2 1 p jk P UijTw (fj,ij ,—), Uij 2 = where pijk = exp(xy /3). For case p = l , fc p (7.11) k namely Poisson distribution, p = 1. 2 79 C2) The random effects Un,..., are positive random variables with Uj m m E(U ) = l, lj and Cov[U , st \a 2 Uij] — 5( s>i) + <J( j)W } • 2 t) Here we make no further parametric assumptions on the random effects. In addition, we do not explicitly assume any first level random effects. Note that the random effects Uu,..., U mJm have the same first and second mo- ments as in conventional Tweedie mixed models; therefore it follows from assumption Cl) that Cov[U ,Yijk], st E[Yij \\J] k and Cov[Y u,Yij ] s k are also the same as in the conventional Tweedie mixed models since the derivation based on the conditioning techniques in Section 3.1 did not involve the Tweedie distributional assumption of Uij. These moments then imply that the orthodox BLUP predictors of the random effects U^, the estimated score function and the sensitivity matrix have exactly the same expressions as in conventional Tweedie mixed models. The adjusted Pearson estimator for p also has the same expression as in conven2 tional Tweedie mixed models. However the adjusted Pearson estimators for LO and 2 a for conventional Tweedie mixed models involve Ui although we do not assume the 2 existence of Ui. On the other hand, in conventional Tweedie mixed models Ui is a linear function of the responses and the adjusted Pearson estimators are moment estimators; therefore the estimating equation induced by the adjusted Pearson estimators for to and a remain valid if we use the linear expression for Ui as an intermediate 2 2 quantity to estimate to and a in the current semiparametric models. The difference 2 2 80 is that the linear expression for Ui does not necessarily have a parametric interpretation beyond the Tweedie family. In summary, the estimating equations associated with the estimated score function and adjusted Pearson estimators given in (7.6), (7.8), (7.9) and (7.10) are still unbiased estimating equations for both the regression and dispersion parameters in these semiparametric models. The regression and dispersion parameter estimators are also consistent asymptotically normal with the inverse of — S(/3) being the asymptotic variance of /3 under the same regularity conditions. As the expression for Var(U — U) in (4.2) c(ij) = Var((7 - U ) = E(U ) ij i3 i3 - Cow(U , Y )Var" (Y )Cov(Y , U ) 1 i3 i i l i3 depends on the first and second moments only, the expression for c(ij) will not be changed. Thus both small dispersion and large sample consistency results Uij in conventional Tweedie mixed models can be extended to these semiparametric models under the same regularity conditions. In conclusion, besides its parametric interpretation with the Tweedie random effects distributions, the orthodox BLUP approach to generalized linear mixed models is robust against misspecification of the random effects distributions. 7.2.2 Log-normal random effects An interesting example of non-Tweedie distributions for U^ in the conventional Tweedie mixed models setting is the log-normal. Let V^s be log-normal random 81 variables logVy ~ A T ( 0 , r ) , 2 with T = log(l + a + to ) and covariance structure 2 2 2 S(.,i) {a + S(t,j)u } 2 cov[v.„vy= , Vij\ Clearly U - = i3 V l + ^ + u ) 2 - 2 1 + ff2 + w 2 satisfies E(£/«) = 1, and {cr + <5(j)o; } . Cov[U , Uij] = 2 st 2 t Assume Cl) in the last section holds, then we would have logE(l^lU) = x; /3 + logf/y = xJ /3 - I l g(l + a + a; ) + log V , fc 2 fc 0 2 I3 where, under log link, the term \ log(l+a +a; ) affects only the value of the intercept. 2 2 Now the conditional expectations of the responses given the random effects can be expressed as a linear combination of regression parameters plus the normal random effects. This example shows the connection between generalized linear models with Tweedie random effects distributions and those with the normal random effects distribution. Note that the covariance structure assumptions are made about the lognormal random effects instead of the normal random effects; therefore, unlike Tweedie random effects cases, the partially observed joint log likelihood for both the response and log-normal random effects is not explicitly given here. 82 7.2.3 Binomial mixed models Binomial mixed models are usually handled by logistic-normal and binomial-beta models because of their analytic tractability. While our orthodox BLUP approach to Tweedie mixed models is not directly applicable to this case since the binomial distribution does not belong to the Tweedie family, we can deal with binomial mixed models via the conventional Poisson-Tweedie models based on the the following wellknown relationship between the binomial and Poisson random variables: Suppose that Y\ and Y are independent with 2 Y ~ k Poisson(fj,), k then Fi + Y ~ Poisson (pi + /i ), 2 2 and Fi|Fx + Y = N ~ binomial(7V, ^/(^ 2 + fi )). 2 (cf. McCullagh and Nelder 1989). Suppose that there are m pairs of observations, (Ri, Ni) i = 1,2,... ,m, Ri being the number of successes in Ni trials. Consider the following paired Poisson-Tweedie models: Yn = Ri\U = u ~ Poisson (un exp[xj(a Yi = Ni — Ri\XJ = 2 where F , Y" ,..., F , Y u 12 m l m2 + (3)f) u ~ Poisson (u expfx^a]j , i2 are independent given U — u. Assume further that the random effects U = ( U \ , U , U n , U j ) m m m T satisfy assumption B2) and B3) for conventional Tweedie mixed models with J, = 2 for all i as usual. Then we have 83 Yn\Yn + Y = Ni,U ~ i2 binomial(Ni,pi), where Pi Un exp[x7(a + /3)] Un -exp[x7(a + /3)] + U expfx^a i2 Hence logit pi = xj(3 + logitU + U U t l n i2 Clearly a = cu = 0 corresponds to binomial models without the random effects and 2 2 a = 0 is equivalent to one level random effects models. 2 The binomial-beta model is the special case with a = 0 and U - ~ Gamma(l,w ) 2 2 i3 j = 1,2. This is because Vi = .^ . is known to follow a symmetric beta distribution if U^ ~ Gamma(l,w ) j = 1,2. Asymmetric beta distributed random effects can u 2 u 2 easily be incorporated by allowing different dispersion parameters to for U^ j = 1,2. 2 Lee and Nelder (1996) also derived the binomial-beta model from a pair of Poissongamma models with the common covariates xj, but different regression parameters (3j for each of the two groups (J = 1,2). Under logit link, they have logit pi = x; (P - p ) + logitUn + Ua' x 2 with (3 — /3 corresponding to our (3 here. But they derived their estimating equa1 2 tions for MHLEs directly from the binomial-beta model. We derive our binomial mixed models based on paired Poisson-Tweedie models with the common regression parameter (e* ,/3 ) , but different covariates for each of T T T the two groups, that is, (x^, x^) for Ri and (xj, 0) for A ^ - i i j . This approach enables 84 us to obtain the fixed effects estimate /3 as well as its standard error and the random effects predictors Vi for binomial mixed models by fitting the paired Poisson-Tweedie models. Clearly J3 is a consistent estimate for (3 as m — > oo since (a,/3) is shown to be a consistent estimate for (a,/3) in Section 5.2.2. There are many different approaches to Poisson mixed models in the literature. The paired Poisson mixed models trick provides a way of handling binomial mixed models via Poisson mixed models. The extension of this trick to handle multinomial mixed models via Poisson mixed models is straightforward. Such multinomial data examples can be found in McCullagh and Nelder (1989). 85 Chapter 8 Data analysis Illustrative examples for the orthodox BLUP approach to analyzing data of different types based on conventional Tweedie mixed models are presented in this chapter. We will concentrate mainly on the analysis of the epilepsy data to illustrate the use of variations of the models. We will also analyze the seed germination data and cake baking data to illustrate the use of orthodox BLUP approach to handle binomial and continuous data. 8.1 Count data In this section, we illustrate the orthodox BLUP approach to count data with analyses of the epilepsy data described in Section 2.4.1. Preliminary analysis indicated that the counts were much lower during the fourth visit so Breslow and Clayton (1993) introduced a linear trend covariate Visit (coded (-0.3,-0.1,0.1,0.3)). To facilitate visual understanding of the data, Breslow (1996) plotted the logarithm of the seizure counts reported at baseline (visit 0) and at each of the four visits as in Figure 8.1. Counts at the four follow-up visits were increased by 0.5 when taking log-transform 86 0 1 2 3 4 visit number Figure 8.1: Profile plot of log transformed epilepsy seizure counts. to avoid infinities. The baseline counts were divided by four since the baseline counts were measured over eight weeks instead of two weeks. The baseline seizure counts are less variable because the period over which they were measured was four times as long as for each of the subsequent visits. Patient 207 had exceptionally high counts at baseline count and all subsequent visits. 8.1.1 Overview of the previous analyses Breslow and Clayton (1993) reanalyzed the epilepsy data using the penalized quasilikelihood approach, whereas Lee and Nelder (1996) also reanalyzed the same data based on hierarchical generalized linear models. To make a systematic comparison of 87 these analyses, we summarize their models as follows. Lee and Nelder (1996) considered a model where the response Y given the random effects (U, V) follows a Poisson distribution and Yij\U = u , V = v where Us and V~ijS ~ Poisson(iJ,ijUiVij) i = 1,..., 59, j — 1,2,3,4, (8.1) are all mutually independent having E(£/j) = E(Vjj) = 1 and Var(£/j) = a and Var(Vjj) = <r. Lee and Nelder assume that Ui and 2 2 follow gamma distributions. Thall and Vail (1990) considered a similar model, their model 22, but with V~ij replaced by Vj with known first two moments only. Breslow and Clayton (1993) also considered models similar to Lee and Nelder's models. Since these approaches yield very similar results, we list the results of Lee and Nelder in Table 8.1 for later comparisons in this section. Model GLM in Table 8.1 is the standard Poisson regression model. Model HGLM 1 has Oj = 0 for all j = 1,2,3,4, namely the model of random subject effects only. Model HGLM 2 has Vij = Vi x Visit,- where Visit, is the jth component of covariate Visit for j = 1,2,3,4. The model HGLM 3 is the one with o\ = cr = cr and 0 4 = 0. 2 3 Lee and Nelder set a = 0 since it is reported to tend to zero when fitting. By the way, 4 the estimated GLM regression coefficient for covariate Visit was reported as -0.29 in their paper which appears to be a typographical error. 88 Table 8.1: HGLM parameter estimates for the epilepsy data. GLM HGLM 1 HGLM 2 HGLM 3 Parameter est.* s.e.* est. s.e. est. s.e. est. s.e. Constant -2.80 0.41 -1.35 1.20 -1.32 1.20 -1.21 1.20 0.95 0.04 0.88 0.13 0.89 0.13 0.89 0.12 -1.34 0.16 -0.89 0.39 -0.86 0.39 -0.86 0.38 Base.Trt 0.56 0.06 0.34 0.20 0.33 0.20 0.32 0.19 Age 0.90 0.12 0.51 0.36 0.49 0.36 0.46 0.35 -0.08 0.10 -0.29 0.10 -0.16 0.16 -0.28 0.14 Base Trt Visit 0.27 0.26 0.22 o\ 0.40 0.15 °l 0 0 a 2 * est. and s.e. represent estimates and standard errors respectively. 8.1.2 Poisson-Tweedie models Our reanalysis of these data is based on conventional two level Poisson-Tweedie models. The fixed effects are logarithm of a quarter of baseline seizure counts (Base), logarithm of age (Age), Trt and Visit. Since there is no repetition for each visit, the third index is thus omitted. We denote yij the seizure count for patient i on visit j; therefore the random effects Ui and ity are at patient and visit levels respectively. The formulation is as follows: Yij\\3 = u ~ Poisson ( / i j j U y ) , whereas t y u , = u„ ~ 89 Tw (l,4^), g (8.2) and Ui,...,U iid 59 Tw (l,a ), 2 r where j — 1, 2,3,4 and i = 1,... ,59. The conditional independences corresponding to assumptions BI) and B2) given in Section 7.1 are assumed to hold. As explained in Chapter 7, this model can also be regarded as a semiparametric model. We consider four models with different assumptions for the dispersion parameters Model(o> = 0), Model(a; ), Model(w ) and Model(a>?) assume a; 2 ujfy 2 2 2 = 0, LO , to and LO in (8.2), respectively. Thus Model(o; = 0) is the same one level 2 2 2 2 model as HGLM 1, but without the gamma distributional assumption for the random effects. Model(a> ) is the model with equal to for all (i, j) corresponding to HGLM 2 2 3. Model(a> ), namely the model with distinct co for different visits, has the same 2 2 covariance structure of the response as that of model 22 proposed by Thall and Vail through crossed random effects design. The last model, Model(cj ), has distinct co 2 2 for each subject. We start with the Model ( C J 2 ), including all possible two way interaction terms among the covariates. We removed all those interaction terms except.the interaction term Base.Trt from the model based on the Wald test and residual analysis. The p-value of Wald test for this last interaction term is slightly larger than 5%, but we retained it in all our models to allow a systematic comparison. That is, we take Base, Age, Trt, Visit and Base.Trt as the covariates. We fit the four models and present the results in Table 8.2. 90 Table 8.2: Parameter estimates for the epilepsy data. Model(cj = 0) 2 Model(cj ) Model(cj ) est. est. 2 2 Model(cj ) 2 Parameter est.* Constant -1.35 1.22 -1.35 1.22 -1.37 1.16 -1.11 0.72 0.88 0.14 0.88 0.14 0.87 0.13 0.94 0.07 -0.89 0.41 -0.89 0.42 -0.91 0.40 -0.49 0.28 Base.Trt 0.34 0.21 0.34 0.21 0.35 0.20 0.03 0.13 Age 0.51 0.36 0.51 0.36 0.52 0.34 0.36 0.21 Visit -0.22 0.10 -0.22 0.22 -0.28 0.23 -0.35 0.18 Base Trt s.e.* s.e. est. 0.19 0.16 0.007 0.36 0.33 0.03 w\ 0.32 0.08 u\ 0.67 0.54 "I 0.25 12.9 a 0.28 s.e. 2 u, 2 s e. * est. and s.e. represent estimates and standard errors respectively. 8.1.3 Model checking We have done model checking through residual analyses. The normal plots for Model(cj ), Model(cj ) and Model(a; ) are displayed in the first, second and third 2 2 2 columns, respectively of Figure 8.2. The plots in the three rows correspond to the level 1, 2 and 3 residuals defined in Section 6.1. As expected, Model (CJ?) exhibits the least curvature since there are many more dispersion parameters estimated. Unlike Model (CJ ), the dispersion parameter estimators are known to be consistent for the 2 Model(cj ) and Model(cjJ). There is not much difference between the normal plots 2 for Model(cj ) and Model(cj ). 2 2 We plot the standardized residuals against fitted values for the level 2 and level 91 Normal plot of laval 1 ratiduali Model(u; ) 2 Notrrnl plot ol IMI 2 raaiduaU Model(o; ) 2 Normal ploi ol laval 3 ratidualt Model(o; ) 2 Model(o; ) 2 Normal ptol ol laval 2 Models ) 2 Normal ptol ol I oval 2 raiiduaJi Model (u; ) 2 Model(o; ) 2 Normal plot ol loval 3 raiiduali Model(o; ) 2 Model(w ) 2 Figure 8.2: Normal plots of level 1, 2 and 3 residuals for epilepsy data. 92 Standardized residuals vs fitted values Model(a; ) 2 Standardized residuals vs fitted values Residuals vs lilted values Residuals vs fitted values Model(u; ) Model(w ) Residuals vs log (fitted values) Residuals vs log (fitted values) 2 2 leg<!B«fv«iM> Model(o; ) 2 Model(a; ) Model(o; ) 2 2 Figure 8.3: Scatter plots of level 2 and 3 residuals for epilepsy data. 3 distributions in Figure 8.3. The level 1 residual plot is omitted since the fitted values are fixed at 1. The residual plots for these three models are arranged in the same fashion as the normal plots. The residuals for the level 3 distribution from the three models are approximately around zero, but display distinct lines due to the discreteness of the Poisson distribution. The residuals of Model(u; ) for the level 2 2 distribution show a slightly megaphone shape, whereas the residuals of Model(o; ) for 2 the level 2 distribution exhibit an upward trend. The residuals of Model (u ) for the 2 level 2 distribution show less pattern than the other two. 93 Based on the residual plots, Model(o; ) appears to fit slightly better than Model ( C J 2 2 ). But the inferences about regression parameters are not significantly different for these two models. The dispersion parameter estimators are shown to be consistent for these two models as well. On the other hand, the normal plots for Model(u; ) look less 2 curved. 8.1.4 Comparison of different approaches Thall and Vail (1990) found their model 22 had the best fit among their fitted models. The numerical results from Breslow and Clayton (1993) were reported as similar to those obtained from this model 22. Analyses from the hierarchical generalized linear models and penalized quasi-likelihood are also very similar except for the intercept estimators. The regression parameter estimates and the standard errors obtained from Model(cj = 0) and Model(cj ) are quite similar to those in their counterpart 2 2 from HGLMs except our standard error estimates are slightly larger. Thus all different approaches yield essentially the same results about the regression parameters although different random effects distributional assumptions were made. The conclusions based on Model(c<; = 0), Model(c<;) and Model(a; ) are similar 2 2 2 to those previous studies. The interaction between the treatment and baseline is retained in the models since the approximate p-values for this interaction is just slightly larger than a = 0.05. As Thall and Vail indicated, this means the predicted mean seizure rate for the treatment group is either higher or lower than that for the placebo group, accordingly as the baseline count does or does not exceed a critical threshold (minus the ratio of the regression parameter for Trt to the regression coefficient for 94 the interaction term). This seems to suggest that the new drug progabide may be contraindicated for patients with high seizure rates. However they pointed out this suggestion should be treated with caution since it is based on a single dataset and a particular family of models. In fact, the interaction between the treatment and baseline is statistically insignificant in Model (of). Based on the Wald test, the p-value of dropping the interaction term from Model (LO ) is as high as 0.8. That is, the contraindication of the new drug 2 progabide for patients with high seizure rates disappears after accounting properly for the marginal heterogeneity among patients. Dropping this interaction term, the p-value for the treatment effect in this model becomes extremely small. That is, the effectiveness of this new drug progabide is statistically significant. Thall and Vail (1990) focused their analyses on the regression parameters in the log-linear models for the marginal seizure rates and the covariance parameters in various patterned covariance matrices. They carefully examined the marginal residuals comparing the observed and fitted counts of each subject at each visit, and identified patients 207, 135, 225, 227 and 112 as 'outliers'. The primary interest of Breslow and Clayton (1993) is systematic identification of patients who had extreme levels, or extreme degrees of change over time, in their seizure rates. They identified patient 135 as having marked improvement over time after an initially high seizure rate, and identified patients 227, 225 and 112 as having the highest overall count levels relative to expectation based on the covariates. They also identified patient 232 as having low and zero counts. However Lee and Nelder were more cautious in claiming outliers and indicated a single outlier, patient 227. 95 The marginal heterogeneity among patients is substantial, even after accounting for the treatment and baseline variables. In all models except Model(w ), the dis2 persion parameters are assumed to remain the same across subjects. That is, the marginal variance only changes with the marginal mean in a systematic way. However, in biological studies, due to the genetic, environmental or other unknown factors, the marginal heterogeneity across subjects often exceeds what could be explained by the mean change. We plot the sample variances s against sample means y~i across 2 subjects in Figure 8.4 to study the mean-to-variance relationship. The boxplots of the responses for different subjects are also displayed in Figure 8.4 to serve this purpose. As indicated by (7.1), the variance is a monotone function of the mean for Tweedie mixed models with equal dispersion parameters. Thus the two plots in Figure 8.4 should reflect this monotone pattern if the dispersion parameters do not vary across subjects. These two plots do show a kind of monotone pattern, but not very closely. To account for the irregular marginal heterogeneity, we considered Model (a; ) 2 where different dispersion parameters u> are allowed across patients. 2 This model seems to have effectively captured this irregular heterogeneity. The minimum, first quartile, third quartile and the maximum of Qjs are displayed in Table 8.2 labeled as to for j = 1,2,3,4. The estimated u; s range from 0.03 to 12.9 corresponding 2 2 to patients 213 and 225 respectively. Model(u; ) also seems to be very sensitive in 2 detecting outliers. The scatter plot of cDs in Figure 8.5 identifies patients 135, 227, 2 112, 207 and 225 with large Q . These patients match outliers reported by Thall and 2 Vail (1990) and Breslow and Clayton (1993). 96 S c a t t e r plot 2 0 4 0 within subject s a m p l e m e a n s B o x p l o t s of within s u b j e c t responses B= i«.«01 • • •—1— —•—•—•—•—•—•—•—•—•—r—c^jcOojcQcQojoocOc^ Figure 8.4: Heterogeneity plots for epilepsy data: (a) scatter plot of within subject sample variances versus sample means; (b) boxplots of within subject responses for all patients. 97 ( Figure 8.5: Heterogeneity plot for epilepsy data. 98 8.1.5 Computational aspects With regards to computing, we ran Splus on a Sparc ultra machine. All parameter estimates for the Model(u; ) updated in monotone directions after five iterations. 2 Those included regression and dispersion parameters as well as standard errors for regression parameters. All estimates at the fifth iteration are close to those in Table 8.2. The same two decimal point precision for all parameter estimates was achieved at the eighteenth iteration. At this iteration, the absolute value of the standardized components of ib((3) were less than 0.001. The absolute distance between regression parameter estimates and their updates were less than 0.00003. We stopped the program after sixty iterations. The maximum of the absolute values of the standardized components of i/>(/3) and the absolute distance between regression parameter estimates and their updates were less than 2E — 12 and 3E — 13, respectively. Each iteration took about 10 seconds of user time and 0.01 seconds of system time on the Sparc ultra machine. Twenty iterations took about four minutes. As to the other two level random effects models, the computing time for Model(a; ) 2 was slightly longer than that for Model(a; ). For Model(c<;), all parameter estimates 2 2 stabilized after eight iterations, but they converged much more slowly although still in a monotone direction. Those estimates approached the values reported in Table 8.2 after about fifty iterations. Each iteration took about 15 seconds on the same machine. 8.2 Binomial data We illustrate the orthodox BLUP approach to binomial data with analysis of the seed germination data described in Section 2.4.2. 99 8.2.1 Paired Poisson-Tweedie models We analyze the seed germination data based on the paired Poisson-Tweedie models. The model notation is similar to that for the epilepsy data. The estimates and corresponding standard errors are displayed in Table 8.3. Table 8.3: Parameter estimates for the seed germination data based on paired Poisson mixed models. Paired GLM Parameter Model(a - 0) Model(o; ) 2 Model^J) 2 est.* s.e.* s.e. s.e. Constant -0.558 0.126 0.349 0.181 0.179 Seed (2) 0.146 0.223 0.512 0.289 0.287 Extract (2) 1.318 0.177 0.475 0.250 0.248 Interaction -0.778 0.306 0.708 0.398 0.395 est. est. a 2 0.265 co 2 s.e. est. 0.229 0.226 0.042 0.055 0.026 co\ * est. and s.e. represent estimates and standard errors respectively. As indicated in Chapter 7, the asymmetric binomial mixed models can be considered based on the paired two level Poisson-Tweedie models with the second level random effects having unequal dispersion parameters co s between the two groups 2 j = 1,2. We started our modeling with this asymmetric model, Model (co ), but there 2 is little evidence for this asymmetry as the difference between u 2 relative to a . 2 and 0J\ is small Thus we proceed to fit the symmetric paired Poisson mixed model with one and two levels of random effects, namely Model(a = 0) and Model(t<; ). 2 100 2 The estimates and standard errors presented in the first and second columns were obtained from the standard paired Poisson generalized linear model, namely, the paired Poisson generalized linear models without random effects. This standard paired Poisson generalized linear model corresponds to the standard binomial generalized linear models. The regression parameter estimates and their corresponding standard errors obtained from these two models are exactly the same. The regression parameter estimates obtained from all above paired Poisson-Tweedie models with various distributional assumptions for random effects are almost identical to those obtained from the standard binomial generalized linear models. Breslow and Clayton (1993) also reported similar results from their simulation study based on some simulated data sets generated from a balanced binomial mixed model design. Their marginal quasi-likelihood (MQL) approach to binomial mixed models led to regression coefficient estimates identical to those obtained from standard logistic regression; however neither the MQL nor PQL regression parameter estimates for this seed germination data were identical to those obtained from standard logistic regression. On the other hand, the standard error estimates obtained from the paired Poisson mixed models are much larger than those obtained from the standard paired Poisson regression models. Furthermore, the standard error estimates clearly vary from one random effects model to another. That is, the standard error estimates depend heavily on the random effects model assumptions. This is different from the results for the epilepsy data. 101 8.2.2 Model checking We performed residual analysis based on the paired Poisson-Tweedie models. Since the procedure is exactly the same as those presented in last section, we present the residual plots for Model(a; ) only. The normal plots of the level 1, 2 and 3 residuals, 2 for Model(u; ) are displayed in Figure 8.6 and the scatter plots of level 2 and 3 resid2 uals are displayed in Figure 8.7. Normal plots at all three levels exhibit moderate curvature, whereas the scatter plots of level 2 and 3 residuals against log-fitted values show moderate upward trends. Some patterns usually exist for a small data set with discrete responses. 8.2.3 Comparison of different approaches Breslow and Clayton (1993) and Lee and Nelder (1996) also analyzed the seed germination data. The estimates and the corresponding standard errors are presented in Table 8.4. Ta Die 8.4: Paramet:er estimates for t ie seed data based on binomial mocels. GLM PQL HGLM Parameter est.* s.e.* est. s.e. est. Constant -0.558 0.126 -0.542 0.190 -0.543 0.187 Seed (2) 0.146 0.223 0.077 0.308 0.080 0.303 Extract (2) 1.318 0.177 1.339 0.270 1.337 0.265 Interaction -0.778 0.306 -0.825 0.430 a 2 0.098 s.e. -0.822 0.423 0.045 * est. and s.e. represent estimates and standard errors respectively. The standard errors obtained from the paired Poisson-Tweedie models with one 102 N o r m a l plot of level 1 q u a n t i l e s of s t a n d a r d residuals n o r m a l N o r m a l p l o t o f l e v e l 2. r e s i d u a l s q u o n t i l o o of s t a n d a r d normal N o r m a l plot of l e v e l 3 qLiantilos o f s t a n d a r d residuals normal Figure 8.6: Normal plots of residuals for seed germination data. 103 Residuals vs fitted values 0.2 0.4 O.e 0.8 L O level 1 r a n d o m effects 1.2 1.4 1.6 Residuals vs log(fitted values) CM I.O 1.5 2.0 2 . 5 log(fitted v a l u e s ) 3.0 3.5 Figure 8.7: Scatter plots of residuals for seed germination data. 104 4 . 0 level of random effects are almost twice those obtained from penalized quasi-likelihood and maximum hierarchical likelihood approaches, whereas the symmetric paired Poisson with two levels of random effects model gives similar, but slightly smaller standard error estimates than those obtained from penalized quasi-likelihood and maximum hierarchical likelihood approaches. The Model(cr = 0) corresponds to binomial-beta 2 model (Lee and Nelder 1996); however the former gives quite different standard error estimates than the latter. This difference is not easy to explain. Further investigation is needed. 8.3 Continuous data We illustrate our approach to clustered continuous data with an analysis of the cake baking data described in Section 2.4.3. Cochran and Cox (1957) analyzed the cake baking data using analysis of variance. The normal plot of residuals for the ANOVA model in Figure 8.8 shows a curvature. Firth and Harris (1991) re-analyzed the cake baking data using multiplicative models and they found strong support for the hypothesis of a constant coefficient of variation. Thus we analyze the cake baking data based on the conventional gamma mixed models, Model(a; ), where the two levels of 2 random effects are at the batch and observation levels. Temperature and recipe are taken as factors. We screen these factors based on the Wald test. We found that the recipe effect and any interaction between temperature and recipe appear negligible, but the temperature effect is highly statistically significant. In fact, Figure 8.9 shows that there is an increasing trend in breaking angles as the temperature increases. Cochran and Cox (1957) and Firth and Harris (1991) reached similar results to 105 ours except they found the recipe effect is significant at the statistical significance level of 0.05. The normal plots of level 1, 2 and 3 residuals in Figure 8.10 do not reveal serious patterns except some curvature for the level 1 residuals. The scatter plot for level 2 residuals in Figure 8.11 does not exhibit any serious pattern, whereas the scatter plot for level 3 residuals shows an upward trend. Noting that the dispersion parameter p is nearly zero, one would expect that the effect of level 3 residuals 2 on the model be small. To have a systematic comparison, we display estimates and corresponding standard errors for ANOVA, gamma and gamma mixed models without interaction in Table 8.5. The parameter estimates from ANOVA model are quite different from those obtained using standard generalized linear model and gamma mixed model; however the regression parameters for the latter two models are again almost identical for this balanced design. After accounting for random effects, the recipe factor becomes less statistically significant, whereas the temperature factor becomes more statistically significant. 106 Table 8.5: Parameter estimates for cake baking data. ANOVA Parameters Model(o; ) 2 s.e. est. s.e. est. s.e. 32.122 0.474 3.466 0.015 3.466 0.030 RI 0.989 0.821 -0.023 0.018 -0.023 0.037 R2 0.819 0.474 -0.008 0.010 -0.008 0.021 Tl 0.598 0.335 0.034 0.026 0.034 0.014 T2 1.092 0.260 0.028 0.015 0.028 0.008 T3 0.647 0.212 0.020 0.010 0.020 0.006 T4 -0.739 0.581 0.033 0.008 0.033 0.005 T5 -0.261 0.335 0.020 0.007 0.020 0.004 Constan a est. GLM 2 0.038 0.019 LO 2 P 0.00025 2 R and T represent recipe and temperature factors. 107 108 175 185 195 205 225 temperature Figure 8.9: Boxplots of cake baking data by temperature. 109 N o r m a l plot oflevel 1 r e s i d u a l s -2 -1 O qtiantileo of s t a n d a r d 1 N o r m a l plot oflevel 2 -3 -2 -1 residuals O quantiloo of s t a n d a r d 2 n o r m a l 1 2 3 2 3 n o r m a l N o r m a l plot oflevel 3 residuals -3 -2 -1 O quantilea of s t a n d a r d 1 n o r m a l Figure 8.10: Normal plots of residuals for cake baking data. 110 L e v e l 2 r e s i d u a l s v s level 1 r a n d o m e f f e c t s o.a 1.2 1.0 1.4 level 1 r a n d o m effects Level 3 residuals v s log-fitted v a l u e s • • • • * * • s - i • ! !: ! • • . • I t • • • • • * •" : . - * • • - t t • • 1 • I * '. • - • • • • i • : : •- : " •• • * • • t• • - • • • • • 3.0 3.2 3.4 log-fitted 3.6 3.8 4.0 values Figure 8.11: Scatter plots of residuals for cake baking data. Ill Chapter 9 Simulation In this chapter, we evaluate the performance of the orthodox BLUP approach based on a simulation study. We will focus on Poisson mixed models since these models play an important role in the analysis of discrete data. Illustrative examples of applications of Poisson mixed models to analyze count data and binomial data were presented in the previous chapter. A possible connection between Poisson mixed models and the random effects Cox proportional hazard models for censored data will be discussed in next chapter. 9.1 Results over 100 simulations To investigate the performance of the orthodox BLUP approach under realistic conditions, we 'replicate' the epilepsy data set 100 times via simulation using Poissongamma model. To generate the random effects and responses via simulation from this model, we take the covariates Constant, Base, Trt, Age, Visit and Base.Trt of epilepsy data as the covariates for simulation. The corresponding regression and dispersion parameter estimates for the epilepsy data are taken as true model parameters. There are listed in Table 9.1 as 'true values' (3 , p , ..., @ , a and LO , respectively. Each 2 0 x 112 5 2 of these 100 data sets is then simulated through the following three steps: • Step 1: Generate 59 samples from Gamma(l, a ), denoted by Vq°\ . . . , u^; 2 • Step 2: Generate 4 samples from Gamma(uf , w « f ) for each i, denoted by 2 j = 1,2,3,4 and i = 1,..., 59; • Step 3: Generate a sample from Poisson (uij exp(x^/3)) for each denoted by y£\ j = 1, 2, 3,4 and i = 1,..., 59, where the design matrix X is taken exactly the same as that of epilepsy data. The regression and dispersion parameter were fixed at the corresponding estimates for epilepsy data, that is, (3 = (1.30,0.88,-0.88,0.50,-0.23,0.34), a = 0.24 and 2 co — 0.44, respectively. The 100 data sets are obtained by repeating this procedure 2 for k = 1,...,100. We analyzed each of these 100 data sets based on the standard Poisson regression model (GLM) and Poisson-Tweedie models with and without degree of freedom correction for small samples, denoted by BLUP.c and BLUP, respectively. The summaries are presented in the next subsection. 113 9.1.1 Summary statistics The averages of regression and dispersion parameter estimates over 100 simulations are displayed in Table 9.1. Table 9.1: Averages of parameter estimates over 100 simulations. f% Pi Ps PA f% true value -1.30 0.88 -0.88 0.50 -0.23 0.34 GLM -1.37 0.88 -0.78 0.51 -0.26 0.29 BLUP -1.36 0.87 -0.89 0.51 -0.25 BLUP.c -1.31 0.87 -0.88 0.50 -0.25 P2 a 2 cu 2 0.24 0.44 0.34 0.17 0.50 0.34 0.23 0.61 The regression parameters are reasonably estimated by both GLM and orthodox BLUP approaches with or without small-sample correction. On the other hand, the dispersion parameters a and cu are underestimated and overestimated respectively. 2 2 This phenomenon is different from simulation studies reported by Breslow and Clayton (1993) where the dispersion parameters were always underestimated. However the possibility of overestimation was previously predicted by Engel and Keen (1996). The small-sample correction alleviates the negative bias of the estimate of a , but 2 worsens the positive bias of the estimate of cu . 2 The standard errors of the estimates over 100 simulations and the averages of the 100 estimated standard errors are termed as simulated and estimated standard errors, respectively (Lin and Breslow 1996b). Table 9.2 displays the simulated and estimated standard errors. The expressions of the simulated and estimated standard errors for Pi are as follows: 1 simulated s.e.(A) = , — 100 Y(p\ - 114 k) 1 — 100 Y B\ ) , k) 2 and 1 1 0 0 estimated s.e.(0i) = ' 100 s.e.(3^). v i Table 9.2: Simulated and estimated standard errors GLM BLUP BLUP.c s.e. . 0o 0i A 0s Sim. 1.65 0.17 0.60 0.48 0.28 0.31 Est. 0.41 0.04 0.15 0.12 0.10 0.06 Sim. 1.42 0.15 0.43 0.41 0.22 0.23 Est. 1.25 0.14 0.42 0.36 0.24 0.22 Sim. 1.42 0.15 0.43 0.41 0.22 0.23 Est. 1.39 0.47 0.41 0.26 0.24 0.16 a 2 ui 0.07 0.07 0.07 0.09 As 2 * Simulated s.e.: s.e. of estimates over 100 simulations. ** Estimated s.e.: average of 100 estimated s.e.s. Clearly the GLM standard error estimates are seriously negatively biased. In contrast, the orthodox BLUP approach with small-sample correction slightly overestimates standard errors for regression parameters, except for the intercept (3 . On the 0 other hand, the orthodox BLUP approach without small-sample correction slightly underestimates the standard errors for the regression parameters, severely so for the intercept. Table 9.3: Mean squared errors over 00 simulations. 02 a 2 f% 0i GLM 2.73 0.03 0.36 0.23 0.08 0.10 BLUP 2.01 0.02 0.18 0.17 0.05 0.05 0.01 0.01 BLUP.c 2.01 0.02 0.19 0.17 0.05 0.05 0.01 0.04 03 04 05 LO 2 To evaluate the overall performance of the estimates over 100 simulations, we dis115 play the mean squared errors (MSEs) in Table 9.3. As expected, the mean squared » errors of regression parameter estimates based on standard Poisson model are larger than those based on the orthodox BLUP approaches. According to mean squared errors, the orthodox BLUP approach without small-sample correction performed slightly better than the orthodox BLUP approach with small-sample correction. All three approaches gave exceptionally large mean squared errors of the intercept ft. 9.1.2 Confidence and prediction intervals Based on the simulated and estimated standard errors, we constructed simulated and estimated 95% confidence intervals for the parameters and 95% prediction intervals for the random effects (u\ and un) assuming normality of the estimators and predictors. We present in Table 9.4 the counts of the number of times the true values are covered by 95% confidence intervals or prediction intervals. The kth confidence intervals for ft are defined as follows: simulated C.I. for ft = - 1.96 simulated s.e.(ft), $[ k) + 1.96 simulated s . e . ( f t ) ) , and estimated C.I. for ft = ($* - 1.96 estimated s.e.0[ ), 0{ - 1.96 estimated s.e.0[ )) . } k) k) k) As expected, the GLM estimated 95% confidence intervals performed very poorly. On the other hand, the estimated 95% confidence intervals based on the orthodox BLUP approaches performed reasonably well, especially those with small-sample correction. The orthodox BLUP approaches also gave reasonable 95% prediction intervals, but with lower coverage counts for u\ and higher coverage counts for u . n 116 Table 9.4: Coverage counts of 95% confidence and prediction intervals. GLM BLUP BLUP.c 9.1.3 Ui A) Ai A> A3 Ai As Sim. 95 96 94 94 94 93 Est. 35 28 38 35 44 31 Sim. 96 97 96 95 94 96 88 100 Est. 90 96 94 91 97 92 84 98 Sim. 96 97 94 95 94 97 92 100 Est. 95 96 95 94 98 97 91 100 Normality of parameter estimates and random effects To check the normality of parameter estimators, we did a normal plot for each of regression and dispersion parameter estimates over the 100 simulations. The results are displayed in Figure 9.1. The departures from normality do not seem serious even for the dispersion parameters. To assess the performance of the orthodox BLUP predictors, we plotted the simulated random effects versus their orthodox BLUP predictors over the 100 simulations. In particular, we plotted u[ ^ versus k and versus u[ i over k = 1,..., 100 in the k upper rows of Figure 9.2. Perfect prediction would lead to diagonal lines in these two plots. The prediction for the second level of random effects seems to be better than that for the first level random effects. Normal plots of the orthodox BLUP predictors of random effects Ui and un over 100 simulations are displayed in the lower row of Figure 9.2. The clear curvatures in these normal plots imply that random effects predictors are not well approximated by normal distributions. The numerical results and plots from this simulation study show that the param- 117 3 : y y Po Pi P2 y ; / s y S Pz PA Ps 3 Figure 9.1: Normal plots of parameters over 100 simulations: x and y axes are quantiles of standard normal distribution and ordered parameter estimates over 100 simulations, respectively. 118 Simulated vs predicted first level random effects over 100 simulations o.o 0.5 1.0 1.5 2.0 Simulated vs predicted second level random effects over 100 simulations 2.5 1 simulated first level random effects simulated second level random effects Normal plot of first level random effects Normal plot of second level random effects -1 - 1 quantiles of standard normal 2 - 1 0 1 2 quantiles of standard normal Figure 9.2: Plots for random effects predictors iii and u u 119 over 100 simulations. eters are reasonably well estimated through orthodox BLUP approach. 9.2 Residual analysis In Chapter 8, we did model checking through residual plots. Serious irregularities exhibited in residual plots may indicate serious violation of model assumptions. To help to set up standards for judging irregularity exhibited by residual plots, we will study the behavior of residuals through simulation in this section. In this section, we make comparison and contrast of residual plots for the epilepsy data and simulated data. The normal plots for the level 1, 2 and 3 residuals for the epilepsy data and simulated data are displayed at row 1, 2, and 3 of Figure 9.3 where the normal plots for the epilepsy data and simulated data are on the left and right hand side, respectively. The histograms and scatter plots for the epilepsy data and simulated data are displayed in Figure 9.4 and 9.5, respectively, in the same fashion as the normal plots. The simulation shows that the residual plots may exhibit some patterns even if the data were generated from a Poisson-gamma model. That is, moderate patterns in residual plots may not indicate violation of the model assumptions. 120 Normal plot of level 1 residuals Normal plot of level 1 residuals (simulated data) quantiles of standard normal quantiles of standard normal Normal plot of level 2 residuals Normal plot of level 2 residuals (simulated data) quantiles of standard normal quantiles of standard normal Normal plot of level 3 residuals Normal plot of level 3 residuals (simulated data) quantiles of standard normal quantiles of standard normal Figure 9.3: Comparison of normal plots of level 1, 2 and 3 residuals for epilepsy data and simulated data. 121 Histogram of level 1 residuals Histogram of level 1 residuals (simulated data) -2 -1 Histogram of level 2 residuals Histogram of level 2 residuals (simulated data) Histogram of level 3 residuals Histogram of level 3 residuals (simulated data) Figure 9.4: Comparison of histograms of level 1, 2 and 3 residuals for epilepsy data and simulated data. 122 Level 2 residuals vs level 1 random effects Level 2 residuals vs level 1 random effects (simulated data) (•vol 1 random affects Level 3 residuals vs log-fitted values Level 3 residuals vs log-fitted values (simulated data) k>g(fltlnd valuas) Figure 9.5: Comparison of scatter plots of level 2 and 3 residuals for epilepsy and simulated data. 123 C h a p t e r 10 Discussion 10.1 Conclusion In this thesis we have introduced a new class of generalized linear mixed models, Tweedie mixed models, and adopted the orthodox BLUP in the fitting algorithm. This approach has several advantages. 1. Tweedie mixed models take both the distributional shape and intra-dependence of clustered data into account. The resulting variance component decomposition of the structure of the covariance matrix of responses is not only very interpretable, but also makes the model fitting much simpler than modal predictor approaches; all expressions for the estimating equations were explicitly derived based on (3.13) so there is no need to invert large matrices. 2. The orthodox approach provides us a common method for computing parameter estimates and random effects predictors for all Tweedie mixed models; therefore we can study Tweedie mixed models as a single class, rather than as a collection of 124 unrelated different models. Furthermore, the conventional Tweedie mixed models are shown to have close connection with well-known Poisson-gamma, binomial-beta and Poisson-lognormal models. The semiparametric interpretation of the conventional Tweedie mixed models allows specification only the first and second moments of the unobserved random effects. 3. The orthodox BLUP for random effects in (4.1), estimated score function for regression parameter in Theorem (5.1) and adjusted Pearson estimators for dispersion parameters in (5.10), (5.15) and (5.17) show that our orthodox BLUP approach relies on only the first and second moment structure of (U, Y). It is important to note that the orthodox BLUP approach requires specification only of the first and second moments of the random effects model. Tweedie mixed models play an interpretable role within this larger family of models just as exponential dispersion models do within generalized linear models. 4. The asymptotic justification of orthodox BLUP approach does not rely on any kind of normality of random effects. In contrast, the asymptotic justifications of penalized quasi-likelihood and maximum hierarchical likelihood approaches largely rely on the approximate normality for random effects or 'the right transformed normality for random effects on the right scale' , respectively. It is important to note that our asymptotic variance of regression parameter estimator is not affected by the variability of estimating dispersion parameters. Our orthodox BLUP approach leads to the same estimating equation for the regression parameters as that obtained by the generalized estimating equation approach. This estimating equation was also reached via the quasi-likelihood function 125 (McCullagh 1983) and multilevel modelling (Goldstein 1995) approaches. All these approaches except the orthodox BLUP approach are marginal modelling instead of conditional modelling approaches, that is, they are useful to make the 'populationaveraged' instead of 'subject-specific' inferences on the mean (Schabenberger 1996; Burton et al. 1998). On the other hand, the orthodox BLUP approach explicitly incorporates random effects and mergers the 'subject-specific' and the 'populationaveraged' inferences using the same model. In addition, the marginal covariance structures for our models are generally different from those considered by the other three approaches. In the development of random effects modelling methodologies, the model checking has generally been ignored (see e.g. Zeger and Karim 1991; Breslow and Clayton 1993) or done via normality checking of residuals (Lee and Nelder 1996); however the justification of the normality of residuals has not been found in the literature. It appears from our simulation study that moderate departures from normality may not necessarily indicate departures from the model assumptions. We illustrated the orthodox BLUP approach to analyses of clustered count, binomial and continuous data using Tweedie mixed models. In addition, our compound Poisson-Tweedie models may be applied to positive continuous data with a positive probability component at zero. One data example is car insurance data, for example, where the zero component corresponds to the case of no claims for a given insurance policy and the positive part of the distribution is the claim-size distribution (J0rgensen and Souze 1994). Another data example is the customers' expenditures in a certain shopping center. 126 In this thesis, we take dispersion parameters as nuisance parameters and their standard errors are not estimated; however it would be of interest to develop tests for the presence of random effects as well as issues such as over-dispersion and heterogeneity. Testing for such hypotheses is complicated because the null values lie on the boundary of the dispersion parameter space. Bootstrap methods may be used to construct such confidence intervals (Lele 1991), but these are controversial as regards their ability to fully reflect the relevant uncertainties. 10.2 Further study In this section, we discuss some further research of the orthodox BLUP approach to generalized linear mixed models. 10.2.1 More than two levels of random effects We have presented the Tweedie mixed model with two levels of random effects in detail. This model is useful to handle three-level hierarchical structure; however, hierarchies involving more levels also occur frequently in practice (Goldstein 1995). One such example is multi-center longitudinal clinical trials with centers further nested within geographical regions such as cities or countries. The extension of the Tweedie mixed models with two levels of random effects is straightforward. The estimation procedure for more than two levels of nested random effects models is much the same as that for two-level models. The estimated score function and sensitivity matrix have the same global matrix expression as in Theorem 5.1. The key derivation of the orthodox BLUP approach is the calculation 127 of the orthodox BLUP prediction of random effects which involves the inverse of the covariance matrix of the response. However this inverse can be obtained by repeatedly applying (3.13) to the covariance matrix of the response, twice for the models with two levels of random effects, three times for three levels of random effects, and so on. We have derived explicit expressions for all quantities of interest. These expressions are useful for theoretical study. In practice, those explicit expressions are not necessary for computing purposes. Calculation of the covariance matrix of the response can be programmed based on (3.13), and all quantities of interest can be evaluated in terms of this inverse. This remark will be especially useful to resolve the computing issue for these extended models. 10.2.2 Crossed designs In our analyses of the epilepsy data, we considered the nested random effects only. In clinical trials, the nature and degree of the variability of such epileptic seizure counts over time may be as important as its average behavior (Thall and Vail 1990). This aspect of the epilepsy data may be analyzed through the following simple crossed factor Tweedie mixed model with p — 1: Yijk\U = u,V where Ui,..., U , V\,...,V m n = v ~ Tw p {p (ui ijk + Vj), p (ui 2 + U j ) 1 _ p } , are mutually independent with E(Ui) = E(Vj) = \ . Replacing the random effects in the partially observed score function by their 128 orthodox B L U P predictors would give us an estimated score function as follows: i=lj=lk=l P = X diag(E(Y))Var (Y)(Y-E(Y)). T _1 This estimating equation is clearly unbiased, but its asymptotic properties have yet to be studied. Ways of handling crossed designs are still rather limited in the literature of generalized linear mixed models. Further work on crossed factor random effects models would be of interest. The development of models for complicated crossed factor designs such as salamander mating data (McCullagh and Nelder 1989) is under way. 10.2.3 Survival data analysis Incorporating random effects into Cox proportional hazard models has gained increasing attention in analyses of epidemiological or event history data. However, these models pose considerable theoretical difficulties in the development of estimation and inference procedures (Clayton 1991). Following Whitehead's (1980) idea of fitting the traditional Cox model using Poisson modelling techniques, we showed elsewhere that the random effects Cox model can also be fitted using random effects Poisson modelling techniques. To be more specific, we consider a Cox model with two levels of random effects as follows: Let the kth individual in the j t h sub-cluster of the ith cluster be indexed by k). Let the hazard function for individual k) at time t be denoted by hijk(t). We assume that, given random effects U = u, the hazard functions for individuals 129 are conditionally independent with h (t) = h (t)Uij ijk 0 exp(x; /3), (10.1) FC where h (t) is the baseline hazard function and Ui s are positive random effects, or 0 3 'frailties', shared by all individuals of the same cluster. Suppose further that • T i , . . . , r g are distinct death times; • m/i is the multiplicity of deaths at time T>,; • the risk set at r h is TZ(r ) = {(i,j,k) h : t i3k > Th}, where t i3k is the observed survival time for individual (z, j, k). Let Yijk,h be 1 if individual (i,j,k) dies at Th and 0 otherwise. Thefittingof random effects Cox model is equivalent to fitting the following auxiliary random effects Poisson model: Y ,h\^J ijk —u ~ Poisson (u exp(a i3 h . + xj f3f) jk (10-2) Given random effects, Peto's version of the conditional partial likelihood (Cox and Oakes 1984) and conditional Poisson likelihood are v£(3- Y I U - u) - TT " ( ^ t t f o ) ^ ( m h l ) and /(«, /3; Y | U = u) = JI " ' " ^^ ' > X p ^ +* 5 f ^ * (10.4) It can be shown that, given random effects, Peto's version of the conditional partial likelihood and conditional Poisson likelihood lead to the same maximum likelihood 130 estimates for regression parameters. Hence the orthodox BLUP approach may be adopted to fit the random effects Cox models via Poisson mixed models. 131 Appendix A D a t a sets Table A.l: Seed germination data (a)seeds: O. Aegyptiaca 75 and 73, (b)root tracts: bean and cucumber 0. Aegyptiaca 75 Bean 0. Aegyptiaca 73 Cucumber Bean Cucumber r n r/n r n r/n r n r/n r n r/n 10 39 0.26 5 6 0.83 8 16 0.50 3 12 0.25 23 62 0.37 53 74 0.72 10 30 0.33 22 41 0.54 23 81 0.28 55 72 0.76 8 28 0.29 15 30 0.50 26 51 0.51 32 51 0.63 23 45 0.51 32 51 0.63 17 39 0.44 46 79 0.58 0 4 0 3 7 0.43 10 13 0.77 132 Table A.2: Epilepsy data: successive two-week seizure counts for 59 epileptics, (a) Trt: treatment 0=placebo, l=progabide); (b) Base: eight-week baseline seizure counts; (c)Age (in years). ID n Base Age Trt 3 3 3 11 31 0 3 5 3 3 11 30 0 2 4 0 5 6 25 0 114 4 4 1 4 8 36 0 116 7 18 9 21 66 22 0 118 5 2 8 7 27 29 0 123 6 4 0 2 12 31 0 126 40 20 23 12 52 42 0 Vi 104 5 106 107 Y 2 Y 3 130 5 6 6 5 23 37 0 135 14 13 6 0 10 28 0 141 26 12 6 22 52 36 0 145 12 6 8 4 33 24 0 201 4 4 6 2 18 23 0 202 7 9 12 14 42 36 0 205 16 24 10 9 87 26 0 206 11 0 0 5 50 26 0 210 0 0 3 3 18 28 0 213 37 29 28 29 111 31 0 215 3 5 2 5 18 32 0 217 3 0 6 7 20 21 0 219 3 4 3 4 12 29 0 220 3 4 3 4 9 21 0 222 2 3 3 5 ' 17 32 0 226 8 12 2 8 28 25 0 227 18 24 76 25 55' 30 0 230 2 1 2 1 9 40 0 234 3 1 4 2 10 19 0 238 13 15 13 12 47 22 0 101 11 14 9 8 76 18 1 102 8 7 9 4 38 32 1 133 Table A.2: continued ID Yi Yz Y 2 Yt Base Age Trt 103 0 4 3 0 19 20 1 108 3 6 1 3 10 30 1 110 2 6 7 4 19 18 1 111 4 3 1 3 24 24 1 112 22 17 19 16 31 30 1 113 5 4 7 4 14 35 1 117 2 4 0 4 11 27 1 121 3 7 7 7 67 20 1 122 4 18 2 5 41 22 1 124 2 1 1 0 7 28 1 128 0 2 4 0 22 23 1 129 5 4 0 3 13 40 1 137 11 14 25 15 46 33 1 139 10 5 3 8 36 21 1 143 19 7 6 7 38 35 1 147 1 1 2 3 7 25 1 203 6 10 8 8 36 26 1 204 2 1 0 0 11 25 1 207 102 65 72 63 151 22 1 208 4 3 2 4 22 32 1 209 8 6 5 7 41 25 1 211 1 3 1 5 32 35 1 214 18 11 28 13 56 21 1 218 6 3 4 0 24 41 1 221 3 5 4 3 16 32 1 225 1 23 19 8 22 26 1 228 2 3 0 1 25 21 1 232 0 0 0 0 13 36 1 236 1 4 3 2 12 37 1 134 Table A.3: Cake baking data: breaking angles (deg Temperature Rep.. 175° 185° 195° 205° 215° 225° 1 42 47 42 47 32 26 28 24 39 47 53 2 46 29 32 32 57 45 39 41 45 45 26 47 26 35 34 3 4 5 6 7 Recipe I Recipe II 26 24 24 24 23 33 27 33 13 14 15 33 28 29 24 26 39 31 28 40 28 33 27 31 29 32 8 9 10 11 12 23 28 27 43 24 37 29 27 32 33 35 33 31 34 30 33 35 37 23 33 30 43 33 25 40 37 31 33 31 28 39 29 40 1 39 46 51 49 55 42 2 3 4 35 34 25 47 42 28 39 35 46 52 42 37 61 35 37 5 6 7 31 24 22 46 30 26 30 29 25 35 29 26 40 24 29 8 9 10 11 12 26 27 21 29 29 26 24 32 24 31 28 27 31 34 27 27 32 37 28 31 35 36 35 36 37 33 30 33 29 30 35 33 31 30 23 21 25 21 22 19 21 28 26 27 20 1 . 2 46 44 43 43 45 43 46 46 48 47 63 58 3 33 24 40 37 41 38 4 38 21 24 41 38 30 36 25 31 30 35 30 33 37 35 23 31 21 24 30 24 21 21 28 28 35 35 28 14 15 5 6 7 20 23 32 23 26 24 27 28 35 13 Recipe III 30 22 35 37 35 31 22 25 33 21 35 8 20 24 9 24 23 18 10 11 26 28 27 26 27 28 24 25 30 26 25 35 38 33 28 29 28 43 28 33 28 37 19 21 22 28 . 27 25 25 25 31 35 25 12 13 14 15 135 25 33 35 Bibliography [1] Aichison, J. and Ho, C H . (1989). The multivariate Poisson log-normal distribution. Biometrika 76, 643-653. [2] Anderson, D.A. and Aitkin, M. (1985) Variance component models with binary responses: interviewer variability. J. Roy. Statist. Soc. Ser. B 47, 203-210. [3] Artes, R. and J0rgensen, B. (1998). Longitudinal data estimating equations for dispersion models. Technical report RT-MAE 9802, Institute of Mathematical Statistics, University of Sao Paulo. [4] Atkinson, A.C. (1981). Two graphical displays for outlying and influential observations in regression. Biometrika 68, 13-20. [5] Breslow, N.E. (1984) Extra-Poisson variation in log-linear models. Applied Statist. 33, 38-44. [6] Breslow, N.E. (1996). Generalized linear models: Checking assumptions and strengthening conclusions. Statistica Applicata 8, 23-41. [7] Breslow, N.E. and Clayton, D.G. (1993). Approximate inference in generalized linear mixed model. J. Amer. Statist. Assoc. 88, 9-25. [8] Breslow, N.E. and Lin, X. (1995). Bias correction in generalized linear mixed models with a single component of dispersion. Biometrika 82, 81-91. 136 [9] Brillinger, D. and Preisler, H.K. (1986) Two examples of quantal data analysis. Proceedings of the 13th International Biometrics Conference Seattle: The Biometrics Society. 95-113. [10] Brockwell, P.J. and Davis, R.A. (1991). Time Series: Theory and Methods 2nd ed. New York: Springer-Verlag. [11] Burton, P., Gurrin, L., and Sly, P. (1998). Tutorial in biastatistics. Extending the simple linear regression model to account for correlated responese: an introduction to generalized estimating equations and multi-level mixed modelling. Statist. Med. 17, 1261-1291. [12] Clayton, D.G. (1991) A Monte Carlo method for Bayesian inference in frailty models. Biometrics 47, 467-485. [13] Clayton, D.G. (1996). Discussion of the paper by Lee Y., and Nelder J.A., Hierarchical generalized linear models. J. Roy. Statist. Soc. Ser. B 58, 667. [14] Cochran, W.G. and Cox, D.M. (1957). Experimental Designs. 2ed. New York: Wiley. [15] Cox, D.R. and Oakes, D. (1984) Analysis of Survival Data. London: Chapman & Hall. [16] Crowder, M.J. (1978). Beta-binomial Anova for proportions. Appl. Statist. 27, 34-37. [17] Crowder, M.J. (1986). On consistency and inconsistency of estimating equations. Econometric Theory 3, 305-330. [18] Crowder, M.J. (1987). On linear and quadratic estimating function. 74, 591-597. 137 Biometrika [19] Dean, C B . and Lawless, J.F. (1989). A mixed Poisson-inverse Gaussian regression model. Canad. J. Statist. 17, 171-181. [20] Dempster, A.P., Laird, N. and Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm (with discussion). J. Roy. Statist. Soc. Ser. B 39, 1-38. [21] Draper, D. (1996). Discussion of the paper by Lee Y., and Nelder J.A., Hierarchical generalized linear models J. Roy. Statist. Soc, Ser. B 58, 662. [22] Engel, B. and Keen, A. (1996). Discussion of the paper by'Lee Y., and Nelder J. A., Hierarchical generalized linear models, J. Roy. Statist. Soc, Ser. B 58, 656. [23] Firth, D. and Harris, I.R. (1991). Quasi-likelihood for multiplicative random effects. Biometrika 78, 545-555. [24] Fitzmaurice, G.M. (1995). A caveat concerning independence estimating equations with multivariate binary data. Biometrics 51, 309-317. [25] George, E.I., Makov, U.E. and Smith, I.G. (1987). Conjugate likelihood distributions. Scand. J. Statist. 20, 147-156. [26] Geyer, C J . and Thompson, E.A. (1992). Constrained Monte Carlo maximum likehihood for dependent data. J. Roy. Statist. Soc, Ser. B 54, 657-699. [27] Gilks, W.R., Wang, C C , Yvonnet, B., and Coursaget, P. (1993). Random-effects models for longitudinal data using Gibbs sampling. Biomatrics 49, 441-453. [28] Gilmour, A.R., Anderson, R. and Rae, A.L. (1985) The analysis of binomial data by a generalized linear mixed model. Biometrika 138 72, 593-599. [29] Glifford, P. (1993). Discussion on the meeting on the Gibbs sampler and other Markov chain Monte Carlo methods. J. Roy. Statist. Soc. Ser. B 55, 53-54. [30] Godambe, V.P. (1976). Conditional likelihood and unconditional optimum estimating equations. Biometrika 63, 277-284. [31] Goldstein, H. (1995). Multilevel Statistical Models. New York: Wiley. [32] Gray, J., Jesson, D., Goldstein, H., Hedger, J. and Rasbash, J. (1995). A multilevel analysis of school improvement: changes in schools' performance over time. School Effectiveness and School Improvement, 6, 97-114. [33] Harvey, A.C. (1981). Time Series Models. Oxford: Allan. [34] Harville, D.A. and Mee, R.W. (1984) A mixed model procedure for analyzing ordered categorical data. Biometrics 40, 393-408. [35] Harville, D. (1977). Maximum likelihood approaches to variance component estimation and related problems. J. Amer. Statist. Assoc. 72, 320-340. [36] Hedeker, D. and Gibbons, R.D. (1994). A random-effects ordinal regression model for multilevel analysis. Biometrics 50, 933-944. [37] Henderson, C.R. (1975). Best linear unbiased estimation and prediction under a selection model. Biomatrics 31, 423-447. [38] Hinde, J. (1982) Compound Poisson regression models. In GLIM 82: Proceedings of International Conference on Generalized Models, ed. R. Gilchrist, Berlin: Springer. 109-121. [39] Hobert, J. and Casella, G. (1996) The effect of improper priors on Gibbs Sampling in hierarchical linear mixed models. J. Amer. Statist. Assoc. 91, 1461-1473. 139 [40] J0rgensen, B. (1987). Exponential dispersion models (with discussion). J. Roy. Statist. Soc. Ser. B 49, 127-162. [41] J0rgensen, B. (1997). The Theory of Dispersion Models. London: Chapman & Hall. [42] J0rgensen, B. and Labouriau, R.S. (1995). Exponential Family and Theoret- ical Inference. Lecture notes. Department of Statistics, University of British Columbia, Vancouver, Canada. [43] J0rgensen, B., Labouriau, R. and Lundbye-Christensen, S. (1996a). Linear growth curve analysis based on exponential dispersion model. J. Roy. Statist. Soc. Ser. B 58, 573-592. [44] J0rgensen, B., Lundbye-Christensen, S., Song, X.-K. and Sun, L. (1996b). A longitudinal study of emergency room visits and air pollution for Prince George, British Columbia. Statist. Med. 15, 823-836. [45] J0rgensen, B., Lundbye-Christensen, S., Song, X.-K. and Sun, L. (1996c). Statespace models for multivariate longitudinal data of mixed types. Canad. J. Statist. 24, 385-402. [46] J0rgensen, B., Martinez, J.R. and Tsao, M. (1994). Asymptotic behaviour of the variance function. Scand. J. Statist. 21, 223-243. [47] J0rgensen, B. and Souza, M.C.P. (1994). Fitting Tweedie's compound Poisson model to insurance claims data. Scand. Actuarial J. 1, 69-93. [48] Karim, M.R. and Zeger, S.L. (1992). Generalized linear models with random effects; Salamander Mating Revisited. Biometrics 48, 631-644. 140 [49] Knudsen, S.J. (1998). Estimating functions and Separate Inference. Cand. Scient. thesis, Department of Theoretical Statistics, Aarhus University. [50] Laird, N.M. (1978). Empirical Bayes methods for two-way contingency tables. 65, 581-590. Biometrika [51] Lambert, D. and Roeder, K. (1995). Overdispersion diagnostics for generalized linear models. J. Amer. Statist. Assoc. 90, 1225-1236. [52] Lawless, J.F. (1987). Negative binomial and mixed Poisson regression. Canad. J. of Stat. 15, 209-225. [53] Lee, Y. and Nelder, J.A. (1996). Hierarchical generalized linear models. J. Roy. Statist. Soc. Ser. B 58, 619-678. [54] Lele, S. (1991). Resampling using estimating equations. Theory of Estimating Equations.(ed. Godambe, V.P.) Oxford: Clarendon Press. [55] Liang, K.-Y. and Zeger, S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73, 13-22. [56] Lin, X. and Breslow, N.E. (1996a). Bias Correction in Generalized Linear Mixed Models with Multiple Components of Dispersion. J. Amer. Statist. Assoc. 9 1 , 1007-1016. [57] Lin, X. and Breslow, N.E. (1996b). Analysis of correlated binomial data in logistic-normal models. J. Statist. Comput. Simul. 55, 130-146. [58] Lipsitz, S.R., Fitzmaurice, G. M., Orav, E. J., and Laird, N. M. (1994). Performance of generalized estimating equations in practical situations. 50, 270-278. 141 Biometrics [59] McCullagh, P. (1983). Quasi-likelihood functions Annals of statistics 11, 59-67. [60] McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models. 2nd ed. London: Chapman & Hall. [61] McCulloch, C.E. (1997). Maximum likelihood algorithms for generalized linear mixed models. J. Amer. Statist. Assoc. 92, 162-170. [62] McGilchrist, C A . (1994). Estimation in generalized mixed models. J. Roy. Statist. Soc. Ser. B 56, 61-69. [63] McGilchrist, C.A.and Yau, K.K.W. (1995). The derivation of BLUP, ML, REML estimation methods for generalised linear mixed models. Commun. Statist., Theory and Methods 24, 2963-2980. [64] McLeish, D.L. and Small, C.G. (1988). The Theory and Applications Inference Functions. of Statistical Lecture Notes in Statistics 44, New York: Springer-Verlag. [65] Morton, R. (1987). A generalized linear model with nested strata of extra-Poisson variation. Biometrics 74, 247-257. [66] Natarajan, R. and McCullogh, C.E. (1995). A note on existence of the posterior distribution for a class of mixed models for binomial responses. Biometrika 82, 639-643. [67] Paik, M . C , Tsai, W Y . and Ottman, R. (1994). Multivariate survival analysis using piecewise gamma frailty. Biometrics 50, 975-988. [68] Rao, C.R. (1973). Linear Statistical Inference and Its Applications. York: Wiley. 142 2nd ed. New [69] Ripley, B.D. and Venables, W.N. (1994). Modern Applied Statistics with Splus. New York: Springer-Verlag. [70] Robinson, C.K. (1991). That BLUP Is a Good Thing: The Estimation of Random Effects. Statist. Sci. 6, 15-32. [71] Schabenberger, O. (1996). Population-averaged and subject-specific approaches for clustered categorical data. J. Statist. Simul. 54, 231-253. [72] Schall, R. (1991). Estimation in generalized linear models with random effects. Biometrika 78, 719-727. [73] Smith, A.F.M. and Roberts, G.O. (1993). Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods. J. Roy. Statist. Soc. Ser. B 55, 3-23. [74] Stiratelli, R., Laird, N. and Ware, J. (1984) Random effects models for serial observations with binary responses. Biometrics 40, 961-971. [75] Thall, P.F. and Vail, S.C. (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46, 657-671. [76] Waclawiw, M.A. and Liang, K.Y. (1993). Prediction of random effects in the generalized linear model. J. Amer. Statist. Assoc. 88, 171-178. [77] Whitehead, J. (1980) Fitting Cox's regression model to survival data using GLIM. Applied Statistics 29, 268-275. [78] Williams, D.A. (1982) Extra-binomial variation in logistic linear models. Applied Statist. 3 1 , 144-148. 143 [79] Wolfinger, R. and O'Connell, M. (1993). Generalized linear mixed models: A pseudo-likelihood approach. J. Statist. Cornput. Simul. 48, 233-243. [80] Zeger, S.L. and Karim, M.R. (1991). Generalized linear models with random effects; A Gibbs sampling approach. J. Amer. Statist. Assoc. 86, 79-86. [81] Zeger, S.L. and Liang, K.-Y. (1986). Longitudinal data analysis for discrete and continuous outcomes. Biometrics 42, 121-130. 144
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- An orthodox blup approach to generalized linear mixed...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
An orthodox blup approach to generalized linear mixed models Ma, Renjun 1999
pdf
Page Metadata
Item Metadata
Title | An orthodox blup approach to generalized linear mixed models |
Creator |
Ma, Renjun |
Date Issued | 1999 |
Description | We introduce a new class of generalized linear mixed models assuming Tweedie exponential dispersion model distributions for both the response and the random effects. This class of models accommodates a wide range of discrete, continuous and mixed data. By letting the random effects enter as weights as well as means in the conditional distributions, the variance matrix may be expressed as a sum of variance components. We consider an orthodox BLUP approach to parameter estimation and random effects prediction for this new class of models based on a predictor of the random effects that is truly best linear and unbiased, in contrast to the conventional BLUP which is the conditional mode. We obtain an optimal estimating equation based on the orthodox BLUP, which is solved by a modified Newton algorithm. This approach facilitates analysis of residuals and allows justification of asymptotic results under realistic conditions through standard estimating equation theory. An important feature of this approach is that the principal results depend only on the first and second moment assumptions of unobserved random effects. The common fitting algorithm based on orthodox BLUP enables us to study this new class of models as a single class, rather than as a collection of unrelated different models. This approach is illustrated with the analyses of seed germination data, epilepsy data and cake baking data. By means of asymptotic justifications, simulations and worked examples, we conclude that the orthodox BLUP approach is of practical value for analysis of clustered non-normal data. |
Extent | 5190557 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0089321 |
URI | http://hdl.handle.net/2429/10153 |
Degree |
Doctor of Philosophy - PhD |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1999-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1999-389340.pdf [ 4.95MB ]
- Metadata
- JSON: 831-1.0089321.json
- JSON-LD: 831-1.0089321-ld.json
- RDF/XML (Pretty): 831-1.0089321-rdf.xml
- RDF/JSON: 831-1.0089321-rdf.json
- Turtle: 831-1.0089321-turtle.txt
- N-Triples: 831-1.0089321-rdf-ntriples.txt
- Original Record: 831-1.0089321-source.json
- Full Text
- 831-1.0089321-fulltext.txt
- Citation
- 831-1.0089321.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0089321/manifest