This site contains examples of the use of Joseph Lang's very useful mph.fit function for some of the analyses in the text that relied on that function. Of course, there are much more elegant ways of creating the model matrices than the cell-by-cell entry shown here! The function and a detailed manual can be obtained from Joseph Lang's home page.

1. mph.fit analysis for finding standardized residuals (called "adjusted residuals" in this output) for analysis in Section 3.5.8 on cumulative logit model for region and religious beliefs.

> y <- scan()
1: 92
2: 352
3: 234
4: 274
5: 399
6: 326
7: 739
8: 536
9: 412
10: 192
11: 423
12: 355
13:
Read 12 items
> y <- matrix(y,12,1,byrow=T)
> z <- c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1)
> Z <- matrix(z,12,4)
> 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
> M <- c(1,0,0,0,0,0,0,0,0,0,0,0,
+ 0,1,1,0,0,0,0,0,0,0,0,0,
+ 1,1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,0,0,0,0,0,0,0,0,0,
+ 0,0,0,1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,1,0,0,0,0,0,0,
+ 0,0,0,1,1,0,0,0,0,0,0,0,
+ 0,0,0,0,0,1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,0,0,0,0,0,
+ 0,0,0,0,0,0,0,1,1,0,0,0,
+ 0,0,0,0,0,0,1,1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,0,0,0,
+ 0,0,0,0,0,0,0,0,0,1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,1,
+ 0,0,0,0,0,0,0,0,0,1,1,0,
+ 0,0,0,0,0,0,0,0,0,0,0,1)
> M <- matrix(M,16,12,byrow=T)
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 1 0 0 0 0 0 0 0 0 0
[3,] 1 1 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 1 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 1 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 1 1 0 0 0 0 0 0
[7,] 0 0 0 1 1 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 1 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 1 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 1 1 0 0 0
[11,] 0 0 0 0 0 0 1 1 0 0 0 0
[12,] 0 0 0 0 0 0 0 0 1 0 0 0
[13,] 0 0 0 0 0 0 0 0 0 1 0 0
[14,] 0 0 0 0 0 0 0 0 0 0 1 1
[15,] 0 0 0 0 0 0 0 0 0 1 1 0
[16,] 0 0 0 0 0 0 0 0 0 0 0 1
> C <- c(1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,-1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1)
> C <- matrix(C,8,16,byrow=T)
> C
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 -1 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 1 -1 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 1 -1 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 1 -1 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 1 -1 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 1 -1 0 0
[7,] 0 0 0 0 0 0 0 0 0 0 0 0 1 -1
[8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[,15] [,16]
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] 0 0
[5,] 0 0
[6,] 0 0
[7,] 0 0
[8,] 1 -1
> X <- c(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,
+ 1,0,0,0,0,
+ 0,1,0,0,0)
> X <- matrix(X,8,5,byrow=T)
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 0 0
[2,] 0 1 1 0 0
[3,] 1 0 0 1 0
[4,] 0 1 0 1 0
[5,] 1 0 0 0 1
[6,] 0 1 0 0 1
[7,] 1 0 0 0 0
[8,] 0 1 0 0 0
> y
[,1]
[1,] 92
[2,] 352
[3,] 234
[4,] 274
[5,] 399
[6,] 326
[7,] 739
[8,] 536
[9,] 412
[10,] 192
[11,] 423
[12,] 355
> L.fct <- function(m){
+ C%*%log(M%*%m)}
> L.fct(y)
[,1]
[1,] -1.8515312
[2,] 0.6405034
[3,] -0.9730435
[4,] 0.7248479
[5,] -0.2490566
[6,] 1.1296781
[7,] -1.3992312
[8,] 0.5495045

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

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 53.60521 norm.score= 29.62958
iter= 2 norm.diff= 8.708879 norm.score= 2.526581
iter= 3 norm.diff= 0.6733886 norm.score= 0.1478496
iter= 4 norm.diff= 0.04016614 norm.score= 0.005596941
iter= 5 norm.diff= 0.0001390173 norm.score= 0.0005161505
iter= 6 norm.diff= 2.666436e-06 norm.score= 4.756527e-05
iter= 7 norm.diff= 2.354554e-07 norm.score= 4.399225e-06

> mph.summary(mph.out,cell.stats=T,model.info=T)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 3 ): Gsq = 98.0238 (p = 0 )
Pearson's Score Stat (df= 3 ): Xsq = 97.50184 (p = 0 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = -1.261818 ( 0.064033 )
beta2 = 0.472851 ( 0.061096 )
beta3 = -0.069839 ( 0.093035 )
beta4 = 0.268777 ( 0.083536 )
beta5 = 0.889678 ( 0.075704 )

link1 = -1.331658 ( 0.075439 )
link2 = 0.403012 ( 0.072547 )
link3 = -0.993041 ( 0.061793 )
link4 = 0.741629 ( 0.060814 )
link5 = -0.372141 ( 0.046974 )
link6 = 1.362529 ( 0.051581 )
link7 = -1.261818 ( 0.064033 )
link8 = 0.472851 ( 0.061096 )

CELL-SPECIFIC STATISTICS...
OBS FV StdErr(FV) PROB StdErr(PROB) ADJUSTED RESIDS
y1 92 141.6242 8.452307 0.2088853 0.012466529 -7.7881339
y2 352 264.7765 5.995397 0.3905258 0.008842768 7.7881335
y3 234 271.5993 11.810687 0.4005889 0.017419892 -7.7881328
y4 274 270.0415 12.176088 0.2703118 0.012188276 0.5667163
y5 399 406.6330 7.727063 0.4070400 0.007734798 -0.5667163
y6 326 322.3255 13.277291 0.3226481 0.013290582 0.5667163
y7 739 688.3363 19.140872 0.4080239 0.011346101 7.9021674
y8 536 654.8032 13.214185 0.3881465 0.007832949 -7.9021673
y9 412 343.8605 14.121521 0.2038296 0.008370789 7.9021673
y10 192 214.0412 10.681338 0.2206611 0.011011689 -3.0356333
y11 423 383.5354 7.929392 0.3953974 0.008174631 3.0356333
y12 355 372.4233 14.017542 0.3839416 0.014451074 -3.0356333

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

