‘Counting chickens and other tales’ Using random effect models and MCMC estimation in applied statistics research.

1. ‘Counting chickens and other tales’ Using random effect models and MCMC estimation in applied statistics research. Dr William J. Browne School of Mathematical Sciences University of Nottingham

2. Outline • Background to my research, random effect and multilevel models and MCMC estimation. • Random effect models for complex data structures including artificial insemination and Danish chicken examples. • Multivariate random effect models and great tit nesting behaviour. • Sample size calculations for complex random effect models. • Multilevel modelling of mastitis incidence. • Other research and future work.

3. Background • 1995-1998 – PhD in Statistics, University of Bath. “Applying MCMC methods to multilevel models.” • 1998-2003 – Postdoctoral research positions at the Centre for Multilevel Modelling at the Institute of Education, London. • 2003-2006 – Lecturer in Statistics at University of Nottingham. • 2006- Associate professor of Statistics at University of Nottingham. Research interests: Multilevel modelling, complex random effect modelling, applied statistics, Bayesian statistics and MCMC estimation.

4. Random effect models • Models that account for the underlying structure in the dataset. • Originally developed for nested structures (multilevel models), for example in education, pupils nested within schools. • An extension of linear modelling with the inclusion of random effects. A typical 2-level model is Here i might index pupils and j index schools. Alternatively in another example i might index cows and j index herds. The important thing is that the model and statistical methods used are the same!

5. Estimation Methods for Multilevel Models Due to additional random effects no simple matrix formulae exist for finding estimates in multilevel models. Two alternative approaches exist: • Iterative algorithms e.g. IGLS, RIGLS, that alternate between estimating fixed and random effects until convergence. Can produce ML and REML estimates. • Simulation-based Bayesian methods e.g. MCMC that attempt to draw samples from the posterior distribution of the model. One possible computer program to use for multilevel models which incorporates both approaches is MLwiN.

6. MLwiN • Software package designed specifically for fitting multilevel models. • Developed by a team led by Harvey Goldstein and Jon Rasbash at the Institute of Education in London over past 15 years or so. Earlier incarnations ML2, ML3, MLN. • Originally contained ‘classical’ estimation methods (IGLS) for fitting models. • MLwiN launched in 1998 also included MCMC estimation. • My role in the team was as developer of the MCMC functionality in MLwiN in my time at Bath and during 4.5 years at the IOE. Note: MLwiN core team relocated to Bristol in 2005.

7. MCMC Algorithm • Consider the 2-level normal response model • MCMC algorithms usually work in a Bayesian framework and so we need to add prior distributions for the unknown parameters. • Here there are 4 sets of unknown parameters: • We will add prior distributions

8. MCMC Algorithm (2) One possible MCMC algorithm for this model then involves simulating in turn from the 4 sets of conditional distributions. Such an algorithm is known as Gibbs Sampling. MLwiN uses Gibbs sampling for all normal response models. Firstly we set starting values for each group of unknown parameters, Then sample from the following conditional distributions, firstly To get .

9. MCMC Algorithm (3) We next sample from to get , then to get , then finally To get . We have then updated all of the unknowns in the model. The process is then simply repeated many times, each time using the previously generated parameter values to generate the next set

10. Burn-in and estimates Burn-in: It is general practice to throw away the first n values to allow the Markov chain to approach its equilibrium distribution namely the joint posterior distribution of interest. These iterations are known as the burn-in. Finding Estimates: We continue generating values at the end of the burn-in for another m iterations. These m values are then averaged to give point estimates of the parameter of interest. Posterior standard deviations and other summary measures can also be obtained from the chains.

11. So why use MCMC? • Often gives better (in terms of bias) estimates for non-normal responses (see Browne and Draper, 2006). • Gives full posterior distribution so interval estimates for derived quantities are easy to produce. • Can easily be extended to more complex problems as we will see next. • Potential downside 1: Prior distributions required for all unknown parameters. • Potential downside 2: MCMC estimation is much slower than the IGLS algorithm. • For more information see my book: MCMC Estimation in MLwiN – Browne (2003).

12. School S1 S2 S3 S4 Pupil P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 Nbhd N1 N2 N3 Extension 1: Cross-classified models For example, schools by neighbourhoods. Schools will draw pupils from many different neighbourhoods and the pupils of a neighbourhood will go to several schools. No pure hierarchy can be found and pupils are said to be contained within a cross-classification of schools by neighbourhoods:

