gov = read.table("http://sites.williams.edu/bklingen/files/2013/12/govsuc.dat", header=TRUE)
require("reshape2") #requires installation of reshape2 package
gov.long=melt(gov)
colnames(gov.long)=c("health", "envi", "count")
gov.long
#matrix for forming marginal cumulative probabilities
# P(Y_t<=j) and 1-P(Y_t<=j), j=1,2, t=1,2
A <- matrix(c(
1,0,0,1,0,0,1,0,0,
0,1,1,0,1,1,0,1,1,
1,1,0,1,1,0,1,1,0,
0,0,1,0,0,1,0,0,1,
1,1,1,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,0,0,0,1,1,1),
byrow=TRUE,ncol=9)
# matrix for forming cumulative logits logit(P(Y_t<=j)):
C <- matrix(c(
1,-1,0,0,0,0,0,0,
0,0,1,-1,0,0,0,0,
0,0,0,0,1,-1,0,0,
0,0,0,0,0,0,1,-1),
byrow=TRUE,ncol=8)

C%*%log(A%*%gov.long$count) #sample cumulative logits

#Design matrix for parameters
# alpha1, alpha2, beta"
X <- matrix(c(
1,0,0,
0,1,0,
1,0,1,
0,1,1),
byrow=TRUE,ncol=3)

#function that takens in multinomial probabilities and turns them into cumulative logits
L.fct <- function(p) C%*%log(A%*%p)
L.fct(gov.long$count/sum(gov.long$count)) #sample cumulative logits
# now, use mph.fit:

source("mph_fit.R")  #have to request function from Prof. Lang at Univ. of Iowa
fit.marginal.ML <- mph.fit(y=gov.long$count,L.fct=L.fct,X=X)
fit.marginal.ML$beta   #ML parameter estimates
fit.marginal.ML$covbeta  #covariance matrix
fit.marginal.ML$Xsq  #Pearson X2 GoF statistic
fit.marginal.ML$df  # degrees offreedom for GoF statistic