2. Using mph.fit for partial proportional odds model analysis in Section 3.6.5.

> y <- scan()
1: 334
2: 99
3: 117
4: 159
5: 30
6: 350
7: 307
8: 345
9: 481
10: 67
11:
Read 10 items
> y <- matrix(y,10,1,byrow=T)
> y
[,1]
[1,] 334
[2,] 99
[3,] 117
[4,] 159
[5,] 30
[6,] 350
[7,] 307
[8,] 345
[9,] 481
[10,] 67
> z <- c(1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1)
> Z <- matrix(z,10,2)
> 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

> M <- matrix(16,10,byrow=T)
> M <- c(1,0,0,0,0,0,0,0,0,0,
+ 0,1,1,1,1,0,0,0,0,0,
+ 1,1,0,0,0,0,0,0,0,0,
+ 0,0,1,1,1,0,0,0,0,0,
+ 1,1,1,0,0,0,0,0,0,0,
+ 0,0,0,1,1,0,0,0,0,0,
+ 1,1,1,1,0,0,0,0,0,0,
+ 0,0,0,0,1,0,0,0,0,0,
+ 0,0,0,0,0,1,0,0,0,0,
+ 0,0,0,0,0,0,1,1,1,1,
+ 0,0,0,0,0,1,1,0,0,0,
+ 0,0,0,0,0,0,0,1,1,1,
+ 0,0,0,0,0,1,1,1,0,0,
+ 0,0,0,0,0,0,0,0,1,1,
+ 0,0,0,0,0,1,1,1,1,0,
+ 0,0,0,0,0,0,0,0,0,1)
> M <- matrix(M,16,10,byrow=T)
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 1 1 1 0 0 0 0 0
[3,] 1 1 0 0 0 0 0 0 0 0
[4,] 0 0 1 1 1 0 0 0 0 0
[5,] 1 1 1 0 0 0 0 0 0 0
[6,] 0 0 0 1 1 0 0 0 0 0
[7,] 1 1 1 1 0 0 0 0 0 0
[8,] 0 0 0 0 1 0 0 0 0 0
[9,] 0 0 0 0 0 1 0 0 0 0
[10,] 0 0 0 0 0 0 1 1 1 1
[11,] 0 0 0 0 0 1 1 0 0 0
[12,] 0 0 0 0 0 0 0 1 1 1
[13,] 0 0 0 0 0 1 1 1 0 0
[14,] 0 0 0 0 0 0 0 0 1 1
[15,] 0 0 0 0 0 1 1 1 1 0
[16,] 0 0 0 0 0 0 0 0 0 1
> C <- c(1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
C+ 0,0,1,-1,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,-1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,-1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,1,-1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1)
> C <- matrix(C,8,16,byrow=T)
> C
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 -1 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 1 -1 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 1 -1 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 1 -1 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 1 -1 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 1 -1 0 0
[7,] 0 0 0 0 0 0 0 0 0 0 0 0 1 -1
[8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[,15] [,16]
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] 0 0
[5,] 0 0
[6,] 0 0
[7,] 0 0
[8,] 1 -1
> L.fct <- function(m){
+ C%*%log(M)%*%m}
> X <- c(1,0,0,0,0,0,
+ 0,1,0,0,0,0,
+ 0,0,1,0,0,0,
+ 0,0,0,1,0,0,
+ 1,0,0,0,1,0,
+ 0,1,0,0,1,1,
+ 0,0,1,0,1,2,
+ 0,0,0,1,1,3)
> X <- matrix(X,8,6,byrow=T)
> X
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 1 0 0 0 1 0
[6,] 0 1 0 0 1 1
[7,] 0 0 1 0 1 2
[8,] 0 0 0 1 1 3

> L.fct <- function(m){
+ C%*%log(M%*%m)}
> L.fct(y)
[,1]
[1,] -0.1927461
[2,] 0.3471526
[3,] 1.0681713
[4,] 3.1626581
[5,] -1.2321437
[6,] -0.3069026
[7,] 0.6034780
[8,] 3.0971297
> mph.out <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 49.48404 norm.score= 1.116227
iter= 2 norm.diff= 0.8209376 norm.score= 0.05379518
iter= 3 norm.diff= 0.1169749 norm.score= 0.001880014
iter= 4 norm.diff= 0.0005789811 norm.score= 0.0001078673
iter= 5 norm.diff= 0.0001584866 norm.score= 6.285434e-06
iter= 6 norm.diff= 2.732279e-06 norm.score= 3.773991e-07

> mph.summary(mph.out)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 2 ): Gsq = 3.42938 (p = 0.18002 )
Pearson's Score Stat (df= 2 ): Xsq = 3.44836 (p = 0.17832 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = -0.200966 ( 0.073764 )
beta2 = 0.386943 ( 0.070829 )
beta3 = 1.040001 ( 0.081528 )
beta4 = 3.203276 ( 0.13751 )
beta5 = -1.0171 ( 0.09443 )
beta6 = 0.297788 ( 0.047339 )

link1 = -0.200966 ( 0.073764 )
link2 = 0.386943 ( 0.070829 )
link3 = 1.040001 ( 0.081528 )
link4 = 3.203276 ( 0.13751 )
link5 = -1.218066 ( 0.059772 )
link6 = -0.332368 ( 0.04931 )
link7 = 0.618477 ( 0.052339 )
link8 = 3.07954 ( 0.109599 )

CONVERGENCE STATISTICS...
iterations = 6
norm.diff = 2.73228e-06
norm.score = 3.77399e-07
Original counts used.

3. mph.fit analysis of mean response model in Section 5.6.2.

