##
#####  Function to exponentiate coefficients and produces CIs
##
expfit <- function( model ){
	coef <- summary( model )$coef
	ci95.lo <- exp( coef[,1] - qnorm(.975) * coef[,2] )
	ci95.hi <- exp( coef[,1] + qnorm(.975) * coef[,2] )
	rslt <- round( cbind( exp( coef[,1] ), ci95.lo, ci95.hi, coef[,3:4] ), 4 )
	colnames( rslt )[1] <- "exp( Est )"
	rslt
	}
##
#####  Function to collapse grouped binary data to binomial
##
collapse <- function( data, outcome ){
	index <- (1:length(names(data)))[ names(data)==outcome ]
	y <- data[,index]
	covnames <- names( data )[-index]
	data <- data[,-index]
	if( is.null( dim( data ) ) ){
		rslt <- aggregate( y, list(data), FUN=length)
		rslt <- as.data.frame( cbind( rslt, aggregate( y, list(data), FUN=sum)[dim(rslt)[2]] ) )
	}
	else{
		rslt <- aggregate( y, data, FUN=length)
		rslt <- as.data.frame( cbind( rslt, aggregate( y, data, FUN=sum)[dim(rslt)[2]] ) )
	}
	names( rslt ) <- c( covnames, "n", paste("n.", outcome, sep="") )
	rslt
}
##
#####	Function to compute deviance (LR) test p-Value
##
lrtest <- function( fit1, fit2 ){
	cat( "\nAssumption: Model 1 nested within Model 2\n\n" )
	rslt <- anova( fit1, fit2 )
	rslt <- cbind( rslt, c("", round( pchisq( rslt[2,4], rslt[2,3], lower.tail=FALSE ), 4 ) ) )
	rslt[,2] <- round( rslt[,2], 3 )
	rslt[,4] <- round( rslt[,4], 3 )
	rslt[1,3:4] <- c( "", "" )
	names( rslt )[5] <- "pValue"
	rslt
}
##
#####	H-L goodness of fit test
##
binary.gof <- function( fit, ngrp=10, print.table=TRUE ){
	y <- fit$y
	phat <- fitted( fit )
	fittedgrps <- cut( phat, quantile( phat, seq(0,1,by=1/ngrp) ), include.lowest=TRUE )
	n <- aggregate( y, list( fittedgrps ), FUN=length )[,2]
	Obs <- aggregate( y, list( fittedgrps ), FUN=sum )[,2]
	Exp <- aggregate( phat, list( fittedgrps ), FUN=sum )[,2]
	if( print.table==TRUE ){
		cat( "\nFitted Probability Table:\n\n" )
		rslt <- as.data.frame( cbind( 1:ngrp, n, Obs, Exp ) )
		names( rslt )[1] <- "group"
		print( rslt )
	}
	chisqstat <- sum( (Obs - Exp)^2 / Exp )
	df <- ngrp - length( fit$coef )
	pVal <- pchisq( chisqstat, df, lower.tail=FALSE )
	cat( "\n Hosmer-Lemeshow GOF Test:\n\n" )
	cbind( chisqstat, df, pVal )
}
##
#####
#####	robust.se() is a function to compute the Huber-White sandwich variance estimator
#####	for the linear regression model
#####	
##
robust.se.lm <- function( model) { 
  s <- summary( model) 
  X <- model.matrix( model )
  eps <- model$residuals
  sandwich.cov <- solve( t(X)%*%X ) %*%( t(X) %*% diag(eps^2) %*% X ) %*% solve( t(X)%*%X )
  sand.se <- sqrt( diag( sandwich.cov )) 
  t <- model$coefficients/sand.se
  p <- 2*pt( -abs( t ), dim(X)[1]-dim(X)[2] ) 
  ci95.lo <- model$coefficients - qt( .975, dim(X)[1]-dim(X)[2] ) * sand.se
  ci95.hi <- model$coefficients + qt( .975, dim(X)[1]-dim(X)[2] ) * sand.se
  rslt <- cbind( model$coefficients, sand.se, t, p, ci95.lo, ci95.hi ) 
  dimnames(rslt)[[2]] <- c( dimnames( s$coefficients )[[2]], "ci95.lo", "ci95.hi" ) 
  dimnames(rslt)[[2]][2] <- "Robust SE"
  rslt 
} 
##
#####	Function to compute robust se for glms
##
robust.se.glm<-function(glm.obj){
	## 	Compute robust (sandwich) variance estimate
	if (is.matrix(glm.obj$x)) 
		xmat<-glm.obj$x
	else {
		mf<-model.frame(glm.obj)
		xmat<-model.matrix(terms(glm.obj),mf)		
	}
	umat <- residuals(glm.obj,"working")*glm.obj$weights*xmat
	modelv<-summary(glm.obj)$cov.unscaled
	robust.cov <- modelv%*%(t(umat)%*%umat)%*%modelv
	
	##	Format the model output with p-values and CIs
	s <- summary( glm.obj) 
	robust.se <- sqrt( diag( robust.cov )) 
	z <- glm.obj$coefficients/robust.se
	p <- 2*pnorm( -abs( z ) ) 
	ci95.lo <- glm.obj$coefficients - qnorm( .975 ) * robust.se
	ci95.hi <- glm.obj$coefficients + qnorm( .975 ) * robust.se
	rslt <- cbind( glm.obj$coefficients, robust.se, z, p, ci95.lo, ci95.hi ) 
	dimnames(rslt)[[2]] <- c( dimnames( s$coefficients )[[2]], "ci95.lo", "ci95.hi" ) 
	dimnames(rslt)[[2]][2] <- "Robust SE"
	rslt 
	}

