LONGITUDINAL ANALYSIS FOR BINARY AND COUNT DATA By Jinko Graham B. Sc. (Mathematics) University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF STATISTICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 1991 © Jinko Graham, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature Department of 3 i A i 15 TICS The University of British Columbia Vancouver, Canada Date DE-6 (2/88) N.C. 33 199i ) Abstract Longitudinal data sets consist of repeated observations of an outcome over time, and a corresponding set of covariates for each of many subjects. In many fields, multivariate analysis-ofvariance is commonly used to analyse longitudinal data. Such an analysis is appropriate when responses for each subject are multivariate Gaussian with a common covariance matrix for all subjects. In many cases, however, the longitudinal response cannot be transformed to satisfy these assumptions. An alternate analysis might rely on specification of an appropriate likelihood to obtain estimates of regression parameters and their standard errors. Unfortunately, the correlation structure in the data for each subject may not be well understood, making such parametric modelling difficult. This thesis discusses two methods for the analysis of longitudinal data which require only minimal assumptions about the true correlation structure in the data for each subject to yield consistent estimates of regression parameters and their standard errors. The first method is based on the use of a "working" likelihood and extends the results of Zeger, Liang and Self (Biometrika, 1985) to the case of time-dependent covariates. The second method, first presented by Liang and Zeger (Biometrika, 1986), is based on quasi-likelihood theory. This method uses generalized estimating equations to arrive at consistent estimates of the regression coefficients and their standard errors, and can be applied to any longitudinal response with univariate marginal distributions for which the quasi-likelihood formulation is sensible. This ii includes Gaussian, Poisson, binomial (bernoulli), gamma, and inverse Gaussian distributions. Both methods are extensively illustrated using the results from an experiment on hummingbird learning. These methods enable much more information to be extracted from the hummingbird data set than a more traditional analysis-of-variance, and therefore provide useful and powerful tools for researchers in this subject area. iii Table of Contents ii Abstract^ List of Tables^ vii List of Figures^ viii Acknowledgement^ 1 2 3 xi Introduction 1 1.1 Background to the Experiment ^ 1 1.2 Description of the Experiment ^ 3 1.3 Description of the Data ^ 5 1.4 Outline of the Thesis ^ 6 Classical Analysis 9 2.1 Univariate Analysis-of-Variance for Repeated Measures ^ 12 2.2 The Statistical Model ^ 13 2.3 Results of the Analysis ^ 14 2.4 Conclusions and Limitations ^ 18 Analysis Based on a Working Likelihood Approach 20 3.1 20 Logistic Regression for Binary Repeated Measures iv ^ 3.2 The Working Likelihood Approach ^ 20 3.3 Incorporating Time-dependent Covariates ^ 28 3.3.1 The Independence Working Model ^ 28 3.3.2 The First-Order Working Model ^ 30 3.4 Exploring the Binary Hummingbird Data ^ 34 3.5 Model Fitting With Sex Excluded ^ 37 3.5.1 The Initial Model ^ 37 3.5.2 An Informal Check of the Working Likelihoods ^ 38 3.6 Model Reductions With Sex Excluded ^ 40 3.6.1 Strategy for Model Reductions ^ 40 3.6.2 Description of the Reductions ^ 41 3.7 Residual Analysis for the Final Model with Sex Excluded ^ 47 3.7.1 Analysis of the Independence Residuals ^ 47 3.7.2 Analysis of the First-Order Residuals ^ 49 3.8 Summary of Analysis with Sex Excluded ^ 51 3.9 The Initial Model With Sex Included ^ 53 3.10 Model Reductions with Sex Included ^ 54 3.10.1 Strategy for Model Reductions ^ 54 3.10.2 Description of the Reductions for the Independence Working Likelihood ^ 55 3.10.3 Description of the Reductions for the First-Order Working Likelihood . ^ 57 3.11 Summary of Analysis with Sex Included ^ 60 4 Analysis Based on the GEE Approach ^ 4.1 Description of the GEE approach ^ 4.1.1 Brief Description of Quasi-Likelihood ^ 63 63 64 4.1.2 Generalized Estimating Equations (GEEs) for Longitudinal Data ^ 66 4.2 GEE Analysis of the Binary Hummingbird Data ^ 70 4.2.1 Description of the Model Reductions ^ 71 4.2.2 Fits of the Final Model Obtained by Other Working Correlations ^ 74 4.3 GEE Analysis of the Binomial Hummingbird Data ^ 76 4.3.1 Exploring the Binomial Hummingbird Data ^ 77 4.3.2 Model Fitting for the Binomial Hummingbird Data ^ 79 4.3.3 Fits of the Final Model Obtained by Other Working Correlations ^ 84 5 Conclusion^ 87 5.1 Discussion of Results ^ 87 5.2 Comparison of Methods ^ 89 Figures^ 92 117 Bibliography^ vi List of Tables Table 2.1. Analysis of Variance for Blocks of 10 ^ 16 Table 2.2. Analysis of Variance for Blocks of 15 ^ 16 Table 3.1. Comparison of Z-scores from the Initial Model ^ 39 Table 3.2. Final Model Fits from the First Reduction Path ^ 43 Table 3.3. Final Model Fits from the Second Reduction Path ^ 46 Table 3.4. Estimated Autocorrelations for the Independence Residuals ^ 49 Table 3.5. Estimates of Transformed Lag One Autocorrelations ^ 50 Table 3.6. Model Reductions for the Independence Working Likelihood with Sex Included ^ 57 Table 3.7. Model Reductions for the First-Order Working Likelihood with Sex Included 59 Table 4.1. Reductions for the Stationary 1-Dependent Working Correlation with Sex Included ^ 74 Table 4.2. Robust z-Scores (Parameter Estimates) for the Binary Final Model ^ 75 Table 4.3. Model Reductions for the Binomial Data ^ ^Table 4.4. Robust z-Score (Parameter Estimate) for the Binomial Final Model vii 81 84 List of Figures Figure 1.1. Feeding Arrays for the Hummingbird Experiment ^ 92 Figure 2.1. Mean Treatment (Transformed) Response versus Measurement for Blocks of 15 ^ 93 Figure 2.2. Mean Treatment (Transformed) Response versus Measurement for Blocks of 10 ^ 93 Figure 2.3. Estimated Variances versus Measurement for Blocks of 15 ^ 94 Figure 2.4. Estimated Variances versus Measurement for Blocks of 10 ^ 94 Figure 3.1. Average Treatment Response in the Logit Scale for Blocks of 3 Trials . . . 95 Figure 3.2 Average Treatment Response in the Logit Scale for Blocks of 5 Trials . . . 96 Figure 3.3. Average Treatment Response in the Logit Scale for Blocks of 10 Trials . . 97 Figure 3.4. Smoothed Logit of Average Treatment Response versus Trial (Window = 10) ^ 98 Figure 3.5. Smoothed Logit of Average Treatment Response versus Trial (Window = 20) ^ 99 Figure 3.6. Smoothed Logit of Average Treatment Response versus Trial (Window = 30) ^ viii 100 Figure 3.7. Smoothed Individual Logits for Treatment 1 ^ 101 Figure 3.8. Smoothed Individual Logits for Treatment 2 ^ 102 Figure 3.9. Smoothed Individual Logits for Treatment 3 ^ 103 Figure 3.10. First Independence Fits without Sex ^ 104 Figure 3.11. Second Independence Fits without Sex ^ 104 Figure 3.12. Fitted Success Probabilities for Treatment 1 from the Independence Working Model ^ Figure 3.13. Independence Pearson Residuals by Bird ^ 105 106 Figure 3.14. Boxplot of Independence Pearson Residuals by Sex ^ 107 Figure 3.15. Smoothed Logit of Average Response versus Trial for Treatment 1 (Window = 60) ^ 108 Figure 3.16. Smoothed Logit of Average Response versus Trial for Treatment 2 (Window = 60) ^ 108 Figure 3.17. Smoothed Logit of Average Response versus Trial for Treatment 3 (Window = 60) ^ 108 Figure 3.18. Binary Independence Fits for Females ^ 109 Figure 3.19. Binary First-Order Fits for Females ^ 109 Figure 3.20. Fitted Success Probabilities for Females in Treatment 1 ^ 110 Figure 4.1. Binary Fits for Females Obtained from the Stationary 1-Dependent Working Correlation ^ 111 Figure 4.2. Average Logit of Response for Treatment 1, Blocks of 10 Trials ^ 112 ix Figure 4.3. Average Logit of Response for Treatment 2, Blocks of 10 Trials ^ 112 Figure 4.4. Average Logit of Response for Treatment 3, Blocks of 10 Trials ^ 112 Figure 4.5. Smoothed Logit of Success Proportion versus Trial (Window = 30) ^113 Figure 4.6. Smoothed Logit of Success Proportion versus Trial (Window = 60) ^113 Figure 4.7. Smoothed Individual Logits (Binomial Data) for Treatment 1 (Window = 60) ^ 114 Figure 4.8. Smoothed Individual Logits (Binomial Data) for Treatment 2 (Window = 60) ^ 114 Figure 4.9. Smoothed Individual Logits (Binomial Data) for Treatment 3 (Window = 60) ^ 114 Figure 4.10. Initial Model Fits for Treatment 1 (Binomial Data) ^ 115 Figure 4.11. Initial Model Fits for Treatment 2 (Binomial Data) ^ 115 Figure 4.12. Initial Model Fits for Treatment 3 (Binomial Data) ^ 115 Figure 4.13. Final Fits for Females from the Stationary 1-Dependent Working Correlation ^ 116 Figure 4.14. Final Fits for Males from the Stationary 1-Dependent Working Correlation116 x Acknowledgement I would like to thank John Petkau for his guidance and supervision throughout the development of this thesis. I would also like to thank Nancy Heckman for her careful reading of the manuscript and subsequent helpful suggestions. Additionally, I am grateful to Gayle Brown for many instructive conversations about hummingbird learning, and, in particular, for Figure 1.1 of the text. Finally, I would like to thank the Statistical Consulting and Research Laboratory for providing the data. xi Chapter 1 Introduction 1.1 Background to the Experiment A desire for a more complete understanding of how organisms relate objects separated in space motivates much research on "spatial learning" in both the biological sciences and psychology. Currently, in the Department of Zoology at U.B.C., two Ph.D. students are studying spatial learning in hummingbirds. One of these students, Gayle Brown, came to the Statistical Consulting and Research Laboratory in the summer of 1990 with questions about how to analyze data she had collected in one of her experiments. The experiment was designed to answer questions about "spatial association learning" in hummingbirds. In "spatial association learning", organisms learn to associate certain cues with reinforcement (such as food, water, or mates) separated in space from the cues. For example, Cori (1989) has shown that bees use the colours of spots on banner petals to select flowers with greater average nectar rewards. Results from previous experiments conducted as part of the ongoing research of the investigator show that hummingbirds learn spatial associations very quickly relative to all other organisms that have been tested in the laboratory (monkeys, rats, pigeons, and children) and even perform perfectly on problems that may be unlearnable by other animals, according to 1 Chapter 1. Introduction^ 2 documentation in the literature (see reviews by MacKintosh, 1983; and Bowe, 1984). One explanation for the success of hummingbirds in spatial association experiments is that, in nature, hummingbirds must frequently learn floral cues in order to efficiently extract nectar hidden within flowers. The fact that many of these floral cues are visibly connected to the nectar of the flower suggests visible guides may facilitate the learning of spatial associations in hummingbirds. To test this possibility, an experiment on the feeding behaviour of eighteen experimentallynaive adult rufous hummingbirds (10 females and 8 males) was carried out. The hummingbirds were assigned randomly by sex to one of three treatments: "No Tape" (NT or treatment 1), in which no visible guides connected cues with the feeders below them; "Partial Tape" (PT or treatment 2), in which fluorescent orange Dymo type (9 mm wide) provided a discontinuous (i.e., broken in two places) connection between each cue and its feeder; and "Full Tape" (FT or treatment 3), in which the visible guide between each cue and its feeder (fluorescent orange Dymo tape) was continuous. The NT group (treatment 1) had 4 females and 2 males, while the PT (treatment 2) and FT (treatment 3) groups each had 3 males and 3 females. The investigator was interested in comparing learning patterns of the birds between these treatments. She expected different learning patterns on each treatment; in particular, she expected hummingbirds in the two treatment groups with visible guides between cues and feeders (2=PT and 3=FT) to become aware of the cue-feeder relationship sooner than hummingbirds in the treatment group with no visible guides (1=NT). From psychological principles of visual perception, she also expected birds in the treatment groups with discontinuous and continous visible guides to become aware of the cue-feeder relationship at roughly the same time. Chapter I. Introduction^ 3 1.2 Description of the Experiment On a given day, three birds (one bird from each treatment group) were tested separately in three randomly assigned rooms, 1.3 x 2.5 x 2.5 metres high. A portion of one end wall in each room contained a horizontal array of six feeders spaced 3 cm apart on a thin metal panel. Figure 1.1 displays the feeding arrays for the NT, PT, and FT groups. Feeders in each feeding array were marked by 19 mm round fluorescent orange Avery labels with central 3 mm holes. Hummingbirds hovered at the feeders and probed their bills through the holes to small food reservoirs. If the feeder was the correct or "profitable" one, a miniature solenoid valve immediately released 2 pl of 20% sucrose solution into the reservoir. Hummingbirds were free to fly and visit feeders at all times during the experiment. Between brief foraging bouts to the feeding array and short non-foraging flights around the room, they spent most of their time on a 1.5 metre high central perch, located 1.8 metres from the feeding array and designed to place birds at eye level with the feeders. A small red light, 4 mm in diameter, protruded slightly through the metal panel 12 cm above each feeder. These small lights served as the spatial cues. During each feeding trial one was lit to signal the feeder below it as profitable. A computer controlled the lights, dispensed the food, and recorded the time and duration of all visits to feeders and perches. To allow the three hummingbirds to become accustomed to the feeding array used in the experiment, on the morning of the day before testing, the birds were placed in training cages identical to their home cages, except that food was available from a wall feeder similar in appearance to the feeding array. The birds quickly learned to feed from these feeders and were Chapter 1. Introduction^ 4 moved from their training cages to the experimental rooms that afternoon. Until experimental training began the following morning, the hummingbirds fed from a standard commercial feeder marked with a fluorescent orange Avery label around the access hole. The feeder was hung in front of a central array feeder on which training would start. The other five feeders were covered. At the beginning of the next day, the standard feeder was removed from each room, exposing a central feeder in the array to the hummingbird in each room. After the birds had fed several times from this feeder, the light above it was lit and training continued. After several more feedings, all other feeders were uncovered and the cue and profitability then reassigned to a different cue-feeder pair. This transfer procedure was repeated until birds had fed several times from each feeder when the light above it was lit. Training lasted approximately 2 hours. Testing started between 10:30-11:00 am, immediately after training was finished. Each hummingbird was tested for 60 feeding trials on each of three successive days. Birds received initial training (as described in the preceding paragraph) on each testing day. On the first and second days, at the end of the day's 60 feeding trials, a commercial feeder was hung in front of a middle feeder in the array, the other feeders in the array were covered, and the birds remained in the experimental rooms overnight. On the third and final day, birds were returned to their home cages immediately after testing was completed. The one profitable feeder and its cue were reassigned randomly among the six feeders each feeding trial, but the same random sequence was used for all birds. Feeding trials began 2 minutes after preceding trials ended, or as soon as preceding trials ended and the birds perched. At the beginning of each feeding trial the following events occurred simultaneously: a soft buzzer sounded for 0.5 seconds, the cue to the profitable feeder was lit, and the profitable Chapter 1. Introduction^ 5 feeder was set to provide 2 ,a1 of solution on each visit. Birds could visit any sequence of feeders and obtain food from the correct feeder up to 11 times, for a total of 22 id each feeding trial. This quantity provided food at a rate slightly more than average ad libitum feeding rates. No limit was imposed on the time for birds to respond in feeding trials. Feeding trials ended only after birds had probed the correct feeder 11 times, or after they had visited at least one feeder and returned to their perch. When a feeding trial ended, the light cue was turned off and the previously profitable feeder delivered no more food. All feeders remained exposed between trials. 1.3 Description of the Data Learning and performance of the hummingbirds was to be evaluated from first visits in feeding trials. If the first visit was to a profitable feeder, the feeding trial was scored as a success (coded as 1); if the first visit was to any one of the five unprofitable feeders, the feeding trial was scored as a failure (coded as 0). The binary data set is comprised of 18 binary series of length 180, corresponding to the successes and failures of each bird on the first feeding foray in a trial, for 3 days of 60 consecutive feeding trials. In addition to the information in these binary series on the first feeding foray of each trial, information on both the number of successful feeding forays and the total number of feeding forays in each trial is available. What we will refer to as the binomial data set consists of the total number of feeding forays (denoted by n) and the number of successful feeding forays (denoted by y) in each trial. Note that the number of successful forays in each trial is not really a binomial random variable because: Chapter 1. Introduction^ 6 • the feeding forays within a trial are not independent, and • even if the forays within a trial were independent, for those trials in which the bird feeds successfully the maximum number of 11 times (roughly 60% of the 3240 trials), the total number of feeding forays made during the trial (denoted by n), rather than the number of successful feeding forays made during the trial (denoted by y), would be the random quantity (a negative binomial random variable) because of the way the experiment has been set up. 1.4 Outline of the Thesis This thesis discusses several approaches to the analysis of binary longitudinal data, the last of which is also applicable to count (and continuous) responses. The first approach discussed is univariate analysis-of-variance for repeated measures using the arcsine square-root variance stabilizing transformation for proportions; this was the approach taken by the investigator in her analysis of the binary hummingbird data. To use this approach, the data for each bird must be blocked over trials before the elements in the resulting response vector are transformed. For example, the response vector could consist of 3 elements, the arcsine square-roots of the proportions of successful trials on each of the 3 days in the experiment. The arcsine square-root transformation may not be entirely effective if the response is a vector of proportions because, although the variances of the proportions are stabilized, covariances remain functions of the success probabilities. Because the investigator expects success probabilities to differ across treatment groups, the covariance matrices are also expected to be different. In this context, univariate analysis-of-variance for repeated measures is of little Chapter I. Introduction^ 7 utility because it requires a common covariance matrix for valid statistical tests. Even if we assume a common covariance matrix is shared by all birds, conservative statistical tests are required (with this data) for checking the effects of measurement and measurement-by-treatment interaction in the analysis-of-variance. Moreover, the sex of each bird is not easily incorporated into the analysis because of the unequal number of males and females in the experiment. An alternate technique for analysis of binary longitudinal data which accounts for correlation over time within subjects and permits flexible modelling of the response would be preferable to this more classical approach. The "working likelihood" approach of Zeger, Liang, and Self (1985), described in Chapter 3, is just such a technique. Here, a working likelihood (not necessarily the true likelihood) is used to generate estimating equations for regression parameters which model the marginal probability of success. Sections 3.1 and 3.2 outline the approach of Zeger, Liang, and Self (1985), while Section 3.3 extends to this approach to time-dependent covariates. This extension permits success probabilities on the first feeding foray to be modelled as functions of time, while taking into account correlation between measurements on each bird. Hence, comparisons of the learning trajectories of hummingbirds across treatment groups are more direct than comparisons achieved by the classical analysis. The remainder of Chapter 3 presents an analysis of the binary hummingbird data using this extension. Note that with only 18 birds in the experiment, efficient estimation of the regression parameters is an important issue. One weakness of this approach is that only one working likelihood is available for explicitly modelling correlation in the data (although, the results could presumably be extended to include working likelihoods modelling correlation in a variety of ways). Therefore, specification of a working likelihood Chapter 1. Introduction^ 8 "close" to the true likelihood (for reasons of efficiency) may be difficult. Further, when modelling correlation explicitly, asymptotic results only hold under direct assumptions on the true correlation structure. Chapter 4 presents the "generalized estimating equation" (GEE) approach of Zeger and Liang (1986) and Liang and Zeger (1986). The GEE approach extends the use of quasi-likelihood to longitudinal data and may be used for any multivariate response with univariate marginal distributions for which the quasi-likelihood formulation is sensible. The GEE approach is more flexible than the working likelihood approach not only because the marginal likelihoods need not be specified, but also because many possible "working correlations" are available for explicitly modelling correlation in the responses. Moreover, no assumptions about the true correlation structure are necessary for valid asymptotic results with any of these working correlations. To illustrate the flexibility of the GEE approach, both the binary and binomial hummingbird data are analyzed using different working correlation structures. In Chapter 5, the thesis concludes with a summary of the important conclusions made resulting from the analyses and a comparison of the methods discussed. Chapter 2 Classical Analysis A standard approach to the analysis of binary longitudinal data such as the hummingbird data (for the first feeding foray) might be to calculate the proportion of correct responses over blocks of an equal number of successive time points for each subject, and then "normalize" these proportions via the arcsine square-root transformation. The resulting subject response vectors within the ith treatment group might then be assumed to be approximately normally distributed with mean vector and covariance matrix Ei. Moreover, if there is a physical basis for doing so, or if there is insufficient evidence in the data to conclude otherwise, the covariance matrices may be assumed to be equal, a sufficient assumption for the use of multivariate analysis-of-variance (MANOVA) of the response vectors. Before using this approach, the equality of the covariance matrices should be tested. Assuming the covariance matrices of birds within each treatment group are the same, the estimated covariance matrix common to each treatment group will be non-singular only if q, the number of elements in each subject's response vector, is fewer than one less the number of subjects in the treatment group. The hummingbird experiment has only 6 birds per treatment group, so q must be 4 or less for these estimated covariance matrices to be non-singular. The assumption of equal covariance matrices cannot be formally checked if q > 5 because the test requires the estimated covariance matrices to be invertible. To check the equality of covariance matrices, 9 Chapter 2. Classical Analysis^ 10 subject data must therefore be blocked over at least 45 trials (leading to q = 4 measurements per subject). Aggregating the data of each bird over 45 trials is awkward, though, because there are 60 trials per day and day effects could be important in describing the performance of the birds. Instead, individual data might be blocked over 60 trials (one day), allowing MANOVA with q = 3 "measurements" on each bird. But with only 3 measurements (corresponding to different days) on each bird, such an analysis is likely to be too crude to capture the true nature of the temporal learning pattern on each treatment. The aggressive blocking necessary to carry out MANOVA with only 3 measurements per subject does not permit detailed inference on how success probabilities within the treatment groups change as a function of time. An alternative method of analysis able to provide higher resolution of the temporal learning patterns would be ideal. In the MANOVA context, a step towards this ideal would be to block over fewer trials. If data were blocked over 30 trials instead of 60, there would be 6 measurements on each subject, and the previous confounding of blocks and days would be eliminated, allowing for a more detailed representation of temporal learning patterns. But the resulting MANOVA tests would not be valid unless all treatment groups shared a common covariance matrix; that is, E i = E, i = 1,2,3. This assumption of common covariance matrices across treatments is inconsistent with the investigators beliefs and contradicted by the data.' Moreover, if the 1 The data for q = 3 measurements per bird (where each measurement is the arcsine square-root of a raw success proportion) indicate that covariance matrices for the transformed success proportions are not equal across treatments. Note that if these covariance matrices are not equal when there are q = 3 measurements per bird, it is unlikely that they will be equal when less aggressive blocking results in q > 3 such measurements. Covariance matrices can be estimated within treatment groups when there are 3 or fewer measurements per bird, but the distribution of the test statistic for their equality is only chi-squared as the number of subjects in each treatment group (6 for this experiment) becomes large. Moreover, the test is very sensitive to non-normality ^ Chapter 2. Classical Analysis^ 11 success probabilities for birds differ across treatment groups, the covariance matrices for the transformed response must differ across treatment groups as well because, even though the arcsine square-root transformation can be expected to stabilize the variances, the covariances will still be functions of the success probabilities. Although this analysis would provide a more detailed picture of temporal learning patterns, its conclusions would be somewhat tenuous because the assumption of a common covariance matrix is supported by neither the investigator's beliefs, the data, nor the nature of the variance-stabilizing transformation. Even though assuming a common covariance matrix may not be correct, it allows blocking over as few as 15 trials. The resulting 12 measures per bird permit representations of temporal response patterns across treatments that are much more detailed than the representations obtained when the equality of covariance matrices can be tested (and we are limited to 3 or fewer measurements per bird). Another benefit of assuming a common covariance matrix is that univariate analysis-of-variance for repeated measures, a simpler method more accessible to researchers familiar with experimental design, may be used as an alternative to MANOVA. The details of such analyses are presented in the remainder of this chapter. in the data (see Morrison 1976, page 252). Hence, both the test results and the estimated covariance matrices should be examined to judge whether underlying covariance matrices appear to be approximately equal. The estimated covariance matrices (common to birds within treatment groups) for the case of q = 3 measurements per bird (where the measurements are transformed success proportions based on blocks of 60 successive trials) are, 2.5 3.7 2.7^2.2^1.1^-0.1^5.5 3.7^1.1 3.7 5.9 4.2 i , .01 [ -0.1 1.1^2.1^0.3 I , .01 [ 3.7 9.8^7.4 I , .01 [ 2.7 4.2 8.1 ^0.3^3.0 1.1^7.4^11.0 for treatments 1, 2, and 3 respectively. These appear quite different, and the same conclusion is suggested by the chi-squared test; the value of the test statistic is 44.1 on 12 degrees of freedom (p < .001, based on the asymptotic chi-squared distribution). It thus seems unlikely that the covariance matrices are equal when there are q > 3 measurements per subject. 12 Chapter 2. Classical Analysis^ 2.1 Univariate Analysis-of-Variance for Repeated Measures Univariate analysis-of-variance for repeated measures is strictly valid only if observations are multinormally distributed with a common covariance structure corresponding to the "type-H" pattern described by Huynh and Feldt (1970); i.e., E = (ak,/), has the pattern defined by ork,1 = ak + ai + A k = 1. a k,1 = ak + cei k 1. If the covariance structure of the observations does not correspond to the type-H pattern, Greenhouse and Geisser (1959) have proposed conservative tests based on adjusting the degrees of freedom for the univariate tests. Even though the univariate approach may sometimes require use of these conservative tests, it still has a strong appeal because it is easy to interpret and commonly used by those familiar with experimental design. Further, unlike multivariate analysis-of-variance, test statistics are still defined if the number of repeated measurements per subject is greater than the within treatment degrees of freedom. The binary hummingbird data may thus be blocked over even fewer than 15 (say 10, for example) trials to perform the univariate analysis, yielding representations of temporal response patterns that are more detailed than those permitted by MANOVA. However, if the binary hummingbird data is blocked over fewer than 15 trials, it will not be possible to check for the type-H pattern in the underlying covariance matrix because the estimated covariance matrix will once again be singular. Even so, an analysis using conservative univariate tests may still be carried out. Because the investigator wishes to understand how success probabilities for birds in different treatment groups compare over time, detailed representations of temporal response patterns are Chapter 2. Classical Analysis^ 13 desirable. Since these are permitted by univariate analysis-of-variance for repeated measures and since this method of analysis is well understood and commonly used by researchers familiar with experimental design, this technique will be used to analyze the binary hummingbird data (on the first feeding foray), even though a key assumption for this technique (i.e., that a common covariance matrix for the transformed response vectors be shared by all subjects) is not supported by either the beliefs of the investigator, the data, or the nature of the arcsine square-root variance stabilizing transformation for binary data. 2.2 The Statistical Model In order to "normalize" the hummingbird data, the proportion of correct responses for blocks of fifteen and then ten successive trials was calculated for each bird and transformed via the arcsine square-root transformation. The resulting 12- and 18-dimensional response vectors for each bird in the i th treatment group were then assumed to be approximately normally distributed with mean pi and common covariance matrix E. The treatment means of these transformed proportions are shown in Figures 2.1 and 2.2. These figures suggest that the (transformed) response grows roughly linearly with time, particularly for birds in treatments 2 (PT) and 3 (FT). The lack of systematic discontinuities in these response curves at measurements corresponding to new days (measurements 7 and 13 in Figure 2.1, and measurements 5 and 9 in Figure 2.2) suggests that modelling a "day" effect is not necessary. Thus, for both blocks of 10 and blocks of 15 observations, a univariate analyis-of-variance for repeated measures was carried out under the following model: 14 Chapter 2. Classical Analysis^ Yijk = arcsine(V _Ask) = It -I- treatment= + measurements + (treatment * measurement)is subject(ok + (measurement * subject)(osk Cijk, for i = 1, 2,3; j = 1, ..,12 or 18 (for blocks of 15 and 10 observations respectively); k = 1, .., 6. Here Pijk is the raw success proportion for the kth subject within the i th treatment at the i th measurement, and ciik is the noise term for the kth subject within the i th treatment at the i th measurement. We assume that (Eijk )j=1,..,12 or 18 (i = 1,2,3; k = 1, ..,6), the subject noise vectors, are independent normally distributed random vectors with mean 0 and common covariance matrix E. "Treatment" and "measurement" are fixed effect factors, and "subject" is a random effect factor nested within "treatment". 2.3 Results of the Analysis Before carrying out univariate analysis-of-variance on the binary hummingbird data (for the first feeding foray) blocked over 15 trials, the underlying covariance matrix was tested for the type-H pattern. Because the value of the appropriate chi-squared statistic (see Morrison 1976, page 251) is 259.17 with 65 degrees of freedom, the covariance matrix does not seem to have the type-H pattern. The estimated covariance matrix for the binary hummingbird data blocked over 10 trials is singular because there are 18 measurements per bird and only 18 birds in the study, so the underlying covariance matrix cannot be tested for the type-H pattern. However, it seems unlikely that the covariance matrix for blocks of 10 trials will have the type-H pattern since the covariance matrix for blocks of 15 trials does not. Univariate analysis-of-variance was Chapter 2. Classical Analysis^ 15 therefore carried out using the conservative tests of Greenhouse and Geisser (1959) for both blocks of 15 and blocks of 10 observations. For the univariate analysis-of-variance tests to be strictly valid, the transformed response vectors for the individual subjects must be normally distributed with a common covariance matrix E; that is, (Yi303=1,..,120, 18 ", N(µi, E). The investigator does not believe that a common covariance matrix is shared by all birds, but separate covariance matrices common to birds within treatment groups cannot be modelled due to the small number of birds within each treatment group. Because these estimated covariance matrices, Ei (i = 1, 2,3), are singular, a formal test for the equality of the underlying covariance matrices is not available. However, elements of the estimated covariance matrices may still be compared visually. Figures 2.3 and 2.4 plot the estimated variances of the transformed responses (the diagonal elements of Ei) against measurement by treatment group for the 12- and 18-measurement analyses, respectively. Neither figure suggests an obvious pattern in the treatment variances over time, but both figures suggest that the variances for birds in treatment 3 are generally greater than those for birds in other treatments, while the variances for birds in treatment 2 are generally lower than those for birds in other treatments (the estimated variances of responses for birds in treatment 3 are anywhere from 1-11 times larger than the estimated variances for birds in treatment 2 for the 12-measurement analysis, and anywhere from 0.4-22 times larger for the 18-measurement analysis). These results suggest that the assumption of a common covariance matrix is not valid for either the 12- or the 18-measurement analysis. Hence, the distributions of test statistics formulated assuming equal covariance matrices can only be expected to roughly approximate the true distributions of the test statistics. In spite of this, the tests should still provide Chapter 2. Classical Analysis^ 16 useful information about the relative magnitude of the components of variability, and strong conclusions from the analyses-of-variance should be consistent with conclusions from analyses which permit unequal covariance matrices. The analyses-of-variance are given in Tables 2.1 (18-measurement analysis) and 2.2 (12measurement analysis). Tables 2.1 and 2.2 confirm one clear indication from Figures 2.1 and 2.2: "measurement" is an important factor. But Tables 2.1 and 2.2 give no strong evidence of "treatment" effects, either through main effects or through interactions with "measurement". On the other hand, Figures 2.1 and 2.2 clearly indicate that the (transformed) responses for treatment 3 are generally greater than those in the other treatments, suggesting that "treatment" main effects do, in fact, exist but are not being detected because there is too much variability between subjects within treatment groups. Table 2.1. Analysis of Variance for Blocks of 10 Source measurement treatment subject (within treatment) measurement *treatment subject*measurement (within treatment) total SS 14.74 .57 5.67 1.26 15.33 37.57 d.f. 17 2 15 34 255 323 MS F conservative d.f. .87 .29 .38 .037 .060 14.20 .75 1,15 p .002 .49 .62 2,15 .55 Table 2.2. Analysis of Variance for Blocks of 15 Source measurement treatment subject (within treatment) measurement*treatment subject*measurement (within treatment) total SS 3.56 1.70 7.22 .958 6.15 19.59 d.f. 11 2 15 22 165 215 MS .32 .85 .48 .044 .037 F 8.68 1.76 conservative d.f. 1,15 p .01 .21 1.17 2,15 .34 The conservative F-tests used in these univariate analyses-of-variance are based on approximate distributions of the sums of squares when the subject response vectors are multivariate Chapter 2. Classical Analysis^ 17 normal with a common covariance matrix E. Greenhouse and Geisser (1959) have shown that when the hypothesis of no measurement effects is true, the F-statistic for measurement has an approximate F-distribution with degrees of freedom (K — 1)c and (K —1)(N — K)c, where K denotes the number of measurements per subject, and N denotes the number of subjects. The scale factor K2(0 _ F.12 E — — 1)(E i E j a — 2K a72 K2Tr 72) is computed from the elements of the (common) covariance matrix E, where o71 is the mean of the p diagonal terms, cr: is the grand mean of all variances and covariances, and c)7 is the mean of the elements in the i th row of E. They have also shown that when the hypothesis of no measurement-by-treatment interaction is true, the F-statistic for this interaction has an approximate F-distribution with degrees freedom (K —1)(p-1)c and (K — 1)(N — p)€, where p denotes the number of treatments. In practice, E is never known and estimating c from the sample covariance matrix is undesirable because more uncertainty is then introduced into the approximate analysis-of-variance. Greenhouse and Geisser have shown that c > 1/(K — 1), implying that the degrees freedom for the approximate tests cannot be less than 1 and (N — K), and (K — 1) and (N — K), respectively. A conservative F-test results from assuming the maximum possible reduction in degrees freedom, and therefore requires larger values of the F-statistic to reject the null hypothesis (i.e, the hypothesis of no measurement effects, or the hypothesis of no treatment-by-measurement interactions) than tests which assume less than the maximum reduction in degrees freedom. These conservative F-tests are used in both the 12- and 18-measurement analyses described Chapter 2. Classical Analysis^ 18 above, but the conclusions from these analyses are the same as the conclusions from the analyses which estimate the scale factor c from the data (e is .33 for the 12-measurement analysis, and .24 for the 18-measurement analysis, so the use of e results in, respectively, degrees of freedom 3.63 and 4.08 times the conservative degrees of freedom). 2.4 Conclusions and Limitations In both univariate analyses, the conservative test for the measurement main effects finds "measurement" to be an important factor. Figures 2.1 and 2.2 also indicate that "measurement" is an important factor. The test (not conservative) for treatment main effects in both univariate analyses-of-variance does not find "treatment" to be important, nor does the conservative test for the interactions of treatment with measurement find "treatment*measurement" to be important. Thus, both analyses-of-variance indicate that treatment group has no effect either directly, or through interaction with measurement, on the (transformed) response pattern. On the other hand, plots of the fitted values (Figures 2.1 and 2.2) clearly indicate that the (transformed) responses for treatment 3 are greater than those in the other treatments, suggesting that a treatment main effect is present but cannot be detected by these methods because variability between subjects within treatments is too large. The same result would therefore be expected in the corresponding MANOVA. Some of the within-treatment variability between subjects may be explained by their sex, but modelling of a sex effect or any interactions of sex with the other factors is difficult in both univariate and multivariate analysis-of-variance because of the unequal number of males and females in treatment 1 (NT). An approach allowing modelling of sex, despite the imbalanced number of males and females across treatment groups, Chapter 2. Classical Analysis^ 19 might reduce the unexplained within-treatment variability, and thereby provide more sensitive tests of all effects. Note that many other "classical" approaches to analyzing the binary hummingbird data (for the first feeding foray) have not been discussed in this section. These include parametric modelling of the "measurement" effects (for example, analysis-of-covariance, which requires that observations within each subject be independent), and analysis-of-variance (either univariate or multivariate) in the logit scale (logit(p) = log( TF_T,), p between 0 and 1) rather than in the arcsine square-root scale to correspond to logistic regression. Although univariate and multivariate analysis-of-variance are standard techniques used to analyze data in which repeated measurements are taken for each subject, both seem to model how the data change over time in indirect ways which lack intuitive appeal. To allow changes in response over time, univariate analysis-of-variance includes "measurement" main effects and interactions in its model. Multivariate analysis-of-variance allows changes in response over time through the different elements of its treatment mean response vectors. More intuitive approaches, presented in the next sections, model treatment response patterns directly as a function of time, while taking into account the correlation of measurements within subjects. Treatment response patterns can then be compared visually by examining the fitted treatment response models. Chapter 3 Analysis Based on a Working Likelihood Approach 3.1 Logistic Regression for Binary Repeated Measures A natural approach to the analysis of the binary hummingbird data is to extend logistic regression to the case where the binary outcome variable is observed repeatedly for each subject. Zeger, Liang, and Self (1985) present the details of one way of implementing such an approach in which, rather than the actual likelihood, a "working likelihood" based on simple assumptions about the time dependence within each subject's data is used to generate estimating equations which lead to consistent estimates of the regression parameters, as well as standard errors for these estimates. Not only does this approach allow for correlation over time within subjects, but it allows flexible modelling of the response and is powerful enough (provided the number of subjects is large) to detect important factors and covariates, while remaining simple and easy to interpret. 3.2 The Working Likelihood Approach To establish notation, let Yi, t (t = 1, .., n t ) be a stationary binary time series, and let Z, be an sxl vector of time-independent covariates for the it def i th subject (i = 1, K) . For = Pr(Yi, t = 1 I Zi) , it is assumed throughout this section that logit in = ZP. The 20 Chapter 3. Analysis Based on a Working Likelihood Approach^ primary objective is to make inference about the vector of parameters, p. 21 There may be additional nuisance parameters in the working likelihood, so let 9 denote all the parameters modelled by the working likelihood. Let S(9) be the working likelihood score function, and let H(0) be the corresponding Hessian matrix. Finally, let 0 be the estimate of 0 obtained by maximizing the working likelihood; in general, we have S(9) = 0. Throughout, score functions and Hessian matrices will be evaluated under the working likelihood, E(.) will denote the expectation under the true (unknown) likelihood and Ew(.) will denote the expectation under the working likelihood. Motivation for the results of Zeger, Liang, and Self (1985) is provided by the Taylor expansion for S(9) about 0, and the Mean Value Theorem which, together with the fact that S(0) = 0, allow us to write S(9) = SA+ H(0.)(0 — 0) = —H(0,0(9 — 9), where 0,, is between 0 and O. If the data for the different subjects is independent, then the working log-likelihood is the sum of the independent log-likelihoods for the different subjects. Thus, the score function and the Hessian are also sums of independent random quantities, and we can write, - 0 =^S(9) = -[- Hi(0.)]-1 - E Si(0), i=1^K i=1 where Si and Hi are the i th subject's score function and Hessian, respectively. If E Si(0) 0, then provided regularity conditions for Lindeberg's Central Limit Theorem hold, under the true model, ^1 EK Si(0) =^N(0, K2 i=1 Chapter 3. Analysis Based on a Working Likelihood Approach^ 22 where /*(0) = E S(0)S'(0). Under regularity conditions ensuring B is consistent for 9, we have K^ 1 K^ 1 urn --E E Hi(0) = lim —E H(0). — EHi 0.) K --+ oo K K" K i=1 K 1=1 ( .4 00 Thus, by Slutsky's Theorem, — 9 N( 0, [E H(0)] -1 I*(0)[E H(0)] -1 ), where all expectations are taken under the true (unknown) likelihood, and convergence is under the true (unknown) distribution. Note that to use this result for inference about 0, both /*(0), the variance-covariance matrix for the score function under the true (unknown) model, and E H(0), the expectation of the Hessian under the true (unknown) model, must be estimated from the data. Provided 0 is consistent for 0, I*(0) = E S(0)S'(0) can be consistently estimated by EL Si(d)S(d). The expected Hessian under the working model is identical to the expected Hessian under the true (unknown) model provided simple assumptions about the time dependence within each subject's data are valid. Therefore, if d is consistent for 0, EwH(0) is a consistent estimate of E H(9). Because H(0) is a sum of independent random quantities, if 0 is consistent for 0 then an alternative consistent estimate of E H(0) is H(0), the observed Hessian evaluated at O. The extension of the Zeger, Liang, and Self approach to the case of time-dependent covariates, presented in the next section, uses the observed Hessian evaluated at d, rather than EwH(0), to (consistently) estimate E H(9). Zeger, Liang, and Self (1985) present two estimators of 0. The first, ijo , is obtained from the working assumption that repeated observations for a subject are independent of one another. Chapter 3. Analysis Based on a Working Likelihood Approach^ 23 This estimator is consistent (as K^oo) for any set of stationary binary time series for which logit ii = 4/3. The second estimator, fji , is obtained from the weaker working assumption that each time series has a stationary 1-dependent correlation structure. A . is consistent for any set of stationary binary time series for which both logit = 40 and the first lag autocorrelation, p = is common to all subjects. For the independence working model with logit 7ri = 40, the working log-likelihood can be written as K ni ^1 0(0) = E E [ yi,t^+ log(1 — ri )] 1=1 t=1 where logit irY = 4[3. Let S(/3) be the score function for the vector /3, and let Sti (13) be its U th element (u = 1, s). The maximum likelihood equations under the working model are K n, 0 = su(i3) = E E(.Yi>i — zi)Ziu, u = 1, ^s. 1=1 t=1 Define A) to be the solution to these estimating equations. Note that the expected values of the scores under the true (unknown) model are 0 provided the mean structure has been correctly specified; that is, E(1 7i, t Zi) = [1 + exp(-40)] -1 , or logit E(Yi, t 1Z1) = Z'i 13. This ensures that ijo is an asymptotically unbiased estimate of /3. The (u, v)th element of the Hessian matrix is given by K n, Iluv(0) = EE i=1 t=1 ziuZivri(1 — which does not depend upon the Yi, t 's. Therefore, the expected Hessian is the same under both the true likelihood and the working likelihood; in other words, the working likelihood Fisher information matrix is the same as — E H(/3) in this case. As K —> oo, we obtain the result that Chapter 3. Analysis Based on a Working Likelihood Approach^ 24 [40 is consistent and asymptotically normal for any set of stationary binary time series such that logit 7i = 40, as stated by Zeger, Liang, and Self (1985) in their: PROPOSITION 2.1. Let Yi, t (t = 1, ..,ni < oo; i = 1,..,K) be stationary binary series such that logit E(Yi, t I Zi) = 43. Then /Jo — is asymptotically multivariate Gaussian with expectation 0 and covariance matrix V0 (/3) = -1O(0) -1 -1;(0)Io(0) -1 , where MO) is the working Fisher information matrix and PAP) = E S(/3)S'(/3). A consistent estimate of /0 (0), the Fisher information matrix under the working model, is given by fo 10(Qo). A consistent estimate of ^the variance-covariance matrix of the score function, is provided by E sioosmc,), where Si(.) is the score function for the i th subject. Combining these results leads to a consistent estimate of V0 (0) given by fo —1 f —1 o The estimator^is obtained similarly, after replacing the working assumption of independence by the working assumption of a common lag one (and no other) autocorrelation. Let p= COTT(Yi,t,Yi,t-1) and 0 = (p, If pi t def E (Yi, t Zi), the true lag one conditional mean, the assumption of a common first-lag autocorrelation implies that ^ Chapter 3. Analysis Based on a Working Likelihood Approach^ 25 Pit =^POri,t-1^= Ew (Yi,t I Yi,t—i, Zi). The working log-likelihood can thus be written ni /AO) = E [yi , 1 4,3 + log(1 — 70+ E[yi, t log pi t + (1 — yi, t )log(1 — NA], i=1^ t=2 where logit ri = 40 as before, and pit =^p(yi,t-i - 70• For this working model, let S0 (0) be the score function for p, and let Su (0) be the score function for i3„ (u = 1, s). The maximum likelihood equations under the working model are K 0 = So(9) = ,Q7 °P ^= ni E E (yi,t -^- ^2=1 ri) t=2^PitPit ^Ziori7T(Yi,t — Pit) 0 = Su (0) =^= E{Ziu(Yi,1 —^+ ( 1 P) upu^i=1^ t=2^PitPit }, where 73 87 = 1 - pi t , and IT = 1 — ri. Define B = (ji„di) to be the solution to these estimating - equations. Note that the expected values of the scores under the true (unknown) model are 0 provided the lag one conditional mean structure has been correctly specified; that is, E (Yi, t Yi,t _ i , Zi) = P(Yi,t-1 — ri). This can be seen by conditioning on Yi,t_i for each Yi, t in the estimating equations. Thus, 0 1 is an asymptotically unbiased estimate of /3. The elements of the Hessian matrix are given by: K n, E E^(yit 2/ ap^ 21 = - pit) ,^2pi d ,OPit )2 , j^I i=1 t=11-PiiPit^(Pit.7)2 K =E i=1 —ri3f7ZiuZiv i + n^(Yi,t t=2 l — [pit1pit + ((Ypitit— ivP)i2t)(1 2piol ( aaPoiu.t )( aaPiLt Pit) ( 02 Pit ) PitPit 0/3,,0/3v and 021 1 00 0. K ni^ { EE \ n21).t (K t — Pit)( " 1-2 i =1 t =1^PitTi 0P 0 /3u ) [ 1 + (Yit Pit) )9 (1^2Pit )1 LPitPit^(PitPit ( )( °Pit )} 0/3u^OP )1 Chapter 3. Analysis Based on a Working Likelihood Approach ^ 26 where °Pit op -= Yi,t-1 — °Pi t = (1 — p)riTr7Ziu, 490u 82 Pit = (1 — p)(1 — 271 - i)ri3t7ZiuZiv, afivaot, and, 0 2 pit ^ = ri7r7ziu. apafi. The assumption that, under the true likelihood, a common lag one autocorrelation is shared by all subjects, that is, that the lag one conditional mean structure has been correctly specified, yields E(Yi,t Zi) = 7ri P(Yi,t-i - ri), without any further assumptions concerning the higher lag autocorrelations. Under this assumption E(Yi,t ^Zi) = Ew(Yi,t and it follows that E H(9) = Ew H(0); in other words, under the assumption that the lag one conditional mean structure has been correctly specified, the working likelihood Fisher information matrix and —E H(0) are identical. The Fisher information matrix under the first-order working model, / 1 (B), can be partitioned as I1=[ Ipp Ina 10P 100 where the (u, v)th element of 100 is + P }, E ziuzivri7Tfi + (ni — 1)(1 W ilri(1 3P) (1 — p)2 71- 07-7 p the u th element of /o p is ( 2 r; — 1) P E(ni (1— p)27rjr7 p, — Chapter 3. Analysis Based on a Working Likelihood Approach^ 27 and , \ K^ Inv ^( 1 ri 7ri ^1)( P)(n 2^(1 — P) 2 roTT p ) As K^oo, we obtain the result that B is consistent and asymptotically normal for any set of stationary binary time series such that logit E(Yi,t j Zi) = 4 13 and corr(li,t, 17i,t--1) = p, as stated by Zeger, Liang, and Self (1985) in their: PROPOSITION 2.2.^Let Yi,t (t = 1,^< oo; i = 1, K) be stationary binary series such that logit E(Yi, t Zi) = Zp and corr(li,t,li,t-1) = p. Let 0 (p, 0). Then 0 — 0 is asymptotically multivariate Gaussian with expectation 0 and covariance matrix V 1 (0) = Ii (0) -1 I1(0)I1 (0) -1 , where 11 (0) is the Fisher information matrix under the first-order working model and IV O) = E S(0)S'(0). A consistent estimate of h (0), the Fisher information matrix under the first-order working model, is given by A consistent estimate of /17(0), the variance-covariance matrix of the score function, is provided by _ E si(o)s:(o), i=i where Si(• is the score function for the i th subject. Combining these results leads to a consistent estimate of V1 (0) given by 1 -1^-1 V1 =I Ii I1 A technical point which might arise during the Newton-Raphson iteration to obtain 0 is that, depending on the values of the ri's, it may be necessary to constrain p to a smaller region Chapter 3. Analysis Based on a Working Likelihood Approach^ 28 than (-1, 1) in order to prevent the pit = Ew (Yi,t I^Zi) = ai +^— ir) = 1, K) from taking on values outside the range [0, 1]. Since the covariates, Zi, are time independent, the success probability ri = Pr(Yi, t = 1 Zi), does not depend on time. Zeger, Liang, and Self (1985) were concerned with short binary time series, in which there is no need for time-dependent success probabilities. However, the binary time series for the first feeding forays of the hummingbirds are long (180 time points), and the success probabilities seem to increase with time. Figures 3.1, 3.2, and 3.3 plot the transformed average response for birds within a treatment group over consecutive blocks of 3, 5, and 10 time points respectively. Note that, except for the vertical scale, Figure 3.3, which plots the average response in the logit scale, and Figure 2.2, which plots the average response in the arcsine square root scale used for the classical analysis, are identical. All the plots clearly indicate that success probabilities generally increase with time. In the next section, the results of this section are extended to the case of time-dependent covariates to allow the modelling of success probabilities which change with time. 3.3 Incorporating Time-dependent Covariates 3.3.1 The Independence Working Model Suppose we start with the working assumption that repeated observations for a subject are independent, but now let Zi, t be an sxl vector of time-dependent covariates for the (i = 1, .., K). If ri,t Pr(Yi,t = 1 I Zi, t ), it is assumed throughout this section that i th subject ^ Chapter 3. Analysis Based on a Working Likelihood Approach^ logit 7ri t = 29 44 0. The estimator )40 is obtained from the working log-likelihood K ni lo(Q) = E E i=1 [yiA,03 + log(1 — ri,t)] where logit 7ri t = 4 4 13. , Let S u ( i3) be the score function for O u (u = 1, ..,^) under this working likelihood. Then the maximum likelihood equations under the working model are K ni , m^ o 0 = Su (13) = = E E(yi t — i=1 t.i Define go to be the solution to these estimating equations. Note that the expected values of the scores under the true (unknown) model are 0 provided the mean structure has been correctly specified; that is, 7r,, t^E(Yi,t Zi,t) = [1 + exp(—Z2/3)] -1 , or logit E(Yi,t I Zi,t) = 4, t ,(3. Therefore, go is asymptotically unbiased for 0. The (u, v)th element of H o (/3), the Hessian matrix for this independence working model, is given by K n, — E E^— i=1 t=1 which does not depend upon the Yi, t 's. Therefore, the expected Hessian is the same under both the true likelihood and the working likelihood; in other words, the working likelihood Fisher information matrix is the same as —E Ho (/3) in this case. Under regularity conditions sufficient to guarantee the consistency of g o , the following generalization of Zeger, Liang, and Self's (1985) Proposition 2.1 is immediate: PROPOSITION 1. Let Yi, t (t = 1 ..,ni < oo; i = 1,..,K) be binary series such that , logitE(Yo I Zi,t) = 4,0. Then &— is asymptotically multivariate Gaussian with expectation Chapter 3. Analysis Based on a Working Likelihood Approach^ 30 0 and covariance matrix vo(0) = [E H0(0)] -1 1O(0) [E 110(0)] -1 , where I8`(0)= E S(0)SV). A consistent estimate of E HO), the expectation of the Hessian under the true (unknown) model, is given by Ho (00 ). A consistent estimate of /8`(0), the variance-covariance matrix of the score function, is provided by K Esi(go)sago), i=i where Si(.) is the score function for the i th subject. Combining these results leads to a consistent estimate of Vo(0) given by Vo = [Ho(go)] - 1 Io [Ho(g0)] -1 • 3.3.2 The First-Order Working Model The first-order working model with time-dependent covariates assumes each binary series has a stationary 1-dependent correlation structure with logit 7ri t = Zii,t 0 and , corr( 170,Yi,t-i.) = P. If Pit 41 E (Yi,t I Yi,t-i, Zi), the assumption of a common lag one autocorrelation implies that Pit^p 71i,t7ri,t ^ (yi,t-i - 70-1 71- 0-1 where lri ,t = ( 1 — ri,t)• Chapter 3. Analysis Based on a Working Likelihood Approach^ 31 Let 0 = (p,/3). Then the working log-likelihood is ni /1( 0 ) = E [yoz:,io + log(1 —^E[yi,tiog Pit + (1 — yi,t)log(1— t=2 pit)]] where i ,t Pit = 711 ,t^P^ — 7ri,t-1)• 7ri,t—i 7ri,t—i For this working model, let So(0) be the score function for p, and let S.(0) be the score function for ou . The maximum likelihood equations under the working model are all K EE ° it (Yi t - Pit) P L'P^i=1 t=2 PitPit^OP = So(0) =^= , and o = su(e = ) where pi t ni vK^ Pit 'Pit na =^—^E Yi't uP.^ i=1^t=2 PitF ift 0,3 14 all , ( ) = 1 - pi t . The partial derivatives appearing in these equations are given by: t^ Ori,t7ri,t^— afitt +22 7ri,tiri,t — 7ri,t-1] [(1 — 27ri,t)Zi,t,u — (1 — 2 7ri,t—i)Zi,t—i,u], and °pit Op Define 0 ^ 7ri,t7ri,t [Yi,t—i — ri,t-1] • 7ri,t-1lri,t-1 =^A.) to be the solution to these estimating equations. Note that the expected values of the scores under the true (unknown) model are 0 provided the lag one conditional mean structure has been correctly specified; that is, E(Yi,t I^Zi,t) = 7r • tri t 7r^ 2,t p^t-1 — ir i ,t). This can be seen by conditioning on Yi,t_i for each 11 4 in the estimating equations. Therefore, B is asymptotically unbiased for 0. The Hessian matrix for the first-order working likelihood, H 1 (0), can be partitioned as Chapter 3. Analysis Based on a Working Likelihood Approach ^ H1(9) = Hpp P[ HapII 00O ], where the (u, v)th element of Ho is 0 2 11 K 00.00v i=1 ni 1 (yi,t^-^\i0Pit 0Pit^(yi,t — Pit) 0 2 Pit I Pit) f + t=2 E^aou +^ --= LPit) PitPit (PitPit)2 oov^Pit ^a Oua Ov [{ — the u th element of Ho p is 921 1 3 an — 4 )] ^2pit (prct) it - Pi2t) (1 1^(Yi,t [ pit.7 ^00u u i=1 t=2 [ °pit °pit^(yi,t - pit) 0 2Pit OP^( 0 Oua P ) and 0 2 /1^K IIP P 0 13 2 ni 1 + „ (1 2p it )] _^ ^ = EE i=1 t=2 [ PitPit^kiittPitr^Op (Yi t Pit)^aPit ) 2^(Yi,t — Pit) 0 2 Pit pOW 0,9 2 The partial derivatives appearing in these expressions are given by: = (1 — + 2^' P 2 tri " t-i ri,t -1 (1 — 27i,t) 1r ,t 1 ri,t ^ [Yi,t -1 where ait^.5[(1 — 27ri,t)Zi,t, u — (1 —^— 27ri,t)Zi, t ,„ — (1 — — 2 [ 7ri,tri,tZi,t,uZi,t 7 v —^Zi,t-LuZi,t-1,v] 02, " ^.5 aduaP — r , —^27ri,t)Zi,t,u — (1 — 27i,t-i)Zi,t-1 04.] 32 Chapter 3. Analysis Based on a Working Likelihood Approach^ 33 and O2Pit op 2^• Under regularity conditions sufficient to guarantee the consistency of 0, the following generalization of Zeger, Liang, and Self's (1985) Proposition 2.2 is immediate: PROPOSITION 2.^Let Yi,t (t = 1, ..,ni < oo; i = 1,..,K) be binary series such that logit E(Yi,t I Zi,t) = 4 4 0 and corr(Yi,t,Yi,t—i) = p. Let 0 = (1)0(3 ). Then 0 - 0 is asymp- totically multivariate Gaussian with expectation 0 and covariance matrix V1 (0) = [E H1(0)] -1^111(9)] where Ip(0) = E S(0),5"(0). A consistent estimate of E H 1 (0), the expectation of the Hessian under the true (unknown) model, is given by H 1 (0). A consistent estimate of IV O), the variance-covariance matrix of the score function, is given by Ii =E Si(0)s:(0), i=i where S1(.) is the score function for the i th subject. Combining these results leads to a consistent estimate of V1 (0) given by Vi =^( 0 ' Ii [H1A -1 . )] Because this extension to the Zeger, Liang, and Self approach incorporates time dependent covariates, it provides the necessary tools for the analysis of the binary hummingbird data. The extension presented in this section permits the success probabilities to be modelled as a function of time, while taking into account the correlation between measurements on each Chapter 3. Analysis Based on a Working Likelihood Approach^ 34 bird, and thereby allows comparisons of temporal response patterns across treatments which are more direct than the comparisons achieved by either univariate analysis-of-variance for repeated measures or MANOVA. 3.4 Exploring the Binary Hummingbird Data Prior to applying the extension of the methods of Zeger, Liang, and Self (1985) to the analysis of the binary hummingbird data, a series of plots was examined to determine appropriate parametric forms to be incorporated into modelling of the temporal response patterns. To check if a day effect was apparent, the average treatment responses on the logit scale (log( 1 1—Y-4-cFc ), c = .5) for successive blocks of 3, 5, and 10 trials were plotted (see Figures 3.13.3). Because there were no systematic discontinuities at those trials corresponding to new days (trials 61 and 121), there seemed to be no need to model a day effect. To get an overall picture of how treatment responses changed over time, the averages over the 6 birds within each treatment were calculated separately at each time point, and then the logits of the results were smoothed using robust locally weighted linear least squares ("lowess"). Robust locally weighted linear regression is a method for smoothing a scatterplot, (x i , yi), i = 1, .., n,, in which the fitted value at xk is the value of a linear fit to the data using weighted least squares, where the weight for (xi, yi) is large if xi is close to xk and small if it is not. "Lowess" centres a window at xk which includes the r nearest neighbours of xk on either side. The 2r + 1 points in this window are weighted symmetrically about xk. If xk is near the endpoints, however, "lowess" cannot centre the window there. Points close to the endpoint, including xk, get larger weights than they would if the window was centred at xk. If the data Chapter 3. Analysis Based on a Working Likelihood Approach^ 35 has extreme values near the endpoints, the resulting fit for xk may not be consistent with fits for points at which the window can be centred. Consequently, "lowess" fits near the endpoints may have peculiar features. Figures 3.4-3.6 show the results of this "lowess" smoothing with windows of 10, 20, and 30 points. In these figures, the fits near the endpoints do not seem peculiar because logits of treatment average responses (over 6 birds) have been smoothed. This averaging reduces the number of extreme values in the data and protects against extreme values near the endpoints. Figure 3.6 seems to best display the basic features of the change in transformed responses over time. The average responses for birds in treatments 2 (PT) and 3 (FT) appear to grow roughly linearly with time in the logit scale, with the responses for birds in treatment 3 (FT) being higher, in general, than the responses for birds in treatments 1 (NT) and 2 (PT). This is contrary to the expectations of the investigator, who used psychological principles of visual perception to predict that birds in the treatment groups with discontinuous (PT) and continous (FT) visible guides should have similar learning trajectories. For birds in treatment 1, improvement in the average response appears to be either linear or quadratic. The responses for birds in treatment 1 seem to be higher, in general, than the responses for birds in treatment 2, but not as high as those for birds in treatment 3. Of all the birds, those in treatment 3 seem to have the highest responses throughout the entire experiment. To study individual-to-individual variation within each treatment, the logits of individual results (log(1 ), c = .5) at each time point were also smoothed by "lowess", using a window of 30 consecutive time points. Figures 3.7-3.9 display the results of this smoothing. Step-like patterns associated with sudden awareness of the relationship between light cue and feeder Chapter 3. Analysis Based on a Working Likelihood Approach^ 36 were expected by the experimenter, but are not present in these plots. In all the figures, the peculiarities of "lowess" fits near the endpoints occur because individual logits, rather than logits of treatment averages are smoothed. The individual logits are more likely to have extreme values at the endpoints, leading to fits which may not accurately describe the patterns in the data of each individual. In Figure 3.7, note that two females in treatment 1 share the same "lowess" fits given by the horizontal solid line. This is because both birds have binary series which are virtually all zeroes; in other words, almost all of their first feeding forays were unsuccessful. In addition, a third female has a binary series which is predominantly zero up to the 110th trial. In the interior of the data sets, the plots do not show any obvious change in response variability over time for any of the treatments. Moreover, at any given time, the variability of responses for birds in treatment 1 (NT) appears to be larger than the variability of responses for birds in treatments 3 (FT). At any given time, the variability of responses for birds in treatment 2 (PT) appears to be the smallest of all three treatments. These plots of the "lowess"-smoothed binary responses for each bird do not suggest an obvious difference between the responses of females and males in treatment 2, but in treatments 1 and 3, the male responses seem to be higher, on average, than the female responses. In fact, when designing the experiment, the investigator did not anticipate a difference in the success probabilities for the first feeding foray of male and female hummingbirds. However, since sex may be incorporated into the model with little extra effort, two analyses of the binary hummingbird data are presented in the following sections; the first analysis does not include sex as a covariate, while the second does. The first analysis corresponds to the classical analysis (univariate analysis-of-variance for repeated measures) discussed in the previous section, where Chapter 3. Analysis Based on a Working Likelihood Approach^ 37 sex was ignored. Note that the unequal number of males and females on one of the treatments in the experiment would complicate including sex as a covariate in the classical analysis. For the binary hummingbird data on the first feeding foray, both the analysis excluding sex and the analysis including sex will be repeated under the independence and the first- order working model. The independence working model assumes that observations within individuals, as well as between individuals, are independent. The first - order working model, assumes that the observations for all individuals share a common first lag autocorrelation. 3.5 Model Fitting With Sex Excluded 3.5.1 The Initial Model For the analysis of the binary data which excludes sex as a covariate, the exploratory plots suggest that the transformed probability of success on the first feeding foray may be modelled, initially, as logit ri,t =^+ 01,it 02,it 2^i = 1,..,3, T = 1, ..,180, t = where p i T—T s.d.(T) (i = 1, 2, 3) may be interpreted as the main effect of the i th treatment. Since birds in different treatment groups are not restricted to have the same linear and quadratic coefficients for the pattern over time, interactions between time and treatment group are permitted by this model. Note that the time covariate has been standardized to orthogonalize the highly collinear linear and quadratic time covariates. Chapter 3. Analysis Based on a Working Likelihood Approach^ 38 3.5.2 An Informal Check of the Working Likelihoods When calculating the standard errors of the parameter estimates, we must take into account that the working likelihood may not be the same as the true likelihood. Section 3.3 provides consistent estimators of the covariance matrix for the parameter estimates from which "robust" standard error estimates can be obtained. Because these "robust" standard error estimates account for the fact that the working likelihood may be incorrectly specified, we expect them to be different from the "naive" standard error estimates calculated assuming that the working likelihood correctly specifies the true likelihood. If the working likelihood is correctly specified, the robust standard error estimates should be the same as the naive standard error estimates as the number of subjects becomes large. Thus, an informal method for checking how well the working likelihood approximates the true likelihood (when the number of subjects is sufficiently large) might compare the naive and robust standard error estimates or the corresponding naive and robust z-scores. Differences between robust and naive standard error estimates could result from two possibilities: (i) although the working likelihood adequately described the true likelihood, the number of subjects was not large enough for accurate estimation of the true standard errors; or, (ii) the working likelihood was not adequate for the data. Table 3.1 presents the naive and robust z-scores resulting from the initial model fits for both the independence and first-order working likelihoods. Notice that the naive z-scores for the independence fit are consistently larger, and typically substantially so, than the corresponding robust z-scores. If we assume that the robust z-scores reflect the associations between the outcome variable and the covariates accurately, we might conclude that the naive z-scores for the independence working likelihood tend to overstate the association, a well-known result Chapter 3. Analysis Based on a Working Likelihood Approach^ 39 when the observations on a subject are positively correlated. More efficient parameter and standard error estimates might then be expected from the first-order working likelihood which models some correlation. Although substantial differences still exist for many parameters, the agreement between the naive and robust z-scores appears slightly improved for the first-order working likelihood. This suggests the first-order working likelihood may better reflect the actual likelihood than the independence working likelihood; a gain in efficiency of estimation would then result from using the first-order working likelihood instead of the independence working likelihood. This conclusion is somewhat tenuous since the assumption that the robust z-scores reflect the associations between the outcome variables and the covariates accurately may not be reasonable with only 18 birds in the experiment. Qualitatively, the conclusions obtained from the two working likelihoods are identical: modelling of the treatment main effects and the linear terms seems worthwhile, while modelling of the quadratic terms for treatments 2 and 3 may be unnecessary. Table 3.1. Comparison of Z-scores from the Initial Model parameter Pi P2 /13 01,1 /31,2 31,3 /32,1 /32,2 /32,3 estimate -.31 -.71 .16 .27 .50 .60 -.21 .04 -.09 Independence robust z-score naive z-score -.97 -3.35 -3.41 -7.23 .42 1.72 2.31 4.14 7.27 3.03 2.70 9.13 -2.37 -2.82 .37 .47 -.68 -1.29 estimate -.32 -.68 .17 .26 .48 .62 -.20 .03 -.09 First-Order robust z-score naive z-score -1.00 -2.95 -3.68 -6.32 .42 1.54 2.37 3.51 3.31 6.42 2.54 7.92 -2.41 -2.23 .41 .35 -.66 -1.02 Chapter 3. Analysis Based on a Working Likelihood Approach^ 40 3.6 Model Reductions With Sex Excluded 3.6.1 Strategy for Model Reductions The primary aim of this analysis is to understand the differences in hummingbird learning across treatment groups. Because birds in different treatments may learn in different ways, the initial model for the success probabilities allows all possible interactions of time and treatment. The purpose of subsequent model reductions is to obtain a parsimonious final model which reliably addresses the investigator's questions about the way hummingbird learning differs across treatments. To focus the analysis on treatment differences, the reduction procedure should eliminate unnecessary terms for negligible interactions between treatment and time first, leaving the examination of treatment main effects for the final steps. The elimination of interaction terms reduces the complexity of the model to be used for interpretation; in particular, a final model with no interactions between treatment and time permits differences in learning patterns for birds in different treatment groups to be easily described by treatment main effects only. One such reduction procedure for the binary hummingbird data (with sex excluded as a covariate) is outlined in the following steps: 1. Fit the initial model. 2. Check the quadradic terms. • Can they be set equal? • If so, can this common term be eliminated? 3. Check the linear terms. Chapter 3. Analysis Based on a Working Likelihood Approach ^ 41 • Can they be set equal? • If so, can this common term be eliminated? 4. Examine the treatment main effects. This procedure is not meant to be followed rigidly. Instead, it is meant as a rough guideline for model reductions, to be used in conjunction with existing scientific theory and exploratory plots of the data. 3.6.2 Description of the Reductions Using the strategy outlined above, two plausible model reduction paths arise for both the independence and first-order working likelihoods. In this section, the two models are compared, and for each working likelihood, the model most consistent with patterns in the data suggested by the exploratory plots is selected to represent the experimental results. To examine how the results of the analysis depend on the choice of the working likelihood, these final models are then compared across working likelihoods. Final models with similar parameter and standard error estimates for different working likelihoods suggest that the efficiency of the two working likelihoods is about the same. Note that with only 18 birds in the experiment, the number of subjects may not be sufficiently large to ensure the adequacy of the asymptotic multivariate normal approximation to the distribution of the parameter estimates. The Wald statistics used to test the model reductions are based on this approximation. The first model reduction path is described by the following steps: 1. Check the quadratic terms. Chapter 3. Analysis Based on a Working Likelihood Approach^ 42 • Can they be set equal? Yes; the p-values for this reduction are .17 and .21 for the independence and first-order working likelihoods, respectively. • Can the common quadratic term be eliminated? Yes; the p-values for this reduction are .17 and .19 for the independence and first-order working likelihoods, respectively. (P-values are calculated using the reduced model from the previous step.) 2. Check the linear terms. • Can they be set equal? Yes; the p-values for this reduction are .29 and .26 for the independence and first-order working likelihoods, respectively. • Can this common linear term be eliminated? No; the p-value for this reduction is less than .001 using both the independence and the first-order working likelihoods. 3. Examine the treatment main effects. • Can they be set equal? Not entirely clear; the p-value for this reduction is .14 using both the independence and first-order working likelihoods. This first plausible reduction path leads to models with coefficients and standard error estimates given in Table 3.2. Note that the predicted success probabilities resulting from the final models are virtually identical. Figure 3.10 displays the final model fit for the independence working likelihood and suggests the success probabilities for birds in treatments 1 and 2 are very similar. Since corr(iii,iii) :::-., 0, the standard errors in Table 3.2 indicate that the model could be further reduced by taking main effects for treatments 1 and 2 to be the same. In fact, the p-values for the reduction are .69 and .71, with the common fitted values for p i and 1/2 being Chapter 3. Analysis Based on a Working Likelihood Approach^ 43 —.60 (s.e.=.18) and —.58 (s.e.=.17), for the independence and first-order working likelihoods, respectively. Table 3.2. Final Model Fits from the First Reduction Path Independence^First-Order Parameters estimate s.e.^estimate s.e. ^.17 ^.03 P -.53^.31^-.52^.30 Pi. -.67^.17^-.65^.15 P2 .07^.33^.08^.34 P3 .45^.11^.45^.10 /3 1,1 = / 3 1,2 ==/3 1,3 /3 2,1 = /3 2,2 = / 32,3 ^ The estimated lag one autocorrelation and standard error provided by the first-order working likelihood suggest the presence of (positive) correlation in the data of each bird. On the other hand, the first-order and independence working likelihoods result in virtually the same estimates of the regression parameters and their standard errors. Both final models indicate the transformed response depends only linearly on time, with a common slope for all three treatment groups; no terms for interaction of treatment and time are included. This lack of interaction in the fitted model allows for a simple interpretation of treatment differences through examination of treatment main effects only. For example, Figure 3.10 suggests the success probabilities of birds in treatment 1 are about the same as the success probabilities of birds in treatment 2, while birds in treatment 3 have higher success probabilities than birds in either treatment 1 or treatment 2 (about .67 units greater in the logit scale at any given time). These final models describe the success probabilities of birds in treatment 1 by a linear trajectory. This seems to disagree with the exploratory lowess plot of the transformed outcomes in Figure 3.6, which suggests the success probabilities of birds in treatment 1 could be better Chapter 3. Analysis Based on a Working Likelihood Approach^ 44 described by a quadratic trajectory. Figure 3.6 also suggests birds in treatment 2 start out with lower success probabilities than birds in treatment 1, but develop higher success probabilities by the end of the experiment. A possible explanation for the discrepancy between the lowess plot and the first final model emerges upon examining the path of model reductions leading to the first final model. Although the p-values for setting the quadratic terms for all three treatments equal are .17 and .21 for the independence and first-order working likelihoods respectively, and the p-values for setting this common quadratic term to be equal to zero are .17 and .19 for the independence and first-order working likelihoods, respectively, the p-values for proceeding directly from the full model to a model with no quadratic terms are a marginal .12 and .14 for the independence and first-order working likelihoods, respectively. This suggests the final model may be neglecting some curvature present in the data. In fact, z-scores for the quadratic term in treatment 1 ( -2.37 and -2.23 for the independence and first-order working likelihoods, respectively; see Table 3.1) lend support to what is suggested by Figure 3.6: much of the curvature present in the data comes from the learning trajectory of birds in treatment 1. Table 3.3 displays the coefficient and standard error estimates resulting from an alternate reduction path which allows curvature in the fitted learning trajectory of birds in treatment 1. This second reduction path is described by the following steps: 1. Check the quadratic terms. • Can they be eliminated? No; the p-values for this reduction are a marginal .12 and .14 for the independence and first-order working likelihoods, respectively. Chapter 3. Analysis Based on a Working Likelihood Approach ^ 45 2. Check the linear terms. • Can they be set equal? Yes; the p-values for this reduction are .32 and .27 for the independence and first-order working likelihoods, respectively. • Can the common linear term be eliminated? No; the p-values for this reduction are less than .001 for both the independence and first-order working likelihoods. 3. Examine the treatment main effects. • Different quadratic terms for different treatment groups are still present in the model, so checking whether the treatment main effects are equal corresponds to checking whether the success probabilities for all three treatment groups are equal at t = 0, or at trial 90.5. The p-values for this reduction are .12 and .13 for the independence and first-order working likelihoods, respectively, indicating that this reduction may not be justified. Note that even if this reduction was justified, it may be of little practical interest because it does not simplify interpretation of the final model. 4. For interpretability of the final model... • Eliminate the quadratic terms for treatments 2 and 3. The p-values for this reduction are .62 and .74 for the independence and first-order working likelihoods, respectively. These p-values justify what is suggested by the lowess plot: the success probabilities of birds in treatments 2 and 3 can be adequately described by a linear trajectory. Chapter 3. Analysis Based on a Working Likelihood Approach^ 46 Table 3.3. Final Model Fits from the Second Reduction Path Parameters Independence estimate s.e. p /2 2 01,1 = 01,1 = 03,1 02,1 = 02,2 = 02,3 -.31 -.67 .07 .46 -.23 .32 .17 .33 .11 .10 First-Order estimate s.e. .17 .03 .30 -.52 -.65 .15 .34 .08 .45 .10 -.22 .10 As was the case for the first reduction path, the estimated first lag autocorrelation and standard error from the first-order working likelihood suggest the presence of (positive) correlation in the data of each bird, but the predicted success probabilities resulting from the independence and first-order fits are virtually identical. Figure 3.11 plots the final model for the second reduction path using the independence working likelihood. Comparison of Figures 3.10 and 3.11 suggests the only noticeable difference between the final models for the two reduction paths is the fit for the learning trajectory of birds in treatment 1. Because the final model from the second reduction path allows curvature in the fitted learning trajectory of treatment 1 and the final model from the first reduction path does not, the former appears to be more consistent with patterns displayed by the exploratory lowess plots, and may thus be considered the more reliable summary of the experimental results. Note that the final fits for birds in treatment 1 resulting from the two different reduction paths do not differ greatly in the probability scale except at the very beginning and the very end of the experiment; see Figure 3.12. The largest discrepancies in the fitted success probabilities occur at the end of the experiment, with the maximum difference in fitted probabilities on the last trial. Here the linear model fits a success probability of .56 (with a standard error of about Chapter 3. Analysis Based on a Working Likelihood Approach^ 47 .09) and the quadratic model fits a probability of .45 (with a standard error of about .10). 3.7 Residual Analysis for the Final Model with Sex Excluded Analysis is based on the Pearson residuals (given by (yi t — fit )/0" ;st (1 — ii70) from both the in-- - dependence and first-order fits. To determine whether the independence working assumption is reasonable (recall that the independence working assumption isn't necessary for asymptotically valid results), the independence residuals are examined to see if they are approximately uncorrelated with mean 0 and variance 1. After analyzing the independence residuals, the first-order residuals are examined to determine whether the assumption of a common lag one autocorrelation for all birds (necessary for asymptotically valid results when using the first-order working likelihood) is reasonable. 3.7.1 Analysis of the Independence Residuals Under the independence working assumption, the 180 Pearson residuals for each bird should be approximately uncorrelated, with mean 0 and variance 1. Figure 3.13 displays scatterplots of the Pearson residuals for each bird and suggests the distributions of residuals for all birds have mean 0 and variance 1, with those for birds in treatments 1 and 2 being positively skewed, and those for birds in treatment 3 being more symmetric. This is explained by the relative number of O's and l's in each treatment. Only 38% of treatment 1 responses (407/1080) and 35% of treatment 2 responses (374/1080) are l's, whereas 52% of the responses in treatment 3 are l's (558/1080). The final model therefore fits the "1" responses in treatment 3 as well as the "0" responses, but for treatments 1 and 2, fits "0" responses better than "1" responses since the Chapter 3. Analysis Based on a Working Likelihood Approach^ 48 majority of responses in these treatments is 0. Note that the gaps around zero in Figure 3.13 arise because the responses are either 0 or 1, and the denominator of the Pearson residuals (Vflt (1 — CO) is essentially constant. To check if Pearson residuals from the independence final model were approximately uncorrelated within birds, estimates of the first, second, and third lag Pearson residual autocorrelations were calculated for each bird (see Table 3.4). Under the independence working assumption, the standard error of each of these autocorrelation estimates is about 1/ 180 = .075. From Table 3.4, we see that 4 out of 18, or 2/9 of the first lag autocorrelations appear significantly different from zero at level a = .05 (birds 5, 6, 14, and 15 have rh's greater than .15). Because birds are independent, we would expect this to happen for only 1 out of every 20 estimated first lag autocorrelations. Therefore, the assumption that observations within birds are independent seems questionable for these data. It is interesting to note that when the three autocorrelations are compared across birds, birds 14 and 15 (both in treatment 3) appear to have higher first, second, and third lag autocorrelations than the others. Checking the adequacy of the final model is limited by the fact that the only explanatory variables modelled are treatment group and time. The patterns of Pearson residuals across treatment groups can be seen in Figure 3.13 and have already been explained by the relative number of 0's and l's in each treatment group. Note that the Pearson residuals for the i th treatment group can be characterized by two functions of time - one (essentially the negative fit, or —17t ) for those responses which are 0; the other (essentially the negative fit plus one, or - - 1 — fii) for those responses which are 1. Thus, plotting residuals against time is not a useful check of model adequacy. Another variable, excluded as a covariate in the model because its Chapter 3. Analysis Based on a Working Likelihood Approach^ 49 role in hummingbird learning was thought to be minimal, is the sex of each bird. Figure 3.14, a boxplot of the Pearson residuals by sex, suggests sex should be included in the model. The residual distribution for males appears to be positively skewed, whereas for females it appears to be negatively skewed. The skewness of the distributions is a consequence of the females having a smaller proportion (607/1800 or 34%) of successful first feeding forays (response=1) than the males (732/1440 or 51%). Table 3.4. Estimated Autocorrelations for the Independence Residuals P1 P2 P3 Bird (Treatment 1) 1^2^3^4^5^6 .07 -.02 .000 -.06^.16 .16 .11 -.09^.09^-.05 .02 .02 .08 -.09^.03^-.04 .10 .08 Bird (Treatment 2) 7^8^9^10^11^12 P1 - .005 .13 -.09 -.07^.10^.02 P2^-.05^.07^.13 -.09^.03^.04 /53 -.009 .10 -.07^.10^-.06 .002 13 /51 P2 143 .01 .14 .02 Bird (Treatment 3) 14^15^16^17 .24 .34 -.07 .02 .16 .28 .15 .03 .28 .37 -.006 -.01 18 .10 .07 .12 3.7.2 Analysis of the First-Order Residuals For the first-order working model, let r2 denote the estimate of the l th subject's lag 1 autocorrelation obtained by regressing the Pearson residuals for that subject on their predecessors. The variance of ri, approximated by 1/ni(1 - lq), where ni and pi are, respectively, the number Chapter 3. Analysis Based on a Working Likelihood Approach^ 50 of time points and the lag one autocorrelation for the i th subject, can be stabilized by taking the inverse hyperbolic sine of ri. If the number of time-points for each subject is large, sinh -1 (ri) is approximately normally distributed with mean sinh - l(pi) and variance l/ni. The standard error of all the transformed individual first-lag autocorrelation estimates is thus approximately 1/ 180 = .075, allowing a check on the assumption of a common first lag autocorrelation for all subjects. Testing the values in Table 3.5 for homogeneity (average= .16, sample variance= .021), leads to a chi-squared statistic of 65.2 on 17 degrees of freedom; the data thus provide quite convincing evidence that a common lag one autocorrelation is not shared by all subjects. Table 3.5. Estimated Lag One Autocorrelations for the First-Order Residuals ri sinh l(ri) - Ti sinh -1 (ri) ri sinh -1 (ri) 1 .16 .16 7 -.005 -.005 13 .04 .04 Bird 2 .04 .04 (Treatment 1) 3 4 5 .27 .04 .33 .27 .04 .33 Bird (Treatment 2) 8^9^10^11 .13 .06 -.07 .17 .13 .06 -.07 .17 Bird 14 .27 .27 (Treatment 3) 15 16 17 .47 .15 .45 .45 .15 .44 6 .19 .19 12 .03 .03 18 .15 .15 The assumption of a common lag one autocorrelation is critical to deriving the asymptotic results for the first-order working model, so the fact that this assumption appears to be violated casts doubt on the validity of such an analysis. However, except for birds 3, 5, 14, 15 and 17, the estimated first-lag autocorrelations are not particularly large. The reason for compelling Chapter 3. Analysis Based on a Working Likelihood Approach^ 51 evidence against homogeneity is that the standard errors of the estimates are quite small due to the large number of responses (ni = 180) for each bird. Although the assumption of a common lag one autocorrelation does not appear to hold for the binary responses of the 18 birds in this experiment, the conclusions drawn from the analysis based on the first-order working likelihood are virtually identical to those from the analysis based on the independence working likelihood. Therefore, the conclusions drawn from the first-order analysis appear to be reasonable. The Pearson residuals obtained from the first-order working likelihood are virtually the same as those obtained by the independence working likelihood because the final models are almost identical. The plot of the first-order residuals against sex has the same features as the plot of the independence residuals displayed in Figure 3.14. Hence, analysis of both types of Pearson residuals indicates sex is an important covariate. 3.8 Summary of Analysis with Sex Excluded The models obtained by the independence and first-order working likelihoods are virtually identical. For both working likelihoods, the final model corresponding to the second reduction path has been selected to represent the experimental results. This final model indicates an increase in success probabilities over successive trials for treatments 2 (PT) and 3 (FT), and an increase, peaking at about trial 110 (near the end of the second day), followed by a decline in success probability for treatment 1 (NT). The final model uses a quadratic trend in the logit scale to describe the pattern of improvement in treatment 1, and a linear trend to describe the patterns of improvement in treatments 2 and 3. Its main features are summarized as follows: Chapter 3. Analysis Based on a Working Likelihood Approach^ 52 (i) birds in treatment 3 (FT) have greater fitted success probabilities than birds in treatments 1 (NT) and 2 (PT) throughout the entire experiment; (ii) the fitted rate of improvement for birds in treatment 2 (PT) is the same as the fitted rate of improvement for birds in treatment 3 (FT); and (iii) although birds in treatment 2 (PT) have lower fitted success probabilities than birds in treatment 1 (NT) for the intermediate trials of the experiment (trials 30-120), sometime during the course of the third day of the experiment (at about the 150t h trial), they develop greater fitted success probabilities than birds in treatment 1. Note that different working assumptions lead to virtually identical fits. This is because both approaches give similar parameter and standard error estimates, suggesting the difference in efficiency between the two approaches is minimal. The estimated lag one auto-correlation and standard error provided by the first-order working likelihood seem to indicate the presence of positive correlation in the data (see Tables 3.2 and 3.3), suggesting the first-order parameter estimates and standard errors may be more efficient. On the other hand, the first-order Pearson residuals suggest a common lag one autocorrelation may not be shared by all subjects. Because this assumption is critical in deriving the asymptotic results for the first-order parameter estimates and standard errors, the independence model fit (which does not require any assumptions about the time-dependence within each subject's responses) may be considered the more reliable summary of the experimental results. Both the exploratory plots and the residual analysis suggest sex may be an important covariate worth including in the model. In the next section, the binary hummingbird data is Chapter 3. Analysis Based on a Working Likelihood Approach^ 53 analyzed with sex included as a covariate. 3.9 The Initial Model With Sex Included The importance of sex as a covariate is clearly indicated by Figures 3.15-3.17 which plot the lowess smoothed logits of the average response by sex for each treatment group. On average, males seem to have a higher success probability than females in all three treatments throughout the entire experiment. For treatments 2 (PT) and 3 (FT), both males and females appear to have success probabilities which increase roughly linearly with time, whereas for treatment 1 (NT), the male and female success probabilities appear to increase over the first two days, but to decrease over at least part of the third day. For all three treatment groups, it appears a model incorporating an additive (on the logit scale) effect for sex, with males having higher transformed success probabilities than females, may be adequate to explain the patterns in the data. Furthermore, sex-by-treatment interaction may also be present, since the figures suggest that a different additive effect for sex in each treatment group may be necessary. These exploratory plots suggest the transformed probability of success on the first feeding foray may be modelled, initially, as T—T logit^=^sexi,s^01,1,st + 132,i,8t 2^i = 1, .., 3, s = 1, 2, T = 1, .., 180, t = ^ s.d.(T) where pi may be interpreted as the main effect of the female-to-male differential sex effect for the the i th treatment group, sex2, 1 is the i th treatment group, and sex2, 2 s-g 0 (for males in i th treatment group). The female-to-male differential sex effect is permitted to vary across treatments so treatment-by-sex interaction is incorporated. Similarly, the three-way interaction of sex, treatment and time is incorporated because different linear and quadratic coefficients Chapter 3. Analysis Based on a Working Likelihood Approach^ 54 for each sex and each treatment group are permitted. Again, the time covariate has been standardized to orthogonalize the highly collinear linear and quadratic time covariates. As in Sections 3.5-3.8, the analysis is repeated under both the independence and the first-order working model. 3.10 Model Reductions with Sex Included 3.10.1 Strategy for Model Reductions The reduction strategy for the initial model which includes sex as a covariate is similar to its counterpart with sex ignored. Because the investigator does not anticipate complicated sex-by-time interactions, the reduction strategy checks for the presence of such interactions before it checks for the interaction of treatment group and time. In addition, quadratic terms common to males and females within a treatment are treated differently from the corresponding linear terms because Figure 3.15 suggests males and females in the treatment 1 have success probabilities (for the first feeding foray) which may be described by a quadratic trend in the logit scale. The general reduction strategy is outlined by the following steps: 1. Fit the initial model. 2. Check the quadradic terms. • Can the coefficients for males and females be set equal within all treatment groups? • If so, can they be elimated? 3. Check the linear terms. • Can the coefficients for males and females be set equal within all treatment groups? Chapter 3. Analysis Based on a Working Likelihood Approach^ 55 • If so, can these coefficients be set equal for all treatment groups? • If so, can this common term be eliminated? 4. Check the differential sex effects. • Can the coefficients be set equal for all treatment groups? • If so, can this common term be eliminated? 5. Examine the treatment main effects. As in Section 3.6, this procedure is meant as a rough guideline for model reductions, to be used in conjunction with existing scientific theory and exploratory plots of the data. 3.10.2 Description of the Reductions for the Independence Working Likelihood The model reduction path for the independence working likelihood including sex as a covariate is described by the following steps: 1. Check the quadratic terms. • Can the coefficients for males and females be set equal within all treatment groups? Yes; the p-value for this reduction is .80. • Can the quadratic coefficients be eliminated? No; the p-value for the simultaneous elimination of all three quadratic coefficients is a marginal .12. Leave the quadratic terms for each treatment in the model for the time being. 2. Check the linear terms. Chapter 3. Analysis Based on a Working Likelihood Approach ^ 56 • Can the linear coefficients for males and females be set equal within each treatment group? Yes; the p-value for this reduction is .91. • Can the resulting linear coefficients be set equal for all treatment groups? Yes; the p-value for this reduction is .38. • Can this common linear term be eliminated? No; the p-value for this reduction is less than .001. 3. Check the differential sex effects. • Can they be set equal for all treatment groups? Yes; the p-value for this reduction is .40. • Can this additive sex effect be eliminated? No; the p-value for this reduction is .007. 4. Examine the treatment main effects. • Different quadratic terms for different treatment groups are still present in the model, so checking whether the treatment main effects are equal corresponds to checking whether the success probabilities for all three treatment groups are equal at t = 0, or at trial 90.5. The p-value for this reduction is .10, indicating this reduction may not be justified. As explained in Section 3.6, a reduction like this is of little practical interest because it does not simplify interpretation of the final model. 5. For interpretability of the final model... • Eliminate the quadratic terms for treatments 2 and 3. The p-value for this reduction is .74, justifying what is suggested by the lowess plots in Figures 3.16 and 3.17: the Chapter 3. Analysis Based on a Working Likelihood Approach^ 57 success probabilities of birds in treatments 2 and 3 can be adequately described by a linear trajectory. Note that quadratic terms for all three treatment groups are permitted to remain in the model until the final step, leading to a final model which is consistent with the patterns observed in the lowess plots. The reduction path for the independence working likelihood is summarized in Table 3.6. Table 3.6. Model Reductions for the Independence Working Likelihood with Sex Included /32,2,2 1 .36(.03) -.59(.31) .53(.75) -1.03(.42) -.25(.40) -.72(.76) .21(.09) .55(.21) .54(.38) .41(.31) .46(.23) .74(.27) -.21(.08) -.04(.05) -.16(.15) -.22(.21) .10(.18) /32,3,2 -.007(.18) Parameters Pi P2 P3 sex1,1 sex2,1 sex3,1 01,1,1 )31,2,1 01,3,1 /31,1,2 /3 1,2,2 /32,1,1 /32,2,1 /32,3,1 /32,1,2 Working Deviance Pearson X 2 degrees freedom 4029.67 3235.89 3222 2 .36(.17) -.53(.26) .60(.71) -1.02(A2) -.37(.30) -.86(.66) .21(.09) .54(.21) .54(.38) .41(.31) .46(.24) .73(.27) -.21(.09) .03(.10) -.09(.12) set=82,1,1 set=82,2,1 Set=/32,3,1 4031.50 3233.99 3225 Model 3^4 .39(.17) .36(.17) -.53(.24) -.53(.24) .58(.67) .60(.69) -1.03(.43) -1.06(.44) -.36(.31) -.36(.30) -.86(.64) -.83(.60) .48(.11) .29(.14) set=,31 , 1 , 1 .50(.16) .63(.24) set=/31,1,1 set=/31,1,1 set=/31,1,1 set=01,2,1 set=/31,1,1 set=/31,3,1 set 4 - 314,1 -.24(.09) -.22(.08) .04(.10) .04(.10) -.10(.14) -.10(.15) =, 2,1, i set=i32,1,1 set^3 set=82,2,1 set=82,2,1 set=g2,3,1 set=02,3,1 4048.42 4035.64 3236.09 3252.95 3228 3230 5 .18(.27) -.36(.26) .54(.48) -.75(.28) set=sexi,i set=sex1, 1 .47(.11) set= i31 , 1 , 1 set= /3 1,1,1 set=fii,i,i set=81,1,1 set=/31,1,1 -.24(.10) .04(.11) -.09(.14) set 7-132,1,1 set=82,2,1 set=f32,3,1 4062.55 3252.44 3232 6 .18(.27) -.32(.20) .45(.41) -.75(.28) set=sex 1 , 1 set=sex1,1 .48(.11) set=th, i , i set=/31,1,1 set=/31,1,1 set=/31,1,1 set=/31,1,1 -.24(.10) set=/32,1,1 - 4064.58 3252.62 3234 3.10.3 Description of the Reductions for the First-Order Working Likelihood The model reduction path for the independence working likelihood including sex as a covariate is described by the following steps: Chapter 3. Analysis Based on a Working Likelihood Approach^ 58 1. Check the quadratic terms. • Can the coefficients for males and females be set equal within all treatment groups? Yes; the p-value for this reduction is .74. • Can the quadratic coefficients be eliminated? Yes; the p-value for this reduction is .20. 2. Check the linear terms. • Can the linear coefficients for males and females be set equal within each treatment group? Yes; the p-value for this reduction is .88. • Can the resulting linear coefficients be set equal for all treatment groups? Yes; the p-value for this reduction is .34. • Can this common linear term be elimated? No; the p-value for this reduction is less than .001. 3. Check the differential sex effects. • Can they be set equal for all treatment groups? Yes; the p-value for this reduction is .33. • Can this additive sex effect be eliminated? No; the p-value for this reduction is .02. 4. Examine the treatment main effects. • Can they be set equal? Not entirely clear; the p-value for this reduction is a marginal .10. Chapter 3. Analysis Based on a Working Likelihood Approach^ 59 The reduction path for the first-order working likelihood is summarized in Table 3.7. Since corr(rii,iii) .c.--, 0, the standard errors for the final model in Table 3.7 indicate that taking the main effects for treatments 1 and 2 to be equal would be a permissible further reduction of the model. In fact, the p-value for this reduction is .35, with the fitted value common to p i and p 2 being —.19 (s.e.= .21). This reduction could be carried out for interpretability of the final model but is not included in this analysis. Table 3.7. Model Reductions for the First Order Working Likelihood with Sex Included - Parameters P iii P2 Pa sexio. sex2,1 sex3,1 01,1,1 01,2,1 01,3,1 01,1,2 /31,2,2 /31,3,2 /32,1,1 /32,2,1 /32,3,1 /32,1,2 /32,2,2 ,32,3,2 Working Deviance Pearson X 2 degrees freedom 1 .13(.03) .36(.04) -.58(.30) .59(.87) -1.03(.40) -.23(.37) -.79(.88) .20(.08) .52(.20) .54(.38) .42(.35) .46(.22) .77(.33) -.19(.09) -.03(.05) -.15(.15) -.21(.24) .09(.18) -.04(.19) 3980.30 3230.33 3221 2 .13(.03) .35(.20) -.53(.24) .65(.81) -1.02(.42) -.34(.27) -.89(.75) .20(.08) .51(.20) .53(.38) .42(.34) .46(.22) .76(.33) -.20(.10) .03(.10) -.10(.13) set=024 , 1 set=/32,2,1 set=02,3,1 3981.21 3230.80 3224 Model 3^4 .13(.03) .13(.03) .15(.25) .16(.29) -.50(.18) -.49(.19) .54(.69) .56(.76) -1.02(.42) -1.02(.41) -.34(.27) -.33(.27) -.87(.69) -.89(.76) .27(.13) .19(.07) .49(.15) .52(.20) .53(.38) .64(.26) .42(.36) .47(.23) .77(.37) - 3988.20 3234.58 3227 set=th,i,i set=fi1,2,1 set=i3 1 ,3, 1 - 5 .13(.03) .17(.26) -.50(.17) .51(.64) -1.06(.42) -.33(.27) -.84(.64) .47(.11) set=i31,1,1 set=g1,1,1 set=/31,1,1 set=-#1,1,1 set=/31,1,1 - 6 .14(.03) -.04(.27) -.30(.20) .46(.42) -.73(.28) set=sesi,i set=ses1,1 .46(.11) set=th,i,i set=th,i,i set=fli,i,i set=/31,1,1 set=/31,1,1 - 3993.00 3229.04 3230 - 4005.18 3234.18 3232 - 4017.33 3228.03 3234 As in the analysis where sex was not included as a covariate, the estimated lag one auto- correlation and standard error provided by the first-order working likelihood suggest the data for each bird are positively correlated. As before, the fits resulting from the independence and Chapter 3. Analysis Based on a Working Likelihood Approach^ 60 first-order working likelihoods are similar. When sex is included in the analysis, the only difference in the fits is the learning trajectory for birds in treatment 1. The independence working likelihood models a quadratic learning trajectory for birds in treatment 1, whereas the firstorder working likelihood models a linear learning trajectory. The fitted learning trajectories resulting from the independence and first-order working likelihoods are virtually identical for both treatments 2 and 3. Figures 3.18 and 3.19 display the fitted female learning trajectories in the logit scale for the independence and first-order working likelihoods respectively; because the effect of sex has been modelled as additive, the male trajectories are parallel to those for females, with the fitted (transformed) response for males being greater than that for females by an additive constant (-75 and -.73 for the independence and first-order working likelihood, respectively). 3.11 Summary of Analysis with Sex Included Except for the fitted learning trajectory for birds in the first treatment group, the final models obtained from the independence and first-order working likelihoods are virtually identical. Both final models describe the effect of sex as additive only; interactions of sex with time or with treatment group are not present in either model. Both models describe the success probabilities of males as being about .74 (s.e.= .28) units higher in the logit scale than those for females throughout the entire experiment; sex is clearly an important covariate for these binary response data. As in the previous models where sex was not included as a covariate, the transformed success probability for birds in treatments 2 (PT) and 3 (FT) appears to improve linearly at the same rate in both final models, with the transformed success probabilities for birds in Chapter 3. Analysis Based on a Working Likelihood Approach^ 61 treatment 3 being about .77 (s.e.= .35) units higher than those for birds in treatment 2, for both the working likelihoods, throughout the entire experiment. For the learning trajectories of birds in the first treatment group, different working likelihoods lead to different final models. The final model resulting from the independence working likelihood describes the learning trajectory for birds in treatment 1 (NT) by a quadratic trend in the logit scale for both males and females, whereas the final model resulting from the first-order working likelihood describes the same trajectory by a linear trend. For birds in treatment 1, the independence final model suggests an increase, peaking at about trial 110 (near the end of the second day), followed by a decline in success probability in the logit scale. The first-order final model, on the other hand, suggests a steady increase throughout the entire experiment (see Figures 3.18 and 3.19). Figure 3.20 displays the final fits resulting from the two working likelihoods for female birds in treatment 1. Although these fits do not differ greatly in the probability scale, they may be interpreted quite differently. For example, the independence fit suggests success probabilities for birds in treatment 2 are higher than those for birds in treatment 1 at the end of the experiment (see Figure 3.18). On the other hand, the first-order fit suggests the success probabilities for birds in treatment 1 are higher than those for birds in treatment 2 (see Figure 3.19). Collecting more data on the first treatment group in the form of additional subjects or further trials might resolve this conflict in interpretation, since the presence or absence of a quadratic effect would presumably become more clear. In this analysis, the final model resulting from the independence working likelihood is viewed as the more reliable summary of the experimental results because the independence fits are more consistent with patterns in the exploratory plots (see Figures Chapter 3. Analysis Based on a Working Likelihood Approach^ 62 3.15-3.18) than the first-order fits. The results of this analysis indicate that sex is important in describing the probabilities of success on the first feeding foray for the hummingbirds in this experiment. One possible explanation for the importance of sex is the higher body-weight to wing-span ratio of males. Because males must work harder to obtain their food, they be more likely to succeed on the first feeding foray than females. Since sex is important in this analysis, it seems that any complete analysis of the hummingbird data (binary or binomial) should include sex as a covariate. In the next chapter both the binary and binomial hummingbird data are analyzed with sex included as a covariate, using the generalized estimating equation (GEE) approach of Liang and Zeger (1986) and Zeger and Liang (1986). The GEE approach extends quasi-likelihood theory (see McCullagh and Nelder, 1989) to arrive at estimating equations for the regression parameters, and is applicable to any longitudinal response with univariate marginal distributions for which the quasi-likelihood formulation is sensible. Unlike the working likelihood approach described in this section, the GEE approach does not require assumptions about the true correlation structure when explicitly modelling correlation in the responses to improve the efficiency of estimation. In these cases, the GEE approach permits many alternatives to the independence "working correlation" structure, unlike the working likelihood approach where only one alternative is available. Chapter 4 Analysis Based on the GEE Approach This chapter presents the generalized estimating equation (GEE) approach described in Liang and Zeger (1986) and Zeger and Liang (1986). This unified approach to the analysis of longitudinal data can be applied to any repeated measures data with univariate marginal distributions for which the quasi-likelihood formulation is sensible. The GEE approach is first used to analyze the binary hummingbird data. Then, to illustrate an application to non-binary longitudinal data, it is used to analyze the binomial hummingbird data.' 4.1 Description of the GEE approach Let /Ti t , i = 1, .., K and t = 1, .., n1, be the outcome variable for the i th subject at the t th measurement, and xi t be the corresponding pxl vector of covariates. Suppose that we are interested in how the outcome variable Yi t depends on the covariates xi t . If there were only one time of observation for each subject, a generalized linear model (see McCullagh and Nelder, 1989) could be used to model a variety of continuous and discrete outcome variables. When dealing with longitudinal data, however, the situation is complicated by the fact that there may 1 A11 the GEE analyses reported here were carried out using a SAS macro received from K.Y. Liang, and written by M. Rezaul Karim. The SAS macro permits modelling of the following "working" correlation matrices: independence, stationary M-dependent, non-stationary M-dependent, AR-M, exchangeable, unspecified, and user-specified. Only working correlations which are common to all subjects may be specified. In addition, if observation times are unequally spaced or there is missing data, only the independence, exchangeable, or user-specified working correlations should be used. 63 Chapter 4. Analysis Based on the GEE Approach^ 64 be correlation among the repeated observations for a given subject. This correlation must be taken into account to obtain a valid analysis. Modelling the correlation directly is often difficult for non-Gaussian longitudinal data because there are few natural models for the joint distribution of the repeated observations. Liang and Zeger (1986) and Zeger and Liang (1986) propose a methodology based on an extension of the quasi-likelihood approach described by Wedderburn (1974) and McCullagh (1983). They model the marginal distribution of the univariate outcome at each time point by specifying a known function of the mean as a linear function of the covariates, and by assuming the (marginal) variance is a known function of the mean. Rather than specify a form for the joint distribution of the outcome vector for each subject, they merely specify a "working" correlation matrix for this outcome vector. This partial specification of the joint distribution gives rise to generalized estimating equations (GEEs) which, under weak assumptions, provide consistent estimates of the regression coefficients and their asymptotic covariance matrix. 4.1.1 Brief Description of Quasi-Likelihood Since the GEE approach is an extension of the quasi-likelihood approach to modelling independent univariate outcomes, a brief description of the pertinent aspects of quasi-likelihood theory is provided to motivate the GEE approach and establish notation. Quasi-likelihood is a unified methodology for the regression analysis of (independent) discrete or continuous outcomes first proposed by Wedderburn (1974) and examined further by McCullagh (1983). A desirable property of quasi-likelihood analysis is its flexibility. Unlike likelihood analysis, the actual form of the outcome distribution does not have to be specified; Chapter 4. Analysis Based on the GEE Approach^ 65 only the relationships between the outcome mean and the covariates and between the mean and the variance are required. To fix notation, let Yi be the univariate response for the i th individual, and let xi be the corresponding p-dimensional vector of covariates for that individual. Define pi and vi to be the expectation and variance of Yi. The expectation is assumed to be a known function of a linear combination of the covariates, i.e., = h(43), where /3 is a px1 vector of regression parameters. In addition, except for an unknown scale parameter 0, the variance is assumed to be a known function of the expectation, i.e., vi = 9(p i )10. Inference for the regression parameters /3, rather than for the nuisance parameter 0, is the primary objective. The quasi-likelihood estimator for /3 is obtained by solving a score-like system of equations where the equation for the Ph regression parameter, 0/, is given by MO ) = r_, i=1 814 ^- iti) = 0, 1 = 1, ..,p. (4.1) These are in fact score equations for /3 when the distribution of the outcome is in the exponential family. The solution to these quasi-score equations is the quasi-likelihood estimator of the regression parameters, and can be obtained by the method of iteratively re-weighted least squares (for example, by the use of GLIM). McCullagh (1983) has shown that, under mild regularity conditions, the quasi-likelihood estimator of the regression parameters is asymptotically Gaussian. Chapter 4. Analysis Based on the GEE Approach^ 66 4.1.2 Generalized Estimating Equations (GEEs) for Longitudinal Data To extend the quasi-likelihood approach to longitudinal data, Zeger and Liang (1986) adopt the structure described above for each of the univariate outcomes Yi t . Their extension involves the specification of a "working" correlation matrix for Yi, the response vector of the i th subject, and the incorporation of this correlation into the estimating equations for 0. The observation times and the working correlation matrices are allowed to differ across subjects in the GEE approach. Let Ri(a) be the nixni "working" correlation matrix for the i th subject. Zeger and Liang (1986) call Ri(ct) the working correlation matrix because it is not necessarily specified correctly. The collection of working correlation matrices is assumed to be fully specified by a, an sx1 vector of unknown parameters. The working covariance matrix for Yi is then = V 2 Ri(ce)A 1 /2 /0, where Ai is an nixni diagonal matrix with g(pi t ) as the t th diagonal element. Liang and Zeger (1986) extend the earlier quasi-likelihood equations in (4.1) and propose generalized estimating equations, K^K si(0, a) = 1 (yi - pi), where Si( i3, a) denotes the p-dimensional quasi-score function for the i th subject, pi = E(Yi) = (z11, ••, = (1(x=10), -,h(x'in, O)Y, and Di = 0140 13. Note that these equations reduce to (4.1) in the case ni=1 for all i. More generally, the quasi-score function for the i th subject, S1(/3, a) = vi -1( y, - pi), is equivalent to the estimating function suggested by Wedderburn (1974), except that here the V's are functions of a as well as 0. Because Di Vi-1 does not 67 Chapter 4. Analysis Based on the GEE Approach ^ depend on the y i 's, as long as E(Y —^= 0 (i.e., as long as the mean is correctly specified), the average of the quasi-score functions will converge to zero, ensuring that the equations have consistent roots. Since the nuisance parameters a (a parameter defined only for the working correlation) and 0 are unknown and the primary goal is inference for the regression parameters /3, an obvious way to proceed is to replace a in the expression for Vi by a consistent estimator, a__(y i , yK, /3, 0), replace 0 in a by a consistent estimator 0 yK , /3), and solve the altered estimating equations for /3. Liang and Zeger (1986) adopt this approach and define /3G, the GEE estimate of /3, as the solution of K YK^= O. Provided a is consistent given (/3, 0), is consistent given /3, and mild regularity conditions hold, Liang and Zeger (1986, Theorem 2, page 16) show that as the number of subjects, K, becomes large, ;JG consistently estimates /3. They also show that under these same conditions, K(/3G — /3) is asymptotically multivariate Gaussian with zero mean and a covariance matrix given by --1 -1 ^VG = liM K (E .1Z Vi-1 Di) [E D: Vi -1 cov(YOK -1 Di]^Vi-1 Di 1^i=1^ bill K--.co K17 -1 V. V -1 , i=1 (4.2) where cov(Yi) is the actual (unknown) covariance matrix for the response vector of the i th individual. Although in general one might expect VG, the asymptotic covariance matrix of ,jc, to depend on which estimators are used for a and 0, the above result indicates that it does not depend upon how a and 0 are estimated, as long as they are consistently estimated. The Chapter 4. Analysis Based on the GEE Approach ^ 68 covariance matrix VG can be consistently estimated by replacing cov(Yi) by (yi — rti)(yi — where^= (h(eii i3), h(eini MY , and a and d by consistent estimates in (4.2). The estimate of the regression coefficient, SG, is obtained by iterating between solving the generalized estimating equations for and consistently estimating a and cb. These two steps are repeated until convergence of all the estimates. To consistently estimate a and 0, Liang and Zeger (1986) recommend the use of the standardized residuals based on the updated estimate of /3, rit = (yit — vi Pit) ([ 17 — 1 ]tt) - 112 . Provided the fourth moments of the Yit's are finite, at a given iteration, 0 can be consistently estimated by the longitudinal analog of the Pearson X 2 statistic K^2 0 -1 = v rit , 2=1 N — p' where N = Ef= c i ni. As in many quasi-likelihood problems (where the mean-variance relationship is specified by a scale parameter times a function of the mean), direct estimation of the parameter cb may not be necessary for estimation of Q and VG. In particular, if the working correlation matrix is common to all subjects, i.e., Ri(a) = R(a), and the elements of this common working correlation matrix are all multiples of the parameter a, an estimate of be calculated. (Because 0 need not 0 appears as a multiplicative constant in the formula for ii, it cancels in the calculation of Vi = A/ 2 R(a),4 211 2 , so is not required for estimation of /3 and VG.) This is the case for many choices of the working correlation matrix which are of practical interest. To estimate a, we need to borrow strength across subjects. The specific estimator depends upon the choice of the Ri(a), but when R 2 (a) = R(a) for i = 1, .., K, the general idea is to estimate Chapter 4. Analysis Based on the GEE Approach^ 69 a by a simple function of IZuv^riuriv p• :=1 N. In practice, the true correlation structure of the response is rarely known. At best, we may only have vague ideas about the actual correlation structure present in the data. Nevertheless, we may wish to use this knowledge in order to obtain more efficient estimates of the regression parameters. The GEE approach allows us to incorporate this knowledge into the estimation process by specifying working correlation matrices we expect might be close to the true correlation structure. The GEE approach has a strong appeal because it does not require the working correlations to be correctly specified for asymptotically valid results. For binary longitudinal data, the GEE approach is more flexible than the working likelihood approach of Section 3 because working correlations other than independence and stationary 1-dependent may be used, permitting increased efficiency of estimation in many different situations. Further, when a working correlation other than the identity matrix is specified in the GEE approach, observations belonging to different subjects need not share the same correlation structure for asymptotically valid results. This is not the case in the working likelihood approach. (Recall that for valid asymptotic results using the first-order working likelihood, all subjects are required to share a common first-order autocorrelation.) However, a trade-off between efficiency and validity exists in the GEE approach, just as in the working likelihood approach. When specifying the working correlation R 2 (a) as something other than the independence working correlation in the GEE approach, the parameter a must be consistently estimated from the data in order to obtain asymptotically valid parameter and standard error estimates. Because a (a parameter defined only for the working correlation) must Chapter 4. Analysis Based on the GEE Approach^ 70 be consistently estimated, it should describe some aspect of the true correlation structure, at least indirectly. An analogous situation arises when applying the working likelihood approach developed for binary data. When using the first-order working likelihood, the stationary 1dependent working correlation is consistently estimated from the data provided all subjects have a common lag one autocorrelation. One final point regarding the GEE approach is that it may be used for any multivariate outcomes with (univariate) marginal distributions for which the quasi-likelihood formulation is sensible. A wide variety of univariate distributions such as Gaussian, Poisson, binomial (binary), gamma, and inverse Gaussian, are included in this family. Liang and Zeger (1986) suggest the approach may even be extended to include multivariate multinomial and ordinal data by using formulations described by McCullagh and Nelder (1989) for these types of variables and choosing an independence working correlation. 4.2 GEE Analysis of the Binary Hummingbird Data For the GEE analysis of the binary hummingbird data, the logit link with the binomial meanvariance relationship (up to an unknown scale factor 0) is used. If the independence working correlation is specified using this choice of link and mean-variance relationship, the scale parameter 0 does not enter into the estimation of 0 and VG. Therefore, the GEE parameter and standard error estimates using the independence working correlation will be identical to those obtained in the earlier analysis based on the independence working likelihood. In contrast, fits for the first-order working likelihood may be different from the GEE fits based on a stationary Chapter 4. Analysis Based on the GEE Approach ^ 71 1-dependent working correlation (the quasi-likelihood analog of the first-order working likelihood) because the working likelihood models the correlation directly, whereas quasi-likelihood treats the correlation as a nuisance parameter, obtaining a method-of-moments estimator after first estimating the regression parameters at each iteration of a process which continues until convergence of all the parameter estimates. In order to compare the fits obtained by the working likelihood approach and the GEE approach, in particular, the fits obtained by the first-order working likelihood and its quasilikelihood analog, the binary hummingbird data was analyzed using the GEE approach with a stationary 1-dependent working correlation. Sex was included as a covariate in the GEE analysis because of the clear indications of its importance in the working likelihood analyses of Chapter 3. To examine the robustness of parameter and standard error estimates across working correlations, the final model obtained from the stationary 1-dependent working correlation was re-fit using several other working correlations. 4.2.1 Description of the Model Reductions For the GEE analysis including sex as a covariate, the same general reduction procedure outlined in Section 3.10.1 was applied to the initial model for the success probabilities on the first feeding foray described in Section 3.9. The model reduction path for the stationary 1-dependent working correlation is described by the following steps: 1. Check the quadratic terms. • Can the coefficients for males and females be taken to be equal within all treatment groups? Yes; the p-value for this reduction is .81. Chapter 4. Analysis Based on the GEE Approach^ 72 • Can the quadratic coefficients be eliminated? No; the p-value for the simultaneous elimination of all three quadratic coefficients is a marginal .12. Leave the quadratic terms for each treatment in the model for the time being. 2. Check the linear terms. • Can the linear coefficients for males and females be set equal within each treatment group? Yes; the p-value for this reduction is .91. • Can the linear coefficients for all treatment groups be set equal? Yes; the p-value for this reduction is .21. (Note that the fitted linear coefficients of Model 3 in Table 4.1 suggest different slopes might be identified if more birds were included in the experiment.) • Can this common linear term be eliminated? No; the p-value for this reduction is less than .001. 3. Check the differential sex effects • Can they be set equal across treatment groups? Yes; the p-value for this reduction is .48. • Can this additive sex effect be eliminated? No; the p-value for this reduction is .006. 4. Examine the treatment main effects. • Different quadratic terms for different treatment groups are still present in the model, so checking whether the treatment main effects can be set equal corresponds to checking whether the success probabilities for all three treatment groups can be set Chapter 4. Analysis Based on the GEE Approach ^ 73 to be equal at t = 0 or at trial 90.5. The p-value for this reduction is .02, so it is not justified by the data. 5. For interpretability of the final model... • Eliminate the quadratic terms for treatments 2 and 3. The p-value for this reduction is .79, justifying what is suggested by the lowess plots in Figures 3.16 and 3.17; the success probabilities of birds in treatments 2 and 3 can be adequately described by a linear trajectory. Interestingly, the reduction path for this GEE analysis with a stationary 1-dependent working correlation matches that for the independence working likelihood (described in Section 3.10.2), rather than that for the first-order working likelihood (described in Section 3.10.3); the path leads to a final model with a quadratic term for the trajectory of birds in treatment 1. Figure 4.1 plots the fitted final model for females (for the males, the curves within each treatment are simply shifted .77 units higher on the logit scale), while Table 4.1 summarizes the model reductions. Although the path of reductions is the same, the final model fits obtained from the independence working likelihood and the stationary one-dependent working correlation are slightly different. In particular, the fitted linear coefficient obtained here (.39, with standard error .10), is somewhat less than that obtained from the independence working likelihood (.48, with standard error .11; see Table 3.6), although these estimates are still reasonably close considering the standard errors. The estimated female-to-male differential sex effect (-.77, with standard error .28) is very similar to that obtained from the independence working likelihood (-.75, with standard error .28). 74 Chapter 4. Analysis Based on the GEE Approach ^ Table 4.1. Reductions for the Stationary 1 Dependent Working Correlation with Sex Included - /31,1,2 /31,2,2 .41(.31) .46(.23) 2 .36(.17) -.52(.24) .60(.71) -1.02(.42) -.37(.30) -.86(.66) .22(.09) .54(.21) .54(.38) .41(.31) .46(.24) /32,1,1 /32,2,1 -.21(.08) -.21(.09) Model 3^4 .37(.16) .35(.16) -.53(.24) -.50(.24) .57(.66) .60(.70) -1.02(.41) -1.05(.42) -.39(.30) -.44(.30) -.81(.59) -.86(.64) .39(.10) .18(.11) .42(.17) set=th,i,i set=/31,1,1 .63(.24) set=/31,1,1 set=/31,1,1 set=/31,2,1 set=/31,1,1 set=/31 , 3 , 1 set=/31,1,1 -.23(.09) -.21(.09) /32,3,1 -.16(.15) /32,1,2 /32,2,2 '3 2,3,2 -.22(.21) -.008(.17) -.09(.12) set=/32,1,1 set=/32,2,1 set =i32,3, i -.10(.15) set=/32,1,1 set=fl2,2,1 set =fl2,3, i /31,2,1 1 .36(.03) -.59(.31) .53(.75) -1.03(.42) -.25(.40) -.72(.76) .22(.09) .55(.22) )31,3,1 .55(.38) Parameters iii P2 /1 3 sexi,i sex 2 , 1 sex 3 , 1 /31,1,1 /31,3,2 .74(.27) - .04(.06) .10(.18) .73(.27) .03(.10) .03(.10) .01(.11) -.10(.14) set=02,1,1 set=#2,2,1 set =02,3, i 5 .19(.26) -.35(.26) .55(.48) -.77(.28) set=sex1,1 set=sexi,i .39(.10) set=th,i,i set=/31,1,1 set=/31,1,1 6 .19(.26) -.34(.20) .46(.41) -.77(.28) set=sex1,1 set=sexi,i .39(.10) set=/31,1,1 set=/31,1,1 set=/31,1,1 set=/31,1,1 -.22(.10) set=/31,1,1 -.22(.10) set=/31,1,1 set=/31,1,1 .01(.11) - -.09(.14) set=/32 ,1, 1 set=/32,2,1 set =/32,3, i - set=/32,1,1 - - Estimates of the lag one (working) autocorrelation at each stage of the model reductions are not included in Table 4.1. These are .12, .12, .13, .13, .13, and .13, for Models 1 through 6, respectively, essentially the same as the fitted lag one autocorrelations for the first-order working likelihood which were roughly .13 (s.e. P..,- .03) for all model fits. 4.2.2 Fits of the Final Model Obtained by Other Working Correlations The final model obtained using the stationary 1-dependent working correlation was re-fit using a variety of working correlation structures; the results are summarized in Table 4.2. Qualitatively, the conclusions obtained with each working correlation are the same; all terms in the final model appear to be important. In general, the robust z-scores do not reflect as much association between the covariates and the response as the corresponding naive z-scores (not presented), which are up to 4 times larger. Results across working correlation structures are Chapter 4. Analysis Based on the GEE Approach ^ 75 remarkably similar, although the exchangeable working correlation fits a larger male-to-female sex differential and a larger (negative) main effect for treatment 2 (as well as larger standard errors for both these estimates) than the other working correlations. This suggests parameter and standard error estimates may be reasonably accurate, even with only 18 subjects. Table 4.2. Robust z-Scores (Parameter Estimates) for the Binary Final Model* Working Corr. Independent P1 P2 P3 sex(11 .68(.18) -1.54(-.32) 1.09(.45) -2.64(-.75 ) 4.31(.48) -2.47(-.24) Stat. 1-dep. Stat. 2-dep. Stat. 5-dep. Stat. 8-dep. .73(.19) .73(.19) .76(.20) .75(.20) -1.71(-.34) -1.71(-.34) -1.74(-.34) -1.69(-.33) 1.12(.46) 1.12(.46) 1.13(.46) 1.14(.45) -2.76(-.77) -2.75(-.77) -2.76(-.77) -2.75(-.76) 3.89(.39) 3.89(.39) 3.98(.40) 4.03(.40) -2.27(-.22) -2.26(-.22) -2.34(-.23) -2.40(-.23) AR-1 AR-2 AR-5 .73(.19) .73(.19) .76(.20) -1.71(-.34) -1.71(-.34) -1.74(-.34) 1.12(.46) 1.12(.46) 1.13(.46) -2.76(-.77) -2.75(-.77) -2.76(-.77) 3.89(.39) 3.90(.40) 4.01(.40) -2.27(-.22) -2.26(-.23) -2.38(-.24) Exchangeable .55(.16) -1.83(-.54) 1.14(.50) -2.69(-.95) 3.51(.40) -2.28(-.23) 13121,1 , 133 133i,1 * Entries in table are z-score (parameter estimate). 1. sex 1 , 1 = sex i^common , 2 = sex 1 , 3 sex, the comn differential sex effect. 2. )31,1,1 — 131,1,2 - 01,2,1 - 131,2,2 - 131,3,1 - 31,3,2, the common linear coefficient. 3 . i(32,1,1 = )32 ,1, 2 , the common quadratic coefficient for birds in treatment 1. / Note that among the working correlations for which results are presented in Table 4.2, the exchangeable structure is scientifically the least plausible candidate for the true structure because it models the correlations for responses distant in time to be the same as those for responses adjacent in time. Oddly enough, of all the working correlations presented, the naive and robust z-scores for the exchangeable structure are the closest together. Although this does not necessarily imply the exchangeable structure best approximates the true correlation structure, it is interesting to note that when the working correlation structure is specified to be stationary 8-dependent, the estimated lag one through eight working correlations are all Chapter 4. Analysis Based on the GEE Approach ^ 76 about .12 (slightly larger than .08, the fitted exchangeable working correlation), suggesting correlations for responses relatively far apart in time are not too different from correlations for responses adjacent in time. 4.3 GEE Analysis of the Binomial Hummingbird Data To illustrate the versatility of the GEE approach, we turn now to an analysis of the "binomial" data. Recall that the maximum number of successful feeding forays on any given trial in the experiment was 11, while the number of unsuccessful feeding attempts was unrestricted. Because the total number of feeding forays in a trial is random, the "binomial" response, i.e., the number of successful feeding forays on each trial, cannot have a binomial distribution. Nevertheless, we assume the mean-variance relationship for the number of successful feeding forays in a trial is that of a binomial random variable and interpret the fitted success probability for a particular trial as a rough average of the success probabilities for all the (possibly dependent) forays within that trial, rather than as a success probability for a binomial experiment (i.e., the trial) in which the total number of (independent) forays is non-random. The binomial data is not as meaningful to the investigator as the binary data because hummingbird researchers believe the first feeding forays are the most accurate measure of a bird's expectations about feeder profitability. Many of the birds which are unsuccessful on the first foray have not yet associated the light cue with the feeder. For these birds, any forays subsequent to the initial unsuccessful attempt will reflect their strategy for finding food after failing the first time. This strategy varies with each individual; so, on a particular trial, a bird's binomial response (which summarizes information on all of its forays during that trial) Chapter 4. Analysis Based on the GEE Approach^ 77 is expected to convey less information about the bird's feeding expectations than its binary response at the first feeding foray. Even though analysis of the binomial hummingbird data may not be as scientifically relevant as the earlier analysis of the binary data, it is worthwhile as an example of the GEE approach applied to non-binary data. 4.3.1 Exploring the Binomial Hummingbird Data A series of plots was used to determine appropriate parametric forms to be incorporated into modelling of the temporal response patterns. To check if a day effect was apparent, the logits of the individual responses (logit(y) = log(-n 2i e - - ), c = .5, where y is the number of successful forays out of n total forays in the trial for a particular bird) were averaged over all birds of the same sex within each treatment group for successive blocks of 10 trials and then plotted (see Figures 4.2-4.4). Because there were no systematic discontinuities at those trials corresponding to new days (trials 61 and 121), there seemed to be no need to model a day effect. To get a better overall picture of the binomial data, the logits of responses within each treatment-sex group were smoothed using robust locally weighted linear least squares ("lowess") with windows of 30 and 60 points. Figures 4.5 and 4.6 show the results of this "lowess" smoothing. Comparing across treatment groups, the success probabilities of birds in treatment 3 (FT) are generally higher than those of birds in the other treatments, suggesting the presence of treatment main effects. Success probabilities of males in treatment 1 (NT) appear to be almost as high as those of males in treatment 3, and roughly the same as those of females in treatment 3. In contrast, success probabilities of females in treatment 1 appear to be lower than those of females in treatment Chapter 4. Analysis Based on the GEE Approach ^ 78 3 and closer to those of females in treatment 2 (PT). In general, males appear to have higher probabilities of success than females, especially for birds in treatment 1. This suggests the presence of sex-by-treatment interaction in this data. The success probabilities of all birds appear to increase roughly linearly over time, with the slope for females in treatment 1 being somewhat less than the slopes for the other birds. This suggests a three-way interaction of sex, treatment group and time may be present in this data. To study individual-to-individual variation within each treatment, the logits of individual results at each time point were smoothed using lowess, with a window of 60 consecutive time points; see Figures 4.7-4.9. Peculiarities of these "lowess" fits near the endpoints occur because only the individual logits, rather than logits for all the individuals of a certain sex within a treatment, have been smoothed. The individual logits are more likely to have clusters of extreme values with few moderate values to counteract their effect at the endpoints, leading to fits for an individual which may not accurately describe the patterns in the individual's data near the endpoints. In the interior of the the data sets, the plots do not show any obvious change in response variability over time for any of the treatments. The plots also suggest that, at any given time, the variability of responses in treatments 1 (NT) and 2 (PT) is roughly the same, but that responses in treatment 3 (FT) are more variable at than responses in treatments 1 and 2. Chapter 4. Analysis Based on the GEE Approach ^ 79 4.3.2 Model Fitting for the Binomial Hummingbird Data These exploratory plots suggest the transformed probability of success may be modelled, initially, as logit^=^sexi,s^ 2,i,st2^= 1, .., 3, s = 1,2, T = 1, ..,180, t=T—T s.d.(T) where pi may be interpreted as the main effect of the i th treatment group, sexi, i is the femaleto-male differential sex effect for the i th treatment group, and sexo s=-St 0 (for males in the i th treatment group). As in the binary analysis, a quadratic learning trajectory (in the logit scale) has been modelled to allow for some departure from linearity in the responses. In addition, the time covariate has been standardized to orthogonalize the highly collinear linear and quadratic time covariates. The female-to-male differential sex effect is permitted to vary across treatments to allow for treatment-by-sex interaction. Similarly, the three-way interaction of sex, treatment and time is incorporated because different linear and quadratic coefficients for each sex and each treatment group are permitted. Because the stationary 1-dependent working correlation models some correlation in the responses, it is expected to provide more efficient estimates than the independence working correlation. Model reductions have therefore been carried out using the stationary 1-dependent working correlation and are summarized in Table 4.3. The path of these reductions is described by the following steps: 1. Check the quadratic terms. • Can the coefficients for males and females be taken to be equal within all treatment groups? Yes; the p-value for this reduction is .65. Set these coefficients to be equal. Chapter 4. Analysis Based on the GEE Approach^ 80 • Can the quadratic coefficients in the resulting model be eliminated? Yes; the p-value for the simultaneous elimination of all three quadratic coefficients is .58. Remove the quadratic terms from the model. 2. Check the linear terms in the resulting model. • Can the linear coefficients for males and females be set equal within each treatment group? No; the p-value for this reduction is less than .001. From the coefficients and standard errors of Model 3 in Table 4.3, it appears that essentially all the evidence against this reduction comes from the responses for birds in treatment 1. • Can the linear coefficients for all treatment groups be set equal within sexes? Yes; the p-value for this reduction is .25. To simplify the model, set two common linear slopes: one for females and one for males. • Can the resulting linear terms for females or males be eliminated? For the females, yes; the p-value for the reduction is .67, so remove this term from the model. For the males, no; the p-value for this reduction is less than .001. 3. Check the differential sex effects. • Can they be taken to be the same for all treatment groups? Yes; the p-value for this reduction is .47. • Can this the common sex effect be eliminated? No; the p-value for this reduction is .05. 4. Examine the treatment main effects. 81 Chapter 4. Analysis Based on the GEE Approach ^ • Can they be taken to be equal? No; the p-value for this reduction is .02. For interpretability of the final model, a further reduction might set the main effects for treatments 1 and 2 to be equal, as suggested by standard errors of these estimates under Model 6 (see Table 4.3). The p-value for this reduction is .63, indicating the reduction is permissible. The resulting model (not included in Table 4.3) leads to a simpler representation of the learning trajectories because, within each sex, birds in treatments 1 and 2 are then modelled as having the same learning trajectory in the logit scale. Table 4.3. Model Reductions for the Binomial Data Parameter Pi P2 P3 sex 1 , 1 sex 2 ,1 sex3,1 th,i,i /31,2,1 /31,3,1 /31,1,2 /3 1,2,2 /3 1,3,2 /32,1,1 /32,2,1 /32,3,1 /32,1,2 /32,2,2 / 1 1.08(.15) 1.08(.14) 2.22(.71) -.35(.49) -.12(.20) -.76(.71) -.06(.09) .23(.09) .10(.24) .31(.02) .39(.18) .48(.30) -.01(.11) -.04(.08) -.14(.04) -.01(.17) .04(.15) .05(.16) 2 1.08(.08) 1.12(.12) 2.36(.66) -.35(.36) -.20(.23) -.95(.63) -.06(.09) .23(.09) .10(.25) .31(.02) .39(.18) .45(.30) -.01(.09) -.01(.08) -.09(.07) set=i3 24 , 1 set--7432 , 2 , 1 set=/32,3,1 Model 3^4 1.08(.02) 1.07(.01) 1.11(.13) 1.12(.12) 2.27(.63) 2.24(.57) -.35(.36) -.38(.38) -.19(.23) -.20(.23) -.93(.57) -.96(.64) -.06(.09) .06(.09) .23(.10) set=/31,1,1 .10(.26) set=/31,1,1 .31(.01) .37(.09) .39(.18) set=th,1,2 .47(.33) set=13 1,1,2 - 5 1.08(.02) 1.12(.13) 2.25(.57) -.37(.37) -.18(.23) -.93(.57) .37(.09) set=/3 1,1,2 set=/31,1,2 - 6 1.12(.18) 1.24(.17) 1.91(.30) -.43(.22) set=sex 1 ,1 set=sexi,i .37(.09) set="31,1,2 set=/31,1,2 - - - - - - - - - The fits for the initial model with the stationary 1-dependent working correlation given in Table 4.3 are plotted in Figures 4.10-4.12; similar fits are obtained with both the independence and the AR-1 working correlations. These fits do not reflect some features of the lowess plots in Figures 4.5 and 4.6. For example, the lowess plots suggest the male-to-female differential sex effect should be largest for birds in treatment 1 (NT), yet the GEE fits for the initial model Chapter 4. Analysis Based on the GEE Approach^ 82 suggest this effect is largest for birds in treatment 3 (FT). In addition, the lowess plots suggest success probabilities of treatment 1 males and treatment 3 females and males are higher than those of birds in treatment 2 (PT) on the logit scale, but in the fits for the initial model, only the success probabilities (on the logit scale) of treatment 3 males are clearly higher than those of birds in treatment 2. Finally, the lowess plots suggest success probabilities for females in treatment 1 increase slightly over time, whereas the fits for the initial model indicate they decrease slightly over time. The initial fit for females in treatment 1 is probably overly influenced by the female with a large dip in success probabilities near trial 140; see Figure 4.7. Discrepancies such as these may occur because, unlike the lowess fits for intermediate trials, the GEE fits at intermediate trials may be influenced by extreme data at the beginning and end of the experiment. Runs of extreme data are more likely to occur for a particular individual at the beginning or end of the experiment and may have high leverage because the number of subjects in the experiment (18) is very small relative to the number of trials (180). Such discrepancies between lowess and initial model fits may disappear once data on more birds are collected under these conditions. A noteable peculiarity of the GEE results in Table 4.3 is the large drop (from .08 to .01) in the standard error of the estimated effect for treatment 1 compared to those for the other treatments (from .12 to .12, and from .66 to .63, for treatments 2 and 3, respectively) when the quadratic terms are removed from Model 2. The corresponding naive standard error estimate stays relatively constant for this model reduction, dropping from .11 to .09. Another jump (from .02 to .18) in the standard error of the estimated effect for treatment 1 occurs when the differential sex effects for the treatment groups are set equal (Model 5 to Model 6 in the table). Chapter 4. Analysis Based on the GEE Approach ^ 83 Again, the corresponding naive standard error estimate remains relatively constant (dropping from .09 to .07). These results suggest the robust standard error estimate for the main effect of treatment 1 may not accurately reflect the true standard error; perhaps there are too few birds (18) in the experiment for the asymptotic approximation to be adequate. If the asymptotic approximation to the covariance matrix of the parameter estimates is inaccurate, the validity of all the model reductions (and hence the final model) is questionable. Somewhat reassuringly, the values of standard error estimates for the other parameters are fairly constant. Figures 4.13 and 4.14 display the final model fits for the binomial hummingbird data obtained with the stationary 1-dependent working correlation. Model reductions are intended to lead to a parsimonious final model which describes essential features of learning trajectories for hummingbirds in general, rather than those for the eighteen hummingbirds in this experiment. Hence, some features evident in the lowess plots are absent in the final model fits, presumably because of the large variability of the responses. For example, the variability of female success probabilities, primarily in treatments 1 and 3 (see Figures 4.7-4.9), permits the fit of a common linear slope for females in Model 4 of Table 4.3. This feature of the fits is not consistent with the lowess plots, which suggest success probabilities for females in treatments 2 and 3 increase more quickly over time on the logit scale than those for females in treatment 1. When averaging all the female curves, those for the four females in treatment 1 dominate, leading to a negligible fitted slope which is eliminated in Model 5 of Table 4.3. Because both a common (negligible) female and a common male slope are modelled, a common differential sex effect is permissible in Model 6. This feature of the fits is also inconsistent with the lowess plots, which suggest a higher male-to-female differential sex effect in treatment 1. Specifying a common differential 84 Chapter 4. Analysis Based on the GEE Approach ^ sex effect affects the placement of the fitted curve for males in treatment 1 relative to the fitted curves for males in other treatments. For this reason, the fitted curve for males in treatment 1 is lower than the curve suggested by the lowess plots. In order to identify which of these patterns more accurately reflects the true response probabilities of hummingbirds under these experimental conditions, more birds need to be included in the experiment. However, one feature of the binomial data suggested by the lowess plots is confirmed by the GEE analysis: the (binomial) success probabilities for birds in treatment 3 are higher, in general, than those for birds in treatments 1 and 2. 4.3.3 Fits of the Final Model Obtained by Other Working Correlations To examine the robustness of parameter and standard error estimates across working correlations, the final model was re-fit using a variety of working correlation structures; results are summarized in Table 4.4. Table 4.4. Robust z-Score (Parameter Estimate) for the Binomial Final Model* Working Corr. Independent /21 6.03(1.08) /12 6.48(1.13) 113 5.76(1.81) sex „1) -1.81(-.40) /3121).,2 3.95(.36) Stat. 1-dep. Stat. 2-dep. Stat. 3-dep. 6.23(1.12) 6.44(1.18) 5.82(1.31) 7.33(1.24) 8.56(1.43) 9.86(1.75) 6.32(1.91) 7.10(2.04) 7.39(2.23) -1.94(-.43) -2.13(-.47) -2.16(-.54) 4.12(.37) 4.06(.37) 3.41(.38) AR-1 AR-2 6.27(1.12) 6.37(1.22) 7.49(1.26) 9.23(1.57) 6.42(1.92) 7.46(2.15) -1.97(-.43) -2.19(-.50) 4.13(.37) 3.98(.38) ( * Entries in table are z-score (parameter estimate). 1. sex1, 1 = sex 1 ,2 = sex1,3 = sex, the common differential sex effect. 2. 01,1,2 = 01,2,2 = 01,3,2, the common linear coefficient for males. Chapter 4. Analysis Based on the GEE Approach ^ 85 The robust z-scores for different working correlations are reasonably close, leading to qualitatively similar conclusions: all terms in the final model are important. For the working correlations used in the analysis of this binomial hummingbird data, the robust z-scores do not reflect as much association between the covariates and the response as the naive z-scores (not presented), which are up to 5 times larger. When naive and robust z-scores are compared across the working correlations used in this analysis, the differences between naive and robust z-scores for the stationary 3-dependent working correlation are the smallest. This, together with the fact that the stationary 3-dependent correlation structure is scientifically plausible, might suggest it approximates the true correlation better than the other working correlations presented in the table. For the stationary 3-dependent working correlation, estimates for the first, second and third lag autocorrelations are .24, .26, and .26, respectively, suggesting non-negligible correlation may indeed be present in the data. Note that estimates of all the effects, and, in particular, the treatment main effects, increase as more correlation structure is incorporated into the working assumptions. In fact, when solving the generalized estimating equations iteratively using the stationary 5-dependent, stationary 8-dependent, AR-3, and exchangeable working correlations (not included in the table), successive estimates of the treatment main effects become very large and do not converge. This is in contrast to the binary data on the first feeding foray analyzed in the previous section, where convergence was achieved with all of these working correlations. For the binomial hummingbird data, lack of convergence with more detailed working correlations may be related to the large number of trials in the experiment (180) relative to the small number of subjects (18). Chapter 4. Analysis Based on the GEE Approach ^ 86 This same feature of the data set (180 trials on only 18 subjects) may explain discrepancies between the lowess fits and the GEE fits for the full model in the binomial analysis: because the number of trials (180) is large relative to the number of subjects (18), the GEE fits at intermediate trials may be overly influenced by a run of extreme data for a particular individual at the beginning or end of the experiment. This problem does not arise in the binary analysis because there are no extreme responses for binary data (binary data takes on only two values). Hence, the GEE and lowess fits for the binary data are reasonably consistent. The binomial analysis seems to illustrate that when a long series of observations is made on relatively few subjects, the GEE approach should be used with caution since it relies on combining strength across subjects to (i) obtain accurate estimates of the regression parameters; and (ii) compensate for making no strong assumptions about the correlation structure of each subject's response. Chapter 5 Conclusion 5.1 Discussion of Results Exploratory plots of the binary responses for each bird on the first feeding foray (Figures 3.7-3.9) do not reveal the step-like pattern associated with sudden awareness of the spatial association between light cue and feeder expected by the experimenter. These exploratory plots suggest a more gradual increase in individual success probabilities over time, but with considerable fluctuation throughout the experiment. Evaluating the performance of birds across treatment groups by the three methods described in this thesis, both the working likelihood and GEE analysis indicate that, as expected, birds in treatment 3 (FT) have the highest success probabilities throughout the entire experiment. The univariate ANOVA, however, does not detect an important difference in success probabilities across treatment groups. Unlike the ANOVA, both the working likelihood and GEE approaches permit modelling of the response directly as a function of time. When the learning trajectories of birds are modelled by a quadratic trend in the logit scale using these approaches, birds in treatments 2 (PT) and 3 (FT) appear to improve linearly at the same rate. However, success probabilities for birds in treatment 2 appear to start out (and remain) lower than those for birds in treatment 3 87 Chapter 5. Conclusion^ 88 throughout the entire experiment. Unlike the linear trajectories (in the logit scale) for birds in treatments 2 and 3, the trajectory for birds in treatment 1 (NT) appears to be quadratic, with success probabilities at the beginning of the experiment closer to those for birds in treatment 3 than those for birds in treatment 1. However, at the end of the experiment, success probabilities of birds in treatment 1 appear to drop closer to and perhaps below those of birds in treatment 2. Contrary to the expectations of the investigator, birds in treatment 2 (PT) appear to have lower success probabilities than the other birds except, perhaps, at the end of the experiment, where their success probabilities may be greater than or equal to those for birds in treatment 1. One possible explanation proposed by the investigator is that the orange Dymo tape for the partially taped feeding array creates three horizontal bands which may interfere with the intended vertical pattern to be perceived by the birds in treatment 2 (see Figure 1.1). Contrary to the accepted wisdom among researchers in hummingbird learning, both the working likelihood and GEE analyses reveal that sex is important in describing the probabilities of success on the first feeding foray for the hummingbirds in this experiment. Since the effect of sex appears to be additive for the binary data, the fitted models indicate the odds of success for males are more than twice the odds of success for females. (For the fit resulting from the independence working likelihood, the odds for males are about 2.1 the odds for females, with a standard error of roughly 0.6.) One possible explanation for the importance of sex is the higher body-weight to wing-span ratio of males. Because males must work harder to obtain their food, they may be more likely to succeed than females. Note that in any ANOVA of the results in this experiment, sex would be difficult to incorporate because of the imbalance in the number of males and females in treatment 1 (NT). Chapter 5. Conclusion^ 89 Although not as scientifically relevant as the analysis of the binary data for the first feeding foray, the analysis of the binomial data for all the feeding forays in a trial confirms that birds in treatment 3 (FT) have higher success probabilities, in general, than birds in the other treatment groups. Again, sex appears to be an important covariate, although its effect on the binomial response (where sex appears to interact with time) seems to be more complicated than its effect on the binary response (where the effect of sex is simply additive). Because both the working likelihood and GEE approaches model success probabilities as a function of time while taking into account the (unknown) correlation present in the data, they enable much more information to be extracted from the hummingbird data set than a more traditional analysis-of-variance. These methodologies therefore provide useful and powerful tools for researchers in this subject area. 5.2 Comparison of Methods Longitudinal data consists of repeated measurements over time for each of many subjects. Hence, the response of each subject may be considered a random vector with possibly correlated elements corresponding to measurements at each time of observation. If the response vector is normally distributed, classical methods permit modelling of the joint distribution as a function of the covariates. For example, when sex is not included as a covariate in the analysis of the binary hummingbird data, MANOVA permits testing of parallelism and equality of treatment means, provided the number of subjects is greater than the number of time points. MANOVA may be carried out when the number of subjects is less than the number of time points if data are blocked over time points, but this leads to a loss of information. Modelling of responses which Chapter 5. Conclusion^ 90 are not normally distributed is more difficult. Transformations which stabilize variance for univariate responses do not necessarily stabilize covariance matrices for multivariate responses. Hence, equality of covariance matrices, a necessary assumption for valid F-tests in analysis-ofvariance, may not hold. Instead of attempting to transform non-normal data into forms more suitable for analysis-ofvariance, a more natural approach might be to model the marginal expectations of the responses as functions of the covariates while taking into account the (unknown) correlation present in the data of each subject. For example, with binary longitudinal data, logistic regression might be extended in the manner described by Zeger, Liang and Self (1985). Rather than the "actual" likelihood, this approach uses a "working" likelihood based on simple assumptions about the correlation within each subject's data to generate estimating equations for the regression parameters. Under mild assumptions about the true correlation structure, these estimating equations lead to consistent estimates of the regression parameters and their standard errors. The approach described by Zeger, Liang, and Self (1985) was developed for short binary series. Because detailed models for the change in response over time are unnecessary in this context, the authors considered time-independent covariates only. In the case of the hummingbird learning experiments, where birds are expected to improve over a long series of measurements, the probability of success should be modelled as a function of time. In this thesis, the results of Zeger, Liang, and Self (1985) have been extended to allow this. The resulting "working likelihood" approach may be applied to any binary series with time-dependent covariates. One weakness of this approach is that when correlation is explicitly modelled (in order to increase the efficiency of estimation), asymptotic results only hold under certain assumptions on the Chapter 5. Conclusion^ 91 true correlation structure. Other than the independence working likelihood, only the firstorder Markov working likelihood, which models stationary one-dependent correlation in the responses, has been considered in this thesis. Unlike the working likelihood approach, the "generalized estimating equation" or GEE approach of Liang and Zeger (1986) and Zeger and Liang (1986) requires no strong assumptions about the true correlation structure when modelling correlation explicitly. The GEE approach extends quasi-likelihood theory to the case where the response is observed repeatedly over time for each subject, and uses a "working" correlation for these responses to arrive at estimating equations for the regression parameters. Hence, the GEE approach is applicable to any longitudinal response with univariate marginal distributions for which the quasi-likelihood formulation is sensible. This includes binary, binomial, count, and continuous longitudinal responses. The GEE approach thus provides a unified methodology for the analysis of longitudinal data. Future work investigating the extension of the GEE approach to the case of ordinal longitudinal data, and to the case of data with censored or missing observations would be valuable. For instance, the approach could be used to analyze the results of an ongoing clinical trial, in which U.B.C. is one of the participating centres, concerning the effects of beta-interferon on the progression of multiple sclerosis. In this clinical trial, several outcome variables are observed on each patient repeatedly over time, including the so-called "Kurtzke score", an ordinal response variable measuring the severity of disease. Figure 1.1. Feeding Arrays for the Hummingbird Experiment lit cue • ooc oo 4-- cues Key 12 cm 000000^feeders 3 cm • 0 0 0 0 0 No Tape Treatment 0 0 0 0 0 0 r--- Partial Tape Treatment Full Tape Treatment 92 Figure 2.1. Mean Treatment (Transformed) Response versus Measurement for Blocks of 15 2 ^ 4 ^ 6 ^ 8 ^ 10 ^ Measurement Figure 2.2. Mean Treatment (Transformed) Response versus Measurement for Blocks of 10 5 ^ 10 ^ Measurement 15 12 Figure 2.3. Estimated Variances versus Measurement for Blocks of 15 2 ^ 4 ^ 6 ^ 8 ^ 10 ^ Measurement Figure 2.4. Estimated Variances versus Measurement for Blocks of 10 O O 5 10 Measurement 15 12 Figure 3.1. Average Treatment Response in the Logit Scale for Blocks of 3 Trials 0 10 20 30 Trial Block of 3 40 50 60 Figure 3.2. Average Treatment Response in the Logit Scale for Blocks of 5 Trials 0 ^ 10 ^ 20 Trial Block of 5 ^ 30 Figure 3.3. Average Treatment Response in the Logit Scale for Blocks of 10 Trials 0 ^ 5 ^ 10 Trial Block of 10 ^ 15 Figure 3.4. Smoothed Logit of Average Treatment Response versus Trial (Window=10) .88 .82 .73 .62 .50 .38 .27 .18 ^ treatment 1 .12 ^ treatment 2 — - — treatment 3 .08 0 50 100 Trial 150 Figure 3.5. Smoothed Logit of Average Treatment Response versus Trial (Window = 20 Trials) 0 treatment 1 ^ treatment 2 treatment 3 — - — cy - iI 0 ^ 50 100 Trial I 150 Figure 3.6. Smoothed Logit of Average Treatment Response versus Trial (Window = 30) 0 ^ 50 ^ 100 Trial ^ 150 Figure 3.7. Smoothed Individual Logits for Treatment 1 0 ^ 50 ^ 100 Trial ^ 150 Figure 3.8. Smoothed Individual Logits for Treatment 2 Trial Figure 3.9. Smoothed Individual Logits for Treatment 3 0 ^ 50 ^ 100 Trial ^ 150 Figure 3.10. First Independence Fits without Sex ^ Figure 3.11. Second Independence Fits without Sex .73 .62 .50 .4 o .38^•-• . .27 0 50 1 ^ Treatment 1 2 ^ Treatment 2 1 ^ Treatment 1 2 ^ Treatment 2 3 — - — Treatment 3 3 — - — Treatment 3 100 trial 150 0 50 100 trial 150 .18 - Figure 3.12. Fitted Success Probabilities for Treatment 1 from the Independence Working Model trial Pearson Residuals -1 aa••••••01M -s 2 1 0 all••••••••• OM N •• MO •• • OW •^• - all• OP OD • 40 • • • co - a • ail.^• •^• • • •^•^• M4041 ••• • • •^•• • •^• • - - •. •. • MO •• • • • 01•••••••^ •••^ Gal•••■••• • afa ••■••••••^ Miala10• • •• • • • • •• 4, • • - a.••••• ••N• WOO. • • • CO - 4••4■•^ •••••■••1••• • •• • • a • • • •••••• a • •• • • OD • • • • • ••1•10111MMILONSM• 0111••••■••^•D • -L. a•1•1•0•1044•••■•••a• 041•00 •• 4•allaa•••••••• • • • • •• • • 01 aNIONIaMeaLl• ••• .. • .0 •••• • • ••••• • 411.• • ^••••••• •••••■• • • • • MO MO NO Ma ^401•0000 • ••■••■•••••1414aMMID 106 Figure 3.14. Boxplot of Independence Pearson Residuals by Sex male ^ female sex Figure 3.15. Smoothed Logit of Average Response vs. Trial for Treatment 1 (Window = 60) Figure 3.16. Smoothed Logit of Average Response vs. Trial for Treatment 2 (Window = 60) 0 trial 50 100 trial Figure 3.17. Smoothed Logit of Average Response vs. Trial for Treatment 3 (Window = 60) 150 trial Figure 3.18. Binary Independence Fits for Females Figure 3.19. Binary First-Order Fits for Females .62 .50 tr) 0 .62 0 .50 .38 .27 .18 .12 0 50 100 trial 150 1^Treatment 1 to 2 ^ Treatment 2 3 --- Treatment 3 .18 .12 cJ- 0 50 100 trial 150 Figure 3.20. Fitted Success Probabilities for Females in Treatment 1 I - Independence Fit F ^ First-Order Markov Fit trial Figure 4.1. Binary Fits for Females Obtained from the Stationary 1-Dependent Working Correlation O. 0 iY/. 1 - Treatment 1 2 ^ Treatment 2 3 --- Treatment 3 0 50 100 trial 150 Figure 4.2. Average Logit of Response by Sex for Treatment 1, Blocks of 10 Trials Figure 4.3. Average Logit of Response by Sex ^Figure 4.4. Average Logit of Response by Sex for Treatment 2, Blocks of 10 Trials^ for Treatment 3, Blocks of 10 Trials .92 — .92 .88 — .88 .82 .73 .82 so .38 5 10 Trial Block of 10 15 10 Trial Block of 10 15 5 10 Trial Block of 10 15 Figure 4.5. Smoothed Logit of Success ^Figure 4.6. Smoothed Logit of Success Proportion versus Trial (Window = 30) Proportion versus Trial (Window = 60) N - T. 0 - 0 F^Treatment 1, Female Treatment 1, Female M^Treatment 1, Male F ^ Treatment 2, Female Treatment 1, Male M ^ Treatment 2, Male Treatment 2, Male F^Treatment 3, Female M^Treatment 3, Male Treatment 3, Female 50 100 Trial Treatment 2, Female Treatment 3, Male 150 0 50 100 Trial 150 Figure 4.7. Smoothed Individual Logits (Binomial Data) for Treatment 1 (Window = 60) Figure 4.8. Smoothed Individual Logits (Binomial Data) for Treatment 2 (Window = 60) .73 Figure 4.9. Smoothed Individual Logits (Binomial Data) for Treatment 3 (Window = 60) .95^r, .95 .88 ea .73 .73 1 50^100^150 Trial 0 ^ 50 .50 .27^'7 .27 .12^°4 .12 50^100^150 Trial Trial Figure 4.10. Initial Model Fits for Treatment 1 (Binomial Data) Figure 4.11. Initial Model Fits for Treatment 2 (Binomial Data) Figure 4.12. Initial Model Fits for Treatment 3 (Binomial Data) 20 ^ F Females ^ M Males co (,) a) 0 rn he . 50 ^ Trial ^ ^ 100 150 ^ 0 ^ 50 ^ Trial ^ ^ 100 150 ^ 0 ^ ^ 50 Trial 100 ^ 150 Figure 4.13. Final Model Fits for Females from the Stationary 1-Dependent Working Correlation Figure 4.14. Final Model Fits for Males from the Stationary 1-Dependent Working Correlation .92 1 ^ Treatment 1 2 ^ Treatment 2 1^Treatment 1 2 ^ Treatment 2 3 — - — Treatment 3 3 — - — Treatment 3 0. .88 01 q .73 2- 1^ .62 0 50 100 Trial 150 0 50 100 Trial 150 Bibliography [1] Bowe, C.A. (1984). Spatial relations in animal learning and behavior. Psychological Record, 34 181-209. [2] Cori, D.F. (1989). Floral color change in Lupinus argenteus (Fabaceae): why should plants advertise the location of unrewarding flowers to pollinators? Evolution, 43 870-881. [3] Greenhouse, S.W., and Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24 95-112. [4] Huynh, H., and Feldt, L.S. (1970). Conditions under which mean square ratios in repeated measurements designs have exact F-distributions, Journal of the American Statistical Association, 65 1582-1589. [5] Liang, Kung-Yee, and Zeger Scott L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73 13-22. [6] Mackintosh, N.J. (1983). Conditioning and Associative Learning. Oxford: Clarendon Press. [7] McCullagh, P. (1983). Quasi-likelihood functions. Annals of Statistics, 11 59-67. [8] McCullagh, P., and Nelder, J.A., (1989). Generalized Linear Models (2nd edition). Chapman and Hall, London. [9] Morrison, Donald F. (1976). Multivariate Statistical Methods. McGraw-Hill, New York. [10] Wedderburn, R.W.M. (1974). Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method. Biometrika, 61 439-447. [11] Zeger, Scott L., and Liang, Kung-Yee (1986). Longitudinal data analysis for discrete and continuous outcomes. Biometrics, 42 121-130. [12] Zeger, Scott L., Liang, Kung-Yee, and Self, Steven G. (1985). The analysis of binary longitudinal data with time-independent covariates, Biometrika, 72 31-78. 117
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Longitudinal analysis for binary and count data
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Longitudinal analysis for binary and count data Graham, Jinko 1991
pdf
Page Metadata
Item Metadata
Title | Longitudinal analysis for binary and count data |
Creator |
Graham, Jinko |
Date | 1991 |
Date Issued | 2008-09-05T16:55:14Z |
Description | Longitudinal data sets consist of repeated observations of an outcome over time, and a corresponding set of covariates for each of many subjects. In many fields, multivariate analysis-of-variance is commonly used to analyse longitudinal data. Such an analysis is appropriate when responses for each subject are multivariate Gaussian with a common covariance matrix for all subjects. In many cases, however, the longitudinal response cannot be transformed to satisfy these assumptions. An alternate analysis might rely on specification of an appropriate likelihood to obtain estimates of regression parameters and their standard errors. Unfortunately, the correlation structure in the data for each subject may not be well understood, making such parametric modelling difficult. This thesis discusses two methods for the analysis of longitudinal data which require only minimal assumptions about the true correlation structure in the data for each subject to yield consistent estimates of regression parameters and their standard errors. The first method is based on the use of a "working" likelihood and extends the results of Zeger, Liang and Self (Biometrika, 1985) to the case of time-dependent covariates. The second method, first presented by Liang and Zeger (Biometrika, 1986), is based on quasi-likelihood theory. This method uses generalized estimating equations to arrive at consistent estimates of the regression coefficients and their standard errors, and can be applied to any longitudinal response with univariate marginal distributions for which the quasi-likelihood formulation is sensible. This includes Gaussian, Poisson, binomial (bernoulli), gamma, and inverse Gaussian distributions. Both methods are extensively illustrated using the results from an experiment on hummingbird learning. These methods enable much more information to be extracted from the hummingbird data set than a more traditional analysis-of-variance, and therefore provide useful and powerful tools for researchers in this subject area. |
Extent | 4759748 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2008-09-05 |
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.0086470 |
Degree |
Master of Science - MSc |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/1647 |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_1992_spring_graham_jinko.pdf [ 4.54MB ]
- Metadata
- JSON: 1.0086470.json
- JSON-LD: 1.0086470+ld.json
- RDF/XML (Pretty): 1.0086470.xml
- RDF/JSON: 1.0086470+rdf.json
- Turtle: 1.0086470+rdf-turtle.txt
- N-Triples: 1.0086470+rdf-ntriples.txt
- Original Record: 1.0086470 +original-record.json
- Full Text
- 1.0086470.txt
- Citation
- 1.0086470.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 25 | 3 |
United States | 6 | 0 |
Russia | 3 | 0 |
Japan | 3 | 0 |
Canada | 2 | 1 |
Thailand | 1 | 0 |
City | Views | Downloads |
---|---|---|
Beijing | 22 | 2 |
Ashburn | 6 | 0 |
Saint Petersburg | 3 | 0 |
Shenzhen | 3 | 1 |
Tokyo | 3 | 0 |
Halifax | 2 | 1 |
Unknown | 1 | 3 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
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-0086470/manifest