> y <- scan()
1: 28
2: 120
3: 89
4: 202
5: 51
6: 37
7: 10
8: 20
9: 79
10: 124
11: 362
12: 120
13: 90
14: 18
15: 3
16: 12
17: 26
18: 128
19: 107
20: 211
21: 47
22: 42
23: 201
24: 136
25: 320
26: 83
27: 63
28: 18
29: 33
30: 87
31: 107
32: 459
33: 123
34: 92
35: 19
36: 5
37: 19
38: 29
39: 177
40: 121
41: 183
42: 52
43:
Read 42 items
> y <- matrix(y,42,1,byrow=T)
> y
[,1]
[1,] 28
[2,] 120
[3,] 89
[4,] 202
[5,] 51
[6,] 37
[7,] 10
[8,] 20
[9,] 79
[10,] 124
[11,] 362
[12,] 120
[13,] 90
[14,] 18
[15,] 3
[16,] 12
[17,] 26
[18,] 128
[19,] 107
[20,] 211
[21,] 47
[22,] 42
[23,] 201
[24,] 136
[25,] 320
[26,] 83
[27,] 63
[28,] 18
[29,] 33
[30,] 87
[31,] 107
[32,] 459
[33,] 123
[34,] 92
[35,] 19
[36,] 5
[37,] 19
[38,] 29
[39,] 177
[40,] 121
[41,] 183
[42,] 52
> z <- c(1,1,1,1,1,1,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,0,0,0,0,0,1,1,1,1,1,1,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,0,0,0,0,0,1,1,1,1,1,1,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,0,0,0,0,0,1,1,1,1,1,1,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,0,0,0,0,0,1,1,1,1,1,1,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,0,0,0,0,0,1,1,1,1,1,1,1)
> Z <- matrix(z,42,6)
> Z
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 1 0 0 0 0 0
[3,] 1 0 0 0 0 0
[4,] 1 0 0 0 0 0
[5,] 1 0 0 0 0 0
[6,] 1 0 0 0 0 0
[7,] 1 0 0 0 0 0
[8,] 0 1 0 0 0 0
[9,] 0 1 0 0 0 0
[10,] 0 1 0 0 0 0
[11,] 0 1 0 0 0 0
[12,] 0 1 0 0 0 0
[13,] 0 1 0 0 0 0
[14,] 0 1 0 0 0 0
[15,] 0 0 1 0 0 0
[16,] 0 0 1 0 0 0
[17,] 0 0 1 0 0 0
[18,] 0 0 1 0 0 0
[19,] 0 0 1 0 0 0
[20,] 0 0 1 0 0 0
[21,] 0 0 1 0 0 0
[22,] 0 0 0 1 0 0
[23,] 0 0 0 1 0 0
[24,] 0 0 0 1 0 0
[25,] 0 0 0 1 0 0
[26,] 0 0 0 1 0 0
[27,] 0 0 0 1 0 0
[28,] 0 0 0 1 0 0
[29,] 0 0 0 0 1 0
[30,] 0 0 0 0 1 0
[31,] 0 0 0 0 1 0
[32,] 0 0 0 0 1 0
[33,] 0 0 0 0 1 0
[34,] 0 0 0 0 1 0
[35,] 0 0 0 0 1 0
[36,] 0 0 0 0 0 1
[37,] 0 0 0 0 0 1
[38,] 0 0 0 0 0 1
[39,] 0 0 0 0 0 1
[40,] 0 0 0 0 0 1
[41,] 0 0 0 0 0 1
[42,] 0 0 0 0 0 1

> M <- c(1,2,3,4,5,6,7,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,0,0,0,0,0,1,2,3,4,5,6,7,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,0,0,0,0,0,1,2,3,4,5,6,7,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,0,0,0,0,0,1,2,3,4,5,6,7,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,0,0,0,0,0,1,2,3,4,5,6,7,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,0,0,0,0,0,1,2,3,4,5,6,7)
> M <- matrix(M,6,42,byrow=T)
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 2 3 4 5 6 7 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 1 2 3 4 5 6 7
[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0 0 0
[3,] 1 2 3 4 5 6 7 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 1 2 3 4 5
[5,] 0 0 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 0 0
[,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0 0 0
[4,] 6 7 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 1 2 3 4 5 6 7 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 1 2 3
[,39] [,40] [,41] [,42]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0
[4,] 0 0 0 0
[5,] 0 0 0 0
[6,] 4 5 6 7
> X <- c(1,0,1,0,1,0,0,1,1,0,0,0,1,1,1,0,1,1,0,1,1,1,0,0)
> X <- matrix(X,6,4,byrow=T)
> X
[,1] [,2] [,3] [,4]
[1,] 1 0 1 0
[2,] 1 0 0 1
[3,] 1 0 0 0
[4,] 1 1 1 0
[5,] 1 1 0 1
[6,] 1 1 0 0
> L.fct <- function(m){
+ p <- diag(c(1/Z%*%t(Z)%*%m))%*%m
+ M%*%p
+ }
> L.fct(y)
[,1]
[1,] 3.519553
[2,] 4.014760
[3,] 5.144195
[4,] 3.535342
[5,] 3.982609
[6,] 4.957338
> source("mph.r")
> mph.out <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 12.97496 norm.score= 0.8194791
iter= 2 norm.diff= 0.05699647 norm.score= 0.02810338
iter= 3 norm.diff= 0.02331553 norm.score= 0.002753730
iter= 4 norm.diff= 0.0004092286 norm.score= 0.0003541458
iter= 5 norm.diff= 0.0002700452 norm.score= 4.625873e-05
iter= 6 norm.diff= 1.070353e-05 norm.score= 7.268501e-06
iter= 7 norm.diff= 5.196447e-06 norm.score= 1.070365e-06

Time Elapsed: 0.432 seconds

> mph.summary(mph.out)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 2 ): Gsq = 4.17616 (p = 0.12392 )
Pearson's Score Stat (df= 2 ): Xsq = 4.16588 (p = 0.12456 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 5.08051 ( 0.042231 )
beta2 = -0.062937 ( 0.039552 )
beta3 = -1.512681 ( 0.052151 )
beta4 = -1.049344 ( 0.047567 )

link1 = 3.567828 ( 0.043844 )
link2 = 4.031165 ( 0.036551 )
link3 = 5.08051 ( 0.042231 )
link4 = 3.504892 ( 0.039599 )
link5 = 3.968229 ( 0.035114 )
link6 = 5.017573 ( 0.041713 )

CONVERGENCE STATISTICS...
iterations = 7
norm.diff = 5.19645e-06
norm.score = 1.07037e-06
Original counts used.

4. mph.fit analysis of global odds ratio model in Section 6.6.2, using some matrices shows that were input in an earlier run.

> y
[,1]
[1,] 13
[2,] 80
[3,] 37
[4,] 1
[5,] 22
[6,] 17
[7,] 0
[8,] 3
[9,] 4
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 1 1 0 1 1
[3,] 0 1 1 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 1 0 0
[5,] 1 1 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 1
[7,] 0 0 1 0 0 0 0 0 0
[8,] 0 0 0 1 1 0 1 1 0
[9,] 1 0 0 1 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 1 1
[11,] 0 1 1 0 1 1 0 0 0
[12,] 0 0 0 0 0 0 1 0 0
[13,] 1 1 0 1 1 0 0 0 0
[14,] 0 0 0 0 0 0 0 0 1
[15,] 0 0 1 0 0 1 0 0 0
[16,] 0 0 0 0 0 0 1 1 0

> C
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 1 -1 -1 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 1 1 -1 -1 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 1 1 -1 -1 0 0
[4,] 0 0 0 0 0 0 0 0 0 0 0 0 1 1
[,15] [,16]
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] -1 -1
> L.fct <- function(m){C%*%log(M%*%m)}
> z <- y+0.000000001
> mph.out <- mph.fit(y=z,Z=Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 3.584269 norm.score= 0.6011723
iter= 2 norm.diff= 0.2841504 norm.score= 0.03586807
iter= 3 norm.diff= 0.007202971 norm.score= 0.001762248
iter= 4 norm.diff= 0.0004832279 norm.score= 6.538237e-06
iter= 5 norm.diff= 2.10661e-06 norm.score= 5.995607e-08

Time Elapsed: 0.02 seconds

> mph.summary(mph.out,cell.stats=T,model.info=T)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 3 ): Gsq = 1.39139 (p = 0.70755 )
Pearson's Score Stat (df= 3 ): Xsq = 1.06792 (p = 0.78482 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 0.809007 ( 0.337416 )
link1 = 0.809007 ( 0.337416 )
link2 = 0.809007 ( 0.337416 )
link3 = 0.809007 ( 0.337416 )
link4 = 0.809007 ( 0.337416 )

