1 / 52

Modelling Count Data: Outline

Modelling Count Data: Outline. Characteristics of count data and the Poisson distribution Applying the Poisson: Flying bomb strikes in South London Deaths by horse-kick : as a single-level model Poisson model fitted in MLwiN

eliza
Download Presentation

Modelling Count Data: Outline

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. Modelling Count Data: Outline • Characteristics of count data and the Poisson distribution • Applying the Poisson: Flying bomb strikes in South London • Deaths by horse-kick: as a single-level model Poisson model fitted in MLwiN • Overdispersion: types and consequences, the unconstrained Poisson, the Negative Binomial • Taking stock: 4 distributions for modeling counts • Number of extramarital affairs: the incidence rate ratio (IRR); handling categorical & continuous predictors; comparing model with DIC • Titanic survivor data; taking account of exposure, the offset • Multilevel Poisson and NBD models; estimation and VPC • Applications: HIV in India and Teenage employment in Glasgow • Spatial models: Lip cancer in Scotland; respiratory cancer in Ohio

  2. Some characteristics of count data • Very common in the social sciences • Number of children; Number of marriages • Number of arrests; Number of traffic accidents • Number of flows; Number of deaths • Counts have particular characteristics • Integers & cannot be negative • Often positively skewed; a ‘floor’ of zero • In practice often rare events which peak at 1,2 or 3 and rare at higher values • Modelled by • Logit regression models the log odds of an underlying propensity of an outcome; • Poisson regression models the log of the underlying rateof occurrence of a count.

  3. Theoretical Poisson distribution I • The Poisson distribution results if the underlying number of random events per unit time or space have a constant mean ()rate of occurrence, and each event is independent Simeon-Denis Poisson (1838) Research on the Probability of Judgments in Criminal and Civil Matters • Applying the Poisson: Flying bomb strikes in South London • Key research question: falling at random or under a guidance system • If random independent events should be distributed spatially as a Poisson distribution • Divide south London into 576 equally sized small areas (0.24km2) • Count the number of bombs in each area and compare to a Poisson • Mean rate =  = [229(0) + 211(1) + 93(2) + 35(3) + 7(4) +1(5)]/576 • = 0.929 hits per unit square • Very close fit; concluded random

  4. Theoretical Poisson distribution II Probability mass function(PMF) for 3 different mean occurrences • When mean =1; very positively skewed • As mean occurrence increases (more common event), distribution approaches Gaussian; • So use Poisson for ‘rareish’ events; mean below 10 • Fundamental property of the Poisson: mean = variance • Simulated 10,000 observations according to Poisson • Mean Variance Skewness • 0.993 0.98 1.00 • 4.001 4.07 0.50 • 10.03 10.38 0.33 • Variance is not a freely estimated parameter as in Gaussian

  5. Death by Horse-kick I: the data • Bortkewicz L von.(1898) The Law of Small Numbers, Leipzig • No of soldiers killed annually by horse-kicks in Prussian cavalry; 10 corps over 20 years (occurrences per unit time) • The full data 200 corps years of observations • As a frequency distribution (grouped data) • Deaths 0 1 2 3 4 5+ • Frequency 109 65 22 3 1 0 Mean Variance Number of obs • 0.61 0.611 200 • Interpretation: mean rate of 0.61 deaths per cohort year (ie rare) • Mean equals variance, therefore a Poisson distribution

  6. Death by Horse-kick II: as a Poisson Again :The Poisson results if underlying number of random events per unit time or space have a constant mean () rate of occurrence, and each event is independent • With a mean (and therefore a variance) of 0.61 • Deaths 0 1 2 3 4 5+ • Frequency 109 65 22 3 1 0 • Theory 109 66 20 4 1 1 • Formula for Poisson PMF • e is base of the natural logarithm (2.7182) •  is the mean (shape parameter): the average number of events in a given time interval • x! is the factorial of x EG: mean rate  of 0.6 accidents per corps year; what is probability of getting 3 accidents in a corps in a year?

  7. Horse-kick III: as a single level Poisson model General form of the single-level model Observed count is distributed as an underlying Poisson with a mean rate of occurrence of ; That is as an underlying mean and level-1 random term of z0 (the Poisson ‘weight’) Mean rate is related to predictors non-linearly as an exponential relationship Model loge to get a linear model (log link) The Poisson weight is the square root of estimated underlying count, re-estimated at each iteration Variance of level-1 residuals constrained to 1, • Modeling on the loge scale, cannot make prediction of a negative count on the raw scale • Level-1 variance is constrained to be an exact Poisson, (variance =mean)

  8. Horse-kick IV: Null single-level Poisson model in MLwiN • The raw ungrouped counts are modeled with a log link and a variance constrained to be equal to the mean • -0.494 is the mean rate of occurrence on the log scale • Exponentiate -0.494 to get the mean rate of 0.61 • 0.61 is interpreted as RATE; the number of events per unit time (or space), ie 0.61 horse-kick deaths per corps-year

  9. Overdispersion I: Types and consequences • So far; equi-dispersion, variances equal to the mean • Overdispersion: variance > mean; long tail, eg LOS (common) • Un-dispersion: variance < mean; data more alike than pure Poisson process; in multilevel possibility of missing level • Consequences of overdispersion • Fixed part SE’s "point estimates are accurate but they are not as precise as you think they are" • In multilevel, mis-estimate higher-level random part • Apparent and true overdispersion: thought experiment • : number of extra-marital affairs: men women with different means • apparent: mis-specified fixed part, not separated out distributions with different means • true: genuine stochastic property of more inherent variability • in practice model fixed part as well as possible, and allow for overdisperion

  10. Overdispersion II: the unconstrained Poisson • Deaths by horse-kick • estimate an over- dispersed Poisson • allow the level-1 variance to be estimated • Not significantly different from 1; • No evidence that this is not a Poisson distribution

  11. Overdispersion III: the Negative binomial • Instead of fitting an overdispersed Poisson, could fit a NBD model • Handles long-tailed distributions; • An explicit model in which variance is greater than the mean • Can even have an over-dispersed NBD • Same log-link but NBD has 2 parameters for the level-1 variance; that is quadratic level-1 variance, v is the overdispersion parameter

  12. Overdispersion IV: the Negative binomial • Horse=kick analysis: • Null single-level NBD model; essentially no change, v is estimated to be 0.00 (see with Stored model; Compared stored model) • Overdispersed negative binomial; • No evidence of overdispersion; deaths are independent

  13. Linking the Binomial and the Poisson: I • First Bernoulli and Binomial • Bernoulli is a distribution for binary discrete events • y is observed outcome; ie 1 or 0 • E(y) = ; underlying propensity/probability for occurrence • Var(y) = /1-  • Binomial is a distribution for discrete events out of a number of trials • y is observed outcome; n is the number of trials, • E(y) = ; underlying propensity/ of occurrence • Var(y) = [/(1- )]/n • Least variation when denominator is large (more reliable), and as underlying probability approaches 0 or 1

  14. Linking the Binomial and the Poisson: II • Poisson is limit of a binomial process in which prob → 0, n→∞ • Poisson describes the probability that a random event will occur in a time or space interval when the probability is very small, the number of trials is very large, and each event is independent • EG : The probability that any automobile is involved in an accident is very small, but there are very many cars on a road in a day so that the distribution (if each crash is independent) follows a Poisson count • If non-independence of crashes (a ‘pile-up’), then over-dispersed Poisson/NBD, latter used for ‘contagious’ processes • In practice, Poisson and NBD used for rare occurrences, less than 10 cases per interval,& hundreds or even thousands for denominator/ trials [Clayton & Hills (1993)Statistical Models in Epidemiology OUP]

  15. Taking stock: 4 distributions for counts • If common rate of occurrence (mean >10) then use raw counts and Gaussian distribution (assess Normality assumption of the residuals) • If rare rate of occurrence, then use over-dispersed Poisson or NBD; the level-1 unconstrained variance estimate will allow assessment of departure from equi-dispersion; improved SE’s, but biased estimates if apparent overdispersion due to model mis-specification • Use the Binomial distribution if count is out of some total and the event is not rare; that is numerator and denominator of the same order • Mean variance relations for 4 different distributions that could be used for counts

  16. Modeling number of extra-marital affairs Single- level Poisson with Single categorical Predictor Extract of raw data (601 individuals from Fair 1978) Fair, R C(1978) A theory of extramarital affairs, Journal of Political Economy, 86(1), 45-61 • Single categorical predictor Children with NoKids, as the base • Understanding customized predictions…………….

  17. Single- level Poisson with Single categorical Predictor Understanding customized predictions • Log scale: • NoKids: -0.092 WithKids-0.092+0.606= 0.514 • First use equation to get underyling log-number of events then exponeniate to get estimated count (since married) • As mean/median counts • NoKids: expo(-0.092)= 0.91211 • Withkids: expo(-0.092 + 0.606) = 1.6720 • Those with children have a higher average rate of affairs (but have they been married longer?)

  18. Modeling number of extra marital affairs: the incidence rate ratio (IRR) So far: mean counts NoKids: expo(-0.092)= 0.91211 Withkids: expo(-0.092 + 0.606) = 1.6720 But also as IRR; comparing the ratio of those with and without kids IRR = 1.6720/0.91211 = 1.8331 That is Withkids have a 83% higher rate BUT can get this directly from the model by exponentiating the estimate for the contrasted category expo(0.606) = 1.8331 • Rules: • exponeniating the estimates for the (constant plus the contrasted category) gives the mean rate for the contrasted category • b) ) exponeniating the contrasted category gives IRR in comparison to base category

  19. Why is the exponentiated coefficient a IRR? • As always, the estimated coefficient is the change in response corresponding to a one unit change in the predictor • Response is underlying logged count • When Xi is 0 (Nokids); log count| Xi is 0 = β0 • But when Xi is 1(Withkids); log count| Xi is 1 = β0 + β1 X1i • Subtracting the first equation from the second gives • (log count| Xi is 1)-(log count| Xi is 0) = β1 • Exponentiating both sides gives (note the division sign) • (count| Xi is 1)/(count| Xi is 0) = exp(β1) • Thus, exp(β1) is a rate ratio corresponding to the ratio of the mean number of affairs for a with-child person to the mean number of affairs without-child person • Incidence: number of new cases • Rate: because it number of events per time or space • Ratio: because its is ratio of two rates

  20. Modeling number of extra-marital affairs: changing the base category • Previously contrasted category: Withkids: +0.606 • Now contrasted category: Nokids : -0.606 • Changing base simply produces a change of sign on the loge scale • Exponentiating the contrasted category: • Before: expo(+0.606) = 1.8331 • Now: expo(-0.606) = 0.5455 • Doubling the rate on loge scale is 0.693; Halving the rate on loge scale is -0.693 • IRR of 0.111 = IRR of 9-fold increase, difficult to appreciate • Advice: choose base category to be have the lowest mean rate; get positive contrasted estimates; always then comparing a larger value to a base of 1

  21. Affairs: modeling a set of categorical predictors • A model with ‘years married’ included with < 4 years as base Customised predictions: mean rate, IRR, & graph with 95% CI’s

  22. Affairs: modeling a continuous predictor Age Age as a 2nd order polynomial centred around 17 years (the youngest person in survey; also lowest rate • To get mean rate as it changes with age • Expo (-0.990 + 0.149*(Age-17) -0.003*(Age-17)2) • To get IRR in comparison for a person aged 33 compared to 17 • Expo( 0.149*16) –(0.003 *162) (drop the constant!) • Easiest interpreted as graphs!

  23. Affairs: a SET of predictors & models • Notice substantial overdispersion • Poisson & extra-Poisson no change in estimates; some change to NBD • Notice larger SE’s when allow for overdispersion; NBD most conservative • In full model, WithKids not significant

  24. NBD model for Marital Affairs • IRR of 1 for Under 4 years married, Very religious, No children, aged 17 • Previous Age effect is really length of marriage • Used comparable vertical axes, range of 4

  25. NBD model for number of Extra-Marital Affairs • With 95% confidence intervals • NB that they are asymmetric on the unlogged scale

  26. Affairs: Evaluating a sequence of models using DIC • Likelihood and hence the Deviance are not available for Poisson and NBD models fitted by quasi-likelihood • DIC criterion available though MCMC; typically needs larger number of iterations than Normal & Binomial (suggested default is 50k not 5k) • Currently MCMC not available in MLwiN for over-dispersed Poisson nor NBD models; so have to use Wald tests in Intervals and tests window

  27. Titanic survivor data: Taking account of exposure • So far, response is observed count, now we want to model a count given exposure: EG only 1 high-class female child survived but only 1 exposed! Here 2 possible measures of exposure a) the number of potential cases; could use a binomial b) the expected number if everyone had the same exposure (i indexes cell) Death rate = Total Deaths/ Total exposed =817/1316 = 0.379179 Expi = Casesi * Survival rate Latter often used to treat the exposure as a nuisance parameter & allows calculation of Standardised Rate’s SRi = Obsi/Expi * 100 Previous examples Horsekick: Exposure removed by design: 200 cohort years Affairs: included length of marriage; theoretically interesting

  28. Modeling SR’s: the use of the OFFSET Model: SRi = (Obsi/Expi) = F(Agei, Genderi, Classi) Where i is a cell, groups with same characteristics Aim: are observed survivors greater or less than expected, and how these ‘differences’ are related to a set of predictor variables? As a non-linear model: E(SRi)= E(Obsi/Expi) = As a linear model (division of raw data is subtraction of a log) Loge (Obsi) - Loge(Expi)) = As a model with an offset, moving Loge(Expi) to the right-hand side, and constraining coefficient to be 1; ie Exp becomes predictor variable Loge (Obsi) = 1.0* Loge(Expi) + NB MLwiN automatically loge transforms the observed response; you have to create the loge of the expected and declare it as an offset Sir John Nelder

  29. Surviving on the Titanic as a log-linear model • Include the offset • As a saturated model; ie Age *Gender*Class, (2*2*3), 12 terms for 12 cells • Make predictions on the loge scale (must include constant); exponentiate all terms to get departures from the expected rate, that is modeled SR’s

  30. Titanic survival: parsimonious model • Remove insignificant terms starting with 3-way interactions for High*women*children Customized predictions: Very low rates of survival for Low and Middle class adult men; large gender gap for adults, but not for children

  31. Titanic survival: parsimonious model • Modeled SR’s and descriptive SR’s • Ordered by worse survival • Estimated SR’s only shown if 95% CI’s do not include 1.0

  32. Two-level multilevel Poisson One new term, the level 2 differential, on the loge scale, is assumed to come from Normal distribution with a variance of • Can also fit Poisson multilevel with offset and NBD multilevel in MLwiN

  33. Estimation of multilevel Poisson and NBD in MLwiN I • Same options as for binary and binomial • Quasi-likelihood and therefore MQL or PQL fitted using IGLS/RIGLS; fast, but no deviance (have to use Wald tests); may be troubled by small number of higher-level units; simulations have shown that MQL tends to overestimate the higher-level variance parameters • MCMC estimates; good quality and can use DIC to compare Poisson models; but currently MCMC is not possible for extra-Poisson nor for NBD • MCMC in MLwiN often produces highly correlated chains (in part due to the fact that the parameters of the model are highly correlated; variance =mean) Therefore requires substantial number of simulations; typically much larger than for Normal or for Binomial

  34. Estimation of multilevel Poisson and NBD in MLwiN II • Possibility to output to WinBUGS and use the univariate AR sampler and Gamerman (1997) method which tends to have less correlated chains, but WinBUGS is considerably slower generally [Gamerman, D. (1997) Sampling from the posterior distribution in generalized linear mixed models. Statistics and Computing 7, 57-68] • Advice: start with IGLS PQL; switch to MCMC, be prepared to make 500,000 simulations (suggest use 1 in 10 thinning to store the chains); use Effective sample size to assess required length of change, eg need ESS of at least 500 for key parameters of interest; compare results and contemplate using PQL and over-dispersed Poisson • Freely available software MIXPREG for multilevel Poisson counts including offsets; uses full information maximum likelihood estimated using quadrature http://tigger.uic.edu/~hedeker/mixpcm.PDF

  35. VPC for Poisson models • Can either use • Simulation method to derive VPC (modify the binomial procedure) • Use exact method: http://people.upei.ca/hstryhn/iccpoisson.ppt (Henrik Stryhn) • VPC for two level random intercepts model (available for other models) Clearly VPC depends on  and

  36. Modeling Counts in MLwiN: HIV in India Aim: investigate the State geography of HIV in terms of riskData: nationally representative sample of 100k individuals in 2005- 2006Response: HIV sero-status from blood samplesStructure: 1720 cells within 28 States; cells are a group of people who share common characteristics [Age-Groups(4), Education(4), Sex(2), Urbanity(2) and State (28)]Rarity: only 467 sero-positives were foundModel: Log count of number of seropositives in a cell related to an offset of Log expected count if national rates appliedPredictors of Age, Sex and Education and UrbanityTwo-level multilevel Poisson, extra-Poisson & NBD

  37. HIV in India: Standardized Morbidity Rates Higher educated females have the lowest risk, across the age-groups

  38. HIV in India: some results Risks for different States relative to living in urban and rural areas nationally.

  39. Modeling proportions as a binomial in MLwiN • exactly the same procedure as for binary models • except that observed y is a proportion (not just 1 and 0, the denominator (n) is variable (not just 1) and extra-dispersion at level 1 is allowed (not just exact binomial) • Reading: • Subramanian S V, Duncan C, Jones K (2001) Multilevel perspectives on modeling census data Environment and Planning A33(3) 399 – 417

  40. Data: teenage employment In Glasgow districts • “Ungrouped” data that is individual data • Model binary outcome of employed or not and two individual predictors

  41. Same data as a multilevel structure: a set of tables for each district GENDER QUALIF MALE FEMALE Postcode UnErate LOW 5 out of 6 3 out of 12 G1A 15% HIGH 2 out of 7 7 out of 9 LOW 5 out of 9 7 out of 11 G1B 12% HIGH 8 out of 8 7 out of 9 LOW 3 out of 3 - G99Z 3% HIGH 2 out of3 out of 5 • Level 1 cell in table • Level 2: Postcode sector • Margins: define the two categorical predictors • Internal cells: the response of 5 out of 6 are employed

  42. Teenage unemployment: some results from a binomial, two-level logit model

  43. A B Person in A Affected by A(SH) and B,C,D (MM) C D E F Person in H Affected by H(SH) and E,I,K,G (MM) G H I J K L M Multiple membership defined by common boundary; weights as function of inverse distance between centroids of areas Spatial Models as a combination of strict hierarchy and multiple membership: counts are commonly used

  44. Scottish Lip Cancer Spatial multiple-membership model • Response: observed counts of male lip cancer for the 56 regions of Scotland (1975-1980) • Predictor: % of workforce working in outdoor occupations (Agric;For; Fish) Expected count based on population size • Structure areas and their neighbours defined as having a common border (up to 11); equal weights for each neighbouring region that sum to 1 • Rate of lip cancer in each region is affected by both the region itself and its nearest neighbours after taking account of outdoor activity • Model Log of the response related to fixed predictor, with an offset, Poisson distribution for counts; • NB Two sets of random effects • 1 area random effects; (ie unstructured; non-spatial variation); • 2 multiple membership set of random effects for the neighbours of each region

  45. MCMC estimation: 50,000 draws Poisson model Fixed effects: Offset and Well-supported + relation Well-supported Residual neighbourhood effect NB: Poisson highly correlated chains

  46. Scottish Lip Cancer: CAR model CAR: CAR one set of random effects, which have an expected value of the average of the surrounding random effects; weights divided by the number of neighbours where ni is the number of neighbours for area i and the weights are typically all 1 • MLwiN: limited capabilities for CAR model; ie at one level only (unlike Bugs)

  47. MCMC estimation: CAR model, 50,000 draws Poisson model Fixed effects: Offset and Well-supported + relation Well-supported Residual neighbourhood effect

  48. NB Scales: shrinkage

  49. Ohio cancer: repeated measures (space and time!) • Response: counts of respiratory cancer deaths in Ohio counties • Aim: Are there hotspot counties with distinctive trends? (small numbers so ‘borrow strength’ from neighbours) • Structure: annual repeated measures (1979-1988) for counties • Classification 3: n’hoods as MM (3-8 n’hoods) • Classification 2: counties (88) • Classification 1: occasion (88*10) • Predictor: Expected deaths; Time • Model Log of the response related to fixed predictor, with an offset, Poisson distribution for counts (C1); • Two sets of random effects • 1 area random effects allowed to vary over time; trend for each county from the Ohio distribution (c2) • 2 multiple membership set of random effects for the neighbours of each region (C3)

  50. MCMC estimation: repeated measures model, 50,000 draws General trend N’hood variance Variance function for between county time trend Default priors

More Related