13. Notation With hierarchical models we use a subscript notation that has one subscript per level and nesting is implied reading from the left. For example, subscript pattern ijk denotes the i’th level 1 unit within the j’th level 2 unit within the k’th level 3 unit. If models become cross-classified we use the term classification instead of level. With notation that has one subscript per classification, that also captures the relationship between classifications, notation can become very cumbersome. We propose an alternative notation introduced in Browne et al. (2001) that only has a single subscript no matter how many classifications are in the model.

14. i nbhd(i) sch(i) 1 1 1 2 2 1 3 1 1 4 2 2 5 1 2 6 2 2 School S1 S2 S3 S4 7 2 3 Pupil P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 8 3 3 Nbhd N1 N2 N3 9 3 4 10 2 4 11 3 4 12 3 4 Single subscript notation We write the model as where classification 2 is neighbourhood and classification 3 is school. Classification 1 always corresponds to the classification at which the response measurements are made, in this case pupils. For pupils 1 and 11 equation (1) becomes:

15. Neighbourhood School Pupil Classification diagrams In the single subscript notation we lose information about the relationship (crossed or nested) between classifications. A useful way of conveying this information is with the classification diagram. Which has one node per classification and nodes linked by arrows have a nested relationship and unlinked nodes have a crossed relationship. School Neighbourhood Pupil Cross-classified structure where pupils from a school come from many neighbourhoods and pupils from a neighbourhood attend several schools. Nested structure where schools are contained within neighbourhoods

16. Donor Donation Woman Cycle Example : Artificial insemination by donor 1901 women 279 donors 1328 donations 12100 ovulatory cycles response is whether conception occurs in a given cycle In terms of a unit diagram: Or a classification diagram:

17. Parameter Description Estimate(se) intercept -4.04(2.30) azoospermia 0.22(0.11) semen quality 0.19(0.03) womens age>35 -0.30(0.14) sperm count 0.20(0.07) sperm motility 0.02(0.06) insemination to early -0.72(0.19) insemination to late -0.27(0.10) women variance 1.02(0.21) donation variance 0.644(0.21) donor variance 0.338(0.07) Model for artificial insemination data We can write the model as Results: Note cross-classified models can be fitted in IGLS but are far easier to fit using MCMC estimation.

18. Extension 2: Multiple membership models • When level 1 units are members of more than one higher level unit we describe a model for such data as a multiple membership model. • For example, •  Pupils change schools/classes and each school/class has an effect on pupil outcomes. • Patients are seen by more than one nurse during the course of their treatment.

19. n1(j=1) n2(j=2) n3(j=3) p1(i=1) 0.5 0 0.5 p2(i=2) 1 0 0 p3(i=3) 0 0.5 0.5 p4(i=4) 0.5 0.5 0 Here patient 1 was seen by nurse 1 and 3 but not nurse 2 and so on. If we substitute the values of w(2)i,j, i and j. from the table into (2) we get the series of equations : Notation Note that nurse(i) now indexes the set of nurses that treat patient i and w(2)i,jis a weighting factor relating patient i to nurse j. For example, with four patients and three nurses, we may have the following weights:

20. nurse Here patients are multiple members of nurses, nurses are nested within hospitals and GP practice is crossed with both nurse and hospital. patient hospital nurse GP practice patient Classification diagrams for multiple membership relationships Double arrows indicate a multiple membership relationship between classifications. We can mix multiple membership, crossed and hierarchical structures in a single model.

21. Farm House Parent flock Child flock Example involving nesting, crossing and multiple membership – Danish chickens Production hierarchy 10,127 child flocks 725 houses 304 farms Breeding hierarchy 10,127 child flocks 200 parent flocks As a unit diagram: As a classification diagram:

22. Parameter Description Estimate(se) intercept -2.322(0.213) 1996 -1.239(0.162) 1997 -1.165(0.187) Results: hatchery 2 -1.733(0.255) hatchery 3 -0.211(0.252) hatchery 4 -1.062(0.388) parent flock variance 0.895(0.179) house variance 0.208(0.108) farm variance 0.927(0.197) Model and results Response is cases of salmonella Note multiple membership models can be fitted in IGLS and this model/dataset represents roughly the most complex model that the method can handle. Such models are far easier to fit using MCMC estimation.

23. Random effect modelling of great tit nesting behaviour • An extension of cross-classified models to multivariate responses. • Collaborative research with Richard Pettifor (Institute of Zoology, London), and Robin McCleery and Ben Sheldon (University of Oxford).