CELL-SPECIFIC STATISTICS...
OBS FV StdErr(FV) PROB StdErr(PROB) ADJUSTED RESIDS
y1 1.3e+01 12.0137090 3.1573537 0.067874062 0.0178381566 0.8895230
y2 8.0e+01 81.8904149 6.3263986 0.462657711 0.0357423652 -0.9475956
y3 3.7e+01 36.0047606 5.2649427 0.203416726 0.0297454387 1.0151442
y4 1.0e+00 1.7781173 0.7070047 0.010045860 0.0039943766 -0.6930922
y5 2.2e+01 20.0982442 3.6730108 0.113549403 0.0207514737 0.9144439
y6 1.7e+01 18.1851070 3.6793401 0.102740718 0.0207872321 -0.7108806
y7 1.0e-09 0.2659720 0.1447422 0.001502666 0.0008177524 -0.5377595
y8 3.0e+00 3.1613416 1.2815885 0.017860687 0.0072406131 -0.1334173
y9 4.0e+00 3.6023335 1.4902023 0.020352166 0.0084192219 0.3476665

CONVERGENCE STATISTICS...
iterations = 5
norm.diff = 2.10661e-06
norm.score = 5.99561e-08
Original counts used.

MODEL INFORMATION...
Linear Predictor Model Link Function L.fct:
function(m){
+ C%*%log(M%*%m)}

5. mph.fit for ML fitting of marginal cumulative logit model for a 4x4 table that is an updating from the 2008 GSS of Table 10.5 on premarital and extramarital sex in the 2002 2nd edition of my book "Categorical Data Analysis". This uses version 3.1 of mph.fit. The model is the one discussed for a 3x3 table on pages 231-232 of "Analysis of Ordinal Categorical Data"

> source("mph.Rcode.txt")
> y <- c(324,6,1,0,95,16,1,0,185,30,17,1,462,120,51,28)
> M <- c(1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,
+ 1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,
+ 1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,
+ 0,1,1,1,0,1,1,1,0,1,1,1,0,1,1,1,
+ 1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,
+ 0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,
+ 1,1,1,0,1,1,1,0,1,1,1,0,1,1,1,0,
+ 0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1)
> M <- matrix(M,12,16,byrow=T)
> C <- c(1,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,-1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,-1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,-1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,-1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,-1)
> C <- matrix(C,6,12,byrow=T)
> X <- c(1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,1,0,1,0,0,1,1)
> X <- matrix(X,6,4,byrow=T)
> L.fct <- function(m){C%*%log(M%*%m)}
> L.fct(y)
[,1]
[1,] -1.11161898
[2,] -0.70213601
[3,] 0.02243924
[4,] 1.36954978
[5,] 2.52613260
[6,] 3.80895870

> mph.out <- mph.fit(y=y,L.fct=L.fct,X=X)

mph.fit, version 3.1, 5/20/09 , running...
Did NOT meet norm score convergence criterion [1e-06]!
Did NOT meet norm diff convergence criterion [1e-06]!
Time Elapsed: 0.809 seconds , Iterations: 100

> mph.summary(mph.out,cell.stats=T,model.info=T)

MODEL GOODNESS OF FIT: Test of Ho: h(p)=0 vs. Ha: not Ho...

Likelihood Ratio Stat (df= 2 ): Gsq = 105.1386 (pval = 0 )
Pearson's Score Stat (df= 2 ): Xsq = 99.84842 (pval = 0 )
Generalized Wald Stat (df= 2 ): Wsq = 78.74009 (pval = 0 )

Adj Resids: -9.413 -8.728 ... 9.951 9.987 , Number |Adj Resid| > 2: 13

