Fisher scoring for some univariate discrete distributions Thomas Yee - PowerPoint PPT Presentation

Fisher scoring for some univariate discrete distributions Thomas Yee University of Auckland 26 August 2010 t.yee@auckland.ac.nz http://www.stat.auckland.ac.nz/~yee Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 1/62 1 / 62

Outline of This Talk Outline of This Talk Introduction to VGLMs and VGAMs 1 VGLMs 2 VGAMs 3 Some Implemented Models 4 Zero-inflated Poisson model 5 Robustness of negbinomial() 6 Yet to do . . . 7 Closing Comments 8 Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 2/62 2 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs I The VGAM package for R implements several large classes of regression models of which vector generalized linear and additive models (VGLMs/VGAMs) are most commonly used. The primary key words are iteratively reweighted least squares (IRLS), maximum likelihood estimation, Fisher scoring, additive models. Other concepts are reduced-rank regression, constrained ordination, vector smoothing. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 3/62 3 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs II Basically . . . VGLMs model each parameter, transformed if necessary, as a linear combination of the explanatory variables. That is, g j ( θ j ) = η j = β T j x = β ( j )1 x 1 + · · · + β ( j ) p x p (1) where g j is a parameter link function . VGAMs extend (1) to g j ( θ j ) = η j = f ( j )1 ( x 1 ) + · · · + f ( j ) p ( x p ) (2) i.e., an additive model for each parameter. Estimated by smoothers, this is a data-driven approach. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 4/62 4 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs III Example: negative binomial Y has a probability function that can be written as � y + k − 1 � � � y � � k µ k P ( Y = y ; µ, k ) = y µ + k k + µ where y = 0 , 1 , 2 , . . . . Parameters µ > 0 and k > 0. VGAM can fit β T log µ = η 1 = 1 x , β T log k = η 2 = 2 x , with vglm(y ∼ x2 + x3 + · · · + xp, negbinomial(zero = NULL)) Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 5/62 5 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs IV The framework extends GLMs and GAMs in three main ways: (i) y not restricted to the exponential family, (ii) multivariate responses y and/or linear/additive predictors η are handled, (iii) η j need not be a function of a mean µ : η j = g j ( θ j ) for any parameter θ j . This formulation is deliberately general so that it encompasses as many distributions and models as possible. We wish to be limited only by the assumption that the regression coefficients enter through a set of linear or additive predictors. Given the covariates, the conditional distribution of the response is intended to be completely general. More general = ⇒ more useful. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 6/62 6 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs V The scope of VGAM is very broad; it potentially covers univariate and multivariate distributions, categorical data analysis, quantile and expectile regression, time series, survival analysis, mixture models, extreme value analysis, nonlinear regression, reduced-rank regression, ordination, . . . . It conveys GLM/GAM-type modelling to a much broader range of models. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 7/62 7 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs VI QRR−VGLM RR−VGLM RR−VGAM VGLM VGAM Generalized RR−VLM Normal errors VLM VAM LM Parametric Nonparametric Figure: Flowchart for different classes of models. Legend: LM = linear model Y = X β + ε , V = vector, G = generalized, A = additive, RR = reduced-rank. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 8/62 8 / 62

Introduction to VGLMs and VGAMs Introduction to VGLMs and VGAMs VII t η = ( η 1 , . . . , η M ) T Model S function Reference 2 x 2 (= B T x ) B T 1 x 1 + B T VGLM vglm() Yee & Hastie (2003) p 1 + p 2 P H k f ∗ B T 1 x 1 + k ( x k ) VGAM vgam() Yee & Wild (1996) k = p 1 +1 B T 1 x 1 + A ν RR-VGLM rrvglm() Yee & Hastie (2003) 0 1 ν T D 1 ν B . C B T 1 x 1 + A ν + . QRR-VGLM cqo() Yee (2004) @ A . ν T D M ν P R B T 1 x 1 + f r ( ν r ) RR-VGAM Yee (2006) cao() r =1 Table: VGAM & its framework. The latent variables ν = C T x 2 , x T = ( x T 1 , x T 2 ). More abbreviations: C = constrained, O = ordination, Q = quadratic. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 9/62 9 / 62

