MPH and HLP Descriptions

Multinomial-Poisson Homogeneous (MPH) and Homogeneous Linear Predictor (HLP) Models: A Brief Description

Author Information: Joseph B. Lang, Department of Statistics and Actuarial Science, University of Iowa, Iowa City, IA 52242 USA < joseph-lang@uiowa.edu, http://www.stat.uiowa.edu/~jblang >

Citation Information: Lang, J.B. "Multinomial-Poisson Homogeneous (MPH) and Homogeneous Linear Predictor (HLP) Models: A Brief Description," Online HTML Document, 03/09/2011 11:47 AM -0600, <mph.hlp.description.htm>.

Primary References:

Lang, J.B (2004). "Multinomial-Poisson Homogeneous Models for Contingency Tables," Annals of Statistics, 32, 340-383 .

Lang, J.B. (2005). "Homogeneous Linear Predictor Models for Contingency Tables," JASA, 100, 121-134

DESCRIPTION of MPH and HLP MODELS (Special-Case):

MPH and HLP models are broad classes of models for contingency tables.

In the following, let y be the vector of c table counts, p be the unknown vector of c table probabilities, s be a vector of c strata identifiers, and F be the set of strata with a priori fixed sample sizes.

Although the R program, mph.fit, can fit, via maximum likelihood, more general models (see below), two important special cases include:
#
#     MPH (Special-Case):    y <- Y ~ MP(nu, p | strata=s, fixed=F),   h(p) = 0.
#
#     HLP (Special-Case):     y <- Y ~ MP(nu, p | strata=s, fixed=F),   L(p) = X beta.
#
Here, h(.) is a smooth constraint function and L(.) is a smooth link function. It is assumed that the constraints in h(p) = 0 are non-redundant so that the Jacobian, d h(p)'/d p,   is full column rank.

The link L(.) is allowed to be many-to-one and row-rank deficient, so this HLP model is quite general. It is only required that the implied constraints, U'L(p) = 0, where null(U') = span(X), are non-redundant. Here, MP stands for the Multinomial-Poisson distribution. The parameters are nu, the vector of expected sample sizes, and p, the vector of table probabilities.

The notation y <- Y ~ MP(nu, p| strata=s, fixed=F) means that y is a realization of random vector Y, which is composed of independent blocks of multinomial and/or Poisson random variables. The strata vector s determines the blocks and F determines which blocks are multinomial and which blocks comprise independent Poisson random variables. More specifically, suppose there are K strata, so s contains K distinct strata identifiers. The components in Y corresponding to s=identifier[k] make up a block. If identifier[k] is included in F, then this block has a multinomial distribution and nu[k] is the a priori known, i.e. fixed, sample size. If identifier[k] is not in F, then this block comprises independent Poisson random variables and nu[k] is an unknown expected sample size.

Note: Given the observed counts y, the pair (strata,fixed)=(s, F) contains the same information as the sampling plan triple (Z,ZF,nF) described in Lang (2004,2005). Specifically,

Z=Z(s), the strata/population matrix, is determined by s.
ZF = ZF(s,F), the sampling constraint matrix, is determined by s and F.
nF = ZF'y is the vector of a priori fixed sample sizes.

Special case MP distributions include...

Full Multinomial:    MP(nu,p| strata=1, fixed="all")
                               A simple random sample of fixed size nu is taken from a single strata (population).

Product Multinomial: MP(nu,p| strata=s, fixed="all")
                              A stratified random sample of fixed sample sizes nu=(nu[1],...,nu[K]) is taken from the K strata determined by s.

Full Poisson:        MP(nu,p| strata=1, fixed="none")
                            A simple random sample is taken from a single strata (population). The sample size is random and follows a Poisson distribution with unknown mean nu.

Product Poisson:     MP(nu,p| strata=s, fixed="none")
                            A stratified random sample is taken from K strata. The sample sizes are all random and distributed as Poissons with unknown means in nu=(nu[1],...,nu[K]).

Specifying the MP distribution in 'mph.fit'...

In older versions of mph.fit (2.1 and lower), the user had to enter the sampling plan triple (Z,ZF,n). In 'mph.fit', version 3.0 and higher, the user need only enter ('strata', 'fixed.strata'), the input variables corresponding to (s, F). Keywords, fixed.strata="all" ["none"] means that all [none] of the strata have a priori fixed sample sizes.

To fit MPH (Special Case), the user must enter the counts 'y', the constraint function 'h.fct' (alias 'constraint'), and the sampling plan variables, 'strata' and 'fixed.strata'. Note: The user can omit the sampling plan variables if the default, multinomial sampling (strata=1,fixed="all"), can be assumed.

To fit HLP (Special Case), the user must enter the counts 'y', the link function 'L.fct' (alias 'link'), the model matrix 'X', and the sampling plan variables, 'strata' and 'fixed.strata'. Note: The user can omit the sampling plan variables if the default, multinomial sampling, can be assumed.

IMPORTANT: When specifying the model and creating the input objects for 'mph.fit', keep in mind that the interpretation of the table probabilities p depends on the sampling plan!
Specifically, if the i^{th} count y[i] is in block k (i.e. corresponds with strata identifier[k]), then the i^{th} table probability p[i] is the conditional probability defined as

