ML FITTING of MPH MODELS using MPH.FIT:  NUMERICAL EXAMPLES

Author:   Joseph B. Lang,   Department of Statistics and Actuarial Science, University of Iowa, Iowa City, IA 52242 USA   (September 27, 2002, last updated: April 15, 2005, February 7, 2007)

The  numerical examples below were fitted using mph.fit, VERSION 1.0.   The current version is VERSION 2.0.    Please read about the changes in the header of the file mph.fit.  (Some of the convergence issues described in these examples have been ironed out in version 2.0.)

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

Lang, J.B. (2002, 2007). "Maximum Likelihood Fitting of Multinomial-Poisson Homogeneous (MPH) Models for Contingency Tables using  MPH.FIT."  online html document.

Lang, J.B. (2002, 2007).  "Multinomial-Poisson Homogeneous and Homogeneous Linear Predictor Models: A Brief Description,"  online html document.

Lang, J.B. (2002, 2007). "ML Fitting of MPH Models using MPH.FIT: Numerical Examples."  online html document.   [This document.]

Numerical Examples:

EXAMPLE  1. Mean Response Model (colds in children)
EXAMPLE  2. Dispersion Linear Trend Model (Danish marital status)
EXAMPLE  3. Marginal Homogeneity Model (eye grade)
EXAMPLE  4. Loglinear Models under Different Sampling Plans (bike helmet use)
EXAMPLE  5. Loglinear Models of Independence (bike helmet use)
EXAMPLE  6. Logit Models under Different Sampling Plans (bike helmet use)
EXAMPLE  7. Sparse Table Example (fine tuning the fitting algorithm)
EXAMPLE  8. Given Marginal Distributions and Odds Ratios, Compute Joint Distribution in a Two-Way Table.
EXAMPLE  9. Model Constraint:  ROC Area Equal to a Constant
EXAMPLE 10. Marginal Homogeneity of Multidimensional Contingency Tables
EXAMPLE 11. Contingency Tables with Given Marginals.
EXAMPLE 12. Testing for (Conditional) Pairwise Independence
EXAMPLE 13. Kappa Regression using Direct and Indirect Smoothing Models.
EXAMPLE 14. Marginal Cumulative Probit Model (in Conjunction with Loglinear Association Model)
EXAMPLE 15. Marginal Cumulative Logit and Linear-by-Linear Loglinear Association Model (Generalized Loglinear and Indirect Smoothing Model Specification)

EXAMPLE 1.  ML fit of mean-response models  

Table 13.1 Colds in Children,  source:  Stokes, Davis, and Koch (2000), "Categorical Data Analysis Using the SAS System."

  
                        Periods with Colds
  Sex    Residence     0       1       2    Total
  ----   ---------    -------------------   -----
    F    Rural         45      64      71    180
    F    Urban         80     104     116    300
    M    Rural         84     124      82    290
    M    Urban        106     117      87    310


Models fitted below:

Random Component:

Condition on row totals--use product multinomial sampling distribution. Specifically,

(Yij1, Yij2, Yij3)  indep ~ mult(nij, pij1, pij2, pij3),  i, j = 1,2.

Using the MP notation of Lang (2002)...

Y ~ MPZ(m|ZF,n),  where Z = ZF = kronecker(diag(4), matrix(1,3,1)) and n = (180, 300, 290, 310).   The expected count vector m = D(Zn)p,  where p = (p111, p112, p113,...,p223).  As with all MP models,  p = D-1(ZZTm)m.

Systematic Components (Linear Predictor Models):

Saturated mean-response model:

Mij = b(0) + b(SEX)i + b(RES)j + b(SEX*RES)ij,

No-interaction mean-response model:

Mij = b(0) + b(SEX)i + b(RES)j.

Here, Mij  = 0*pij1 + 1*pij2 + 2*pij3 = mean number of periods for SEX=i, RES=j.

These mean-response models have the generic form F(p) = Xb.  These are examples of  0-order HLP models. 

IMPORTANT:  mph.fit requires the HLP link function to be defined in terms of the expected counts m.  Therefore, you must define the link as L(m) = F(p(m)) for this example.  Recall  p(m) = D-1(ZZTm)m  [see primary references]. 

y <- scan()
45 64 71
80 104 116
84 124 82
106 117 87

Z <-  ZF <- kronecker(diag(4),matrix(1,3,1))

A <- kronecker(diag(4),matrix(c(0,1,2),1,3))

L.fct <- function(m) {
  p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
  A%*%p
}

X <- scan()
1 1 1 1
1 1 0 0
1 0 1 0
1 0 0 0

X <- matrix(X,4,4,byrow=T)


# SATURATED MEAN-RESPONSE MODEL....

a <- mph.fit(y,Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 0.0004368583  norm.score= 2.184184e-06 
  iter= 2  norm.diff= 3.384649e-08  norm.score= 6.852226e-14 

 Time Elapsed: 0.03 seconds


> mph.summary(a)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 0 ):  Gsq =  0
    Pearson's Score Stat  (df= 0 ):  Xsq =  0
    Generalized Wald Stat (df= 0 ):  Wsq =  0


LINEAR PREDICTOR MODEL RESULTS...
             BETA StdErr(BETA)    Z-ratio     p-value
beta1  0.93870968   0.04467893 21.0101193 0.000000000
beta2  0.18129032   0.06423383  2.8223497 0.004767317
beta3  0.05439377   0.06300701  0.8632971 0.387974137
beta4 -0.02994933   0.09779569 -0.3062438 0.759418994

       OBS LINK   ML LINK  StdErr(L) LINK RESID
link1 1.1444444 1.1444444 0.05885860          0
link2 1.1200000 1.1200000 0.04614952          0
link3 0.9931034 0.9931034 0.04442608          0
link4 0.9387097 0.9387097 0.04467893          0


CONVERGENCE STATISTICS...
    iterations = 2
    norm.diff  = 3.38465e-08
    norm.score = 6.85223e-14
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 


# NO-INTERACTION MEAN-RESPONSE MODEL...

 b <- mph.fit(y,Z,ZF,L.fct=L.fct,X=X[,1:3])

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 29.7406  norm.score= 0.01901283 
  iter= 2  norm.diff= 0.02408615  norm.score= 0.0001054862 
  iter= 3  norm.diff= 0.001056702  norm.score= 7.894742e-07 
  iter= 4  norm.diff= 3.624985e-06  norm.score= 9.750705e-09 

 Time Elapsed: 0.14 seconds

> mph.summary(b,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 1 ):  Gsq =  0.09383 (p =  0.75936 )
    Pearson's Score Stat  (df= 1 ):  Xsq =  0.09386 (p =  0.75933 )
    Generalized Wald Stat (df= 1 ):  Wsq =  0.09379 (p =  0.75942 )


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)    Z-ratio      p-value
beta1 0.94497081   0.03975182 23.7717642 0.0000000000
beta2 0.16834269   0.04842903  3.4760696 0.0005088201
beta3 0.04194803   0.04817685  0.8707093 0.3839129012

       OBS LINK   ML LINK  StdErr(L) LINK RESID
link1 1.1444444 1.1552615 0.04695725 -0.3063624
link2 1.1200000 1.1133135 0.04070944  0.3063624
link3 0.9931034 0.9869188 0.03957190  0.3063624
link4 0.9387097 0.9449708 0.03975182 -0.3063624


CELL-SPECIFIC STATISTICS...
    OBS        FV StdErr(FV)      PROB StdErr(PROB) ADJUSTED RESIDS
y1   45  44.11273   4.991427 0.2450707   0.02773015       0.3063624
y2   64  63.82746   6.393533 0.3545970   0.03551963       0.3063624
y3   71  72.05981   5.589724 0.4003323   0.03105402      -0.3063624
y4   80  80.94135   7.047109 0.2698045   0.02349036      -0.3063624
y5  104 104.12325   8.235446 0.3470775   0.02745149      -0.3063624
y6  116 114.93540   7.669822 0.3831180   0.02556607       0.3063624
y7   84  84.90553   7.163143 0.2927777   0.02470049      -0.3063624
y8  124 123.98247   8.424577 0.4275258   0.02905027       0.3063624
y9   82  81.11200   7.072743 0.2796965   0.02438877       0.3063624
y10 106 104.99696   7.662589 0.3386999   0.02471803       0.3063624
y11 117 117.06512   8.533036 0.3776294   0.02752592      -0.3063624
y12  87  87.93791   7.322569 0.2836707   0.02362119      -0.3063624


CONVERGENCE STATISTICS...
    iterations = 4
    norm.diff  = 3.62498e-06
    norm.score = 9.7507e-09
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
  A%*%p
}

Derivative of Transpose Link Function derLt.fct: 
[1] "Numerical derivatives used."

Linear Predictor Model Design Matrix  X: 
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    0
[3,]    1    0    1
[4,]    1    0    0

U = Orthogonal Complement of X: 
           [,1]
[1,]  0.1052927
[2,] -0.1052927
[3,] -0.1052927
[4,]  0.1052927

Constraint Function  h.fct:
function(m) {
          t(U)%*%L.fct(m)
       }
<environment: 01AC7678>

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
      [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    1    0    0    0
 [3,]    1    0    0    0
 [4,]    0    1    0    0
 [5,]    0    1    0    0
 [6,]    0    1    0    0
 [7,]    0    0    1    0
 [8,]    0    0    1    0
 [9,]    0    0    1    0
[10,]    0    0    0    1
[11,]    0    0    0    1
[12,]    0    0    0    1

Sampling Constraint Matrix ZF:
      [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    1    0    0    0
 [3,]    1    0    0    0
 [4,]    0    1    0    0
 [5,]    0    1    0    0
 [6,]    0    1    0    0
 [7,]    0    0    1    0
 [8,]    0    0    1    0
 [9,]    0    0    1    0
[10,]    0    0    0    1
[11,]    0    0    0    1
[12,]    0    0    0    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

EXAMPLE 2.  Dispersion Linear Trend Model. 

Marital Status of Danes.  Source:  Lloyd, C. (1999), "Statistical Analysis of Categorical Data," p. 72.
 
   
   Age       Single    Married   Divorced    Total
   17-21         17         1         0        18
   21-25         16         8         0        24
   25-30          8        17         1        26
   30-40          6        22         4        32
   40-50          5        21         6        32
   50-60          3        17         8        28
   60-70          2         8         6        16
    70+           1         3         5         9

Models Fitted Below:

Random Component:  Conditional on row totals,

(Yi1, Yi2, Yi3) indep ~ mult(ni, pi1, pi2, pi3), i=1,...,8.

That is,  Y ~ MPZ(m|ZF,n),  where Z = ZF = kronecker(diag(8),matrix(1,3,1)),  n = (18, 24, 26, ...,16,9),  m = D(Zn)p, where p = (p11, p12, p13, p21, ..., p83).  As always for MP models, p = D-1(ZZTm)m.

Systematic Components:

Linear-in-Age Model:  Gi = b(0) + b(AGE)*xi

Saturated Model:       Gi  = bi,

where  Gi  = 1 - (pi12 + pi22 + pi32) = Gini dispersion for age category i and x = (x1,...x8) = (19,23,27,35,45,55,65,80) is the vector of AGE scores [midpoints of the age intervals].

IMPORTANT: These are HLP models of the generic form F(p) = Xb.  As in example 1, the user is reminded that the HLP model link must be defined in terms of the expected counts m.  Therefore, use L(m) = F(p(m)),  where as always p(m) = D-1(ZZTm)m.

y <- scan()
17 1 0
16 8 0
8 17 1
6 22 4
5 21 6
3 17 8
2 8 6
1 3 5

y <- matrix(y,24,1)
Z <- ZF <- kronecker(diag(8),matrix(1,3,1))


L.fct <- function(m) {
  p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
  L <- c()
  for (i in 1:8) {
      L <- rbind(L,1-(p[1+3*(i-1)]^2 + p[2+3*(i-1)]^2 + p[3+3*(i-1)]^2))
  }
  L
}


# LINEAR-IN-AGE MODEL...
X <- scan()
1 19
1 23
1 27
1 35
1 45
1 55
1 65
1 80

X <- matrix(X,8,2,byrow=T)


a <- mph.fit(y,Z,ZF,L.fct=L.fct,X=X,norm.diff.conv=2,norm.score.conv=1e-7)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 18.22503  norm.score= 14.18274 
  iter= 2  norm.diff= 17.10626  norm.score= 8.686708 
  iter= 3  norm.diff= 17.00170  norm.score= 2.545349 
  iter= 4  norm.diff= 3.405157  norm.score= 0.4992686 
  iter= 5  norm.diff= 1.492782  norm.score= 0.08570673 
  iter= 6  norm.diff= 1.411508  norm.score= 0.0669449 
  iter= 7  norm.diff= 1.425170  norm.score= 0.05527805 
   .                        .
   .                        .  
  iter= 63  norm.diff= 1.270491  norm.score= 2.540270e-07 
  iter= 64  norm.diff= 1.270491  norm.score= 1.918706e-07 
  iter= 65  norm.diff= 1.270491  norm.score= 1.443599e-07 
  iter= 66  norm.diff= 1.270491  norm.score= 1.089043e-07 
  iter= 67  norm.diff= 1.270491  norm.score= 8.184554e-08 

 Time Elapsed: 3.94 seconds
REMARK:  The fact that norm.diff does not converge to 0 usually means that the ML estimates do not exist owing to 0 counts; indeed, for this example, a look at the fitted values indicates that the iterate estimates approach 0 for the 6th cell.   N.B. To get convergence, I set norm.diff.conv = 2.0.  The  resulting estimates can be regarded as "extended" ML estimates, which allow for fitted values equal to 0, the boundary value.


> mph.summary(a,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 6 ):  Gsq =  6.61522 (p =  0.3579 )
    Pearson's Score Stat  (df= 6 ):  Xsq =  4.96987 (p =  0.54768 )
    Generalized Wald Stat (df= 6 ):  Wsq =  10.09688 (p =  0.12063 )


LINEAR PREDICTOR MODEL RESULTS...
             BETA StdErr(BETA)  Z-ratio      p-value
beta1 0.331544223  0.072438461 4.576909 4.718973e-06
beta2 0.003625855  0.001398731 2.592245 9.535181e-03

       OBS LINK   ML LINK  StdErr(L) LINK RESID
link1 0.1049383 0.4004355 0.04914697 -2.7342475
link2 0.4444444 0.4149389 0.04471273  0.4734820
link3 0.4763314 0.4294423 0.04056645  0.6290227
link4 0.4765625 0.4584491 0.03356001  0.2198400
link5 0.5097656 0.4947077 0.02879638  0.1907121
link6 0.5382653 0.5309662 0.03038880  0.1102546
link7 0.5937500 0.5672248 0.03753687  0.4409538
link8 0.5679012 0.6216126 0.05358163 -1.0999550


CELL-SPECIFIC STATISTICS...
    OBS           FV   StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1   17 1.342063e+01 7.917228e-01 7.455907e-01 4.398460e-02    2.143869e+00
y2    1 3.642503e+00 1.177341e+00 2.023613e-01 6.540785e-02   -2.143869e+00
y3    0 9.368645e-01 8.349467e-01 5.204803e-02 4.638593e-02   -2.143869e+00
y4   16 1.694951e+01 1.300864e+00 7.062294e-01 5.420266e-02   -5.237165e-01
y5    8 7.050495e+00 1.300864e+00 2.937706e-01 5.420266e-02    5.237165e-01
y6    0 4.166042e-40 2.041088e-20 1.735851e-41 8.504535e-22   -4.166042e-40
y7    8 6.835339e+00 1.495050e+00 2.628976e-01 5.750194e-02    6.956249e-01
y8   17 1.839519e+01 1.165229e+00 7.075074e-01 4.481652e-02   -6.956249e-01
y9    1 7.694699e-01 7.980423e-01 2.959500e-02 3.069394e-02    6.956249e-01
y10   6 5.696848e+00 1.693282e+00 1.780265e-01 5.291507e-02    2.249922e-01
y11  22 2.253683e+01 9.857207e-01 7.042760e-01 3.080377e-02   -2.249922e-01
y12   4 3.766318e+00 1.498098e+00 1.176974e-01 4.681558e-02    2.249922e-01
y13   5 4.768420e+00 1.627600e+00 1.490131e-01 5.086250e-02    1.951100e-01
y14  21 2.148669e+01 9.149380e-01 6.714590e-01 2.859181e-02   -1.951100e-01
y15   6 5.744893e+00 1.733194e+00 1.795279e-01 5.416232e-02    1.951100e-01
y16   3 2.894242e+00 1.305963e+00 1.033658e-01 4.664152e-02    1.121326e-01
y17  17 1.725375e+01 1.225208e+00 6.162052e-01 4.375743e-02   -1.121326e-01
y18   8 7.852012e+00 1.976947e+00 2.804290e-01 7.060526e-02    1.121326e-01
y19   2 1.607440e+00 8.987154e-01 1.004650e-01 5.616971e-02    4.913694e-01
y20   8 8.718400e+00 1.352845e+00 5.449000e-01 8.455284e-02   -4.913694e-01
y21   6 5.674160e+00 1.795040e+00 3.546350e-01 1.121900e-01    4.913694e-01
y22   1 1.577521e+00 9.291630e-01 1.752801e-01 1.032403e-01   -8.729608e-01
y23   3 3.157069e+00 1.420296e+00 3.507855e-01 1.578107e-01   -8.729608e-01
y24   5 4.265410e+00 1.239264e+00 4.739345e-01 1.376960e-01    8.729608e-01


CONVERGENCE STATISTICS...
    iterations = 67
    norm.diff  = 1.27049
    norm.score = 8.18455e-08
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
  L <- c()
  for (i in 1:8) {
      L <- rbind(L,1-(p[1+3*(i-1)]^2 + p[2+3*(i-1)]^2 + p[3+3*(i-1)]^2))
  }
  L
}

Derivative of Transpose Link Function derLt.fct: 
[1] "Numerical derivatives used."

Linear Predictor Model Design Matrix  X: 
     [,1] [,2]
[1,]    1   19
[2,]    1   23
[3,]    1   27
[4,]    1   35
[5,]    1   45
[6,]    1   55
[7,]    1   65
[8,]    1   80

U = Orthogonal Complement of X: 
           [,1]        [,2]       [,3]       [,4]        [,5]       [,6]
[1,]  2.2249920  3.14366770 -0.8965420 -3.5567993  0.37289294  3.5889480
[2,]  0.2449733 -3.80315112 -3.3172856  2.5217949 -3.98805150 -4.5907728
[3,] -0.8528412 -0.91969089  0.7624972 -2.2155426  2.52343416 -2.8872172
[4,] -2.3658711  1.48479685  3.4409134  1.2757091  0.34933676  3.0953646
[5,]  2.4948821  1.05102490  1.9742120  3.8240440  2.30338745 -0.2170431
[6,] -3.2969343 -1.57599767 -1.0975353  1.9734021 -1.06265310  1.8826530
[7,] -0.1947147  0.63718086  0.4574562 -3.3291532  0.07053789  1.1339674
[8,]  1.7455140 -0.01783064 -1.3237160 -0.4934551 -0.56888460 -2.0058999

Constraint Function  h.fct:
function(m) {
          t(U)%*%L.fct(m)
       }
<environment: 00F10E84>

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]    1    0    0    0    0    0    0    0
 [2,]    1    0    0    0    0    0    0    0
 [3,]    1    0    0    0    0    0    0    0
 [4,]    0    1    0    0    0    0    0    0
 [5,]    0    1    0    0    0    0    0    0
 [6,]    0    1    0    0    0    0    0    0
 [7,]    0    0    1    0    0    0    0    0
 [8,]    0    0    1    0    0    0    0    0
 [9,]    0    0    1    0    0    0    0    0
