ML FITTING of MPH MODELS using
MPH.FIT: NUMERICAL EXAMPLESAuthor: 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.)
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 1. ML fit of mean-response models
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)
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:
γ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.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 =
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
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
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 | |