24. Wytham woods great tit dataset • A longitudinal study of great tits nesting in Wytham Woods, Oxfordshire. • 6 responses : 3 continuous & 3 binary. • Clutch size, lay date and mean nestling mass. • Nest success, male and female survival. • Data: 4165 nesting attempts over a period of 34 years. • There are 4 higher-level classifications of the data: female parent, male parent, nestbox and year.

25. Source Number of IDs Median #obs Mean #obs Year 34 104 122.5 Nestbox 968 4 4.30 Male parent 2986 1 1.39 Female parent 2944 1 1.41 Data background The data structure can be summarised as follows: Note there is very little information on each individual male and female bird but we can get some estimates of variability via a random effects model.

26. Diagrammatic representation of the dataset.

27. Univariate cross-classified random effect modelling • For each of the 6 responses we will firstly fit a univariate model, normal responses for the continuous variables and probit regression for the binary variables. For example using notation of Browne et al. (2001) and letting response yibe clutch size:

28. Estimation • We use MCMC estimation in MLwiN and choose ‘diffuse’ priors for all parameters. • We run 3 MCMC chains from different starting points for 250k iterations each (500k for binary responses) and use the Gelman-Rubin diagnostic to decide burn-in length. • We compared results with the equivalent classical model using the Genstat software package and got broadly similar results. • We fit all four higher classifications and do not consider model comparison.

29. Clutch Size Here we see that the average clutch size is just below 9 eggs with large variability between female birds and some variability between years. Male birds and nest boxes have less impact.

30. Lay Date (days after April 1st) Here we see that the mean lay date is around the end of April/beginning of May. The biggest driver of lay date is the year which is probably indicating weather differences. There is some variability due to female birds but little impact of nest box and male bird.

31. Nestling Mass Here the response is the average mass of the chicks in a brood at 10 days old. Note here lots of the variability is unexplained and both parents are equally important.

32. Human example Helena Jayne Browne Born 22nd May 2006 Birth Weight 8lb 0oz Sarah Victoria Browne Born 20th July 2004 Birth Weight 6lb 6oz Father’s birth weight 9lb 13oz, Mother’s birth weight 6lb 8oz

33. Nest Success Here we define nest success as one of the ringed nestlings captured in later years. The value 0.01 for β corresponds to around a 50% success rate. Most of the variability is explained by the Binomial assumption with the bulk of the over-dispersion mainly due to yearly differences.

34. Male Survival Here male survival is defined as being observed breeding in later years. The average probability is 0.334 and there is very little over-dispersion with differences between years being the main factor. Note the actual response is being observed breeding in later years and so the real probability is higher!

35. Female survival Here female survival is defined as being observed breeding in later years. The average probability is 0.381 and again there isn’t much over-dispersion with differences between nestboxes and years being the main factors.

36. Multivariate modelling of the great tit dataset • We now wish to combine the six univariate models into one big model that will also account for the correlations between the responses. • We choose a MV Normal model and use latent variables (Chib and Greenburg, 1998) for the 3 binary responses that take positive values if the response is 1 and negative values if the response is 0. • We are then left with a 6-vector for each observation consisting of the 3 continuous responses and 3 latent variables. The latent variables are estimated as an additional step in the MCMC algorithm and for identifiability the elements of the level 1 variance matrix that correspond to their variances are constrained to equal 1.

37. Multivariate Model Here the vector valued response is decomposed into a mean vector plus random effects for each classification. Inverse Wishart priors are used for each of the classification variance matrices. The values are based on considering overall variability in each response and assuming an equal split for the 5 classifications.

38. Use of the multivariate model • The multivariate model was fitted using an MCMC algorithm programmed into the MLwiN package which consists of Gibbs sampling steps for all parameters apart from the level 1 variance matrix which requires Metropolis sampling (see Browne 2006). • The multivariate model will give variance estimates in line with the 6 univariate models. • In addition the covariances/correlations at each level can be assessed to look at how correlations are partitioned.

39. Partitioning of covariances

40. Correlations from a 1-level model • If we ignore the structure of the data and consider it as 4165 independent observations we get the following correlations: Note correlations in bold are statistically significant i.e. 95% credible interval doesn’t contain 0.

41. Correlations in full model Key: Blue +ve, Red –ve: Y – year, N – nestbox, F – female, O - observation