[10,]    0    0    0    1    0    0    0    0
[11,]    0    0    0    1    0    0    0    0
[12,]    0    0    0    1    0    0    0    0
[13,]    0    0    0    0    1    0    0    0
[14,]    0    0    0    0    1    0    0    0
[15,]    0    0    0    0    1    0    0    0
[16,]    0    0    0    0    0    1    0    0
[17,]    0    0    0    0    0    1    0    0
[18,]    0    0    0    0    0    1    0    0
[19,]    0    0    0    0    0    0    1    0
[20,]    0    0    0    0    0    0    1    0
[21,]    0    0    0    0    0    0    1    0
[22,]    0    0    0    0    0    0    0    1
[23,]    0    0    0    0    0    0    0    1
[24,]    0    0    0    0    0    0    0    1

Sampling Constraint Matrix ZF:
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 [1,]    1    0    0    0    0    0    0    0
 [2,]    1    0    0    0    0    0    0    0
 [3,]    1    0    0    0    0    0    0    0
 [4,]    0    1    0    0    0    0    0    0
 [5,]    0    1    0    0    0    0    0    0
 [6,]    0    1    0    0    0    0    0    0
 [7,]    0    0    1    0    0    0    0    0
 [8,]    0    0    1    0    0    0    0    0
 [9,]    0    0    1    0    0    0    0    0
[10,]    0    0    0    1    0    0    0    0
[11,]    0    0    0    1    0    0    0    0
[12,]    0    0    0    1    0    0    0    0
[13,]    0    0    0    0    1    0    0    0
[14,]    0    0    0    0    1    0    0    0
[15,]    0    0    0    0    1    0    0    0
[16,]    0    0    0    0    0    1    0    0
[17,]    0    0    0    0    0    1    0    0
[18,]    0    0    0    0    0    1    0    0
[19,]    0    0    0    0    0    0    1    0
[20,]    0    0    0    0    0    0    1    0
[21,]    0    0    0    0    0    0    1    0
[22,]    0    0    0    0    0    0    0    1
[23,]    0    0    0    0    0    0    0    1
[24,]    0    0    0    0    0    0    0    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

#SATURATED MODEL...

> b <- mph.fit(y,Z,ZF,L.fct=L.fct,X=diag(8),norm.diff.conv=2,norm.score.conv=1e-7)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1.414345  norm.score= 0.005203488 
  iter= 2  norm.diff= 1.414214  norm.score= 0.00191393 
  iter= 3  norm.diff= 1.414214  norm.score= 0.0007040955 
  iter= 4  norm.diff= 1.414214  norm.score= 0.0002590222 
  iter= 5  norm.diff= 1.414214  norm.score= 9.528896e-05 
  iter= 6  norm.diff= 1.414214  norm.score= 3.505485e-05 
  iter= 7  norm.diff= 1.414214  norm.score= 1.289596e-05 
  iter= 8  norm.diff= 1.414214  norm.score= 4.744158e-06 
  iter= 9  norm.diff= 1.414214  norm.score= 1.745278e-06 
  iter= 10  norm.diff= 1.414214  norm.score= 6.42052e-07 
  iter= 11  norm.diff= 1.414214  norm.score= 2.361977e-07 
  iter= 12  norm.diff= 1.414214  norm.score= 8.689228e-08 

 Time Elapsed: 0.1 seconds


> mph.summary(b,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 0 ):  Gsq =  0
    Pearson's Score Stat  (df= 0 ):  Xsq =  0
    Generalized Wald Stat (df= 0 ):  Wsq =  0


LINEAR PREDICTOR MODEL RESULTS...
           BETA StdErr(BETA)  Z-ratio      p-value
beta1 0.1049383   0.09598275 1.093304 2.742606e-01
beta2 0.4444444   0.06415003 6.928203 4.262146e-12
beta3 0.4763314   0.07284080 6.539348 6.178769e-11
beta4 0.4765625   0.08624766 5.525512 3.285263e-08
beta5 0.5097656   0.08116345 6.280729 3.369891e-10
beta6 0.5382653   0.07087326 7.594759 3.086420e-14
beta7 0.5937500   0.06051536 9.811558 0.000000e+00
beta8 0.5679012   0.10147184 5.596639 2.185471e-08

       OBS LINK   ML LINK  StdErr(L) LINK RESID
link1 0.1049383 0.1049383 0.09598275          0
link2 0.4444444 0.4444444 0.06415003          0
link3 0.4763314 0.4763314 0.07284080          0
link4 0.4765625 0.4765625 0.08624766          0
link5 0.5097656 0.5097656 0.08116345          0
link6 0.5382653 0.5382653 0.07087326          0
link7 0.5937500 0.5937500 0.06051536          0
link8 0.5679012 0.5679012 0.10147184          0


CELL-SPECIFIC STATISTICS...
    OBS           FV   StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1   17 1.700000e+01 0.9718253158 9.444444e-01 5.399030e-02               0
y2    1 1.000000e+00 0.9718253158 5.555556e-02 5.399030e-02               0
y3    0 6.144212e-08 0.0002478752 3.413451e-09 1.377085e-05               0
y4   16 1.600000e+01 2.3094010768 6.666667e-01 9.622504e-02               0
y5    8 8.000000e+00 2.3094010768 3.333333e-01 9.622504e-02               0
y6    0 6.144212e-08 0.0002478752 2.560088e-09 1.032813e-05               0
y7    8 8.000000e+00 2.3533936217 3.076923e-01 9.051514e-02               0
y8   17 1.700000e+01 2.4258226202 6.538462e-01 9.330087e-02               0
y9    1 1.000000e+00 0.9805806757 3.846154e-02 3.771464e-02               0
y10   6 6.000000e+00 2.2079402166 1.875000e-01 6.899813e-02               0
y11  22 2.200000e+01 2.6220221204 6.875000e-01 8.193819e-02               0
y12   4 4.000000e+00 1.8708286934 1.250000e-01 5.846340e-02               0
y13   5 5.000000e+00 2.0539595906 1.562500e-01 6.418624e-02               0
y14  21 2.100000e+01 2.6867731575 6.562500e-01 8.396166e-02               0
y15   6 6.000000e+00 2.2079402166 1.875000e-01 6.899813e-02               0
y16   3 3.000000e+00 1.6366341768 1.071429e-01 5.845122e-02               0
y17  17 1.700000e+01 2.5842932164 6.071429e-01 9.229619e-02               0
y18   8 8.000000e+00 2.3904572187 2.857143e-01 8.537347e-02               0
y19   2 2.000000e+00 1.3228756555 1.250000e-01 8.267973e-02               0
y20   8 8.000000e+00 2.0000000000 5.000000e-01 1.250000e-01               0
y21   6 6.000000e+00 1.9364916731 3.750000e-01 1.210307e-01               0
y22   1 1.000000e+00 0.9428090416 1.111111e-01 1.047566e-01               0
y23   3 3.000000e+00 1.4142135624 3.333333e-01 1.571348e-01               0
y24   5 5.000000e+00 1.4907119850 5.555556e-01 1.656347e-01               0


CONVERGENCE STATISTICS...
    iterations = 12
    norm.diff  = 1.41421
    norm.score = 8.68923e-08
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

EXAMPLE 3.  ML Fit of Marginal Homogeneity Model

Display of unaided distance vision data from 30-39 year old FEMALE employees of United Kingdom Royal Ordnance factories during 1943-1946 (Kendall and Stuart 1961, "The Advanced  Theory of Statistics, Vol. 2: Inference and Relationship," p564 and p 586. See also, Stokes, Davis, and Koch (2000), "Categorical Data Analysis using the SAS System," p. 443)

  
                      Left Grade
   Right Grade     1     2      3    4  Total
   -----------   ---------------------  -----
       1         1520   266   124   66   1976
       2          234  1512   432   78   2256
       3          117   362  1772  205   2456
       4           36    82   179  492    789
                 ---------------------  -----
    Total        1907  2222  2507  841   7477
 

Model of Marginal Homogeneity:  

mi+ = m+i,  i=1,2,3.  That is, h(m) = [m1+ - m+1, m2+- m+2, m3+ - m+3]T = 0.


y <- scan()
1520 266 124 66
234 1512 432 78
117 362 1772 205
36 82 179 492

y <- matrix(y,16,1)
Z <- ZF <- matrix(1,16,1)
Mrow <- kronecker(diag(4),matrix(1,1,4))[-4,]
Mcol <-  kronecker(matrix(1,1,4),diag(4))[-4,]
h.fct <- function(m) {
   Mrow%*%m - Mcol%*%m
}


 a <- mph.fit(y,Z,ZF,h.fct=h.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 0.3885707  norm.score= 1.850956 
  iter= 2  norm.diff= 0.02150224  norm.score= 0.2007438 
  iter= 3  norm.diff= 0.001983644  norm.score= 0.01447043 
  iter= 4  norm.diff= 0.0002421110  norm.score= 0.001959174 
  iter= 5  norm.diff= 3.183782e-05  norm.score= 0.0002554976 
  iter= 6  norm.diff= 4.726952e-06  norm.score= 3.758008e-05 
  iter= 7  norm.diff= 6.464685e-07  norm.score= 5.183805e-06 

 Time Elapsed: 0.05 seconds

# Alternative to numerically computing the derivative of h(), input the derivative explicitly...

derht.fct <- function(m) {
   t(Mrow - Mcol)
}

 
> b <- mph.fit(y,Z,ZF,h.fct=h.fct,derht.fct=derht.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 0.3885707  norm.score= 1.850956 
  iter= 2  norm.diff= 0.02150224  norm.score= 0.2007438 
  iter= 3  norm.diff= 0.001983644  norm.score= 0.01447042 
  iter= 4  norm.diff= 0.000242111  norm.score= 0.00195917 
  iter= 5  norm.diff= 3.183773e-05  norm.score= 0.0002554976 
  iter= 6  norm.diff= 4.726968e-06  norm.score= 3.757911e-05 
  iter= 7  norm.diff= 6.464497e-07  norm.score= 5.186721e-06 

 Time Elapsed: 0.03 seconds


> mph.summary(b,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 3 ):  Gsq =  11.9872 (p =  0.0074271 )
    Pearson's Score Stat  (df= 3 ):  Xsq =  11.9698 (p =  0.0074873 )
    Generalized Wald Stat (df= 3 ):  Wsq =  11.95657 (p =  0.0075334 )


CELL-SPECIFIC STATISTICS...
     OBS         FV StdErr(FV)        PROB StdErr(PROB) ADJUSTED RESIDS
y1  1520 1520.00000  34.799412 0.203290090 0.0046541944   -4.547474e-13
y2   266  252.48209  12.609806 0.033767833 0.0016864793    1.466663e+00
y3   124  111.84293   9.507314 0.014958263 0.0012715412    2.733414e+00
y4    66   56.96589   6.961095 0.007618816 0.0009310011    3.179168e+00
y5   234  247.23710  12.554104 0.033066350 0.0016790295   -1.466663e+00
y6  1512 1512.00000  34.731011 0.202220142 0.0046450463    6.821210e-13
y7   432  409.41752  15.228459 0.054756923 0.0020367071    1.813324e+00
y8    78   70.58517   7.752804 0.009440307 0.0010368870    2.367027e+00
y9   117  131.26859  10.085383 0.017556318 0.0013488542   -2.733414e+00
y10  362  383.13268  15.089140 0.051241497 0.0020180742   -1.813324e+00
y11 1772 1772.00000  36.770200 0.236993447 0.0049177745   -4.547474e-13
y12  205  195.25849  11.211717 0.026114550 0.0014994941    1.213367e+00
y13   36   42.78522   6.163218 0.005722244 0.0008242902   -3.179165e+00
y14   82   91.62502   8.600436 0.012254249 0.0011502523   -2.367027e+00
y15  179  188.39931  11.119551 0.025197179 0.0014871674   -1.213367e+00
y16  492  492.00000  21.438879 0.065801792 0.0028673102    1.136868e-13


CONVERGENCE STATISTICS...
    iterations = 7
    norm.diff  = 6.4645e-07
    norm.score = 5.18672e-06
    Original counts used.

MODEL INFORMATION...

Constraint Function  h.fct:
function(m) {
   Mrow%*%m - Mcol%*%m
}

Derivative of Transpose Constraint Function derht.fct:
function(m) {
   t(Mrow - Mcol)
}

Population Matrix Z: 
      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    1
 [6,]    1
 [7,]    1
 [8,]    1
 [9,]    1
[10,]    1
[11,]    1
[12,]    1
[13,]    1
[14,]    1
[15,]    1
[16,]    1

Sampling Constraint Matrix ZF:
      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    1
 [6,]    1
 [7,]    1
 [8,]    1
 [9,]    1
[10,]    1
[11,]    1
[12,]    1
[13,]    1
[14,]    1
[15,]    1
[16,]    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

EXAMPLE 4.  Loglinear models under different sampling plans

Bicycle Data  (Table 16.4, p. 557, Stokes, Davis, and Koch (2000), "Categorical Data Analysis using the SAS System")

                  Wearing Helmet?
Bicycle Type      Yes       No
   Mountain       34        32
   Other          10        24

y <- scan()
34 32
10 24

y <- matrix(y,4,1)
Z.full.mult <- ZF.full.mult <- matrix(1,4,1)
Z.full.Poisson <- matrix(1,4,1); ZF.full.Poisson <- 0
Z.prod.mult <- ZF.prod.mult <- kronecker(diag(2),matrix(1,2,1))

L.fct <- function(m) {
  log(m)
}

X <- scan()
1 1 1 1
1 1 0 0
1 0 1 0
1 0 0 0

X <- matrix(X,4,4,byrow=T)
a1 <-
mph.fit(y,Z.full.mult,ZF.full.mult,L.fct=L.fct,X=X)
a2 <-
mph.fit(y,Z.full.Poisson,ZF.full.Poisson,L.fct=L.fct,X=X)
a3 <-
mph.fit(y,Z.prod.mult,ZF.prod.mult,L.fct=L.fct,X=X)


> cbind(a1$beta,a2$beta,a3$beta)

            BETA       BETA       BETA
beta1  3.1780538  3.1780538  3.1780538
beta2  0.2876821  0.2876821  0.2876821
beta3 -0.8754687 -0.8754687 -0.8754687
beta4  0.9360934  0.9360934  0.9360934


> cbind(sqrt(diag(a1$covbeta)),sqrt(diag(a2$covbeta)),sqrt(diag(a3$covbeta)))

           [,1]      [,2]      [,3]
beta1 0.1779513 0.2041241 0.1107019
beta2 0.2700309 0.2700309 0.1683846
beta3 0.3763863 0.3763863 0.3763863
beta4 0.4498093 0.4498093 0.4498093


> cbind(a1$p,sqrt(diag(a1$covp)),a2$p,sqrt(diag(a2$covp)),a3$p,sqrt(diag(a3$covp)))

   PROB            PROB                 PROB           
p1 0.34 0.04737088 0.34 0.04737088 0.5151515 0.06151748
p2 0.32 0.04664762 0.32 0.04664762 0.4848485 0.06151748
p3 0.10 0.03000000 0.10 0.03000000 0.2941176 0.07814249
p4 0.24 0.04270831 0.24 0.04270831 0.7058824 0.07814249

 
> cbind(a1$m,sqrt(diag(a1$covm)),a2$m,sqrt(diag(a2$covm)),a3$m,sqrt(diag(a3$covm)))

   FV          FV          FV         
m1 34 4.737088 34 5.830952 34 4.060154
m2 32 4.664762 32 5.656854 32 4.060154
m3 10 3.000000 10 3.162278 10 2.656845
m4 24 4.270831 24 4.898979 24 2.656845

# THE SUMMARIES ...
> mph.summary(a1,F,T)
OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 0 ):  Gsq =  0
    Pearson's Score Stat  (df= 0 ):  Xsq =  0
    Generalized Wald Stat (df= 0 ):  Wsq =  0


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)   Z-ratio    p-value
beta1  3.1780538    0.1779513 17.859121 0.00000000
beta2  0.2876821    0.2700309  1.065367 0.28670971
beta3 -0.8754687    0.3763863 -2.325984 0.02001938
beta4  0.9360934    0.4498093  2.081089 0.03742574

      OBS LINK  ML LINK StdErr(L) LINK RESID