p[i] = probability of a Type i outcome GIVEN that the outcome is one of the types in stratum k.

For example, in an IxJ table with row variable A and column variable B, if row-stratified sampling is used, the table probabilities have the interpretation, p[i,j] = prob of a Type (i,j) outcome GIVEN that the outcome is one of the types in stratum i (i.e. one of (i,1),...,(i,J)). That is, p[i,j] = P(A=i,B=j|A=i) = P(B=j|A=i). For column-stratified sampling, p[i,j] = P(A=i|B=j). And for non-stratified sampling, p[i,j] = P(A=i,B=j).

Log-Linear Models: Log-Linear models specified as log(p) = X beta, are special-case HLP models.

As with any HLP model, log(p) = X beta can be restated as a collection of constraints; specifically, log(p) = X beta is equivalent to h(p) = U'log(p) = 0,where null(U') = span(X). Noting that Z'p = 1, we see that to avoid redundant constraints, span(X) should contain span(Z). Loosely, fixed-by-sampling-design parameters should be included.

Log-Linear models of the form log(p) = X beta are simple to fit using 'mph.fit'. For example, mph.fit(y, link="logp",X=model.matrix(~A+B)) or, equivalently, mph.fit(y, link=function(p){log(p)},X=model.matrix(~A+B))

DESCRIPTION (continued): MORE GENERAL MPH and HLP MODELS...

Instead of (nu, p), the MP distribution can alternatively be parameterized in terms of the vector of expected table counts, m = E(Y). Formally, (nu, p) and m are in one-to-one correspondence and satisfy:

m = D(Z nu)p and nu= Z'm, p = D^{-1}(ZZ'm)m.

HHere, Z = Z(s) is the cxK strata/population matrix determined by strata vector s. Specifically, Z[i,k] = 1(s[i] = identifier[k]).

The MPH (Special-Case) Model given above is a special case because it constrains the expected counts m only through the table probabilities p. Similarly, the HLP (Special-Case) Model given above is a special case because it uses a link function that depends on m only through the table probabilities p.

More generally, 'mph.fit' computes maximum likelihood estimates and fit statistics for MPH and HLP models of the form...
#
# MPH. y <- Y ~ MP(nu, p | strata=s, fixed=F), h(m) = 0.
# HLP: y <- Y ~ MP(nu, p | strata=s, fixed=F), L(m) = X beta.
#
Here, h(.) is a smooth constraint function that must also be Z (i.e. strata s) homogeneous. L(.) is a smooth link function that must also satisfy the HLP conditions with respect to Z (i.e. strata s) and X. That is,

(1) L(.) has HLP link status wrt Z AND
(2) The implied constraint function h(m)=U'L(m) is Z homogeneous. Here, null(U') = span(X).

L(.) has HLP link status wrt Z if, for m = D(Z nu)p...

(1)(a) L(m) = a(nu) + L(p), where a(nu1/nu2)-a(1) = a(nu1) - a(nu2),   i.e. a(nu) has the form C log nu + constant.
                       OR

(1)(b) L(m) = G(nu)L(p), where G(nu) is a diagonal matrix with diagonal elements that are powers of the nu elements, i.e. L(.) is Z homogeneous (see Lang 2004).
                     OR

(1)(c) The components of L(.) are a mixture of types (a) and (b):    L[j](m) = a[j](nu) + L[j](p) or L[j](m) = G[j](nu)L[j](p),          j=1,...,l.

Lang (2005) described HLP models that satisfied (1)(a) and (2), but the definition of HLP models can be broadened to include those models satisfying (1) and (2). That is, HLP models can be defined so they also include models that satisfy (1)(b) and (2) or (1)(c) and (2).     Version 3.1+ of mph.fit uses this broader definition of HLP Model.

Note:    The input variable 'h.mean' must be set to TRUE to fit this more general MPH model. Similarly, the input variable 'L.mean' must be set to TRUE to fit this more general HLP model. (An exception: If the link function is specified using the keyword "logm" then 'L.mean' is automatically set to TRUE.)
Note:   'mph.fit', version 3.0 and above, carries out "necessary-condition" checks of Z homogeneity of h(.) and HLP link status of L(.) for these general models.

Log-Linear Models: Log-Linear models of the form log(m) = X beta are HLP models provided the span(X) contains the span(Z). Loosely, provided fixed-by-design parameters are included, the log-linear model is a special case HLP model.

Log-Linear models of the form log(m) = X beta are simple to fit using 'mph.fit'. For example, mph.fit(y, link="ogm",X=model.matrix(~A+B)) or, equivalently, mph.fit(y, link=function(m){log(m)}, L.mean=T, X=model.matrix(~A+B))

Note: Most reasonable generalized loglinear models, which have the form L(m) = C log Mm = X beta, are also HLP models (see Lang 2005).

QUESTIONS: Is there documentation for mph.fit? Yes, click on mph.fit documentation. Are the examples showing how to use mph.fit? Yes, click on MPH and HLP model fitting examples.

Page Last Updated: 03/09/2011 11:47 AM -0600, Joseph B. Lang