Home> mph.fit documentation

MAXIMUM LIKELIHOOD FITTING of MULTINOMIAL-POISSON
HOMOGENEOUS (MPH) MODELS
for CONTINGENCY TABLES using *MPH.FIT*

Author: Joseph B. Lang, Department of Statistics and Actuarial Science, University of Iowa, Iowa City, IA 52242 USA

**Document Description: ** Herein, you will
find supporting documentation for the R program
**
mph.fit**.
The program
**
mph.fit**
computes maximum likelihood estimates and model assessment statistics for the
broad class of multinomial-Poisson homogeneous (MPH) contingency table models.
This document includes, or provides links to, the following information:

primary references | software availability | implementation instructions | program description |

**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

**SOFTWARE AVAILABILITY:**

To receive a free copy of the R programs referred to in this document, contact me (Joseph B. Lang) at joseph-lang@uiowa.edu. [

Note: The statistical package R is freeware that is available athttp://cran.r-project.org;its syntax is very similar to that of the commercial softwareS-plus.]

**IMPLEMENTATION of
**
**
mph.fit and accompanying R programs
**

In response to your email request, you should receive the file, "mph.Rcode.txt", which can be directly sourced into R, e.g. R> source("mph.Rcode.txt").

Sourcing "mph.Rcode.txt" results in the creation of 9 functions:

mph.fit (version 3.1)

mph.summary

num.deriv.fct

create.U

create.Z.ZF

check.homog

check.HLP

block.fct

M.fct

The main maximum-likelihood fitting program, mph.fit, uses num.deriv.fct, create.U, create.Z.ZF, check.homog, and check.HLP.

The program mph.summary is used to create a useful summary of the object created by mph.fit.

The remaining functions, block.fct and M.fct, are useful for MPH and HLP modeling, but are not essential for using mph.fit and mph.summary.

**PROGRAM DESCRIPTION: **

The following lines are the comments found at the top of the file mph.Rcode.txt . These comments currently serve as user documentation. See also the basic numerical examples, numerical examples , and the primary references above.

# # Sourcing this R file (> source("mph.Rcode.txt") ) results in the creation # of the following 9 functions: # # mph.fit, mph.summary, num.deriv.fct, create.U, create.Z.ZF, check.homog, check.HLP # block.fct, M.fct # # ###################### Begin mph.fit ########################################## mph.fit <- function(y, h.fct=constraint,constraint=NULL, h.mean=F, L.fct=link, link=NULL, L.mean=F, X=NULL, strata=rep(1,length(y)),fixed.strata="all", check.homog.tol=1e-9, check.HLP.tol=1e-9, print.iter=F, maxiter=100, step=1, change.step.after=0.25*maxiter, y.eps=0, iter.orig=5, m.initial=y, norm.diff.conv=1e-6,norm.score.conv=1e-6, max.score.diff.iter=10, derht.fct=NULL,derLt.fct=NULL) { # mph.fit, version 3.1 # # Originally Created: August 25, 2000 (last revised: 2/12/07, 4/23/09, 5/20/09) # # Author: Joseph B. Lang, # Dept. of Statistics and Actuarial Science # Univ. of Iowa, Iowa City, IA 52242 # joseph-lang@uiowa.edu # # 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: # # This program computes maximum likelihood estimates and fit statistics for # multinomial-Poisson homogeneous (MPH) and homogeneous linear # predictor (HLP) models for contingency tables (see Lang 2004, 2005). # # 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 mph.fit can fit 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)). # = 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 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. # # Here, 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 # # (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 # # (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="logm",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). # # SUBROUTINES and COMPLEMENTARY FUNCTIONS: # # Subroutines: 'mph.fit' uses 5 other short subroutines: # 'num.deriv.fct', 'create.U', 'create.Z.ZF', # 'check.homog', and 'check.HLP'. # # Complementary Function: # 'mph.summary' summarizes and displays the results from 'mph.fit'. # # # INPUT: # # Required Input: # y = vector (not matrix) of table counts # # Optional Input: # h.fct = function object that defines the constraint function h(.). # It must return a column vector. # h.fct can also be set to 0, in which case h(.) is viewed as the # 0 function, so no constraints are imposed. # # By default, h(.) is viewed as a function of the # table probabilities p. If h.mean is set to h.mean=T, # then h(.) is viewed as a function of the expected counts m. # Default: h.fct = NULL (If both h.fct=NULL and L.fct = NULL, # then h.fct is set to 0 and no constraints are imposed. # If h.fct = NULL and L.fct is not= NULL then # h.fct will be computed as U'L.fct.) # # constraint is an alias for the argument 'h.fct'. # Argument 'constraint' is secondary to the primary argument 'h.fct' # in the following senses: # If 'constraint' and 'h.fct' are not equal, 'h.fct' is used. # # h.mean = logical argument. If h.mean=F (the default), h.fct is treated as # a function of p. If h.mean=T then h.fct is treated as a function of m. # # L.fct = function object that defines the link L(.) in the HLP model; # it must return a column vector. # # Or... L.fct = keyword, where # candidate keywords include "logp" and "logm". # # Entering L.fct="logp" tells the program to create the # function object as L.fct <- function(p){log(p)}. # L.fct="logm" tells the program to (i) create the # function object as L.fct <- function(m){log(m)} and # (ii) set L.mean = T. # # By default, L.fct is treated as a function of the # table probabilities p (even if the argument in the # L.fct function object is "m"). If L.mean is set to L.mean=T, # then L.fct is treated as a function of the expected counts m. # Default: L.fct = NULL means no constraints of the form L(p)=X beta # are imposed. # # link is an alias for the argument 'L.fct'. # Argument 'link' is secondary to the primary argument 'L.fct' # in the following senses: # If 'link' and 'L.fct' are not equal, 'L.fct' is used. # # L.mean = logical argument. If L.mean=F (the default), L.fct is treated # as a function of p. If L.mean=T, L.fct is treated # as a function of m. # # # X = HLP model design matrix, assumed to be full rank. # (default: X = NULL; i.e., it is left unspecified and unused.) # # strata = vector of the same length as y that gives the stratum # membership identifier. Default: strata = rep(1,length(y)); # i.e. the default is the single stratum (non-stratified) setting. # Examples, strata=A or strata=c(1,1,1,2,2,2,3,3,3) or # strata=paste(sep="","A=",A,",B=",B) # # fixed.strata = the object that gives information on which stratum have # fixed sample sizes. It can equal one of the keywords, # fixed.strata = "all" or fixed.strata = "none", or it can # be a vector of stratum membership identifiers, e.g. # fixed.strata=c(1,3) or fixed.strata=c("pop1","pop5"). # (default: fixed.strata="all") # # check.homog.tol = round-off tolerance for Z homogeneity check # (default: check.homog.tol=1e-9) # check.HLP.tol= round-off tolerance for HLP Link status check # (default: check.HLP.tol=1e-9) # # print.iter = F, Do not print out iteration information. If print.iter=T # then iteration information is printed out. # # maxiter = maximum number of iterations (default: maxiter=100) # step = step-size value (default: step=1) # change.step.after: If the score value increases for more than # change.step.after iterations in a row, then the # initial step size is halved. (default: # change.step.after = 0.25*maxiter) # y.eps = non-negative constant to be temporarily added # to the original counts in y (default: y.eps=0) # iter.orig = iteration at which the original counts will # be used (default: iter.orig=5). # m.initial = initial estimate of m (default: m.initial=y) # See Input Note 7 below. # norm.diff.conv = convergence criteria value; see norm.diff # in the Output section # (default: norm.diff.conv=1e-6) # norm.score.conv = convergence criteria value; see norm.score in # the Output section # (default: norm.score.conv = 1e-6) # max.score.diff.iter, The variable score.diff.iter # keeps track of how long norm.score is # smaller than norm.score.conv, but norm.diff is # greater than norm.diff.conv. If this is the # case too long (i.e. score.diff.iter >= max.score.diff.iter), # then stop the iterations because the solution likely # includes at least one ML fitted value of 0. # (default: max.score.diff.iter = 10) # derht.fct = NULL, function object that computes analytic derivative of # the transpose of the constraint function h(.) with # respect to m. # If h(.) maps from Rp to Rq (i.e. there are q # constraints), then derht.fct returns a # pxq matrix of partial derivatives. # (default: derht.fct=NULL. This means that the # derivative is calculated numerically.) # derLt.fct = NULL, function object that computes analytic derivative of # the transpose of the link function L(.) with # respect to m. # If L(.) maps from Rp to Rq (i.e. there are q # link components), then derLt.fct returns a # pxq matrix of partial derivatives. # (default: derLt.fct =NULL, i.e. by default this # derivative is calculated numerically.) # # Input Notes: # # 1. CONSTRAINT FUNCTION. # # 'constraint' is an alias for 'h.fct'. # # If h.fct is a function object, it must return a column vector. # # By default, h.fct is treated as a function of the table # probabilities p. To treat h.fct as a function of the expected # counts m, you must set h.mean=T (by default, h.mean=F). # # The fitting algorithm will fail if the constraints in h.fct=0 # are redundant. (The redundancy results in a rank-deficient # Jacobian, d h'/d m.) # # 2. MODEL WITH NO CONSTRAINTS. # # The model with no constraints can be specifed using # h.fct = 0. The no-constraint model is the default # when neither h.fct nor L.fct are input (i.e. when # h.fct=NULL and L.fct=NULL). # # 3. HLP MODEL SPECIFICATION. # # 'link' is an alias for 'L.fct'. # # For HLP models, both L.fct and X must be specified. # # The design matrix X must be of full column rank. # # 'mph.fit' recognizes two keywords for link specification, # 'logp' and 'logm'. These are convenient for log-linear modeling. # # If L.fct is a function object, it must return a column vector. # # By default, L.fct is treated as a function of the table # probabilities p. To treat L.fct as a function of the expected # counts m, you must set L.mean=T (by default, L.mean=F). # # The constraint function h.fct is typically left # unspecified for HLP models, but it need not be. # # If h.fct is left unspecified, it is created within the program as # h.fct(m) <- function(m){ t(U)%*%L.fct(m)}, where matrix # U is an orthogonal complement of X. # # If h.fct is specified, the constraints implied by # L.fct(p) =X beta, or L.fct(m)=X beta, are ignored. # # Note: Although the HLP constraints are ignored when # h.fct is specified, estimates of beta and the link are # computed under the model with constraints # h.fct(p)=0 or h.fct(m)=0. # # The fitting algorithm will fail to converge if the implied # constraints, U'L.fct = 0, include redundancies. # # 4. EXTENDED ML ESTIMATES. # When ML estimates are non-existent owing to # zero counts, norm.diff will not converge to # zero, instead it tends to level off at some constant # positive value. This is because at least one ML # fitted value is 0, which on the log scale is # log(0)=-Inf, and the log-scale iterates slowly # move toward -Inf. One solution to this problem is # to set the convergence value norm.diff.conv # to some large number so only the score convergence # criterion is used. In this case, the algorithm # often converges to a solution that can be viewed # as an extended ML estimate, for which 0 estimates # are allowed. Versions 2.0 and higher automates the # detection of such problems. See the description of the # input variable 'max.score.diff.iter' above and the # MISC COMPUTATIONAL NOTES below. # # 5. CONVERGENCE PROBLEMS / FINE TUNING ITERATIONS. # # First check to make sure that the model is correctly # specified and redundant constraints are avoided. # # When ML estimates exist, but there are # non-convergence problems (perhaps caused by zero # counts), a modification to the tuning # parameters `step', 'change.step.after', `y.eps', # and/or `iter.orig' will often lead to convergence. # # With zero counts, it might help to set # y.eps = 0.1 (or some other positive number) # and iter.orig=5 (the default). This tells # the program to initially use y+y.eps rather than # the original counts y. At iteration 5=iter.orig, # after the algorithm has had time to move toward a # non-boundary solution, the original counts are again # used. # # To further mitigate non-convergence problems, # the parameter `step' can be set to a smaller value # (default: step=1) so the iterates do not change as much. # # 6. The initial estimate of m is actually # m.initial + y.eps + 0.01*((m.initial+y.eps)==0) # The program defaults are m.initial=y and y.eps=0. # Note: If m.initial > 0 and y.eps=0, then the # initial estimate of m is simply m.initial. # # # OUTPUT # y = vector of counts used in the algorithm # for ML estimation. Usually, this vector # is identical to the observed table counts # m = vector of ML fitted values # covm = approximate covariance matrix of fitted values # p = vector of cell probability ML estimates # covp = approximate covariance matrix of cell probability estimators # lambda = vector of Lagrange multiplier ML estimates # covlambda = approximate covariance matrix of multiplier estimators # resid = vector of raw residuals (i.e. observed minus fitted counts) # presid = vector of Pearson residuals # adjresid = vector of adjusted residuals # covresid = approximate covariance matrix of raw residuals # Gsq = likelihood ratio stat for testing Ho: h(m)=0 vs. Ha:not Ho # Xsq = Pearson's score stat (same as Lagrange multiplier stat) # for testing Ho: h(m)=0 vs. Ha: not H0. # Wsq = Generalized Wald stat for testing Ho: h(m)=0 vs. # Ha: not H0. # df = degrees of freedom associated with Gsq and Xsq (df= dim(h)) # beta = vector of HLP model parameter ML estimates # covbeta = approximate covariance matrix of HLP # model parameter estimators # L = vector of HLP model link ML estimates # Lobs = vector of HLP model observed link values, L(y). # covL = approximate covariance matrix of HLP model # link estimators # Lresid = vector of adjusted link residuals of the form # (L(obs) – L(fitted))/ase(L(obs)-L(fitted)). # iter = number of update iterations performed # norm.diff = distance between last and second last `theta' iterates # (`theta' is the vector of log fitted values and # Lagrange multipliers) # norm.score= distance between the score vector and zero # h.fct = constraint function used in algorithm # h.input.fct = constraint function as originally input # h.mean = logical variable (h.mean=F, h.fct treated as function of p, # h.mean=T, h.fct treated as function of m) # derht.fct = analytic function used in algorithm that computes # derivative of transpose of constraint function h # L.fct = link function used in algorithm # L.input.fct = link function as originally input # L.mean = logical variable (L.mean=F, L.fct treated as function of p, # L.mean=T, L.fct treated as function of m) # derLt.fct = analytic function used in algorithm that computes # derivative of transpose of link function L # X = HLP model design matrix used in algorithm # U = orthogonal complement of design matrix X # triple = a list object containing the sampling plan triple (Z,ZF,n), # where Z is the population (or strata) matrix, # ZF is the sampling constraint matrix, # and n is the collection of fixed sample sizes. # strata = strata variable used as input # fixed.strata = strata corresponding to fixed sample sizes # version = `mph.fit' version number # warn.message = message stating whether or not the original # counts were used. # warn.message.score = message stating whether or not the norm score # convergence criterion was met. # warn.message.diff = message stating whether or not the norm diff # convergence criterion was met. # # Output Notes # 1. ITERATION HISTORY. # An iteration history is printed out when 'print.iter' is set # equal to TRUE. A single line of the history looks like the following: # "iter= 18[0] norm.diff= 3.574936e-08 norm.score= 1.901705e-15" # Here, iter is the number of update iterations performed. # The number in [] gives the number of step size searches # required within each iteration. # norm.diff and norm.score are defined above. # Finally, the time elapsed is output. # Note: For the model with no restrictions (h.fct = 0) # there are no step size changes. # # 2. STORING and VIEWING RESULTS. # To store the results of mph.fit, issue a command # like the following example # > results <- mph.fit(y, h.fct=h.fct) # Use program "mph.summary" to view the results # of mph.fit. Specifically, if the results of mph.fit # are saved in object "results", submit the command # mph.summary(results) [or mph.summary(results,T) or # mph.summary(results,T,T) depending on how much of # the output you need to see.] # # 3. The output objects beta, covbeta, L, covL, and # Lresid will be set to "NA" unless an HLP model is # specified (i.e. L.fct and X are input) # # MISC COMPUTATIONAL NOTES: # # For computational efficiency (in particular, to speed up # multiplications involving diagonal matrices), the expression # (A*c(x)) is used in place of diag(c(x))%*%A # # Version 1.2 Major Updates: # 1. "<-" replaces "_" to reflect changes in R, version 1.8.1 # 2. Input variable "chscore.criterion" added. This variable # is used to restrict the updating move so that the # norm of the score changes by no more than a specified # amount. This restriction is used to fine tune the # iterative algorithm so that convergence is more likely. # Version 1.3 Update: # 1. Allows user to input initial estimates of m. See input # variate m.initial. # 2. 1.3.1. Fixed bug regarding warning message about original # counts NOT being used... # if ((iter < iter.orig) & (y.eps != 0)) changed to # if ((iter <= iter.orig) & (y.eps != 0)) # Version 2.0 Update: # 1. Z = matrix(1,length(y),1) is now the default; i.e. the # "single strata setting" is the default. # 2. Iterations are fine-tuned a bit so the algorithm is # more robust, especially in sparse count settings. # The 'chscore.criterion' has been replaced. # 3. h.fct, L.fct, and X now have NULL as the default value. # is.null checks are now made in the code. # 4. Added max.score.diff.iter criterion to better handle ML # fitted values of 0. score.diff.iter # keeps track of how long norm.score is # smaller than norm.score.conv, but norm.diff is # greater than norm.diff.conv. If this is the # case too long (i.e. score.diff.iter >= max.score.diff.iter), # then stop the iterations because the solution likely # includes at least one ML fitted value of 0. # 5. Added warn.message.score and warn.message.diff to the # output. These messages state whether or not the # norm score and norm diff convergence criteria were met. # 6. Added change.step.after option. Default: # change.step.after=0.25*maxiter. # If score value increases for more than change.step.after # iterations in a row, the initial step size is halved. # # Version 3.0: Substantial changes were made in version 3.0. # # 1. The sampling plan matrices Z and ZF are no longer entered # directly. With version 3.0, the user specifies # a strata variable that gives the stratum membership identifiers # associated with each count in y. The user can specify the stratum # membership identifiers that correspond to fixed sample # sizes using the 'fixed.strata' input variable. # # 2. Specification of model constraints is simplified. # By default constraint function h(.) is viewed as a function # of table probabilities p, rather than expected counts m. # Similarly, by default, the link function L(.) is viewed as a function # of table probabilities rather than expected counts. # In this case, Z homogeneity is guaranteed. That is, h(p)= 0 # and L(p) = X beta are both MPH models. # # 3. The user can override the default of the previous item. # That is, constraints [links] can be specified in terms of # expected counts, rather than table probabilities. The override # is accomplished by changing the default h.mean=F [L.mean=F] to h.mean=T # [L.mean=T]. # # When h.mean=T...Whereas, h(p) is Z homogeneous, there is no # guarantee that h(m) is Z homogeneous. Version 3.0 now runs # a necessary-condition check to make sure the constraint function # h(m) is Z homogeneous. If it is not, a message to that effect # is returned and the model is not fitted. # # When L.mean=T...Whereas, L(p) satisfies the conditions for # being an HLP link, there is no guarantee that L(m) does. # Version 3.0 now runs a necessary-condition check to make # sure L(m) is an HLP link. If it is not, a message to that effect # is returned and the model is not fitted. # # 4. The output from mph.fit has been expanded to include # more information about the sampling plan. When available, # dimension labels are used. Correspondingly, mph.summary # has been modified. # # 5. Aliases for h.fct and L.fct are available. The user can # use 'constraint' instead of 'h.fct' and 'link' instead of 'L.fct'. # # 6. Keywords for log-linear modeling are available: Rather # than creating a function object, the user can enter L.fct="logp" # or L.fct="logm" # # 7. By default (printer.iter=F), iteration information is not printed. # # 8. The iteration process is initiated using the counts y+y.eps, rather # than y. In version 3.0, when y.eps != 0, convergence will not # occur until after the counts are reset to their original values. # This avoids the problem of convergence based on the adjusted # counts y+y.eps. (If the user wants convergence based on # adjusted counts, they will need to enter the adjusted counts # directly into the y input variable.) # On a related note, if iter.orig >= maxiter, it is reset in the # program as: iter.orig <- maxiter - 1. # # Version 3.1: # # 1. The definition of HLP model has been broadened to include # models with links L(.) that are Z homogeneous of any order. # See the relevant discussion above in the section entitled, # "DESCRIPTION (cont'd): MORE GENERAL MPH and HLP MODELS..." # The HLP link status check has been modified accordingly. # # #

Please do not distribute this documentation or the R program files without express permission of author, Joseph B. Lang ( joseph-lang@uiowa.edu)

*Page Last Updated:
07/20/2009 10:47 PM -0500,
Joseph B. Lang*