link1 3.526361 3.526361 0.1393261          0
link2 3.465736 3.465736 0.1457738          0
link3 2.302585 2.302585 0.3000000          0
link4 3.178054 3.178054 0.1779513          0


CONVERGENCE STATISTICS...
    iterations = 2
    norm.diff  = 5.10992e-07
    norm.score = 1.24962e-12
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  log(m)
}

Derivative of Transpose Link Function derLt.fct: 
[1] "Numerical derivatives used."

Linear Predictor Model Design Matrix  X: 
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    0    0
[3,]    1    0    1    0
[4,]    1    0    0    0

U = Orthogonal Complement of X: 

   0 

Constraint Function  h.fct:
[1] 0

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1

Sampling Constraint Matrix ZF:
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

> mph.summary(a2,F,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 0 ):  Gsq =  0
    Pearson's Score Stat  (df= 0 ):  Xsq =  0
    Generalized Wald Stat (df= 0 ):  Wsq =  0


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)   Z-ratio    p-value
beta1  3.1780538    0.2041241 15.569221 0.00000000
beta2  0.2876821    0.2700309  1.065367 0.28670971
beta3 -0.8754687    0.3763863 -2.325984 0.02001938
beta4  0.9360934    0.4498093  2.081089 0.03742574

      OBS LINK  ML LINK StdErr(L) LINK RESID
link1 3.526361 3.526361 0.1714986          0
link2 3.465736 3.465736 0.1767767          0
link3 2.302585 2.302585 0.3162278          0
link4 3.178054 3.178054 0.2041241          0


CONVERGENCE STATISTICS...
    iterations = 2
    norm.diff  = 5.10992e-07
    norm.score = 1.24962e-12
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  log(m)
}

Derivative of Transpose Link Function derLt.fct: 
[1] "Numerical derivatives used."

Linear Predictor Model Design Matrix  X: 
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    0    0
[3,]    1    0    1    0
[4,]    1    0    0    0

U = Orthogonal Complement of X: 

   0 

Constraint Function  h.fct:
[1] 0

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1

Sampling Constraint Matrix ZF:
     [,1]
[1,]    0


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

 
> mph.summary(a3,F,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 0 ):  Gsq =  0
    Pearson's Score Stat  (df= 0 ):  Xsq =  0
    Generalized Wald Stat (df= 0 ):  Wsq =  0


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)   Z-ratio    p-value
beta1  3.1780538    0.1107019 28.708224 0.00000000
beta2  0.2876821    0.1683846  1.708482 0.08754700
beta3 -0.8754687    0.3763863 -2.325984 0.02001938
beta4  0.9360934    0.4498093  2.081089 0.03742574

      OBS LINK  ML LINK StdErr(L) LINK RESID
link1 3.526361 3.526361 0.1194163          0
link2 3.465736 3.465736 0.1268798          0
link3 2.302585 2.302585 0.2656845          0
link4 3.178054 3.178054 0.1107019          0


CONVERGENCE STATISTICS...
    iterations = 2
    norm.diff  = 5.10992e-07
    norm.score = 1.24962e-12
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  log(m)
}

Derivative of Transpose Link Function derLt.fct: 
[1] "Numerical derivatives used."

Linear Predictor Model Design Matrix  X: 
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    0    0
[3,]    1    0    1    0
[4,]    1    0    0    0

U = Orthogonal Complement of X: 

   0 

Constraint Function  h.fct:
[1] 0

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
     [,1] [,2]
[1,]    1    0
[2,]    1    0
[3,]    0    1
[4,]    0    1

Sampling Constraint Matrix ZF:
     [,1] [,2]
[1,]    1    0
[2,]    1    0
[3,]    0    1
[4,]    0    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

EXAMPLE 5.  Loglinear models of independence

Bicycle Data  (Table 16.4, p. 557, Stokes, Davis, and Koch (2000), "Categorical Data Analysis using the SAS System")

                  Wearing Helmet?
Bicycle Type      Yes       No
   Mountain       34        32
   Other          10        24
 

Loglinear Model of Independence:

log(mij) = b(0) + b(row)i  + b(col)j ;   generically  log(m) = Xb. 

or equivalently,

h(m) = UTlog(m) = 0,   where   U is an orthogonal complement of X.

 

y <- scan()
34 32
10 24

y <- matrix(y,4,1)
Z <- ZF <- matrix(1,4,1)

L.fct <- function(m) {
  log(m)
}

derLt.fct <- function(m) {
  diag(c(1/m))
}

X <- scan()
1 1 1 1
1 1 0 0
1 0 1 0
1 0 0 0

X <- matrix(X,4,4,byrow=T)
X.indep <- X[,-4]

 
> a <- mph.fit(y,Z,ZF,L.fct=L.fct, derLt.fct=derLt.fct, X=X.indep)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 3.653846  norm.score= 1.402688 
  iter= 2  norm.diff= 0.2669929  norm.score= 0.03072015 
  iter= 3  norm.diff= 0.003951403  norm.score= 1.864493e-05 
  iter= 4  norm.diff= 2.620018e-06  norm.score= 3.047489e-10 

 Time Elapsed: 0.02 seconds
Note:  The explicit derivative of L, derLt.fct, is NOT a necessary argument.  If it were omitted, the program would simply use a numerical approximation to the derivative.

 

# Re-fit the independence model using the constraint specification...
U <- create.U(X.indep)
h.fct <- function(m) {
   t(U)%*%L.fct(m)
}

> b <- mph.fit(y,Z,ZF,h.fct=h.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1.814933  norm.score= 1.402688 
  iter= 2  norm.diff= 0.1430537  norm.score= 0.03072015 
  iter= 3  norm.diff= 0.002470891  norm.score= 1.864461e-05 
  iter= 4  norm.diff= 1.584134e-06  norm.score= 4.770012e-10 

 Time Elapsed: 0.02 seconds


> mph.summary(a,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 1 ):  Gsq =  4.55692 (p =  0.032786 )
    Pearson's Score Stat  (df= 1 ):  Xsq =  4.44938 (p =  0.034914 )
    Generalized Wald Stat (df= 1 ):  Wsq =  4.33093 (p =  0.037426 )


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)   Z-ratio     p-value
beta1  2.9465420    0.1651330 17.843448 0.000000000
beta2  0.6632942    0.2111002  3.142083 0.001677505
beta3 -0.2411621    0.2014557 -1.197097 0.231268761

      OBS LINK  ML LINK StdErr(L) LINK RESID
link1 3.526361 3.368674 0.1337116   1.947417
link2 3.465736 3.609836 0.1140555  -2.264984
link3 2.302585 2.705380 0.1792736  -2.562617
link4 3.178054 2.946542 0.1651330   1.874599


CELL-SPECIFIC STATISTICS...
   OBS    FV StdErr(FV)   PROB StdErr(PROB) ADJUSTED RESIDS
y1  34 29.04   3.882984 0.2904   0.03882984        2.109356
y2  32 36.96   4.215491 0.3696   0.04215491       -2.109356
y3  10 14.96   2.681934 0.1496   0.02681934       -2.109356
y4  24 19.04   3.144132 0.1904   0.03144132        2.109356


CONVERGENCE STATISTICS...
    iterations = 4
    norm.diff  = 2.62002e-06
    norm.score = 3.04749e-10
    Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function  L.fct:
function(m) {
  log(m)
}

Derivative of Transpose Link Function derLt.fct: 
function(m) {
  diag(c(1/m)) 
}

Linear Predictor Model Design Matrix  X: 
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    0
[3,]    1    0    1
[4,]    1    0    0

U = Orthogonal Complement of X: 
          [,1]
[1,]  1.280240
[2,] -1.280240
[3,] -1.280240
[4,]  1.280240

Constraint Function  h.fct:
function(m) {
          t(U)%*%L.fct(m)
       }
<environment: 0104FB34>

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1

Sampling Constraint Matrix ZF:
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

> mph.summary(b,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 1 ):  Gsq =  4.55692 (p =  0.032786 )
    Pearson's Score Stat  (df= 1 ):  Xsq =  4.44938 (p =  0.034914 )
    Generalized Wald Stat (df= 1 ):  Wsq =  4.33093 (p =  0.037426 )


CELL-SPECIFIC STATISTICS...
   OBS    FV StdErr(FV)   PROB StdErr(PROB) ADJUSTED RESIDS
y1  34 29.04   3.882984 0.2904   0.03882984        2.109356
y2  32 36.96   4.215491 0.3696   0.04215491       -2.109356
y3  10 14.96   2.681934 0.1496   0.02681934       -2.109356
y4  24 19.04   3.144132 0.1904   0.03144132        2.109356


CONVERGENCE STATISTICS...
    iterations = 4
    norm.diff  = 1.58413e-06
    norm.score = 4.77001e-10
    Original counts used.

MODEL INFORMATION...

Constraint Function  h.fct:
function(m) {
   t(U)%*%L.fct(m)
}

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1

Sampling Constraint Matrix ZF:
     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

EXAMPLE 6.  Logit models under different sampling plans

Bicycle Data  (Table 16.4, p. 557, Stokes, Davis, and Koch (2000), "Categorical Data Analysis using the SAS System")

                  Wearing Helmet?
Bicycle Type      Yes       No
   Mountain       34        32
   Other          10        24

 

Logit Model:

log(mi1/mi2)  = b(0) + b(row)i.

Note:

If the data are the result of a cross-sectional sample from a single population, then mij = γpij,  where pij = P(row var=i, col var = j)  and γ is the expected sample size, which may [full-multinomial sampling plan] or may not  [full-Poisson sampling plan] be fixed a priori.

If the data are the result of two row-stratified random samples, then mij = γi pij,  where pij = P(col var = j | row var = i)  and the  γi's are the expected sample sizes, which may or may not be fixed a priori.

y <- scan()
34 32
10 24

y <- matrix(y,4,1)
Z.full.mult <- ZF.full.mult <- matrix(1,4,1)
Z.full.Poisson <- matrix(1,4,1); ZF.full.Poisson <- 0
Z.prod.mult <- ZF.prod.mult <- kronecker(diag(2),matrix(1,2,1))

L.fct <- function(m) {
  rbind(log(m[1]/m[2]),log(m[3]/m[4]))
}

X <- scan()
1 1
1 0

X <- matrix(X,2,2,byrow=T)
a1 <-
mph.fit(y,Z.full.mult,ZF.full.mult,L.fct=L.fct,X=X)
a2 <-
mph.fit(y,Z.full.Poisson,ZF.full.Poisson,L.fct=L.fct,X=X)
a3 <-
mph.fit(y,Z.prod.mult,ZF.prod.mult,L.fct=L.fct,X=X)

> cbind(a1$beta,a2$beta,a3$beta)

            BETA       BETA       BETA
beta1 -0.8754687 -0.8754687 -0.8754687
beta2  0.9360934  0.9360934  0.9360934

 
> cbind(sqrt(diag(a1$covbeta)),sqrt(diag(a2$covbeta)),sqrt(diag(a3$covbeta)))

           [,1]      [,2]      [,3]
beta1 0.3763863 0.3763863 0.3763863
beta2 0.4498093 0.4498093 0.4498093

 
> cbind(a1$p,sqrt(diag(a1$covp)),a2$p,sqrt(diag(a2$covp)),a3$p,sqrt(diag(a3$covp)))

   PROB            PROB                 PROB           
p1 0.34 0.04737088 0.34 0.04737088 0.5151515 0.06151748
p2 0.32 0.04664762 0.32 0.04664762 0.4848485 0.06151748
p3 0.10 0.03000000 0.10 0.03000000 0.2941176 0.07814249
p4 0.24 0.04270831 0.24 0.04270831 0.7058824 0.07814249

 
> cbind(a1$m,sqrt(diag(a1$covm)),a2$m,sqrt(diag(a2$covm)),a3$m,sqrt(diag(a3$covm)))

   FV          FV          FV         
m1 34 4.737088 34 5.830952 34 4.060154
m2 32 4.664762 32 5.656854 32 4.060154
m3 10 3.000000 10 3.162278 10 2.656845
m4 24 4.270831 24 4.898979 24 2.656845

EXAMPLE 7. Very sparse table:  fine-tuning the fitting algorithm

y <- c(rep(0,99),1)

# > y
# [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
#

Model Fitted:  m1 = m2 = ... = m100 . (This is the model that states all the cell probabilities are 0.01.)


Z <- ZF <- matrix(1,100,1)
C <- cbind(1,-diag(99))
h.fct <- function(m) {
  C%*%m
}
 