Introduction to VGLMs and VGAMs The VGAM Package I See VGAMrefcard.pdf for a summary. How does VGAM differ from other packages? its breadth, its similarity to glm() and gam() , e.g., > n = 20 > y = rchisq(n, df = exp(2)) > fit = vglm(y ~ 1, family = chisq) > fitted(fit) > summary(fit) its size, its room for future extension. Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 10/62 10 / 62

VGLMs VGLMs I Data ( x i , y i ) , i = 1 , . . . , n , from independent“individuals” . Conditional distribution of y given x is Definition f ( y | x ; β ) = h ( y , η 1 , . . . , η M , φ ) , where for j = 1 , . . . , M , η j ( x ) = β T η j = j x , (3) ( β ( j )1 , . . . , β ( j ) p ) T , = β j ( β T 1 , . . . , β T M ) T , = β φ = a vector of scale factors . Often g j ( θ j ) = η j for parameters θ j and link functions g j . Nb. −∞ < η j < ∞ . Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 11/62 11 / 62

VGLMs VGLM Examples I 1 Negative binomial distribution . For y = 0 , 1 , 2 , . . . , � y + k − 1 � � � y � � k µ k f ( y ; µ, k ) = , µ, k > 0 , y µ + k k + µ η 1 = log µ, η 2 = log k are good choices . If k > 1 then could have η 2 = log log k . Use negbinomial(lk = loglog) . Later (Slide ?? ): will be able to fit (i) Var ( Y ) = s 1 µ , (ii) Var ( Y ) = µ + µ s 2 , easily regardless whether s 1 and s 2 are known or unknown. [cf. Var ( Y ) = µ + µ 2 / k ] Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 12/62 12 / 62

VGLMs VGLM Examples II Bivariate (logistic) odds-ratio model 2 Data: ( Y 1 , Y 2 ) where Y j = 0 or 1. Examples ◮ Y 1 = 1 if left eye is blind, Y 2 = 1 if right eye is blind. ◮ Y j = 1 / 0 if Species j is present/absent. ◮ Y 1 = 1 / 0 if person has/hasn’t cancer, Y 2 = 1 / 0 if person has/hasn’t diabetes. p j = P ( Y j = 1) , marginal probabilities, p rs = P ( Y 1 = r , Y 2 = s ) , r , s = 0 , 1 , joint probabilities, ψ = p 00 p 11 / ( p 01 p 10 ) (Odds ratio). Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 13/62 13 / 62

VGLMs VGLM Examples III Model: logit p j ( x ) = η j ( x ) , j = 1 , 2 , log ψ ( x ) = η 3 ( x ) . Recover p rs ’s from p 1 , p 2 and ψ . Q: why not allow a probit or complementary log-log link? > binom2.or("cloglog", exchangeable = TRUE, zero = NULL) Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 14/62 14 / 62

VGLMs VGLM Examples IV Exchangeable data = ⇒ constrain η 1 = η 2 , (e.g., ears and eyes), i.e., cloglog p j ( x ) = η 1 ( x ) , j = 1 , 2 , log ψ ( x ) = η 3 ( x ) . Note:       � � β ∗ η 1 1 0 p p � � (1) k β ∗     =   (1) k β ∗ η 2 1 0 x k =   x k . (4) (1) k β ∗ (2) k η 3 0 1 β ∗ k =1 k =1 (2) k General formula is p � η ( x ) = B T x = H k β ∗ k x k , (5) k =1 � � T where β ∗ k = β ∗ (1) k , . . . , β ∗ . ( r k ) k Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 15/62 15 / 62

VGLMs VGLM Examples V Nonproportional odds model for a categorical response i.e., 3 Y ∈ { 1 , 2 , . . . , M + 1 } logit P ( Y ≤ j | x ) = η j ( x ) , j = 1 , . . . , M . Proportional odds model : constrain η j ( x ) = α j + η ( x ) (aka the parallelism assumption, which stops the probabilities from becoming negative or greater than 1). Have H 1 = I M ; H k = 1 M for k = 2 , . . . , p . Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 16/62 16 / 62

VGLMs VGLM Examples VI Other links (for 0 < p < 1): Φ − 1 ( p ), probit log {− log(1 − p ) } , cloglog tan( π ( p − 1 2 )). cauchit > vglm(y ~ x2, cumulative(link = probit)) Thomas Yee (Auckland University) Fisher scoring univariate discrete distributions 26 August 2010 17/62 17 / 62