## This program explores a problem encountered while doing simulations
## it shows the existance of situations where a diagonal element of omega_i %*% vm is greater than 1


### include some programs
library(yags)
### input modifications to yags, see Splus Software section on web site
source("myags.ssc")
source("myags.pmat.fit.ssc")
## read in original data
## see http://www.phs.wfubmc.edu/data/uipreiss.html
pr<-read.table("preisser.txt",header=T)
pr[1,]
y<-pr$bothered
n<-length(y)
female<-pr$female
age<-pr$age
dayacc<-pr$dayacc
severe<-pr$severe
toilet<-pr$toilet
pract.id<-pr$pract.id
x<-matrix(c(rep(1,n),female,age,dayacc,severe,toilet),n,6)

## define expit
expit<-function(eta){ exp(eta)/(1+exp(eta)) }


## fit the model without ith column of design matrix
beta3<-yags(y~female+dayacc+severe+toilet,id=pract.id,family=binomial,corstr="independence")$coef
fit3<-expit(x[,-3] %*% matrix(beta3,5,1) )

### simulation data from that fit
.Random.seed<-c(3,54,26,16,37,2,45,49,4,37,0,3)
ysim<-rbinom(n,rep(1,n),fit3)
### now do call myags

mout<-myags(ysim~female+age+ dayacc+severe+toilet,id=pract.id,family=binomial,corstr="exchangeable")

### show that for the 18th cluster we can get a value of omega- %*% vm > 1

o18<-mout$omega[18,,]

vminv<-apply(mout$omega,c(2,3),sum)
vm<-solve(vminv)
round(vm,4)
round(vminv,4)
round(o18,4)

diag( o18 %*% vm )





### is there anything special about the 18th cluster?
### I couldn't see anything
dat1<-data.frame(yorig=y,ysim=ysim,female=female,age=age,dayacc=dayacc,
severe=severe,toilet=toilet,pract.id=pract.id)


id18<-unique(pract.id)[18]
dat1[pract.id==id18,]