SAMPLING PLAN INFORMATION...
Number of strata: 1
Strata identifiers: 1
Strata with fixed sample sizes: all
Observed strata sample sizes: 1337

LINEAR PREDICTOR MODEL RESULTS...
BETA StdErr(BETA) Z-ratio p-value
beta1 -1.3474 0.0627 -21.5025 0.00000
beta2 -0.6959 0.0565 -12.3259 0.00000
beta3 0.0695 0.0544 1.2775 0.20143
beta4 2.7802 0.0789 35.2577 0.00000

OBS LINK ML LINK StdErr(L) LINK RESID
link1 -1.1116 -1.3474 0.0627 9.3100
link2 -0.7021 -0.6959 0.0565 -0.4603
link3 0.0224 0.0695 0.0544 -7.7315
link4 1.3695 1.4328 0.0690 -9.8013
link5 2.5261 2.0843 0.0770 10.7799
link6 3.8090 2.8497 0.0852 11.3024

CELL-SPECIFIC STATISTICS...
strata OBS FV StdErr.FV PROB StdErr.PROB ADJUSTED.RESIDS
y1 1 324 271.6084 13.7373 0.2031 0.0103 9.9513
y2 1 6 3.3336 1.8028 0.0025 0.0013 9.7243
y3 1 1 0.8761 0.9325 0.0007 0.0007 1.5921
y4 1 0 0.0000 0.0000 0.0000 0.0000 0.0000
y5 1 95 154.1412 9.5108 0.1153 0.0071 -8.7279
y6 1 16 13.0795 3.5869 0.0098 0.0027 9.9872
y7 1 1 1.7706 1.3138 0.0013 0.0010 -3.7581
y8 1 0 0.0267 0.1633 0.0000 0.0001 -9.4134
y9 1 185 203.7103 11.8945 0.1524 0.0089 -3.3501
y10 1 30 19.8004 4.1973 0.0148 0.0031 7.4189
y11 1 17 19.8454 4.2484 0.0148 0.0032 -2.3216
y12 1 1 3.5234 1.8357 0.0026 0.0014 -6.6431
y13 1 462 449.9595 17.1724 0.3365 0.0128 6.3149
y14 1 120 73.4517 6.7884 0.0549 0.0051 9.6363
y15 1 51 52.2954 5.3369 0.0391 0.0040 -0.2777
y16 1 28 69.5778 5.9870 0.0520 0.0045 -7.5767

CONVERGENCE INFORMATION...
Original counts used.
iterations = 100 , time elapsed = 0.809
norm.diff = 0.0710777 = dist between last and second last iterates.
Did NOT meet norm diff convergence criterion [1e-06]!
norm.score = 0.000250401 = norm of score at last iteration.
Did NOT meet norm score convergence criterion [1e-06]!

MODEL INFORMATION...
Linear Predictor Model Link Function L.fct (L.mean= FALSE ):
function(m) {
p <- m*c(1/Z%*%t(Z)%*%m)
as.matrix(L.input.fct(p))
}


Link information as originally input, L.input.fct:
function(m){C%*%log(M%*%m)}

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

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

6. mph.fit for ML fitting of marginal cumulative logit model for crossover data in Section 8.4.3.

> y <- c(6,4,5,3,13,10,1,8,14,2,3,2,1,3,1,2,1,2,1,0,2,0,0,0,1,1,0)
> Z <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
> Z <- matrix(Z,27,1)
> M <- c(1,1,1,1,1,1,1,1,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,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,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,1,1,1,1,1,1,1,1,1,
+ 1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,
+ 0,0,0,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,
+ 1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,
+ 0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,
+ 1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,
+ 0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,
+ 1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,
+ 0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1)
> M <- matrix(M,12,27,byrow=T)
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 1 1 1 1 1 1 1 1 1 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 1 1 1 1
[3,] 1 1 1 1 1 1 1 1 1 1 1 1 1
[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[5,] 1 1 1 0 0 0 0 0 0 1 1 1 0
[6,] 0 0 0 1 1 1 1 1 1 0 0 0 1
[7,] 1 1 1 1 1 1 0 0 0 1 1 1 1
[8,] 0 0 0 0 0 0 1 1 1 0 0 0 0
[9,] 1 0 0 1 0 0 1 0 0 1 0 0 1
[10,] 0 1 1 0 1 1 0 1 1 0 1 1 0
[11,] 1 1 0 1 1 0 1 1 0 1 1 0 1
[12,] 0 0 1 0 0 1 0 0 1 0 0 1 0
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 1 1 1 1 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 1 1 1 1 1 1 1
[5,] 0 0 0 0 0 1 1 1 0 0 0 0
[6,] 1 1 1 1 1 0 0 0 1 1 1 1
[7,] 1 1 0 0 0 1 1 1 1 1 1 0
[8,] 0 0 1 1 1 0 0 0 0 0 0 1
[9,] 0 0 1 0 0 1 0 0 1 0 0 1
[10,] 1 1 0 1 1 0 1 1 0 1 1 0
[11,] 1 0 1 1 0 1 1 0 1 1 0 1
[12,] 0 1 0 0 1 0 0 1 0 0 1 0
[,26] [,27]
[1,] 0 0
[2,] 1 1
[3,] 0 0
[4,] 1 1
[5,] 0 0
[6,] 1 1
[7,] 0 0
[8,] 1 1
[9,] 0 0
[10,] 1 1
[11,] 1 0
[12,] 0 1
> C <- c(1,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,-1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,-1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,-1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,-1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,-1)
> C <- matrix(C,6,12,byrow=T)
> X <- c(1,0,1,0,0,1,1,0,1,0,0,1,0,1,0,1,1,0,0,0,0,1,0,0)
> X <- matrix(X,6,4,byrow=T)
> X1 <- c(1,0,0,0,0,1,0,0,1,0,1,0,0,1,1,0,1,0,0,1,0,1,0,1)
> X1 <- matrix(X1,6,4,byrow=T)
> L.fct <- function(m){C%*%log(M%*%m)}
> L.fct(y)
[,1]
[1,] 1.0678406
[2,] 2.7850112
[3,] -0.8919980
[4,] 0.6241543
[5,] -1.4008932
[6,] 0.3285041
> source("mph.r")
> mph.out <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 6.601511 norm.score= 0.1235732

iter= 100 norm.diff= 2.005297 norm.score= 3.409773e-10

Time Elapsed: 0.638 seconds

> mph.summary(mph.out,cell.stats=T,model.info=T)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 2 ): Gsq = 0.51744 (p = 0.77204 )
Pearson's Score Stat (df= 2 ): Xsq = 0.51601 (p = 0.77259 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = -1.337457 ( 0.235975 )
beta2 = 0.294815 ( 0.210123 )
beta3 = 2.430313 ( 0.371633 )
beta4 = 0.39203 ( 0.252098 )