> a <- mph.fit(y,Z,ZF,h.fct=h.fct)
 

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 9.950874  norm.score= 0.4514260 
  iter= 2  norm.diff= 9.880235  norm.score= 0.5146033 
  iter= 3  norm.diff= 10.38192  norm.score= 0.5841354 
  iter= 4  norm.diff= 15.03971  norm.score= 0.6132212 
  iter= 5  norm.diff= 34.8195  norm.score= 0.6226252 
  iter= 6  norm.diff= 92.95185  norm.score= 0.616593 
  iter= 7  norm.diff= 241.1516  norm.score= 0.4759282 
  iter= 8  norm.diff= 310.4932  norm.score= 3399782977 
  iter= 9  norm.diff= 691.0164  norm.score= 2568961 
  iter= 10  norm.diff= 10  norm.score= 945067.6 
  iter= 11  norm.diff= 10  norm.score= 347670.6 
  iter= 12  norm.diff= 10  norm.score= 127900.6 
  iter= 13  norm.diff= 9.999999  norm.score= 47051.72 
  iter= 14  norm.diff= 9.999997  norm.score= 17309.08 
  iter= 15  norm.diff= 9.999992  norm.score= 6367.371 
  iter= 16  norm.diff= 9.999978  norm.score= 2342.142 
  iter= 17  norm.diff= 9.99994  norm.score= 861.3436 
  iter= 18  norm.diff= 9.999836  norm.score= 316.5883 
  iter= 19  norm.diff= 9.999554  norm.score= 116.1845 
  iter= 20  norm.diff= 9.99879  norm.score= 42.4614 
  iter= 21  norm.diff= 9.996726  norm.score= 15.34378 
  iter= 22  norm.diff= 9.991203  norm.score= 5.378327 
  iter= 23  norm.diff= 9.976845  norm.score= 1.747325 
  iter= 24  norm.diff= 9.942698  norm.score= 0.5609135 
  iter= 25  norm.diff= 9.88638  norm.score= 0.4740398 
  iter= 26  norm.diff= 10.00577  norm.score= 0.5619719 
  iter= 27  norm.diff= 12.10580  norm.score= 0.6050876 
  iter= 28  norm.diff= 23.83653  norm.score= 0.6220261 
  iter= 29  norm.diff= 62.42647  norm.score= 0.6283873 
  iter= 30  norm.diff= 170.5125  norm.score= 0.630743 
  iter= 31  norm.diff= 465.4956  norm.score= 0.631546 
  iter= 32  norm.diff= 1266.612  norm.score= 619654.1 
Error in qr(a, tol = tol) : NA/NaN/Inf in foreign function call (arg 1)
REMARK:  The MLe's exist for this model of equal probabilities.  However, the many zero counts cause fitting problems. One way to mitigate fitting problems with zeroes is to begin the iterations using modified (non-zero) counts such as y+0.1. 

The next run does just this.  By default, the modified counts y+y.eps are used until the 5th iteration, at which time the original counts y are used.

> a <- mph.fit(y,Z,ZF,h.fct=h.fct,y.eps=0.1)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1.290908  norm.score= 0.4262671 
  iter= 2  norm.diff= 1.557817  norm.score= 0.4444104 
  iter= 3  norm.diff= 2.663936  norm.score= 0.3644358 
  iter= 4  norm.diff= 3.061799  norm.score= 0.1285554 
  iter= 5  norm.diff= 1.183960  norm.score= 0.003548943 
  iter= 6  norm.diff= 0.03237637  norm.score= 1.000010 
  iter= 7  norm.diff= 9.09092  norm.score= 0.6861117 
  iter= 8  norm.diff= 15.48088  norm.score= 0.5463486 
  iter= 9  norm.diff= 26.74028  norm.score= 0.398467 
  iter= 10  norm.diff= 32.49469  norm.score= 0.1676192 
  iter= 11  norm.diff= 16.43828  norm.score= 0.01897850 
  iter= 12  norm.diff= 1.897397  norm.score= 0.0001868761 
  iter= 13  norm.diff= 0.01868667  norm.score= 1.745609e-08 
  iter= 14  norm.diff= 1.745521e-06  norm.score= 1.580232e-14 

 Time Elapsed: 1.31 seconds
REMARK:  To speed up the fitting, we use the analytic derivative of the transpose of the constraint function in the next run...

derht.fct <- function(m) {
  t(C)
}

> b <- mph.fit(y,Z,ZF,h.fct=h.fct,derht.fct=derht.fct,y.eps=0.1)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1.290908  norm.score= 0.4262671 
  iter= 2  norm.diff= 1.557817  norm.score= 0.4444104 
  iter= 3  norm.diff= 2.663936  norm.score= 0.3644358 
  iter= 4  norm.diff= 3.061799  norm.score= 0.1285554 
  iter= 5  norm.diff= 1.183960  norm.score= 0.003548943 
  iter= 6  norm.diff= 0.03237637  norm.score= 1.000010 
  iter= 7  norm.diff= 9.09092  norm.score= 0.6861117 
  iter= 8  norm.diff= 15.48088  norm.score= 0.5463486 
  iter= 9  norm.diff= 26.74028  norm.score= 0.398467 
  iter= 10  norm.diff= 32.49469  norm.score= 0.1676192 
  iter= 11  norm.diff= 16.43828  norm.score= 0.01897850 
  iter= 12  norm.diff= 1.897397  norm.score= 0.0001868761 
  iter= 13  norm.diff= 0.01868667  norm.score= 1.745579e-08 
  iter= 14  norm.diff= 1.745491e-06  norm.score= 1.863317e-14 

 Time Elapsed: 0.78 seconds

> c(b$y)
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
> c(b$m)
 [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [16] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [31] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [46] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [61] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [76] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
 [91] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01

> mph.summary(b)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 99 ):  Gsq =  9.21034 (p =  1 )
    Pearson's Score Stat  (df= 99 ):  Xsq =  99 (p =  0.4811 )
    Generalized Wald Stat (df= 99 ):  Wsq =  0.90826 (p =  1 )

    WARNING: 100% of expected counts are less than 5. 
             Chi-square approximation may be questionable.

CONVERGENCE STATISTICS...
    iterations = 14
    norm.diff  = 1.74549e-06
    norm.score = 1.86332e-14
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 
REMARK. Notice that the CONVERGENCE STATISTICS output tells you that the original counts y (not the modified counts y+y.eps) were used.


EXAMPLE 8.   Given row marginal probs = (0.2, 0.3, 0.5),   col marginal probs = (0.4, 0.1, 0.5),  and  local odds ratios  or11 = 3, or12=2, or21=1, or22=4,   compute the joint distribution.

IMPORTANT:  mph.fit requires the constraint function h() to be defined as a function of the expected counts m. 

y <- matrix(1,9,1)  # any reasonable positive numbers will work here
Z <- ZF <- matrix(1,9,1)  
Mrow <- t(kronecker(diag(3),matrix(1,3,1)))[-3,]
Mcol <- t(kronecker(matrix(1,3,1),diag(3)))[-3,]
h.fct <- function(m) {
   p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
   or11 <- m[1]*m[5]/(m[2]*m[4])
   or12 <- m[2]*m[6]/(m[3]*m[5])
   or21 <- m[4]*m[8]/(m[5]*m[7])
   or22 <- m[5]*m[9]/(m[6]*m[8])
   rbind(Mrow%*%p - c(0.2,0.3),Mcol%*%p-c(0.4,0.1),
         or11-3,or12-2,or21-1,or22-4)
}

