CalcMDI <- function(DAT,S='site',R='Dx',n.boot=100,print.i=FALSE,seed=1,rd=3,ci=c(0.025,0.975)){
              # Create temporary data set, in case manipulations are needed.
              tmp <- DAT  

              # Identify columns of interest
              SiteColumn       <- which(names(tmp)==S)
              ResponseColumns  <- which(unlist(lapply(strsplit(names(tmp),R),length))>1)
              CovariateColumns <- which(1:ncol(tmp)%in%c(SiteColumn,ResponseColumns)==FALSE)

              if(!is.factor(tmp[,SiteColumn])) tmp[,SiteColumn] <- factor(tmp[,SiteColumn])
              s.lv <- levels(tmp[,SiteColumn])
              s.nc <- nchar(names(tmp)[SiteColumn])
              
              # Define all 2-way combinations of S
              bi.comp <- matrix(paste(names(tmp)[SiteColumn],combn(s.lv,2),sep=''),nrow=2)
              bi.name <- apply(bi.comp,2,paste,collapse=' : ')
              
              # Generate bootrapped sample of marginal coefficients              
              Coef <- list()
              for (b in 1:ncol(bi.comp)){ Coef[[b]] <- matrix(NA,ncol=length(ResponseColumns),nrow=n.boot)}
              
              set.seed(seed)
              for (l in 1:n.boot){ if(print.i) {print(l)}
                  rows     <- sample(c(1:nrow(tmp)), nrow(tmp), replace=TRUE)
                  Dat.boot <- tmp[rows,]
              
                  CoefNew <- NULL
                  for (k in 1:length(ResponseColumns)){
                   
                      Dat.boot$Y <- Dat.boot[,ResponseColumns[k]]    
                      mod <- glm( as.formula(paste( "Y~", paste(names(tmp[,c(CovariateColumns, SiteColumn)]), collapse=" + "))), 
                                  data=Dat.boot, family=binomial)
               
                      CoefNew  <- cbind(CoefNew, mod$coef[substr(names(mod$coef), 1, s.nc) == names(DAT)[SiteColumn]]) }
               
               # Add reference line for linear combinations
               CoefNew <- rbind(0,CoefNew)
               rownames(CoefNew)[1] <- paste(names(tmp)[SiteColumn],s.lv,sep='')[1]
    
               for (b in 1:ncol(bi.comp)) { Coef[[b]][l,1:ncol(CoefNew)] <- diff(CoefNew[match(bi.comp[,b],rownames(CoefNew)),]) } }

               ## Stage 1 Results
               est <- NULL
               for (i in 1:length(ResponseColumns)) {
                   tmp$Y <- tmp[,ResponseColumns[i]]
                   est <- cbind(est, glm( as.formula(paste( "Y~", paste(names(tmp[,c(CovariateColumns, SiteColumn)]), collapse=" + "))), 
                                          data=tmp, family=binomial)$coef[substr(names(mod$coef), 1, s.nc) == names(tmp)[SiteColumn]])}
               # Add reference line for linear combinations
               est <- rbind(0,est)
               rownames(est)[1] <- paste(names(tmp)[SiteColumn],s.lv,sep='')[1]

               est.lc <- list()
               for (b in 1:ncol(bi.comp)) {  
                   est.lc[[b]] <- diff(est[match(bi.comp[,b],rownames(est)),]) }

               covs <- lapply(Coef,cov)

               ## Stage 2 Results
               outtab <- matrix('---',ncol=length(s.lv),nrow=length(s.lv))
               MDI    <- NULL
               for (i in 1:length(est.lc)) {
                    inv.cov <- solve(covs[[i]])
                    tr.inv.cov <- sum(diag(inv.cov))
                    MDI <- c(MDI, sqrt((1/tr.inv.cov)*est.lc[[i]]%*%inv.cov%*%t(est.lc[[i]]))) }
               MDI <- format(round(MDI,rd),nsmall=rd)

               outtab[lower.tri(outtab)] <- MDI
               outtab[upper.tri(outtab)] <- MDI
               outtab <- data.frame(outtab)
               colnames(outtab) <- rownames(outtab) <- s.lv
 
               bi.est <- matrix(unlist(est.lc),ncol=length(bi.name))
               bi.lci  <- matrix(unlist(lapply(Coef,function(Z){apply(Z,2,quantile,ci[1])})),ncol=length(bi.name))
               bi.uci  <- matrix(unlist(lapply(Coef,function(Z){apply(Z,2,quantile,ci[2])})),ncol=length(bi.name))

               colnames(bi.est) <- colnames(bi.lci) <- colnames(bi.uci) <- bi.name
               rownames(bi.est) <- rownames(bi.lci) <- rownames(bi.uci) <- names(tmp)[ResponseColumns]

               # Output object
               out <- list(MDI=outtab,est= bi.est,lci=bi.lci,uci=bi.uci,vcov=covs)
               out}