###If model0 is a logistic regression submodel of model, then
###LRtest compares them with a likelihood ratio test and returns
###the p-value.


LRtest=function(model0,model){
  1-pchisq(model0$deviance-model$deviance,
           df=model0$df.residual-model$df.residual)}


###logit calculates the logit of a number p between 0 and 1.


logit=function(p){
  log(p/(1-p))}


###Given a vector of numerical values x, a response vector Y, and a 
###number of groups k, groupplot returns a list.
###x=mean of x groups
###Pi=mean of Y values in each group (with Wilson's adjustment)
###g=logit of Pi


groupplot=function(x,Y,k){
  sortframe=sort(x,index=TRUE)
  x=sortframe$x
  Y=Y[sortframe$ix]
  xmeans=1:k
  Pi=1:k
  s=floor(length(x)/k)  #groupsize
  
  for(i in 1:k){
    index=(((i-1)*s+1):(i*s))
    xmeans[i]=mean(x[index])
    Pi[i]=(sum(Y[index])+2)/(s+4)}
  
  g=logit(Pi)
  return(list(x=xmeans,g=g,Pi=Pi))}