link1 = 1.092856 ( 0.242898 )
link2 = 2.725128 ( 0.32268 )
link3 = -0.945427 ( 0.227661 )
link4 = 0.686845 ( 0.210785 )
link5 = -1.337457 ( 0.235975 )
link6 = 0.294815 ( 0.210123 )

CELL-SPECIFIC STATISTICS...
OBS FV StdErr(FV) PROB StdErr(PROB) ADJUSTED RESIDS
y1 6 6.056901e+00 2.333435e+00 7.042908e-02 2.713296e-02 -1.321493e-01
y2 4 3.563018e+00 1.723877e+00 4.143044e-02 2.004508e-02 6.560601e-01
y3 5 4.836730e+00 2.121688e+00 5.624105e-02 2.467079e-02 6.497257e-01
y4 3 3.533466e+00 1.669715e+00 4.108681e-02 1.941529e-02 -6.885068e-01
y5 13 1.325093e+01 3.140635e+00 1.540806e-01 3.651901e-02 -2.163199e-01
y6 10 1.120856e+01 2.624792e+00 1.303321e-01 3.052084e-02 -7.148621e-01
y7 1 1.023391e+00 1.002567e+00 1.189989e-02 1.165775e-02 -3.001664e-01
y8 8 7.212522e+00 2.258043e+00 8.386653e-02 2.625631e-02 6.410799e-01
y9 14 1.372153e+01 3.360837e+00 1.595526e-01 3.907950e-02 5.720208e-01
y10 2 1.948014e+00 1.318407e+00 2.265133e-02 1.533031e-02 1.277123e-01
y11 3 2.589053e+00 1.469696e+00 3.010527e-02 1.708949e-02 6.935361e-01
y12 2 1.869443e+00 1.307144e+00 2.173771e-02 1.519935e-02 3.766020e-01
y13 1 1.129809e+00 1.015935e+00 1.313731e-02 1.181319e-02 -4.509988e-01
y14 3 2.949435e+00 1.679153e+00 3.429576e-02 1.952504e-02 2.983342e-01
y15 1 1.077289e+00 1.012870e+00 1.252661e-02 1.177756e-02 -3.970701e-01
y16 2 1.973895e+00 1.328021e+00 2.295227e-02 1.544210e-02 6.427481e-02
y17 1 8.731618e-01 9.113024e-01 1.015304e-02 1.059654e-02 6.896597e-01
y18 2 1.893266e+00 1.319437e+00 2.201472e-02 1.534229e-02 3.208369e-01
y19 1 1.099928e+00 9.993917e-01 1.278986e-02 1.162083e-02 -3.386396e-01
y20 0 6.103868e-48 2.470601e-24 7.097521e-50 2.872792e-26 -6.103868e-48
y21 2 2.100177e+00 1.293406e+00 2.442066e-02 1.503961e-02 -1.633726e-01
y22 0 4.677488e-36 2.162750e-18 5.438940e-38 2.514826e-20 -4.677488e-36
y23 0 8.574940e-42 2.928300e-21 9.970861e-44 3.405000e-23 -8.574940e-42
y24 0 6.224554e-38 2.494906e-19 7.237853e-40 2.901053e-21 -6.224554e-38
y25 1 1.116459e+00 1.002743e+00 1.298209e-02 1.165980e-02 -3.749501e-01
y26 1 9.730213e-01 9.063653e-01 1.131420e-02 1.053913e-02 7.197146e-02
y27 0 1.701186e-43 4.124544e-22 1.978124e-45 4.795981e-24 -1.701186e-43

CONVERGENCE STATISTICS...
iterations = 100
norm.diff = 2.0053
norm.score = 3.40977e-10
Original counts used.

> mph.out2 <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X1)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 2.89777 norm.score= 0.1234315

iter= 100 norm.diff= 2.005297 norm.score= 6.56294e-10

Time Elapsed: 0.544 seconds

> mph.summary(mph.out2,cell.stats=T,model.info=T)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 2 ): Gsq = 0.51744 (p = 0.77204 )
Pearson's Score Stat (df= 2 ): Xsq = 0.51601 (p = 0.77259 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 1.092856 ( 0.242898 )
beta2 = 2.725128 ( 0.32268 )
beta3 = -2.038283 ( 0.360318 )
beta4 = -2.430313 ( 0.371633 )

link1 = 1.092856 ( 0.242898 )
link2 = 2.725128 ( 0.32268 )
link3 = -0.945427 ( 0.227661 )
link4 = 0.686845 ( 0.210785 )
link5 = -1.337457 ( 0.235975 )
link6 = 0.294815 ( 0.210123 )

7. mph.fit for ML fitting of Coke taste testing data treating responses as dependent, as discussed in Section 8.6.4.

