demo/mph_inequality.r In hmmm: Hierarchical Multinomial Marginal Models

library(hmmm)

# see Section 4.2, pg. 350,
# "Multinomial-Poisson homogeneous models for contingency tables",
# Lang, J.B.
# The Annals of Statistics, (2004)
#
# Table 2 - 1999 statistics journals citation pattern counts (n_3=225)
# Citing x Cited statistics journals: JASA, BMCS, ANNS

y <- matrix(c(104,24,65,76,146,30,50,9,166),9,1)

# population matrix: 3 strata with 3 observations each
Zmat <- kronecker(diag(3),matrix(1,3,1))

# the 3rd stratum sample size is fixed
ZFmat <- kronecker(diag(3),matrix(1,3,1))[,3]

########################################################################
# Let (i,j) be a cross-citation, where i is the citing journal and j is
# the cited journal. Let m_ij be the expected counts of cross-citations.
# The Gini concentrations of citations for each of the journals are:
# G_i = sum_j=1_3 (m_ij/m_i+)^2  for i=1,2,3.

Gini<-function(m) {
A<-matrix(m,3,3,byrow=TRUE)
GNum<-rowSums(A^2)
GDen<-rowSums(A)^2
G<-GNum/GDen
c(G,G)-c(G,G)
}

####################################################################

# Example of MPH model subject to equality constraints:

##
## --> h = c(G1,G3)-c(G2,G1) = 0
## HYPOTHESIS: G1 = G2 = G3
##

mod_eq <- mphineq.fit(y,Z=Zmat,ZF=ZFmat,h.fct=Gini)

print(mod_eq)

# Example of MPH model subject to inequality constraints:

##
## --> d = c(G1,G3)-c(G2,G1) >= 0
## HYPOTHESIS: G1 > G2, G3 > G1
##

mod_ineq <- mphineq.fit(y,Z=Zmat,ZF=ZFmat,d.fct=Gini)

# Reference model (sat_mod):
# --> model (MPH model subject to inequality constraints) without inequalities
# ==============
# NB sat_mod --> in this case saturated model is properly "the" saturated model
# (Gsq=0), model with no constraints

mod_sat <-mphineq.fit(y,Z=Zmat,ZF=ZFmat)

# HYPOTHESES TESTED:
# NB: testA --> H0=(mod_eq) vs H1=(mod_ineq model)
#     testB --> H0=(mod_ineq model) vs H1=(sat_mod model)

hmmm.chibar(nullfit=mod_eq,disfit=mod_ineq,satfit=mod_sat)

Try the hmmm package in your browser

Any scripts or data that you put into this service are public.

hmmm documentation built on May 2, 2019, 12:27 p.m.