a <- mph.fit(y,Z,ZF,h.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 15.53185  norm.score= 69.5067 
  iter= 2  norm.diff= 57.34751  norm.score= 13.85793 
  iter= 3  norm.diff= 64.35078  norm.score= 3.126103 
  iter= 4  norm.diff= 15.00884  norm.score= 0.4788744 
  iter= 5  norm.diff= 0.8019941  norm.score= 0.05002032 
  iter= 6  norm.diff= 0.07469393  norm.score= 0.000667536 
  iter= 7  norm.diff= 0.0006000141  norm.score= 1.547494e-07 
  iter= 8  norm.diff= 1.910239e-07  norm.score= 1.007794e-10 

 Time Elapsed: 0.11 seconds
 a$p
         PROB
p1 0.15913038
p2 0.01804733
p3 0.02282230
p4 0.13631718
p5 0.04638010
p6 0.11730272
p7 0.10455244
p8 0.03557257
p9 0.35987499


EXAMPLE 9.   Model constraint: area under ROC manifest "curve" (AUC)  equal to a constant 

Neonatal Radiograph Data (source: see Table 1, Lang and Aspelund (2001), Statistical Modelling--An International Journal)
 
              Video Image Rating
               1   2   3   4   5   Total
              -------------------  -----
   Normals     6  17   4   5   1    33
   Abnormals   4   5   5  15  38    67

 
Models Fitted Below...

Random Component:  Product Multinomial

(Yi1, ..., Yi5)  indep ~ mult(ni, pi1, ... , pi5),  i=1,2.

That is, using the MP notation of Lang (2002) [see primary references], 

y = (6,17,4,...,15,38) <- Y ~ MPZ(m|ZF,n),      where

Z = ZF = kronecker(diag(2),matrix(1,5,1)),  n=(33, 67),  m = D(Zn)p.  Here  p = (p11, p12, ..., p15, p21, ..., p25).  As always for MP models, p = D-1(ZZTm)m.

Systematic Component:

Constraint: In words, "Area under the manifest ROC `curve' equals 0.9." In symbols, h(m) = AUC(m) - 0.9 = 0,   where

AUC(m) = the area under the manifest ROC "curve."  Specifically,

AUC(m) = p11*(p22 + p23 + p24 + p25) + p12*(p23  + p24 + p25) + p13*(p24 + p25) + p14*p25 + 0.5*(p11*p21 + p12*p22 + ... + p15*p25) = P(B > A) + 0.5*P(B=A),  where A ~ (p11, ..., p15)  and B ~ (p21, ..., p25).

IMPORTANT: mph.fit requires the constraint function h() to be defined as a function of the expected counts m. 

y <- scan()
6 17 4 5 1
4 5 5 15 38

y <- matrix(y,10,1)
Z <- ZF <- kronecker(diag(2),matrix(1,5,1))


h.fct <- function(m) {
   p <- diag(c(1/(Z%*%t(Z)%*%m)))%*%m
   p[1]*(p[7]+p[8]+p[9]+p[10]) +
   p[2]*(p[8]+p[9]+p[10]) +
   p[3]*(p[9]+p[10]) +
   p[4]*(p[10]) +
   0.5*(p[1]*p[6]+p[2]*p[7]+p[3]*p[8]+p[4]*p[9]+p[5]*p[10]) - 0.90
}

a <-
mph.fit(y,Z,ZF,h.fct=h.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 32.68234  norm.score= 0.7046965 
  iter= 2  norm.diff= 12.64723  norm.score= 0.1868399 
  iter= 3  norm.diff= 2.524178  norm.score= 0.01500246 
  iter= 4  norm.diff= 0.0926745  norm.score= 0.01007404 
  iter= 5  norm.diff= 0.03407739  norm.score= 0.007368128 
  iter= 6  norm.diff= 0.03375312  norm.score= 0.005641723 
  iter= 7  norm.diff= 0.02908642  norm.score= 0.004373105 
  iter= 8  norm.diff= 0.02434840  norm.score= 0.003437024 
  iter= 9  norm.diff= 0.01962449  norm.score= 0.002699434 
  iter= 10  norm.diff= 0.01580134  norm.score= 0.002132258 
   .                        .
   .                        .
  iter= 39  norm.diff= 1.741832e-05  norm.score= 2.306213e-06 
  iter= 40  norm.diff= 1.376869e-05  norm.score= 1.822186e-06 
  iter= 41  norm.diff= 1.087766e-05  norm.score= 1.440207e-06 
  iter= 42  norm.diff= 8.594076e-06  norm.score= 1.138069e-06 

 Time Elapsed: 0.53 seconds

> mph.summary(a,T,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 1 ):  Gsq =  1.94502 (p =  0.16312 )
    Pearson's Score Stat  (df= 1 ):  Xsq =  2.19474 (p =  0.13848 )
    Generalized Wald Stat (df= 1 ):  Wsq =  1.52065 (p =  0.21752 )


CELL-SPECIFIC STATISTICS...
    OBS         FV StdErr(FV)       PROB StdErr(PROB) ADJUSTED RESIDS
y1    6  6.7996921   2.259917 0.20605128   0.06848233       -1.481464
y2   17 17.8647238   2.802296 0.54135527   0.08491805       -1.481464
y3    4  3.8316273   1.836796 0.11610992   0.05566049        1.481465
y4    5  3.9681014   1.733721 0.12024550   0.05253699        1.481465
y5    1  0.5358555   0.654978 0.01623804   0.01984782        1.481462
y6    4  2.5477153   1.220591 0.03802560   0.01821778        1.481465
y7    5  3.8380537   1.732926 0.05728438   0.02586456        1.481464
y8    5  4.6833212   2.076117 0.06990032   0.03098682        1.481464
y9   15 15.2579804   3.428256 0.22773105   0.05116800       -1.481466
y10  38 40.6729294   3.567459 0.60705865   0.05324566       -1.481465


CONVERGENCE STATISTICS...
    iterations = 42
    norm.diff  = 8.59408e-06
    norm.score = 1.13807e-06
    Original counts used.

MODEL INFORMATION...

Constraint Function  h.fct:
function(m) {
   p <- diag(c(1/(Z%*%t(Z)%*%m)))%*%m
   p[1]*(p[7]+p[8]+p[9]+p[10]) + 
   p[2]*(p[8]+p[9]+p[10]) + 
   p[3]*(p[9]+p[10]) +
   p[4]*(p[10]) +
   0.5*(p[1]*p[6]+p[2]*p[7]+p[3]*p[8]+p[4]*p[9]+p[5]*p[10]) - 0.90
}

Derivative of Transpose Constraint Function derht.fct:
[1] "Numerical derivatives used."

Population Matrix Z: 
      [,1] [,2]
 [1,]    1    0
 [2,]    1    0
 [3,]    1    0
 [4,]    1    0
 [5,]    1    0
 [6,]    0    1
 [7,]    0    1
 [8,]    0    1
 [9,]    0    1
[10,]    0    1

Sampling Constraint Matrix ZF:
      [,1] [,2]
 [1,]    1    0
 [2,]    1    0
 [3,]    1    0
 [4,]    1    0
 [5,]    1    0
 [6,]    0    1
 [7,]    0    1
 [8,]    0    1
 [9,]    0    1
[10,]    0    1


FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

EXAMPLE 10.  Marginal homogeneity of multidimensional contingency tables

3x3x3 Table.  See S. Kullback, Annals of Mathematical Statistics, Vol. 42, No. 2. (Apr., 1971), pp. 594-606.

IMPORTANT: mph.fit requires the constraint function h() to be defined as a function of the expected counts m.   

y <- scan()
223 24 6 40 42 2 19 4 12
28 6 9 25 218 6 3 13 9
26 3 18 18 30 24 12 16 164

Z <- matrix(1,27,1)
ZF <- Z

M1 <-
Marg.fct(1,c(3,3,3))
M2 <-
Marg.fct(2,c(3,3,3))
M3 <-
Marg.fct(3,c(3,3,3))

#  Alternatively, the M matrices can be computed as follows:
#    M1 <- kronecker(diag(3),matrix(1,1,9)))
#    M2 <- kronecker(matrix(1,1,3),kronecker(diag(3),matrix(1,1,3)))
#    M3 <- kronecker(matrix(1,1,9),diag(3))

C <- scan()
1 0 0 -1 0 0 0 0 0
1 0 0 0 0 0 -1 0 0
0 1 0 0 -1 0 0 0 0
0 1 0 0 0 0 0 -1 0

C <- matrix(C,4,9,byrow=T)


h.fct <- function(m) {
   marg <- rbind(M1%*%m,  M2%*%m,  M3%*%m)
   C%*%marg
}


a <-
mph.fit(y=y,Z=Z,ZF=ZF,h.fct=h.fct)

mph.summary(a)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 4 ):  Gsq =  50.39723 (p =  2.983e-10 )
    Pearson's Score Stat  (df= 4 ):  Xsq =  48.9876 (p =  5.8737e-10 )
    Generalized Wald Stat (df= 4 ):  Wsq =  47.50884 (p =  1.1946e-09 )


CONVERGENCE STATISTICS...
    iterations = 32
    norm.diff  = 1.00205e-06
    norm.score = 7.81206e-06
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

EXAMPLE 11. Contingency tables with given marginals. 

3x4 Table.  See C. T. Ireland, S. Kullback, Biometrika, Vol. 55, No. 1. (Mar., 1968), pp. 179-188.

Model Fitted Below:

h(m) = [p1+, p2+, p+1, p+2, p+3] - [0.7837288, 0.14831812, 0.07827901, 0.46148631, 0.29658409] = 0.

For the cross-sectional sampling plan assumed,  p = p(m) = m/sum(m). 

IMPORTANT: mph.fit requires the constraint function h() to be defined as a function of the expected counts m. 


y <- scan()
783 7426 4709 2145
517 928 622 703
207 373 337 425

Z <- matrix(1,12,1)
ZF <- Z
M1 <-
Marg.fct(1,c(3,4))
M2 <-
Marg.fct(2,c(3,4))

# Alternatively, the M matrices can be computed as
#   M1 <- kronecker(diag(3),matrix(1,1,4))
#   M2 <- kronecker(matrix(1,1,3),diag(4))

marg <- scan()
15028 2844 1303 1501 8849 5687 3138

margprob <- marg/19175

A <- scan()
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0

A <- matrix(A,5,7,byrow=T)

h.fct <- function(m) {
    p <- m/sum(m)
    obsmp <- rbind(M1%*%p,M2%*%p)
    A%*%(obsmp - margprob)
}

a <-
mph.fit(y=y,Z=Z,ZF=ZF,h.fct=h.fct)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1886.546  norm.score= 5.037459 
  iter= 2  norm.diff= 57.67494  norm.score= 0.05186937 
  iter= 3  norm.diff= 0.6376254  norm.score= 0.0005032949 
  iter= 4  norm.diff= 0.003960461  norm.score= 1.449125e-05 
  iter= 5  norm.diff= 8.15081e-05  norm.score= 5.632072e-07 
  iter= 6  norm.diff= 2.889924e-06  norm.score= 5.794013e-08 

 Time Elapsed: 0.05 seconds

mph.summary(a,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 5 ):  Gsq =  11.30702 (p =  0.045621 )
    Pearson's Score Stat  (df= 5 ):  Xsq =  11.37322 (p =  0.044462 )
    Generalized Wald Stat (df= 5 ):  Wsq =  11.17887 (p =  0.047947 )


CELL-SPECIFIC STATISTICS...
     OBS        FV StdErr(FV)       PROB StdErr(PROB) ADJUSTED RESIDS
y1   783  771.3651   18.12330 0.04022764 0.0009451528    5.732957e-01
y2  7426 7503.4583   26.65348 0.39131464 0.0013900118   -1.247247e+00
y3  4709 4709.0003   24.36141 0.24558020 0.0012704775   -5.946086e-06
y4  2145 2044.1763   23.32243 0.10660633 0.0012162937    2.815559e+00
y5   517  528.7655   17.05606 0.02757578 0.0008894947   -7.873926e-01
y6   928  974.4394   23.27053 0.05081822 0.0012135872   -2.371696e+00
y7   622  646.1227   20.76431 0.03369610 0.0010828847   -1.735519e+00
y8   703  694.6724   20.34763 0.03622802 0.0010611540    5.210075e-01
y9   207  200.8694   12.16505 0.01047559 0.0006344225    8.603393e-01
y10  373  371.1023   15.77822 0.01935345 0.0008228539    1.769853e-01
y11  337  331.8769   15.17238 0.01730779 0.0007912586    5.230557e-01
y12  425  399.1513   15.69299 0.02081624 0.0008184086    2.149785e+00


CONVERGENCE STATISTICS...
    iterations = 6
    norm.diff  = 2.88992e-06
    norm.score = 5.79401e-08
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

cbind(c(M1%*%a$p,M2%*%a$p), margprob)

                  margprob
[1,] 0.78372881 0.78372881
[2,] 0.14831812 0.14831812
[3,] 0.06795306 0.06795306
[4,] 0.07827901 0.07827901
[5,] 0.46148631 0.46148631
[6,] 0.29658409 0.29658409
[7,] 0.16365059 0.16365059
Note: The fitted marginal probabilities do match the given marginal probabilities


EXAMPLE 12. Testing for (conditional) pairwise independence when response includes "unknown" category.

3x3x3 Table.  For a description of the model, see Michael Haber's Example 2, p. 433, in  Biometrics (in Shorter Communications), Vol. 42, No. 2. (Jun., 1986), pp. 429-435.

IMPORTANT: mph.fit requires the constraint function h() to be defined as a function of the expected counts m. 


y <- scan()
201 28 37 21 8 7 12 0 0 27 9 5
14 4 4 2 0 0 142 15 0 27 12 0 0 0 0

Z <- ZF <- matrix(1,27,1)

M12 <- Marg.fct(c(1,2),c(3,3,3))
M13 <- Marg.fct(c(1,3),c(3,3,3))
M23 <- Marg.fct(c(2,3),c(3,3,3))

# Alternatively,
#   M12 <- kronecker(diag(9),matrix(1,1,3))
#   M13 <- kronecker(diag(3),kronecker(matrix(1,1,3),diag(3)))
#   M23 <- kronecker(matrix(1,1,3),diag(9))

M <- rbind(M12,M13,M23)  
Mr <- M[-c(3,6,7,8,9,12,15,16,17,18,21,24,25,26,27),]


# Mr*p gives only those marginal probabilities that do NOT involve the last unknown category.

C <- scan()
1 -1 -1 1 0 0 0 0 0 0 0 0
0 0 0 0 1 -1 -1 1 0 0 0 0
0 0 0 0 0 0 0 0 1 -1 -1 1

C <- matrix(C,3,12,byrow=T)

L.fct <- function(m) {
   p <- m/sum(m)
   C%*%log(Mr%*%p)
}

X <- scan()
1 1 0
1 0 1
1 0 0

X <- matrix(X,3,3,byrow=T)

h.fct <- function(m) {
   L.fct(m)
}
a <-
mph.fit(y=y,Z=Z,ZF=ZF,h.fct=h.fct,step=0.1,maxiter=1000)

  iter= 997  norm.diff= 0.3193339  norm.score= 3.279595e-09 
  iter= 998  norm.diff= 0.3193339  norm.score= 3.922673e-09 
  iter= 999  norm.diff= 0.3193339  norm.score= 3.517499e-09 
  iter= 1000  norm.diff= 0.3193339  norm.score= 2.772561e-09 

 Time Elapsed: 7.36 seconds

mph.summary(a,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 3 ):  Gsq =  33.25029 (p =  2.852e-07 )
    Pearson's Score Stat  (df= 3 ):  Xsq =  38.11701 (p =  2.6698e-08 )
    Generalized Wald Stat (df= 3 ):  Wsq =  33.77626 (p =  2.2088e-07 )


CELL-SPECIFIC STATISTICS...
    OBS           FV   StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1  201 1.826391e+02 1.069712e+01 3.176332e-01 1.860368e-02    5.749431e+00
y2   28 3.712722e+01 5.169310e+00 6.456908e-02 8.990104e-03   -3.225310e+00
y3   37 3.531513e+01 5.745554e+00 6.141762e-02 9.992268e-03    4.589575e+00
y4   21 3.395460e+01 5.143368e+00 5.905148e-02 8.944987e-03   -5.526221e+00
y5    8 5.326702e+00 1.879169e+00 9.263829e-03 3.268120e-03    2.023091e+00
y6    7 9.337218e+00 2.987685e+00 1.623864e-02 5.195974e-03   -4.589575e+00
y7   12 1.174490e+01 3.389471e+00 2.042592e-02 5.894732e-03    1.986529e+00
y8    0 4.556436e-41 6.750138e-21 7.924237e-44 1.173937e-23   -4.556436e-41
y9    0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y10  27 3.866533e+01 5.286542e+00 6.724404e-02 9.193986e-03   -4.094285e+00
y11   9 7.836827e+00 1.643760e+00 1.362926e-02 2.858713e-03    5.187327e-01
y12   5 6.490560e+00 2.512333e+00 1.128793e-02 4.369275e-03   -4.589575e+00
y13  14 7.183925e+00 1.831879e+00 1.249378e-02 3.185877e-03    3.525268e+00
y14   4 1.102329e+00 9.212980e-01 1.917094e-03 1.602257e-03    5.778889e+00
y15   4 1.814156e+00 1.257619e+00 3.155055e-03 2.187164e-03    4.589575e+00
y16   2 2.230031e+00 1.485925e+00 3.878315e-03 2.584218e-03   -1.986529e+00
y17   0 3.526641e-70 1.877935e-35 6.133289e-73 3.265974e-38   -3.526641e-70
y18   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y19 142 1.377167e+02 1.016841e+01 2.395074e-01 1.768419e-02    3.705642e+00
y20  15 1.821994e+01 4.109449e+00 3.168686e-02 7.146868e-03   -3.705642e+00
y21   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y22  27 3.187476e+01 5.327032e+00 5.543437e-02 9.264404e-03   -3.705642e+00
y23  12 6.420549e+00 2.020355e+00 1.116617e-02 3.513660e-03    3.705642e+00
y24   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y25   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y26   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46
y27   0 3.720076e-46 1.928750e-23 6.469697e-49 3.354348e-26   -3.720076e-46


CONVERGENCE STATISTICS...
    iterations = 1000
    norm.diff  = 0.319334
    norm.score = 2.77256e-09
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

EXAMPLE 13. Kappa regression models.

Analysis of Multiple Sclerosis Data. See  Landis and Koch (Biometrics 1977) for data source and a WLS analysis.  See Lang (HLP 2002) for a more detailed description of the models and the corresponding ML analyses outlined below.

Model Systematic Components:  

DIRECT SMOOTHING MODEL:  K(m) = Xb,  where K(m) is the vector of eight weighted kappa statistics.  The design matrix X forces certain kappa coefficients to be equal.

INDIRECT SMOOTHING MODEL:  L(m) = [L1(m), L2(m)],  where  L1(m) = log(m),  L2(m) = K(m), the eight kappa coefficients.   The model constrains the data parameters through

L1(m) = Xb,  L2(m)  = κ.

The model L1(m) =  Xb has (k,i,j) component of the form

log(mkij) =    bk + b(A)ki + b(B)kj + b(A*B)k*i*j + b(agree)k*I(i=j),  where k indexes the two populations, i and j are the responses by neurologist A and B, resp. 

The models have the generic form L(m) = Xb.  They can be seen to be HLP models.

IMPORTANT:  mph.fit requires the HLP link function to be defined in terms of the expected counts m.  Recall  p(m) = D-1(ZZTm)m  [see primary references]. 

#CREATE Weighted Kappa function...
kappa.weighted.fct <- function(w,p,nr) {
#
#  Computes Weighted Kappa Statistic
#  See Landis and Koch 1977.
#  Author:  Joseph B. Lang
#  Created: June 7, 2002 
#
#  Input: w = vectorized version of two-way 
#                  table of weights, w = (w11,w12,...,wRR)
#           p = vectorized version of two-way
#                 probabilities, p = (p11,p12,...,pRR)
#          nr = number of rows (and columns) in the table.
#  Output:    Weighted kappa.
#
  p <- matrix(p,nr,nr,byrow=T)
  prow <- apply(p,1,sum)
  pcol <- apply(p,2,sum)
  pind <- prow%*%t(pcol)
  pind <- matrix(t(pind),nr*nr,1)
  p <- as.matrix(c(t(p)))
  w <- as.matrix(w)
  obsa <- t(w)%*%p
  expa <- t(w)%*%pind
  c((obsa - expa)/(1-expa))
}

w1 <- scan()
1 0 0 0
0 1 0 0 
0 0 1 0 
0 0 0 1

w2 <- scan()
1 1 0 0
1 1 0 0
0 0 1 0
0 0 0 1

w3 <- scan()
1 1 0 0
1 1 0 0
0 0 1 1
0 0 1 1

w4 <- scan()
1 1 0 0
1 1 1 0
0 1 1 1
0 0 1 1

y <- scan()
38  5  0  1 
33 11  3  0 
10 14  5  6  
3  7  3 10   
5 3 0 0  
3 11 4 0 
2 13 3 4 
1 2 4 14


Z <- kronecker(diag(2),matrix(1,16,1))
ZF <- Z
# Fit the DIRECT SMOOTHING MODEL 
L.fct <- function(m) {
   m1 <- m[1:16]
   m2 <- m[17:32]
   p1 <- m1/sum(m1)
   p2 <- m2/sum(m2)
   rbind(
   kappa.weighted.fct(w1,p1 ,4),
   kappa.weighted.fct(w2,p1 ,4),
   kappa.weighted.fct(w3,p1 ,4),
   kappa.weighted.fct(w4,p1 ,4),
   kappa.weighted.fct(w1,p2 ,4),
   kappa.weighted.fct(w2,p2 ,4),
   kappa.weighted.fct(w3,p2 ,4),
   kappa.weighted.fct(w4,p2 ,4))
}
   
XK <- scan()
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0 
0 0 0 1 0 
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 0 1

XK <- matrix(XK,8,5,byrow=T)
a <- mph.fit(y,Z,ZF,L.fct=L.fct,X=XK,norm.diff.conv=10,norm.score.conv=1e-8)
 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 90.69429  norm.score= 0.4710426 
  iter= 2  norm.diff= 6.863225  norm.score= 0.0586231 
  iter= 3  norm.diff= 2.437451  norm.score= 0.005285491 
  iter= 4  norm.diff= 2.260665  norm.score= 0.001090416 
  iter= 5  norm.diff= 2.260652  norm.score= 0.0002428371 
  iter= 6  norm.diff= 2.260610  norm.score= 7.249635e-05 
  iter= 7  norm.diff= 2.260611  norm.score= 2.431151e-05 
  iter= 8  norm.diff= 2.260611  norm.score= 8.772316e-06 
  iter= 9  norm.diff= 2.260611  norm.score= 3.26839e-06 
  iter= 10  norm.diff= 2.260611  norm.score= 1.237528e-06 
  iter= 11  norm.diff= 2.260611  norm.score= 4.727712e-07 
  iter= 12  norm.diff= 2.260611  norm.score= 1.819133e-07 
  iter= 13  norm.diff= 2.260611  norm.score= 7.02036e-08 
  iter= 14  norm.diff= 2.260611  norm.score= 2.741952e-08 
  iter= 15  norm.diff= 2.260611  norm.score= 1.096503e-08 
  iter= 16  norm.diff= 2.260611  norm.score= 5.495152e-09 

 Time Elapsed: 13.28 seconds
 
> mph.summary(a,T)
OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 3 ):  Gsq =  2.16255 (p =  0.53936 )
    Pearson's Score Stat  (df= 3 ):  Xsq =  2.11354 (p =  0.54918 )
    Generalized Wald Stat (df= 3 ):  Wsq =  2.26792 (p =  0.51869 )


LINEAR PREDICTOR MODEL RESULTS...
           BETA StdErr(BETA)  Z-ratio      p-value
beta1 0.2349553   0.04241687 5.539195 3.038648e-08
beta2 0.3150272   0.04935842 6.382442 1.742859e-10
beta3 0.3888738   0.05743424 6.770766 1.281020e-11
beta4 0.5831200   0.06909629 8.439236 0.000000e+00
beta5 0.7934268   0.08080265 9.819316 0.000000e+00

       OBS LINK   ML LINK  StdErr(L)  LINK RESID
link1 0.2079425 0.2349553 0.04241687 -0.97835070
link2 0.3275399 0.3150272 0.04935842  0.31584009
link3 0.4081121 0.3888738 0.05743424  0.44185351
link4 0.5964663 0.5831200 0.06909629  0.41284558
link5 0.2965166 0.2349553 0.04241687  0.94238307
link6 0.3324734 0.3150272 0.04935842  0.26087642
link7 0.3864188 0.3888738 0.05743424 -0.02982431
link8 0.7893773 0.7934268 0.08080265 -0.12157636


CELL-SPECIFIC STATISTICS...
    OBS           FV   StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1   38 4.019436e+01 5.156887e+00 2.697608e-01 3.460998e-02   -1.321327e+00