> y <- scan()
1: 0
2: 1
3: 0
4: 1
5: 0
6: 0
7: 1
8: 4
9: 0
10: 1
11: 0
12: 0
13: 0
14: 0
15: 0
16: 0
17: 0
18: 0
19: 1
20: 2
21: 0
22: 0
23: 0
24: 1
25: 1
26: 1
27: 0
28: 0
29: 0
30: 3
31: 1
32: 0
33: 2
34: 1
35: 2
36: 0
37: 2
38: 1
39: 2
40: 0
41: 0
42: 0
43: 1
44: 2
45: 1
46: 0
47: 0
48: 0
49: 0
50: 0
51: 0
52: 1
53: 0
54: 0
55: 0
56: 1
57: 0
58: 2
59: 0
60: 0
61: 0
62: 0
63: 1
64: 1
65: 0
66: 1
67: 1
68: 0
69: 1
70: 0
71: 0
72: 1
73: 0
74: 1
75: 0
76: 1
77: 0
78: 0
79: 2
80: 0
81: 0
82: 1
83: 0
84: 0
85: 1
86: 0
87: 1
88: 1
89: 2
90: 1
91: 0
92: 2
93: 0
94: 1
95: 0
96: 0
97: 1
98: 0
99: 0
100: 0
101: 0
102: 0
103: 1
104: 0
105: 0
106: 1
107: 0
108: 0
109: 0
110: 0
111: 0
112: 0
113: 0
114: 0
115: 0
116: 0
117: 0
118: 0
119: 0
120: 0
121: 1
122: 0
123: 1
124: 0
125: 0
126:
Read 125 items

> y
[,1]
[1,] 0
[2,] 1
[3,] 0
[4,] 1
[5,] 0
[6,] 0
[7,] 1
[8,] 4
[9,] 0
[10,] 1
[11,] 0
[12,] 0
[13,] 0
[14,] 0
[15,] 0
[16,] 0
[17,] 0
[18,] 0
[19,] 1
[20,] 2
[21,] 0
[22,] 0
[23,] 0
[24,] 1
[25,] 1
[26,] 1
[27,] 0
[28,] 0
[29,] 0
[30,] 3
[31,] 1
[32,] 0
[33,] 2
[34,] 1
[35,] 2
[36,] 0
[37,] 2
[38,] 1
[39,] 2
[40,] 0
[41,] 0
[42,] 0
[43,] 1
[44,] 2
[45,] 1
[46,] 0
[47,] 0
[48,] 0
[49,] 0
[50,] 0
[51,] 0
[52,] 1
[53,] 0
[54,] 0
[55,] 0
[56,] 1
[57,] 0
[58,] 2
[59,] 0
[60,] 0
[61,] 0
[62,] 0
[63,] 1
[64,] 1
[65,] 0
[66,] 1
[67,] 1
[68,] 0
[69,] 1
[70,] 0
[71,] 0
[72,] 1
[73,] 0
[74,] 1
[75,] 0
[76,] 1
[77,] 0
[78,] 0
[79,] 2
[80,] 0
[81,] 0
[82,] 1
[83,] 0
[84,] 0
[85,] 1
[86,] 0
[87,] 1
[88,] 1
[89,] 2
[90,] 1
[91,] 0
[92,] 2
[93,] 0
[94,] 1
[95,] 0
[96,] 0
[97,] 1
[98,] 0
[99,] 0
[100,] 0
[101,] 0
[102,] 0
[103,] 1
[104,] 0
[105,] 0
[106,] 1
[107,] 0
[108,] 0
[109,] 0
[110,] 0
[111,] 0
[112,] 0
[113,] 0
[114,] 0
[115,] 0
[116,] 0
[117,] 0
[118,] 0
[119,] 0
[120,] 0
[121,] 1
[122,] 0
[123,] 1
[124,] 0
[125,] 0

> Z <- matrix(1,nrow=125,ncol=1)

> M <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,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,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,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,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,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,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,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,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,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,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,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,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,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,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,1,1,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,
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,1,1,1,1,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,
1,1,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,1,1,1,1,
1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1,1,1,1,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,1,1,
1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,
1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,
0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,
0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,
0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,
0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,
0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,
0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,
0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,
0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,
0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,
0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,
1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,
0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,
0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,
1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1)

> M <- matrix(M,15,125,byrow=T)
> dim(M)
[1] 15 125
> C <- c(1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,1,0,0,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,0,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,1,-1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0,
+ 0,0,0,0,0,0,1,0,0,-1,0,0,0,0,0,
+ 0,0,0,0,0,0,0,1,0,-1,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,-1,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,0,0,0,-1,
+ 0,0,0,0,0,0,0,0,0,0,0,1,0,0,-1,
+ 0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1)
> C <- matrix(C,12,15,byrow=T)
> C
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 1 0 0 0 -1 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 -1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 -1 0 0 0 0 0 0 0 0
[4,] 0 0 0 1 -1 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 1 0 0 0 -1 0 0 0
[6,] 0 0 0 0 0 0 1 0 0 -1 0 0 0
[7,] 0 0 0 0 0 0 0 1 0 -1 0 0 0
[8,] 0 0 0 0 0 0 0 0 1 -1 0 0 0
[9,] 0 0 0 0 0 0 0 0 0 0 1 0 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 1 0
[11,] 0 0 0 0 0 0 0 0 0 0 0 0 1
[12,] 0 0 0 0 0 0 0 0 0 0 0 0 0
[,14] [,15]
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] 0 0
[5,] 0 0
[6,] 0 0
[7,] 0 0
[8,] 0 0
[9,] 0 -1
[10,] 0 -1
[11,] 0 -1
[12,] 1 -1
> X <- c(0,0,-4,4,1,0,-3,3,0,1,-2,-2,1,0,-1,1,0,0,-4,0,1,0,-3,0,
+ 0,1,-2,0,1,0,-1,0,0,0,0,-4,1,0,0,-3,0,1,0,-2,1,0,0,-1)