##
#####	Function to summarize the multinomial fit
##
summ.mfit <- function( model ){
	s <- summary( model )
	for( i in 1:length(model$coef) ){
		cat( "\nLevel ", model$lev[i+1],  "vs. Level ", model$lev[1], "\n" )
		coef <- s$coefficients[i,]
		rrr <- exp( coef )
		se <- s$standard.errors[i,]
		zStat <- coef / se
		pVal <- 2*pnorm( abs(zStat), lower.tail=FALSE )
		ci95.lo <- exp( coef - qnorm(.975)*se )
		ci95.hi <- exp( coef + qnorm(.975)*se )
		rslt <- cbind( rrr, se, zStat, pVal, ci95.lo, ci95.hi )
		print( round( rslt, 3 ) )
	}
}

##
#####	Function to test linear contrasts
##
LinContr.mfit <- function( contr.names, contr.coef, model ){ 
	beta.hat <- as.vector( t( summary( model )$coefficients ) )
	se <- as.vector( t( summary( model )$standard.errors ) )
	corr <- summary( model )$correlation
	cov.beta <- diag( se ) %*% corr %*% diag( se )

	contr.index <- is.element( dimnames( corr )[[1]], contr.names )
	beta.hat <- beta.hat[ contr.index ]
	cov.beta <- cov.beta[ contr.index,contr.index ]
	est <- contr.coef %*% beta.hat
	rrr.est <- exp( est )
	se.est <- sqrt( contr.coef %*% cov.beta %*% contr.coef )
	zStat <- est / se.est
	pVal <- 2*pnorm( abs(zStat), lower.tail=FALSE )
	ci95.lo <- exp( est - qnorm(.975)*se.est )
	ci95.hi <- exp( est + qnorm(.975)*se.est )
	cat( "\nTest of H_0: " )
	for( i in 1:(length( contr.names )-1) ){
		cat( contr.coef[i], "*", contr.names[i], "+ " )
	}
	cat( contr.coef[i+1], "*", contr.names[i+1], "= 0 :\n\n" )
	rslt <- data.frame( rrr.est, se.est, zStat, pVal, ci95.lo, ci95.hi )
	round( rslt, 3 )
}

##
#####	Function to summarize the proportional odds fit
##
exp.olr <- function( model ){
		coef <- summary( model )$coef[ 1:length(model$coef), ]
		zStat <- coef[,3]
		pVal <- 2*pnorm( abs(zStat),lower.tail=FALSE )
		ci95.lo <- exp( coef[,1] - qnorm(.975) * coef[,2] )
		ci95.hi <- exp( coef[,1] + qnorm(.975) * coef[,2] )
		rslt <- round( cbind( exp( coef[,1] ), zStat, pVal, ci95.lo, ci95.hi), 4 )
		colnames( rslt )[1] <- "exp( Est )"
		rslt
		}