y2    5 4.273782e+00 1.971580e+00 2.868310e-02 1.323208e-02    1.413213e+00
y3    0 1.545839e-09 3.931716e-05 1.037476e-11 2.638736e-07   -1.545839e-09
y4    1 1.022546e+00 1.002682e+00 6.862723e-03 6.729409e-03   -2.237055e-01
y5   33 2.950032e+01 4.184771e+00 1.979887e-01 2.808571e-02    1.411519e+00
y6   11 1.231683e+01 3.222158e+00 8.266331e-02 2.162522e-02   -1.375604e+00
y7    3 3.228929e+00 1.741882e+00 2.167067e-02 1.169048e-02   -6.480198e-01
y8    0 3.745687e-09 6.120202e-05 2.513884e-11 4.107518e-07   -3.745687e-09
y9   10 1.032682e+01 3.018681e+00 6.930754e-02 2.025960e-02   -4.628183e-01
y10  14 1.446523e+01 3.475636e+00 9.708206e-02 2.332642e-02   -4.697435e-01
y11   5 4.904264e+00 2.097584e+00 3.291453e-02 1.407775e-02    1.634700e-01
y12   6 5.758698e+00 2.192534e+00 3.864898e-02 1.471499e-02    2.826312e-01
y13   3 3.093796e+00 1.727300e+00 2.076373e-02 1.159261e-02   -4.373569e-01
y14   7 7.222684e+00 2.573183e+00 4.847439e-02 1.726968e-02   -4.442144e-01
y15   3 2.867554e+00 1.608975e+00 1.924533e-02 1.079849e-02    2.801152e-01
y16  10 9.824181e+00 2.769387e+00 6.593410e-02 1.858649e-02    1.432252e-01
y17   5 3.996490e+00 1.767690e+00 5.792015e-02 2.561869e-02    1.254107e+00
y18   3 4.124129e+00 1.792444e+00 5.976999e-02 2.597744e-02   -1.378728e+00
y19   0 3.378332e-10 1.838024e-05 4.896133e-12 2.663803e-07   -3.378332e-10
y20   0 3.220382e-10 1.794542e-05 4.667220e-12 2.600786e-07   -3.220382e-10
y21   3 4.361039e+00 1.769190e+00 6.320346e-02 2.564044e-02   -1.392465e+00
y22  11 9.959295e+00 2.708989e+00 1.443376e-01 3.926071e-02    9.567614e-01
y23   4 4.072041e+00 1.797114e+00 5.901508e-02 2.604512e-02   -9.284069e-02
y24   0 1.427984e-09 3.778868e-05 2.069543e-11 5.476620e-07   -1.427984e-09
y25   2 1.892490e+00 1.332577e+00 2.742740e-02 1.931271e-02    4.222613e-01
y26  13 1.295768e+01 2.739007e+00 1.877925e-01 3.969575e-02    2.434344e-02
y27   3 2.785997e+00 1.582277e+00 4.037677e-02 2.293155e-02    5.191749e-01
y28   4 4.336488e+00 1.748584e+00 6.284766e-02 2.534180e-02   -3.354161e-01
y29   1 9.580754e-01 9.635969e-01 1.388515e-02 1.396517e-02    3.288490e-01
y30   2 2.019760e+00 1.373853e+00 2.927189e-02 1.991091e-02   -7.305229e-02
y31   4 4.413645e+00 1.684687e+00 6.396586e-02 2.441575e-02   -3.637495e-01
y32  14 1.312287e+01 2.756657e+00 1.901865e-01 3.995155e-02    5.040715e-01


CONVERGENCE STATISTICS...
    iterations = 16
    norm.diff  = 2.26061
    norm.score = 5.49515e-09
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02
# Fit the INDIRECT SMOOTHING MODEL 

 
L.fct <- function(m) {
   m1 <- m[1:16]
   m2 <- m[17:32]
   p1 <- m1/sum(m1)
   p2 <- m2/sum(m2)
   rbind(log(m),
   kappa.weighted.fct(w1,p1 ,4),
   kappa.weighted.fct(w2,p1 ,4),
   kappa.weighted.fct(w3,p1 ,4),
   kappa.weighted.fct(w4,p1 ,4),
   kappa.weighted.fct(w1,p2 ,4),
   kappa.weighted.fct(w2,p2 ,4),
   kappa.weighted.fct(w3,p2 ,4),
   kappa.weighted.fct(w4,p2 ,4))
}
  
A <- factor(gl(4,4))
B <- factor(gl(4,1,16))
Ascore <- c(A)
Bscore <- c(B)

agree <- scan()
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1


XA1 <-  model.matrix(~A + B + Ascore*Bscore - Ascore-Bscore + agree)
XA2 <- XA1

X <-
block.fct(XA1,XA2,diag(8))
b <-
mph.fit(y,Z,ZF,L.fct=L.fct,X=X)

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 11.22854  norm.score= 3.157116 
  iter= 2  norm.diff= 2.684903  norm.score= 0.4388314 
  iter= 3  norm.diff= 0.7834525  norm.score= 0.02531604 
  iter= 4  norm.diff= 0.04914612  norm.score= 9.326271e-05 
  iter= 5  norm.diff= 0.0001894680  norm.score= 2.261607e-09 
  iter= 6  norm.diff= 3.73313e-09  norm.score= 1.116094e-09 

 Time Elapsed: 6.08 seconds


> mph.summary(b,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 14 ):  Gsq =  18.25286 (p =  0.19551 )
    Pearson's Score Stat  (df= 14 ):  Xsq =  22.48127 (p =  0.069252 )
    Generalized Wald Stat (df= 14 ):  Wsq =  13.13214 (p =  0.51615 )


LINEAR PREDICTOR MODEL RESULTS...
              BETA StdErr(BETA)     Z-ratio      p-value
beta1   2.82198935   0.21789484 12.95115269 0.000000e+00
beta2  -0.98252487   0.34280421 -2.86614005 4.155104e-03
beta3  -2.63574417   0.59030845 -4.46502870 8.005802e-06
beta4  -5.07053731   1.00283130 -5.05622162 4.276444e-07
beta5  -2.50209975   0.36908649 -6.77916906 1.208700e-11
beta6  -5.95129014   0.85057662 -6.99677137 2.619238e-12
beta7  -8.20966143   1.38875396 -5.91153053 3.389434e-09
beta8  -0.02783032   0.24285556 -0.11459617 9.087652e-01
beta9   0.80384748   0.15516362  5.18064414 2.211210e-07
beta10  0.40753569   0.40135119  1.01540919 3.099108e-01
beta11 -0.93504420   0.59981172 -1.55889618 1.190210e-01
beta12 -3.00645110   1.21010157 -2.48446177 1.297474e-02
beta13 -6.18631028   2.08586141 -2.96582997 3.018673e-03
beta14 -1.26080589   0.64291838 -1.96106681 4.987123e-02
beta15 -5.23293410   1.44644377 -3.61779298 2.971259e-04
beta16 -8.36055037   2.47380527 -3.37963156 7.258306e-04
beta17  0.02769186   0.34868656  0.07941764 9.367004e-01
beta18  1.04115569   0.29708388  3.50458495 4.573196e-04
beta19  0.20794246   0.05113013  4.06692644 4.763727e-05
beta20  0.33160222   0.05237560  6.33123511 2.432063e-10
beta21  0.41743362   0.06030132  6.92246221 4.438672e-12
beta22  0.55622648   0.07116497  7.81601533 5.551115e-15
beta23  0.29651657   0.07818321  3.79258641 1.490863e-04
beta24  0.38035962   0.07129014  5.33537446 9.534758e-08
beta25  0.47547451   0.07625832  6.23505065 4.516318e-10
beta26  0.73473413   0.09436958  7.78570910 6.883383e-15

        OBS LINK     ML LINK  StdErr(L)    LINK RESID
link1  3.6375862  3.59800651 0.13749315  9.427814e-01
link2  1.6094379  1.92758456 0.28963410 -1.357818e+00
link3       -Inf -0.71775836 0.54095671          -Inf
link4  0.0000000 -2.17228217 0.87175123  7.674599e-01
link5  3.4965076  3.44715943 0.15132358  1.046012e+00
link6  2.3978953  2.52492432 0.24826770 -1.173576e+00
link7  1.0986123  0.71125919 0.37191660  6.585275e-01
link8       -Inf  0.06058286 0.51585603          -Inf
link9  2.3025851  2.59778761 0.20273772 -1.809184e+00
link10 2.6390573  2.50723029 0.21627385  7.877068e-01
link11 1.6094379  1.44175201 0.36104443  5.317376e-01
link12 1.7917595  1.62275346 0.31556472  5.600586e-01
link13 1.0986123  0.96684194 0.44794027  3.168800e-01
link14 1.9459101  1.68013209 0.27927295  8.336371e-01
link15 1.0986123  1.44633160 0.36606789 -1.129845e+00
link16 2.3025851  2.37551990 0.26455789 -5.719039e-01
link17 1.6094379  1.47638325 0.40026306  5.738569e-01
link18 1.0986123  1.22904119 0.42768254 -4.227939e-01
link19      -Inf -1.70193133 0.76941354          -Inf
link20      -Inf -3.78839190 1.42386504          -Inf
link21 1.0986123  1.55480288 0.38587064 -2.085649e+00
link22 2.3978953  2.40400024 0.25548478 -5.932316e-02
link23 1.3862944  0.48649156 0.45024654  1.427064e+00
link24      -Inf -0.55881333 0.71595805          -Inf
link25 0.6931472  0.52455168 0.48406780  2.878691e-01
link26 2.5649494  2.38721287 0.24432611  1.336046e+00
link27 1.0986123  1.56624361 0.38466892 -2.171738e+00
link28 1.3862944  1.53440256 0.35137338 -5.315800e-01
link29 0.0000000 -1.61415181 1.00711518  8.075958e-01
link30 0.6931472  1.28966509 0.40136378 -1.888459e+00
link31 1.3862944  1.48215965 0.38094578 -3.688923e-01
link32 2.6390573  2.54685802 0.23696240  1.051840e+00
link33 0.2079425  0.20794246 0.05113013 -5.652506e-12
link34 0.3275399  0.33160222 0.05237560 -1.162164e-01
link35 0.4081121  0.41743362 0.06030132 -2.387584e-01
link36 0.5964663  0.55622648 0.07116497  1.160200e+00
link37 0.2965166  0.29651657 0.07818321  6.381951e-12
link38 0.3324734  0.38035962 0.07129014 -1.198556e+00
link39 0.3864188  0.47547451 0.07625832 -1.556305e+00
link40 0.7893773  0.73473413 0.09436958  1.884445e+00


CELL-SPECIFIC STATISTICS...
    OBS          FV StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1   38 36.52534895 5.02198531 0.2451365702 0.0337045994      0.96168750
y2    5  6.87288910 1.99062303 0.0461267725 0.0133598861     -1.16301937
y3    0  0.48784460 0.26390281 0.0032741248 0.0017711598     -0.75582368
y4    1  0.11391734 0.09930758 0.0007645459 0.0006664938      2.74804519
y5   33 31.41104031 4.75323093 0.2108123511 0.0319008787      1.07225144
y6   11 12.48995000 3.10085112 0.0838251678 0.0208110814     -1.10209548
y7    3  2.03655406 0.75742826 0.0136681481 0.0050834111      0.80426336
y8    0  1.06245563 0.54807414 0.0071305747 0.0036783500     -1.22315861
y9   10 13.43398393 2.72357521 0.0901609660 0.0182790282     -1.56659247
y10  14 12.27089612 2.65387394 0.0823550075 0.0178112345      0.84198605
y11   5  4.22809698 1.52653085 0.0283764898 0.0102451735      0.57892006
y12   6  5.06702297 1.59897370 0.0340068656 0.0107313671      0.61016795
y13   3  2.62962681 1.17791573 0.0176485021 0.0079054747      0.33870574
y14   7  5.36626477 1.49865259 0.0360151998 0.0100580711      0.95492110
y15   3  4.24750436 1.55487498 0.0285067406 0.0104354025     -0.95432934
y16  10 10.75660406 2.84574444 0.0721919736 0.0190989560     -0.55154593
y17   5  4.37708618 1.75198589 0.0634360317 0.0253910999      0.61378516
y18   3  3.41795081 1.46179790 0.0495355189 0.0211854768     -0.39638227
y19   0  0.18233104 0.14028797 0.0026424789 0.0020331590     -0.45276968
y20   0  0.02263197 0.03222487 0.0003279995 0.0004670270     -0.15404060
y21   3  4.73415323 1.82677076 0.0686109164 0.0264749385     -1.67471149
y22  11 11.06736010 2.82754207 0.1603965232 0.0409788706     -0.05914244
y23   4  1.62659937 0.73237074 0.0235739039 0.0106140687      2.31412319
y24   0  0.57188731 0.40944732 0.0082882218 0.0059340191     -0.90479860
y25   2  1.68970115 0.81792992 0.0244884224 0.0118540568      0.31355899
y26  13 10.88311901 2.65903009 0.1577263624 0.0385366679      1.46213690
y27   3  4.78862642 1.84203576 0.0694003829 0.0266961704     -1.73465268
y28   4  4.63855343 1.62986422 0.0672254120 0.0236212205     -0.49408792
y29   1  0.19905944 0.20047578 0.0028849194 0.0029054461      2.01310856
y30   2  3.63157009 1.45758069 0.0526314506 0.0211243579     -1.42231343
y31   4  4.40244317 1.67709214 0.0638035242 0.0243056832     -0.35176211
y32  14 12.76692730 3.02528177 0.1850279319 0.0438446634      1.10185439


CONVERGENCE STATISTICS...
    iterations = 6
    norm.diff  = 3.73313e-09
    norm.score = 1.11609e-09
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02
Note: The last eight rows in the link estimate matrix of the LINEAR PREDICTOR MODEL RESULTS section correspond to the estimates of the eight kappa coefficients.


EXAMPLE 14.  Marginal cumulative probit model (in conjunction with loglinear association model)

Marijuana use among U.S. male teens.  Data Source: U.S. National Youth Survey (Elliott et al. 1989). Table Source: Lang, J.B. (2002), "Homogeneous Linear Predictor Models for Contingency Tables," Univ of Iowa working paper.

    Use at Age 17
x1 x2 x3 Total
x1 56 13 7 76
Use at Age 15 x2 6 4 10 20
  x3 1 5 15 21
Total 63 22 32 117

The scores x1 < x2 < x3  correspond to the responses, "never used marijuana in the past year," "used no more than once per month in the past year," and "used more than once per month in the past year."

The marginal probit models that are fitted below have the following form...

L(m) = [L1(m), L2(m)] = [X1b, X2a],  where  L2(m) = log(m)  and

L1(m) = Φ-1[p2++ p3+, p3+, p+2 + p+3, p+3]; 

Here pij = mij/m++ and Φ(.) is the standard normal pdf.

The linear predictor for the marginal probit model is defined as

X1b  = [b(cut1), b(cut2), b(cut1)+b(AGE), b(cut2)+b(AGE)],

and the loglinear association model L2(m) = log(m) = X2a has (i,j) component

log(mij) =  a + a(A)i + a(B)j + a(A*B)ij   [saturated association model]  or

log(mij) = a + a(A)i + a(B)j   [independence association model].

The models, which have the generic form L(m) = Xb, can be seen to be HLP models.

IMPORTANT:  mph.fit requires the HLP link function to be defined in terms of the expected counts m.   Recall  p(m) = D-1(ZZTm)m  [see primary references]. 

y <- scan()
56 13  7
 6  4 10
 1  5 15

Z <- ZF <- matrix(1,9,1)

Lprobit.fct <- function(m) {
  p <- matrix(m,3,3,byrow=T)/sum(m)
  pr <- apply(p,1,sum)
  pc <- apply(p,2,sum)
  rbind(
    qnorm(pr[2]+pr[3]),
    qnorm(pr[3]),
    qnorm(pc[2]+pc[3]),
    qnorm(pc[3]),
    log(m)
  )
}

Xprobit <- scan()
1 0 0
0 1 0
1 0 1
0 1 1

Xprobit <- matrix(Xprobit,4,3,byrow=T)

A <- factor(gl(3,3))
B <- factor(gl(3,1,9))
Xindep <- model.matrix(~A+B)
Xsat <- model.matrix(~A*B)
 

# FIT USING THE SATURATED ASSOCIATION MODEL

X <- block.fct(Xprobit,Xsat)
a <-
mph.fit(y,Z,ZF,L.fct=Lprobit.fct,X=X )  

 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 1.356676  norm.score= 0.01463237 
  iter= 2  norm.diff= 0.005404818  norm.score= 0.0002264007 
  iter= 3  norm.diff= 0.0005536982  norm.score= 1.521958e-06 
  iter= 4  norm.diff= 2.897933e-06  norm.score= 4.267208e-08 

 Time Elapsed: 0.17 seconds

mph.summary(a,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 1 ):  Gsq =  0.03427 (p =  0.85314 )
    Pearson's Score Stat  (df= 1 ):  Xsq =  0.03429 (p =  0.8531 )
    Generalized Wald Stat (df= 1 ):  Wsq =  0.03418 (p =  0.85332 )


LINEAR PREDICTOR MODEL RESULTS...
             BETA StdErr(BETA)   Z-ratio      p-value