> L.fct <- function(m){
+ C%*%log(M%*%m)}
L.fct <- function(m){
+ + C%*%log(M%*%m)}
> L.fct(y)
[,1]
[1,] 1.1786550
[2,] 1.5581446
[3,] 1.0116009
[4,] 1.2527630
[5,] 0.4519851
[6,] 0.9444616
[7,] 0.5389965
[8,] 0.6190392
[9,] -0.5389965
[10,] 0.0000000
[11,] 0.1541507
[12,] 0.2876821
> source("mph.r")
> X1 <- c(0,0,4,0,1,0,3,0,0,1,2,0,1,0,1,0,0,0,0,4,1,0,0,3,
+ 0,1,0,2,1,0,0,1,0,0,-4,4,1,0,-3,3,0,1,-2,2,1,0,-1,1)
> X1 <- matrix(X1,12,4,byrow=T)
> X1
[,1] [,2] [,3] [,4]
[1,] 0 0 4 0
[2,] 1 0 3 0
[3,] 0 1 2 0
[4,] 1 0 1 0
[5,] 0 0 0 4
[6,] 1 0 0 3
[7,] 0 1 0 2
[8,] 1 0 0 1
[9,] 0 0 -4 4
[10,] 1 0 -3 3
[11,] 0 1 -2 2
[12,] 1 0 -1 1
> mph.out2 <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X1)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 9.462715 norm.score= 0.6672981

iter= 100 norm.diff= 9.236443 norm.score= 4.270515e-10

Time Elapsed: 3.851 seconds

> mph.summary(mph.out2)

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 8 ): Gsq = 2.4871 (p = 0.96233 )
Pearson's Score Stat (df= 8 ): Xsq = 2.46806 (p = 0.9632 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 0.564523 ( 0.174241 )
beta2 = 0.390272 ( 0.238318 )
beta3 = 0.233658 ( 0.094088 )
beta4 = 0.110771 ( 0.078931 )

link1 = 0.934632 ( 0.376352 )
link2 = 1.265497 ( 0.370829 )
link3 = 0.857588 ( 0.333646 )
link4 = 0.798181 ( 0.219938 )
link5 = 0.443085 ( 0.315726 )
link6 = 0.896836 ( 0.309806 )
link7 = 0.611814 ( 0.30125 )
link8 = 0.675294 ( 0.199431 )
link9 = -0.491547 ( 0.360964 )
link10 = 0.195862 ( 0.292776 )
link11 = 0.144498 ( 0.281587 )
link12 = 0.441636 ( 0.180349 )

> X <- c(0,0,-4,4,1,0,-3,3,0,1,-2,2,1,0,-1,1,0,0,-4,0,1,0,-3,0,
+
+ 0,1,-2,0,1,0,-1,0,0,0,0,-4,1,0,0,-3,0,1,0,-2,1,0,0,-1)

+ + 0,1,-2,0,1,0,-1,0,0,0,0,-4,1,0,0,-3,0,1,0,-2,1,0,0,-1)
> X <- matrix(X,12,4,byrow=T)
> X
[,1] [,2] [,3] [,4]
[1,] 0 0 -4 4
[2,] 1 0 -3 3
[3,] 0 1 -2 2
[4,] 1 0 -1 1
[5,] 0 0 -4 0
[6,] 1 0 -3 0
[7,] 0 1 -2 0
[8,] 1 0 -1 0
[9,] 0 0 0 -4
[10,] 1 0 0 -3
[11,] 0 1 0 -2
[12,] 1 0 0 -1
> mph.out <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 19.25573 norm.score= 0.5606361

iter= 100 norm.diff= 9.236443 norm.score= 3.685312e-10

Time Elapsed: 3.782 seconds

TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 8 ): Gsq = 2.4871 (p = 0.96233 )
Pearson's Score Stat (df= 8 ): Xsq = 2.46806 (p = 0.9632 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 0.564523 ( 0.174241 )
beta2 = 0.390272 ( 0.238318 )
beta3 = -0.110771 ( 0.078931 )
beta4 = 0.122887 ( 0.090241 )

link1 = 0.934632 ( 0.376352 )
link2 = 1.265497 ( 0.370829 )
link3 = 0.857588 ( 0.333646 )
link4 = 0.798181 ( 0.219938 )
link5 = 0.443085 ( 0.315726 )
link6 = 0.896836 ( 0.309806 )
link7 = 0.611814 ( 0.30125 )
link8 = 0.675294 ( 0.199431 )
link9 = -0.491547 ( 0.360964 )
link10 = 0.195862 ( 0.292776 )
link11 = 0.144498 ( 0.281587 )
link12 = 0.441636 ( 0.180349 )

> X2 <- c(0,0,1,0,00,1,1,0,0,0,1,0,0,1,1,0,0,0,1,0,0,1,1,0)
> X2 <- matrix(X2,12,2,byrow=T)
> X2
[,1] [,2]
[1,] 0 0
[2,] 1 0
[3,] 0 1
[4,] 1 0
[5,] 0 0
[6,] 1 0
[7,] 0 1
[8,] 1 0
[9,] 0 0
[10,] 1 0
[11,] 0 1
[12,] 1 0
> mph.out2 <- mph.fit(y=y,Z=Z,L.fct=L.fct,X=X2)

mph.fit, version 1.0, 8/19/01 , running...
iter= 1 norm.diff= 10.52249 norm.score= 2.735496

iter= 100 norm.diff= 9.31842 norm.score= 0.0003783386

Time Elapsed: 3.835 seconds
TEST of Ho: h(m)=0 vs. Ha: not Ho...
Likelihood Ratio Stat (df= 10 ): Gsq = 8.8734 (p = 0.54416 )
Pearson's Score Stat (df= 10 ): Xsq = 7.57836 (p = 0.66995 )

LINEAR PREDICTOR MODEL RESULTS...
beta1 = 0.502986 ( 0.174951 )
beta2 = 0.299645 ( 0.232652 )

link1 = -1e-06 ( 0 )
link2 = 0.502986 ( 0.174951 )
link3 = 0.299644 ( 0.232652 )
link4 = 0.502986 ( 0.174951 )
link5 = 0 ( 0 )
link6 = 0.502986 ( 0.174951 )
link7 = 0.299645 ( 0.232652 )
link8 = 0.502986 ( 0.174951 )
link9 = 0 ( 0 )
link10 = 0.502987 ( 0.174951 )
link11 = 0.299645 ( 0.232652 )
link12 = 0.502986 ( 0.174951 )

This page is maintained by Alan Agresti. Last Updated: November 30, 2009.

(Our Department Home Page)