42. Pairs of antagonistic covariances at different classifications There are 3 pairs of antagonistic correlations i.e. correlations with different signs at different classifications: LD & NM : Female 0.20 Observation -0.19 Interpretation: Females who generally lay late, lay heavier chicks but the later a particular bird lays the lighter its chicks. CS & FS : Female 0.48 Observation -0.20 Interpretation: Birds that lay larger clutches are more likely to survive but a particular bird has less chance of surviving if it lays more eggs. LD & FS : Female -0.67 Observation 0.11 Interpretation: Birds that lay early are more likely to survive but for a particular bird the later they lay the better!

43. Sample size calculations in random effect models • A current ESRC grant (2006-2009) that funds a postdoc (Mousa Golalizadeh). • The grant will consider the problem of deciding on how much data to collect for a research question taking account of the likely structure in the collected data (See later slides). • The grant will also investigate how various MCMC algorithm developments perform in practice when applied to real datasets. • Finally we will investigate when complex models are identifiable in the presence of ‘sparse’ data.

44. Background • Many quantitative research questions are of the form of a hypothesis – A has a significant effect on B. • To answer such a question data is collected that allows the researcher to (hopefully) test whether statistically A has a significant effect on B. (In fact we aim to reject the hypothesis that A doesn’t significantly affect B). • A test is performed and either the researcher is happy and A indeed has a significant effect on B or is left wondering why the data collected do not back up their hypothesis. Is the hypothesis false or was the data not sufficient? • The sufficiency of the data is the motivation for sample size calculations.

45. Example • Suppose I have the research question ‘Are Welshmen on average taller than 175 cms?’ • I now need to get hold of a random sample of n Welshmen and measure each of their heights. • I make some statistical assumption about the distribution of the heights of Welshmen e.g. that they come from a Normal distribution. • I might like to check this assumption by plotting a histogram of the data. • I can then form a statistical hypothesis test and test whether indeed Welshmen are taller than 175cms. • I need to decide how big to make n, my sample of Welshmen.

46. Hypothesis Testing • Let us assume our null hypothesis is that the average height of Welshmen (μ) is 175cm. • So we test H0:μ=175 vs HA:μ>175 (or alternatively H0:θ=0 vs HA:θ>0 where θ=μ-175) • In practice we calculate from our sample its mean ( ) and standard deviation (s2) and use these along with n to form a test statistic which we can compare with the distribution assumed under H0

47. Type I and Type II errors • No hypothesis test is perfect and there is always the possibility of errors • P(Type I error) = α = significance level or size • P(Type II error) = β, 1-β is the power of the test. • In general we fix α to some value e.g. 0.05, 0.01 then 1-β depends on our sample size.

48. Example hypothesis test • Let us assume that in reality our sample mean is 180cms and the population standard deviation (sd) is 5cms (known). • We can then form a test statistic as follows: • Note here that for small n and unknown sd we should use a student-t distribution rather than Normal. • For a 1-sided Z test we wish Z= > 1.645 and so we need our sample to be of size 3 to reject H0, using a student-t distribution increases this to 5. (Here α=0.05) • However if the sample mean had been only 176cms then we would need n > (1.645*5)2 = 68 Welshmen to reject H0

49. Power calculations • Our last slide in some sense is backwards as we cannot get from a given sample mean to choosing a sample size! • What we do instead is use different terminology and play God! • We will choose an ‘effect size’, γ which will represent a guess at the increase in the sample mean for Welshmen. • There then exists an (approximate) formula that links four quantities, size (α), power (1-β), effect size (γ) and sample size (n) • Note that the standard error (SE) of γ is a function of n and σ the population sd which is assumed known. • We can now evaluate one of these quantities conditional on the others e.g. what sample size is required given α,1-β and γ? Here RHS is sum of cases H0 true and H0 false.

50. Welsh height example Here we have looked at two examples with effect sizes 5 and 1 respectively. Assume σ takes the value 5 and so let us suppose we take a sample of size 25 Welshmen. Then Case 1: 5/(5/√25)=1.645+z1-β,z1-β=3.355 β=0.9996 Case 2: 1/(5/ √25)=1.645+z1-β,z1-β=-0.645 β=0.25946 So here a sample of 25 Welshmen from a population with mean 180cms would almost always result in rejecting H0, but if the population mean is 176cms then only 26% of such samples would be rejected. We can plot curves of how power increases with sample size as shown in the next slide.