beta1  -0.3897969   0.11517082 -3.384511 7.130513e-04
beta2  -0.9081314   0.12565068 -7.227429 4.922729e-13
beta3   0.2979799   0.09780294  3.046738 2.313392e-03
beta4   4.0243358   0.09642694 41.734560 0.000000e+00
beta5  -2.2620263   0.40646139 -5.565169 2.618985e-08
beta6  -4.0137515   1.00172815 -4.006827 6.153987e-05
beta7  -1.4331903   0.26785693 -5.350581 8.767210e-08
beta8  -2.0847836   0.40100510 -5.198895 2.004763e-07
beta9   1.0541642   0.71778220  1.468641 1.419303e-01
beta10  3.0701622   1.12907074  2.719194 6.544125e-03
beta11  2.5904165   0.66092097  3.919404 8.876808e-05
beta12  4.7874296   1.10337541  4.338895 1.432012e-05

          OBS LINK     ML LINK  StdErr(L) LINK RESID
link1  -0.38416697 -0.38979694 0.11517082  0.1849613
link2  -0.91732119 -0.90813140 0.12565068 -0.1859397
link3  -0.09655862 -0.09181700 0.11318722 -0.1852044
link4  -0.60224878 -0.61015147 0.11644810  0.1847193
link5   4.02535169  4.02433585 0.09642694  0.1850694
link6   2.56494936  2.59114552 0.21653665 -0.1875993
link7   1.94591015  1.93955228 0.36610760  0.1845754
link8   1.79175947  1.76230956 0.37019818  0.1824503
link9   1.38629436  1.38328343 0.49187564  0.1848848
link10  2.30258509  2.26794253 0.24235739  0.1819747
link11  0.00000000  0.01058434 0.98878275 -0.1861451
link12  1.60943791  1.64755617 0.37838262 -0.1887149
link13  2.70805020  2.71323036 0.23873958 -0.1856434


CELL-SPECIFIC STATISTICS...
   OBS        FV StdErr(FV)        PROB StdErr(PROB) ADJUSTED RESIDS
y1  56 55.943142   5.394426 0.478146509   0.04610620       0.1851634
y2  13 13.345050   2.889692 0.114060255   0.02469823      -0.1851634
y3   7  6.955636   2.546511 0.059449881   0.02176505       0.1851634
y4   6  5.825877   2.156729 0.049793821   0.01843358       0.1851634
y5   4  3.987974   1.961587 0.034085251   0.01676570       0.1851634
y6  10  9.659506   2.341053 0.082559882   0.02000900       0.1851634
y7   1  1.010641   0.999304 0.008637953   0.00854106      -0.1851634
y8   5  5.194270   1.965422 0.044395473   0.01679848      -0.1851634
y9  15 15.077904   3.599692 0.128870974   0.03076660      -0.1851634


CONVERGENCE STATISTICS...
    iterations = 4
    norm.diff  = 2.89793e-06
    norm.score = 4.26721e-08
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02

# TEST WHETHER b(AGE) = 0 (i.e. Test Marginal Homogeneity) vs. b(AGE) > 0 using the Wald statistic Z = bhat(AGE)/ase(bhat(AGE)).
 

a$beta[3]/a$covbeta[3,3]**0.5

   beta3 
3.046738 
1-pnorm(3.046738)
[1] 0.001156696
 
# USE THE INDEPENDENCE ASSOCIATION MODEL
X <- block.fct(Xprobit,Xindep)
b <-
mph.fit(y,Z,ZF,L.fct=Lprobit.fct,X=X )
 mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 4.412758  norm.score= 21.29174 
  iter= 2  norm.diff= 7.327984  norm.score= 3.255189 
  iter= 3  norm.diff= 2.434905  norm.score= 0.2206088 
  iter= 4  norm.diff= 0.1555902  norm.score= 0.002838257 
  iter= 5  norm.diff= 0.003965627  norm.score= 1.855961e-05 
  iter= 6  norm.diff= 4.199287e-05  norm.score= 2.400974e-08 
  iter= 7  norm.diff= 7.596187e-08  norm.score= 3.81935e-09 

 Time Elapsed: 0.23 seconds
mph.summary(b)
OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 5 ):  Gsq =  49.31471 (p =  1.9137e-09 )
    Pearson's Score Stat  (df= 5 ):  Xsq =  45.3922 (p =  1.2075e-08 )
    Generalized Wald Stat (df= 5 ):  Wsq =  28.59303 (p =  2.7864e-05 )


LINEAR PREDICTOR MODEL RESULTS...
            BETA StdErr(BETA)   Z-ratio      p-value
beta1 -0.3890030    0.1162506 -3.346244 8.191435e-04
beta2 -0.9073211    0.1239864 -7.317909 2.517986e-13
beta3  0.2975494    0.1575182  1.888983 5.889404e-02
beta4  3.7106739    0.1092253 33.972660 0.000000e+00
beta5 -1.3639609    0.1994296 -6.839310 7.957635e-12
beta6 -1.2744095    0.2368939 -5.379665 7.462463e-08
beta7 -1.0245378    0.1973052 -5.192654 2.073169e-07
beta8 -0.6828006    0.2159201 -3.162283 1.565374e-03

          OBS LINK    ML LINK StdErr(L)  LINK RESID
link1  -0.38416697 -0.3890030 0.1162506   0.1863585
link2  -0.91732119 -0.9073211 0.1239864  -0.1873838
link3  -0.09655862 -0.0914536 0.1127734  -0.1865775
link4  -0.60224878 -0.6097718 0.1172767   0.1861074
link5   4.02535169  3.7106739 0.1092253   4.9855616
link6   2.56494936  2.6861361 0.1435634  -0.6137610
link7   1.94591015  3.0278733 0.1623681  -9.3092681
link8   1.79175947  2.3467130 0.1540198  -2.2037578
link9   1.38629436  1.3221752 0.2799146   0.1512751
link10  2.30258509  1.6639124 0.1702881   1.6389618
link11  0.00000000  2.4362643 0.2061525 -12.7622752
link12  1.60943791  1.4117266 0.1811976   0.4395295
link13  2.70805020  1.7534638 0.2461471   2.9595106


CONVERGENCE STATISTICS...
    iterations = 7
    norm.diff  = 7.59619e-08
    norm.score = 3.81935e-09
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 
Remark:  As expected, the independence/probit model appears to be untenable for these repeated measures data.  Although the independence-model estimator of the marginal probit parameter b(AGE) is consistent, the reported standard error is not.  In fact, it can be anticipated that the actual standard error of this independence-model estimator is significantly larger than the  saturated-model estimator which correctly accommodates the association between the two responses, AGE15 and AGE17 marijuana use responses.

EXAMPLE 15.   Marginal Cumulative Logit and Linear-by-Linear Loglinear Association Model (Generalized Loglinear and Indirect Smoothing Model Specification).

Opinions about Government Spending (data sources: General Social Survey 1989 and  Table 1 of Lang and Agresti, JASA94).

Cities   1     2     3  
Law Enforcement 1 2 3 1 2 3 1 2 3
Environment Health                  
1 1 62 17 5 90 42 3 74 31 11
  2 11 7 0 22 18 1 19 14 3
  3 2 3 1 2 0 1 1 3 1
2 1 11 3 0 21 13 2 20 8 3
  2 1 4 0 6 9 0 6 5 2
  3 1 0 1 2 1 1 4 3 1
3 1 3 0 0 2 1 0 9 2 1
  2 1 0 0 2 1 0 4 2 0
  3 1 0 0 0 0 0 1 2 3

When writing x = (xij, i=1,2,...I, j=1,2,...,J) etc.  use the convention that the rightmost subscripts move faster.   Thus, x = (x11, x12, ..., xIJ).

Write the vector of observed table counts as y = (yijkl,  i,j,k,l = 1,2,3) and the cell probabilities as p = (pijkl, i,j,k,l = 1,2,3).  Here, the subscripts correspond to the responses as follows:

i = Environment, j=Health, k=Cities, l=Law Enforcement.

Data Model

Data Sampling Model:  We model the  counts as

y ← Y ~ mult(n, p)  or, using the MP notation of Lang (2002),  Y ~ MPZ(m|ZF,n), 

where Z = ZF = 181 and m = np is the vector of expected table counts.

Write the vector of expected table counts as m = (mijkl,  i,j,k,l = 1,2,3).

Data Parameter Model:  Consider the generalized loglinear model  used in Lang and Agresti (JASA94):  L*(m) = [L1(m), L2(m)] = [X1a, X2b] = X*δ,   where the components are defined as follows:

L1(m) = log(m) = log((mijkl, i,j,k,l=1,2,3)).    

L2(m) = vector of oneway cumulative logits = log (cth, t=1,2,3,4, h=1,2)  and cth = odds(Rt <= h).  Here, Rt is the tth response variable (t=1=Environment, 2=Health, 3=Cities, 4=Law Enforcement).

X1a  has (ijkl)th element of the form: 

a(0) + a(1)i + a(2)j+ a(3)k + a(4)l + a(1,2)*ij + a(1,3)*ik + a(1,4)*il + a(2,3)*jk + a(2,4)*jl + a(3,4)*kl

X2b has (t,h)th element of the form:   b(level)h + b(response)t

Note that L*(m) = X*δ  has the generalized loglinear model form ClogAm = X*δ.  This model was denoted J(LxL)∩M(PO) in Lang and Agresti (JASA94), because the association structure is modeled via a  linear-by-linear loglinear model and the one-way marginal structure is modeled via a proportional odds cumulative logit model.

For convenience, we will augment this model so that estimates of the oneway marginals along with their corresponding adjusted residuals will be output.  Specifically,  we will reexpress the  model in the equivalent  form

L(m) = [L1(m), L2(m), L3(m)] = [X1a, X2b, θ] = Xβ,   where  L3(m) = (P(Rt=h), t=1,2,3,4, h=1,2,3) is the vector of oneway marginal probabilities and θ = (θth, t=1,2,3,4,  h=1,2,3).

In this model, L3(m) = θ, is an unrestricted model.  The marginal probabilities L3(m) are indirectly smoothed via the constraints implied by [L1(m), L2(m)] = [X1a, X2b].   See Lang (HLP2002) for a discussion of these "indirect smoothing models." 

Note that L(m) = Xβ and L*(m) = X*δ  are identical models.  However, the latter is of the generalized loglinear model form and the former is not.

The data model (data sampling model along with the data parameter model),  y ←Y ~ MPZ(m|ZF,n),  where Z = ZF = 181  and   L(m) = Xβ,  is a homogeneous linear predictor model  (see Lang HLP2002).  Therefore the program mph.fit can be used to easily fit it. 

IMPORTANT:  mph.fit requires the HLP link function L(m) to be defined in terms of the expected counts m.   Recall  p(m) = D-1(ZZTm)m  [see primary references]. 

y <- scan()
62 17 5 90 42 3 74 31 11
11 7 0 22 18 1 19 14 3
2 3 1 2 0 1 1 3 1
11 3 0 21 13 2 20 8 3
1 4 0 6 9 0 6 5 2
1 0 1 2 1 1 4 3 1
3 0 0 2 1 0 9 2 1
1 0 0 2 1 0 4 2 0
1 0 0 0 0 0 1 2 3

# Write y = (yijkl: i,j,k,l=1,2,3), where the order follows the rule that
# the rightmost subscripts move faster. The subscripts on yijkl,  correspond
# to the four responses as follows:
#        i=Environment, j=Health, k=Cities,  and l= Law Enforcement.
# For example, y1212 = 7 is the 11th observation in the vector y; it means
# that 7 of those sampled responded (Environ=1, Health=2, Cities=1, Law=2).
# See Table 1 of Lang and Agresti (JASA94).

Z <- ZF <- matrix(1,81,1)

## CREATE THE LINK FUNCTION L(m). ##

M1 <- Marg.fct(1,c(3,3,3,3))
M2 <-
Marg.fct(2,c(3,3,3,3))
M3 <-
Marg.fct(3,c(3,3,3,3))
M4 <-
Marg.fct(4,c(3,3,3,3))
M <- rbind(M1,M2,M3,M4)

CUM0 <- scan()
1 0 0
0 1 1
1 1 0
0 0 1

CUM0 <- matrix(CUM0,4,3,byrow=T)
CUM <- kronecker(diag(4),CUM0)
AMarg <- CUM%*%M

C0 <- scan()
1 -1 0 0
0 0 1 -1

C0 <- matrix(C0,2,4,byrow=T)
CMarg <- kronecker(diag(4),C0)


L.fct <- function(m) {
   logm <- log(m)
   onewayclogit <- CMarg%*%log(AMarg%*%m)
   onewayprob <- M%*%m/sum(m)
   rbind(logm,onewayclogit,onewayprob)
}

##CREATE DESIGN MATRIX  X = block[XAssoc, XMarg, I_12] ##

Environ <- factor(gl(3,27,81))
Health  <- factor(gl(3,9,81))
Cities <- factor(gl(3,3,81))
Law <- factor(gl(3,1,81))
Escore <- c(Environ)
Hscore <- c(Health)
Cscore <- c(Cities)
Lscore <- c(Law)

XAssoc <- model.matrix(~Environ+Health+Cities+Law+Escore*Hscore + Escore*Cscore
                  + Escore*Lscore + Hscore*Cscore + Hscore*Lscore + Cscore*Lscore
                  - Escore - Hscore - Cscore - Lscore)

XMarg <- scan()
1 0 0 0 0
0 1 0 0 0
1 0 1 0 0
0 1 1 0 0
1 0 0 1 0
0 1 0 1 0
1 0 0 0 1
0 1 0 0 1

XMarg <- matrix(XMarg, 8,5,byrow=T)

X <-
block.fct(XAssoc, XMarg, diag(12))

## FIT THE MODEL AND SUMMARIZE ##

a <- mph.fit(y,Z,ZF,L.fct=L.fct,X=X)

mph.fit, version 1.0, 6/5/02 , running...
  iter= 1  norm.diff= 31.76308  norm.score= 8.722637 
  iter= 2  norm.diff= 4.729462  norm.score= 1.364037 
  iter= 3  norm.diff= 7.65191  norm.score= 0.3431408 
  iter= 4  norm.diff= 0.5504504  norm.score= 0.04565881 
  iter= 5  norm.diff= 0.04405597  norm.score= 0.004250447 
  iter= 6  norm.diff= 0.00172583  norm.score= 0.0003981413 
  iter= 7  norm.diff= 0.0001241868  norm.score= 4.046984e-05 
  iter= 8  norm.diff= 1.248086e-05  norm.score= 4.241045e-06 
  iter= 9  norm.diff= 1.396062e-06  norm.score= 4.561802e-07 

 Time Elapsed: 1.62 seconds

mph.summary(a,T)

OVERALL GOODNESS OF FIT: TEST of   Ho: h(m)=0 vs. Ha: not Ho...
    Likelihood Ratio Stat (df= 69 ):  Gsq =  71.535 (p =  0.39365 )
    Pearson's Score Stat  (df= 69 ):  Xsq =  64.33728 (p =  0.6365 )
    Generalized Wald Stat (df= 69 ):  Wsq =  43.12471 (p =  0.99382 )


LINEAR PREDICTOR MODEL RESULTS...
               BETA StdErr(BETA)      Z-ratio      p-value
beta1   2.552042346  0.178483120  14.29850818 0.000000e+00
beta2  -2.588205009  0.308823989  -8.38084184 0.000000e+00
beta3  -5.493086727  0.655267977  -8.38296227 0.000000e+00
beta4  -2.595621287  0.297230071  -8.73270082 0.000000e+00
beta5  -5.610605985  0.637689942  -8.79832912 0.000000e+00
beta6  -0.034223012  0.204332191  -0.16748713 8.669868e-01
beta7  -0.937638791  0.413278261  -2.26878323 2.328151e-02
beta8  -1.826510954  0.263654376  -6.92767168 4.278133e-12
beta9  -4.197613889  0.551070672  -7.61719703 2.597922e-14
beta10  0.498689174  0.112274534   4.44169447 8.925323e-06     
beta11  0.314319281  0.103927942   3.02439627 2.491299e-03     
beta12 -0.003494674  0.112011470  -0.03119925 9.751106e-01      
beta13  0.051558694  0.099563584   0.51784691 6.045651e-01     
beta14  0.454699693  0.103442225   4.39568749 1.104227e-05    
beta15  0.198569484  0.089578780   2.21670225 2.664344e-02 Log(m) Model Coef

