1、Multilevel Modeling of Categorical DataSophia Rabe-HeskethGraduate School of Education & Graduate Group in BiostatisticsUniversity of California, BerkeleyInstitute of Education, University of LondonJoint work with Anders SkrondalDivision of EpidemiologyNorwegian Institute of Public HealthMethods of
2、Analysis Program in the Social SciencesStanford University, November 2008GLLAMM p.1Outlinetrianglerightsld I. Ordered categorical responses Latent response and generalized linear model formulations Two-level random intercept models Conditional versus marginal effects Multilevel random coefficient mo
3、dels Estimation methods Example: Cluster randomized trial of sex education in Norwaytrianglerightsld II. Unordered categorical responses Random utility formulation Multilevel random utility models Example: Abuse of antibiotics in ChinaGLLAMM p.2I. Ordered categorical responsestrianglerightsld Small
4、number of mutually exclusive categories, yi=1,.,Strianglerightsld Categories are ordered, so yiis ordinaltrianglerightsld Examples: Severity of symptom (e.g. pain): none, moderate, severe Frequency of behavior: never, occasionally, daily Agreement (Likert scale):disagree strongly, disagree, agree, a
5、gree strongly Educational attainment:high school, college degree, graduate or professional degreeGLLAMM p.3Latent response modelstrianglerightsld Latent (unobserved) continuous response yiunderlies observedordinal response yitrianglerightsld Threshold model determines observed response:yi=1 if yi 12
6、 if 1|xi)0.521xiyiGLLAMM p.5Cumulative versus individual response probabilitiesPr(yi1) = Pr(yi=2)+Pr(yi=3)Pr(yi2) = Pr(yi=3)0.00.20.40.60.81.0Pr(yi1|xi)Pr(yi2|xi)xiCumulativeprobabilitiesPr(yi=3) = Pr(yi2)Pr(yi=2) = Pr(yi1) Pr(yi2)Pr(yi=1) = 1 Pr(yi1)0.00.20.40.60.81.0Pr(yi=1|xi)Pr(yi=2|xi)Pr(yi=3|x
7、i)xiIndividualresponseprobabilitiesGLLAMM p.6From latent response togeneralized linear model formulationstrianglerightsld Cumulative probabilities:Pr(yis|xi)=Pr(yis|xi)=Pr(xprimei + epsilon1is|xi)= Pr(epsilon1i xprimei s|xi)=F(xprimei s) F() is the cumulative density function of epsilon1itrianglerig
8、htsld This is a generalized linear model where F() is inverse link function and xprimei sis linear predictorDistribution of epsilon1i Var(epsilon1i) Link g = F1ModelLogistic 2/3 Logit Prop. oddsStandard normal 1 Probit Ordinal probitGumbel 2/6 Compl. log-log Compl. log-logAsymmetric (non-equivalent
9、model if scale reversed)GLLAMM p.7Logit link and odds ratiostrianglerightsld With a logistic distribution for epsilon1i,alogit link (one covariate xi):lnbraceleftbiggPr(yis|xi)1 Pr(yis|xi)bracerightbiggbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightodds(yis|xi)= xi slnodds(yis|x
10、i= a +1)lnodds(yis|xi= a) = = lnbraceleftbiggodds(yis|xi= a +1)odds(yis|xi= a)bracerightbiggtrianglerightsld Coefficient represents increase in log-odds associated with unitincrease in xitrianglerightsld Exponentiated coefficient exp() represents odds ratio associatedwith unit increase in xiGLLAMM p
11、.8Two-level random intercept modelstrianglerightsld Subjects i nested in clusters j (e.g. schools)trianglerightsld Include random intercept jfor clustersyij= xprimeij + j+ epsilon1ij,j|xij N(0,) jindependent of epsilon1ij jand epsilon1ijindependent across clusters epsilon1ijindependent across units
12、within clustersGLLAMM p.9Intraclass correlation of latent responsestrianglerightsld Total residual ij= j+ epsilon1ijhas varianceVar(ij|xij)= +1 for probit models + 2/3 for logit modelstrianglerightsld Covariance between total residuals ijand iprimejof two subjects insame cluster is and intraclass co
13、rrelation isCor(ij,iprimej|xij,xiprimej)=/( +1) for probit models/( + 2/3) for logit modelsGLLAMM p.10Conditional versus marginal effectstrianglerightsld For probit model, conditional or cluster-specific probabilities arePr(yijs|j,xij)=(xprimeij + j s)trianglerightsld Marginal or population-averaged
14、 response probabilities arePr(yijs|xij)=Pr(yijs)=Pr(xprimeij + ijs)= Pr(ij xprimeij s)= Prparenleftbiggij +1xprimeij s +1parenrightbigg=parenleftbiggxprimeij s +1parenrightbiggtrianglerightsld Therefore, marginal effects / +1attenuated compared withconditional effects GLLAMM p.11Illustration:Conditi
15、onal versus marginal relationship0 20406080100.00.20.40.60.81.0xProbabilitycluster-specific (random sample)medianmarginal or population-averagedGLLAMM p.12Multilevel random coefficient modelstrianglerightsld Consider clustered longitudinal data with occasions i (level 1) nestedin student j (level 2)
16、 in schools k (level 3)trianglerightsld Example of three-level random coefficient model:yijk=(2)0jk+ (3)0k+1+ (2)1jk+ (3)1kx1ijk+ 2x2ijk+ epsilon1ijk (2)0jkand (3)0kare random intercepts at levels 2 and 3 (2)1jkand (3)1kare random coefficients of x1ijkat levels 2 and 3 Random effects at the same lev
17、el, e.g. (2)0jk,(2)1jk), bivariatenormal with zero means, independent across units jktrianglerightsld General three-level random coefficient modelyijk= xprimeijk + z(2)primeijk(2)jk+ z(3)primeijk(3)k+ epsilon1ijkGLLAMM p.13Estimation: Approximate methodstrianglerightsld Estimation of multilevel mode
18、ls with categorical responses difficultbecause likelihood does not generally have closed formtrianglerightsld Penalized Quasilikelihood (PQL) Two versions: First and second order (PQL-1,PQL-2), the latterbeing betterdiamondmath PQL-1 in GLIMMIX in SAS, MLwiN and HLMdiamondmath PQL-2 in MLwiNdiamondm
19、ath Even PQL-2 produces biased estimates for small clusters andhigh intraclass correlationstrianglerightsld Sixth order Laplace in HLMtrianglerightsld H-likelihood in Genstattrianglerightsld Methods do not provide a likelihoodGLLAMM p.14Estimation: Maximum likelihoodtrianglerightsld Numerical integr
20、ation Gauss-Hermite (ordinary) quadrature used in MIXOR/MIXNO(two-level only), aML, NLMIXED in SAS (two-level only) andgllamm in Stata Adaptive quadrature superior to ordinary quadrature, particularlyfor large clusters and high intraclass correlations and available inNLMIXED in SAS (two-level only),
21、 gllamm in Stata, glme inS-PLUStrianglerightsld Monte Carlo integration Simulated maximum likelihood in nlogit Monte Carlo EM - no software?trianglerightsld MCMC with vague priors approximates maximum likelihood andavailable in MLwiN and WinBUGSGLLAMM p.15Example:Cluster randomized trial of sex educ
22、ation in Norwaytrianglerightsld Schools randomized to receive special sex education or nottrianglerightsld Assessments before, 6 and 18 months after randomizationtrianglerightsld One outcome is question relating to Contraceptive self-efficacy: “If my partner and I were about to have intercourse with
23、out eitherof us having mentioned contraception, it would be easy for me toproduce a condom (if I brought one)” Question answered in terms of five ordinal categories:not at all true of me, slightly true of me, somewhat true of me,mostly true of me, completely true of metrianglerightsld Multilevel dat
24、a with responses at occasions (level 1) from 1184students (level 2) in 46 schools (level 3)trianglerightsld Only 570 students always responded, 400 responded on someoccasions and 114 never respondedGLLAMM p.16Models and estimation using gllammtrianglerightsld Occasions t, students j, schools ktriang
25、lerightsld Covariates x1tTime (0, 1, 3) x2jkTreat (yes=1,no=0) x3tjkTreat Timetrianglerightsld Model the probability of exceeding a category s, s =1, 2, 3, 4logitPr(ytjks|(2)jk,(3)k,xtjk) = 1x1t+2x2jk+3x3tjk+(2)jk+(3)kstrianglerightsld Estimation using adaptive quadrature in gllamm:gllamm use treat
26、time treat_time, i(id school) /link(ologit) family(binom) adaptGLLAMM p.17EstimatesSingle-level model Two-level model Three-level modelParameter Est (SE) Est (SE) Est (SE)1Time -0.12 (0.06) -0.13 (0.06) -0.13 (0.06)2Treat -0.05 (0.14) -0.02 (0.19) -0.02 (0.19)3TimeTreat 0.17 (0.08) 0.17 (0.09) 0.17
27、(0.09)Var(2)jk) 2.03 (0.31) 2.03 (0.31)Var(3)k) 0.00 (0.00)1-3.54 (0.17) -4.41 (0.23) -4.41 (0.23)2-2.43 (0.13) -3.15 (0.19) -3.15 (0.19)3-1.18 (0.12) -1.58 (0.16) -1.58 (0.16)40.16 (0.12) 0.25 (0.15) 0.25 (0.15)Log-likelihood -2531 -2471 -2471GLLAMM p.18Predicted probabilitiestrianglerightsld Condi
28、tional probabilities:Pr(ytjk 2|(2)jk=0,(3)k=0,xtjk)=exp(xprimetjkhatwide +0+0hatwide2)1+exp(xprimetjkhatwide +0+0hatwide2)gen u1 = 0gen u2 = 0gllapred p_cond, mu us(u) above(2)trianglerightsld Marginal probabilities:Pr(ytjk 2|xtjk)=integraldisplayintegraldisplayexp(xprimetjkhatwide + (2)jk+ (3)khatw
29、ide2)1+exp(xprimetjkhatwide + (2)jk+ (3)khatwide2)g(2)jk)g(3)k)d(2)jkd(3)kgllapred p_marg, mu marg above(2)GLLAMM p.19Predicted conditional probabilitiestrianglerightsld Probability of responding at least mostly true of me (category 3)trianglerightsld Relationship between probability and occasion fo
30、r two treatmentgroups, when (3)k=0 and (2)jkis -1, 0 and 1:intervention group, control groupTime in monthsAt least mostly true0.60.70.80.906 18GLLAMM p.20Observed proportionsand predicted marginal probabilitiesintervention group, control group, predicted, observedPr(y1) Pr(y2)Time in monthsAt least
31、slightly true0.880.920.961.0006 18Time in monthsAt least somewhat true0.840.880.920.9606 18Pr(y3) Pr(y4)Time in monthsAt least mostly true0.680.720.760.8006 18Time in monthsCompletely true0.360.400.440.4806 18GLLAMM p.21II. Unordered categorical responses or discrete choicetrianglerightsld Small num
32、ber of mutually exclusive categories a, a =1,.,A.trianglerightsld Categories cannot be ordered a prioritrianglerightsld Examples: Candidate voted for: Obama, McCain, Nader, other Brand of cola preferred: Coca Cola, Pepsi, other Method of birth control used: pill, condom, etc. Diagnosis: Autism, Aspe
33、rgers syndrome, pervasivedevelopmental disorder, othertrianglerightsld Responses often correspond to discrete choices among alternatives(categories)GLLAMM p.22Random utility modelstrianglerightsld Unobserved utility Uaiassociated with each alternative a=1, ., A forunit i=1, ., Ntrianglerightsld Rand
34、om utility models composed asUai= Vai+ epsilon1ai Vaiis linear predictor epsilon1aiis residual term (independent over i and a)trianglerightsld Alternative f chosen ifUfiUgifor all g negationslash= fepsilon1aiindependent Gumbel distributedarrowdblbothvPr(fi)=exp(Vfi)summationtextAa=1exp(Vai)multinomi
35、al logitGLLAMM p.23Covariate effects on utilitiesLinear predictor for unit i and alternative a:Vai= ma+ xprimeiga+ xaprimeibtrianglerightsld Covariates and parameters: maalternative-specific constants xivaries over subjects (but not alternatives) and has fixedalternative-specific effects gaExample:
36、Age of subject xaivaries over alternatives (and possibly subjects) and has fixedeffects b, constant across alternativesExample: Cost of treatment alternativeGLLAMM p.24Identificationtrianglerightsld Probability of choosing alternative 1 among alternatives 1, 2 and 3can be expressed in terms of utili
37、ty differencesPr(U1U2 0 and U1U3 0)trianglerightsld Therefore location of Vaiarbitrary:exp(V1i)summationtextaexp(Vai)=exp(V1i+ ci)summationtextaexp(Vai+ ci)trianglerightsld Solution: Last alternative S serves as reference alternative,set mS=0 and gS=0GLLAMM p.25Multilevel random utility modelstriang
38、lerightsld Consider three-level data with patients i (level 1) treated by doctors j(level 2) working in hospitals k (level 3)trianglerightsld Can include random effects in linear predictor:Vaijk=ma+ a(2)0jk+ a(3)0k+xprimeijkga+ a(2)jk+ a(3)k+xaprimeijkb + (2)jk+ (3)ktrianglerightsld Random intercept
39、s: a(2)0jkand a(3)0ktrianglerightsld Random Coefficients I: a(2)jkand a(3)kare alternative-specificrandom coefficients for subject-specific covariates xijktrianglerightsld Random Coefficients II: (2)jkand (3)kare random coefficients foralternative-specific covariates xaijkGLLAMM p.26Abuse of antibio
40、tics in Chinatrianglerightsld Acute respiratory tract infection (ARI) can lead to pneumonia anddeath if not properly treatedtrianglerightsld Inappropriate frequent use of antibiotics common in China in 1990s,leading to drug resistancetrianglerightsld In the 1990s WHO introduced program of case manag
41、ement forchildren under 5 with ARI in Chinatrianglerightsld Data collected on 855 children i (level 1) treated by 134 doctors j(level 2) in 36 hospitals k (level 3) in two counties (one of which wasin WHO program)trianglerightsld Response variable: Abuse defined as prescription of antibioticswhen th
42、ere were no clinical indications based on medical files1. Abuse of several antibiotics2. Abuse of one antibiotic3. Correct use of antibiotics (reference category)GLLAMM p.27Covariatestrianglerightsld 7 covariates xijk Patient level idiamondmath Age Age in years (0-4)diamondmath Temp Body temperature
43、, centered at 36Cdiamondmath Paymed Pay for medication (yes=1, no=0)diamondmath Selfmed Self medication (yes=1, no=0)diamondmath Wrdiag Failure to diagnose ARI early (yes=1, no=0) Doctor level jdiamondmath DRed Doctors education(6 categories from self-taught to medical school) Hospital level kdiamon
44、dmath WHO Hospital in WHO program (yes=1, no=0)Reference: Min Yang (2001). Multinomial Regression. In Goldstein and Leyland (Eds).Multilevel Modelling of Health Statistics, pages 107-123.GLLAMM p.28Modeltrianglerightsld No alternative-specific covariatesVaijk=ma+ a(2)0jk+ a(3)0k+xprimeijkgatriangler
45、ightsld Data:doc child alt choice11 1 011 2 111 3 012 1 012 2 012 3 1trianglerightsld Estimation in gllammgen categ1 = alt = 1gen categ2 = alt = 2eq c1: categ1eq c2: categ2gllamm alt age temp . , i(doc hosp) nrf(2 2) eqs(c1 c2 c1 c2) /link(mlogit) expanded(child choice m) basecat(3) adaptGLLAMM p.29
46、Maximum likelihood estimatesAbuse several Abuse oneParameter Est (SE) Est (SE)ga0Cons -5.72 (0.99) -0.23 (0.55)ga1Age 0.07 (0.09) 0.17 (0.08)ga2Temp -0.27 (0.13) -0.96 (0.12)ga3Paymed 0.92 (0.40) 0.12 (0.32)ga4Selfmed -0.86 (0.29) -0.49 (0.24)ga5Wrdiag 1.85 (0.26) 2.08 (0.23)ga6DRed -0.62 (0.17) 0.0
47、8 (0.11)ga7WHO -2.40 (0.62) -0.88 (0.33)Doctor-level variancesVar(a(2)0jk) 0.46 (0.28) 0.43 (0.22)Cov(1(2)0jk,2(2)0jk) -0.44 (0.13)Hospital-level variancesVar(a(3)0k) 0.88 (0.45) 0.11 (0.12)Cov(1(3)0k,2(3)0k) 0.31 (0.20)GLLAMM p.30Referencestrianglerightsld “Multilevel logistic regression for polyto
48、mous data and rankings”. Psychometrika 68,267-287, 2003. (A. Skrondal & S. Rabe-Hesketh)trianglerightsld “Generalized linear mixed effects models”. In G. Fitzmaurice et al. (Eds.)Longitudinal Data Analysis. Chapman & Hall/CRC, pp. 79-106, 2008.(S. Rabe-Hesketh & A. Skrondal)trianglerightsld “Maximum likelihood estimation of limited and discrete dependent variable models withnested random effects”. Journal of Econometrics 128, 301-323, 2005.(S. Rabe-Hesketh, A. Skrondal & A.