beta16  0.994717362  0.090538025  10.98673576 0.000000e+00
beta17  2.852299490  0.112086258  25.44736117 0.000000e+00
beta18 -0.080955234  0.114848443  -0.70488752 4.808802e-01
beta19 -2.336871532  0.117010601 -19.97145140 0.000000e+00
beta20 -0.462380024  0.119923508  -3.85562456 1.154345e-04 Oneway C-Logit Model Coef

beta21  0.730018682  0.017844267  40.91054531 0.000000e+00
beta22  0.215418743  0.013701431  15.72235329 0.000000e+00
beta23  0.054562575  0.005782026   9.43658464 0.000000e+00
beta24  0.713769393  0.018116177  39.39955875 0.000000e+00
beta25  0.227338142  0.013759216  16.52260873 0.000000e+00
beta26  0.058892465  0.006166318   9.55066933 0.000000e+00
beta27  0.207156029  0.014156292  14.63349476 0.000000e+00
beta28  0.418922021  0.014035304  29.84773476 0.000000e+00
beta29  0.373921950  0.018727908  19.96602865 0.000000e+00
beta30  0.630028094  0.019252284  32.72484997 0.000000e+00
beta31  0.286027282  0.014261072  20.05650644 0.000000e+00
beta32  0.083944625  0.007931189  10.58411540 0.000000e+00 Oneway Prob Coef

           OBS LINK     ML LINK   StdErr(L)   LINK RESID
link1    4.12713439  4.06638400 0.095988109  0.766717419
link2    2.83321334  2.88964755 0.121653595 -0.285222242
link3    1.60943791  1.16831911 0.263546568  0.900833881
link4    4.49980967  4.59660845 0.055480543 -1.322115281
link5    3.73766962  3.61844148 0.083716619  0.884542801
link6    1.09861229  2.09568253 0.146242603 -3.153771089
link7    4.30406509  4.25764013 0.091828124  0.727183836
link8    3.43398720  3.47804264 0.099761757 -0.317381907
link9    2.39789527  2.15385318 0.194390397  0.881753508
link10   2.39789527  2.47571027 0.147114161 -0.315547214
link11   1.94591015  1.75367351 0.151960064  0.499025735
link12         -Inf  0.48704477 0.258577659         -Inf
link13   3.09104245  3.05749341 0.103701811  0.180350700
link14   2.89037176  2.53402614 0.116467519  1.407186242
link15   0.00000000  1.46596688 0.143423648 -3.209435230
link16   2.94443898  2.77008379 0.136629080  0.847325365
link17   2.63905733  2.44518600 0.125889659  0.736911795
link18   1.09861229  1.57569623 0.183270760 -1.151598003
link19   0.69314718  0.46567314 0.314590662  0.313318682
link20   1.09861229  0.19833607 0.279577645  1.046349721
link21   0.00000000 -0.61359298 0.399062615  0.472532093
link22   0.69314718  1.09901497 0.215228332 -0.759955507
link23         -Inf  1.03024739 0.161798207         -Inf
link24   0.00000000  0.41688783 0.265133656 -0.544056636
link25   0.00000000  0.86316404 0.282673601 -1.479719122
link26   1.09861229  0.99296595 0.217235400  0.186282078
link27   0.00000000  0.57817587 0.309194307 -0.849091506
link28   2.39789527  2.28769277 0.159071976  0.403617920
link29   1.09861229  1.10746164 0.189103423 -0.016348805
link30         -Inf -0.61736146 0.330757520         -Inf
link31   3.04452244  3.13223650 0.101436265 -0.492780525
link32   2.56494936  2.15057486 0.130274768  1.325035371
link33   0.69314718  0.62432124 0.207945843  0.098249025
link34   2.99573227  3.10758746 0.120324284 -0.661606694
link35   2.07944154  2.32449530 0.132070516 -0.873285229
link36   1.09861229  0.99681117 0.225586008  0.180947843
link37   0.00000000  1.19570822 0.184460807 -2.314829149
link38   1.38629436  0.47017679 0.191118931  1.196014784
link39         -Inf -0.79994663 0.292969106         -Inf
link40   1.79175947  2.09181064 0.132339573 -0.929070437
link41   2.19722458  1.56484869 0.145148254  1.464695897
link42         -Inf  0.49329476 0.168305244         -Inf
link43   1.79175947  2.11872030 0.142729601 -1.043553823
link44   1.60943791  1.79032783 0.137112068 -0.472669313
link45   0.69314718  0.91734339 0.179786308 -0.370783213
link46   0.00000000 -0.31563974 0.313328219  0.279940179
link47         -Inf -0.58647148 0.280566244         -Inf
link48   0.00000000 -1.40189520 0.394972610  0.709400648
link49   0.69314718  0.63202137 0.193018063  0.087090687
link50   0.00000000  0.55975912 0.141283077 -0.754958660
link51   0.00000000 -0.05709512 0.238341326  0.057086486
link52   1.38629436  0.71048972 0.236612344  1.026103316
link53   1.09861229  0.83679695 0.168630327  0.412416610
link54   0.00000000  0.41851220 0.258680888 -0.545105132
link55   1.09861229  0.19232483 0.340140679  1.077314961
link56         -Inf -0.99140097 0.370492122         -Inf
link57         -Inf -2.71971875 0.552626047         -Inf
link58   0.69314718  1.35118784 0.193879184 -1.403917966
link59   0.00000000  0.36603153 0.221722073 -0.456586987
link60         -Inf -1.16371677 0.398403024         -Inf
link61   2.19722458  1.64085808 0.227096973  1.483811861
link62   0.69314718  0.85427125 0.232411843 -0.264911221
link63   0.00000000 -0.47690756 0.400463523  0.396177783
link64   0.00000000 -0.40097055 0.327427053  0.340783809
link65         -Inf -1.12999665 0.326104251         -Inf
link66         -Inf -2.40361474 0.475427279         -Inf
link67   0.69314718  0.80945116 0.179375461 -0.181353587
link68   0.00000000  0.27899454 0.167189288 -0.327223978
link69         -Inf -0.79605407 0.310174094         -Inf
link70   1.38629436  1.15068010 0.205559438  0.451339032
link71   0.69314718  0.81879296 0.171184372 -0.196222851
link72         -Inf -0.05768616 0.308626546         -Inf
link73   0.00000000 -1.41362933 0.468598686  0.716773918
link74         -Inf -1.68795575 0.431739287         -Inf
link75         -Inf -2.50687414 0.560422554         -Inf
link76         -Inf -0.15164894 0.314340323         -Inf
link77         -Inf -0.22740586 0.255312335         -Inf
link78         -Inf -0.84775478 0.380583675         -Inf
link79   0.00000000  0.24113870 0.340998266 -0.295081481
link80   0.69314718  0.36395125 0.265725370  0.417182713
link81   1.09861229 -0.05782817 0.380972463  1.210454413  Log(m)

link82   1.00207436  0.99471736 0.090538025  0.578628110
link83   2.79379093  2.85229949 0.112086258 -0.420360670
link84   0.91975294  0.91376213 0.088673239  0.422757445
link85   2.79379093  2.77134426 0.111256896  0.170431235
link86  -1.26125559 -1.34215417 0.086191456  1.586060116
link87   0.47321734  0.51542796 0.079998129 -1.671748374
link88   0.50117219  0.53233734 0.082594968 -1.987642815
link89   2.65147985  2.38991947 0.103139188  2.518485179  Oneway C-Logits

link90   0.73146623  0.73001868 0.017844267  0.577647982
link91   0.21087315  0.21541874 0.013701431 -0.477265116
link92   0.05766063  0.05456258 0.005782026  0.431482851
link93   0.71499176  0.71376939 0.018116177  0.422215472
link94   0.22734761  0.22733814 0.013759216  0.000946622
link95   0.05766063  0.05889247 0.006166318 -0.168753257
link96   0.22075783  0.20715603 0.014156292  1.623639811
link97   0.39538715  0.41892202 0.014035304 -1.647600352
link98   0.38385502  0.37392195 0.018727908  1.680441903
link99   0.62273476  0.63002809 0.019252284 -1.995568028
link100  0.31136738  0.28602728 0.014261072  2.196801000
link101  0.06589786  0.08394462 0.007931189 -2.259707404  Oneway Probs


CELL-SPECIFIC STATISTICS...
    OBS          FV StdErr(FV)         PROB StdErr(PROB) ADJUSTED RESIDS
y1   62 58.34560285 5.60048409 0.0961212568 9.226498e-03      0.79048547
y2   17 17.98696892 2.18817944 0.0296325682 3.604908e-03     -0.27732338
y3    5  3.21658139 0.84771899 0.0052991456 1.396572e-03      1.13226452
y4   90 99.14748066 5.50075610 0.1633401658 9.062201e-03     -1.26014140
y5   42 37.27942176 3.12090715 0.0614158513 5.141527e-03      0.93943366
y6    3  8.13098873 1.18909696 0.0133953686 1.958974e-03     -1.99600725
y7   74 70.64307787 6.48702129 0.1163806884 1.068702e-02      0.74432785
y8   31 32.39624900 3.23190673 0.0533710857 5.324393e-03     -0.31049225
y9   11  8.61800122 1.67525668 0.0141976956 2.759896e-03      0.99865951
y10  11 11.89014935 1.74920935 0.0195883844 2.881729e-03     -0.30358241
y11   7  5.77578116 0.87768807 0.0095152902 1.445944e-03      0.55021841
y12   0  1.62749948 0.42083501 0.0026812182 6.933031e-04     -1.35341708
y13  22 21.27416474 2.20616942 0.0350480473 3.634546e-03      0.18341011
y14  18 12.60415019 1.46797411 0.0207646626 2.418409e-03      1.69054385
y15   1  4.33172950 0.62127245 0.0071362924 1.023513e-03     -1.68388668
y16  19 15.95997123 2.18059618 0.0262931981 3.592415e-03      0.92568008
y17  14 11.53269448 1.45184697 0.0189994967 2.391840e-03      0.81319381
y18   3  4.83410609 0.88595030 0.0079639310 1.459556e-03     -0.91582898
y19   2  1.59308619 0.50117004 0.0026245242 8.256508e-04      0.35181763
y20   3  1.21937212 0.34090919 0.0020088503 5.616296e-04      1.69721942
y21   1  0.54140213 0.21605335 0.0008919310 3.559363e-04      0.65232348
y22   2  3.00120829 0.64594506 0.0049443300 1.064160e-03     -0.62464300
y23   0  2.80175888 0.45331956 0.0046157477 7.468197e-04     -1.74317496
y24   1  1.51723231 0.40226935 0.0024995590 6.627172e-04     -0.44489596
y25   1  2.37064968 0.67012008 0.0039055184 1.103987e-03     -0.99116289
y26   3  2.69922838 0.58636796 0.0044468342 9.660098e-04      0.19647796
y27   1  1.78278343 0.55122649 0.0029370402 9.081161e-04     -0.64481862
y28  11  9.85218019 1.56720577 0.0162309394 2.581888e-03      0.42669775
y29   3  3.02666588 0.57235288 0.0049862700 9.429207e-04     -0.01627668
y30   0  0.53936570 0.17839926 0.0008885761 2.939032e-04     -0.75744934
y31  21 22.92519441 2.32544609 0.0377680303 3.831048e-03     -0.47178691
y32  13  8.58979489 1.11903353 0.0141512272 1.843548e-03      1.64176293
y33   2  1.86697829 0.38823037 0.0030757468 6.395888e-04      0.10170899
y34  20 22.36701794 2.69129542 0.0368484645 4.433765e-03     -0.62594652
y35   8 10.22152002 1.34996143 0.0168394070 2.223989e-03     -0.77451438
y36   3  2.70962749 0.61125405 0.0044639662 1.007008e-03      0.19047885
y37   1  3.30589824 0.60980866 0.0054462903 1.004627e-03     -1.35034389
y38   4  1.60027707 0.30584324 0.0026363708 5.038604e-04      1.95772303
y39   0  0.44935295 0.13164653 0.0007402849 2.168806e-04     -0.68391457
y40   6  8.09956730 1.07189328 0.0133436035 1.765887e-03     -0.80264092
y41   9  4.78195134 0.69409189 0.0078780088 1.143479e-03      2.04304813
y42   0  1.63770318 0.27563403 0.0026980283 4.540923e-04     -1.31234369
y43   6  8.32048293 1.18757921 0.0137075501 1.956473e-03     -0.89012085
y44   5  5.99141633 0.82149548 0.0098705376 1.353370e-03     -0.43238399
y45   2  2.50263303 0.44993915 0.0041229539 7.412507e-04     -0.33215871
y46   1  0.72932215 0.22851721 0.0012015192 3.764699e-04      0.32915987
y47   0  0.55628669 0.15607527 0.0009164525 2.571256e-04     -0.76309883
y48   1  0.24613005 0.09721463 0.0004054861 1.601559e-04      1.54991475
y49   2  1.88140977 0.36314607 0.0030995219 5.982637e-04      0.08980750
y50   1  1.75025085 0.24728083 0.0028834446 4.073819e-04     -0.57813358
y51   1  0.94450422 0.22511439 0.0015560201 3.708639e-04      0.05874763
y52   4  2.03498760 0.48150318 0.0033525331 7.932507e-04      1.46613328
y53   3  2.30895942 0.38936058 0.0038038870 6.414507e-04      0.47144206
y54   1  1.51969887 0.39311705 0.0025036225 6.476393e-04     -0.44541658
y55   3  1.21206417 0.41227233 0.0019968108 6.791966e-04      1.75348896
y56   0  0.37105649 0.13747351 0.0006112957 2.264802e-04     -0.62547716
y57   0  0.06589328 0.03641435 0.0001085557 5.999068e-05     -0.25933371
y58   2  3.86201026 0.74876340 0.0063624551 1.233548e-03     -1.02862636
y59   1  1.44200070 0.31972339 0.0023756190 5.267272e-04     -0.38235126
y60   0  0.31232319 0.12443050 0.0005145357 2.049926e-04     -0.57340362
y61   9  5.15959498 1.17172840 0.0085001565 1.930360e-03      1.98508583
y62   2  2.34966145 0.54608915 0.0038709414 8.996526e-04     -0.24467094
y63   1  0.62069991 0.24856767 0.0010225699 4.095019e-04      0.50764158
y64   1  0.66966979 0.21926800 0.0011032451 3.612323e-04      0.41923166
y65   0  0.32303434 0.10534287 0.0005321818 1.735467e-04     -0.57854187
y66   0  0.09039062 0.04297417 0.0001489137 7.079764e-05     -0.30379261
y67   2  2.24667458 0.40299829 0.0037012761 6.639181e-04     -0.17120475
y68   1  1.32180012 0.22099082 0.0021775949 3.640705e-04     -0.28554190
y69   0  0.45110549 0.13992124 0.0007431721 2.305127e-04     -0.68697767
y70   4  3.16034151 0.64963802 0.0052064934 1.070244e-03      0.50894397
y71   2  2.26776090 0.38820523 0.0037360147 6.395473e-04     -0.18439604
y72   0  0.94394615 0.29132684 0.0015551007 4.799454e-04     -1.01930250
y73   1  0.24325881 0.11399076 0.0004007559 1.877937e-04      1.57734044
y74   0  0.18489711 0.07982735 0.0003046081 1.315113e-04     -0.43767273
y75   0  0.08152267 0.04568714 0.0001343042 7.526712e-05     -0.28926867
y76   0  0.85928989 0.27010946 0.0014156341 4.449909e-04     -0.96978006
y77   0  0.79659741 0.20338114 0.0013123516 3.350595e-04     -0.91727393
y78   0  0.42837565 0.16303278 0.0007057260 2.685878e-04     -0.67606031
y79   1  1.27269754 0.43398766 0.0020967011 7.149714e-04     -0.26219902
y80   2  1.43900406 0.38237989 0.0023706821 6.299504e-04      0.49404845
y81   3  0.94381210 0.35956642 0.0015548799 5.923664e-04      2.28035475


CONVERGENCE STATISTICS...
    iterations = 9
    norm.diff  = 1.39606e-06
    norm.score = 4.5618e-07
    Original counts used.

FITTING PROGRAM USED:  mph.fit, version 1.0, 6/5/02 

 Document Last Updated:  02/08/2007 12:00 AM -0600, Joseph B. Lang