.packageName <- "Matchit"
diagnose <-function(formula, match.matrix, pscore, in.sample, data,
                    exact = FALSE, mahvars = NULL, subclass = 0,
                    psclass = NULL, nearest = TRUE, q.cut = NULL,
                    counter = TRUE){
  if(counter){cat("Calculating summary statistics...")}
  treata <- model.frame(formula,data)[,1,drop=FALSE]
  treat <- as.vector(treata[,1])
  names(treat) <- row.names(treata)
  n <- length(treat)
  n0 <- length(treat[treat==0])
  n1 <- length(treat[treat==1])

  if(is.null(names(treat))){names(treat) <- seq(1,n)}
  labels <- names(treat)
  clabels <- labels[treat==0]
  tlabels <- labels[treat==1]

  # setting up weights
  if(nearest==TRUE){
    match.matrix <- match.matrix[in.sample[treat==1],,drop=F]
    num.matches <- dim(match.matrix)[2]-apply(as.matrix(match.matrix), 1, function(x) { sum(is.na(x)) })
    names(num.matches) <- tlabels[in.sample[treat==1]]

    t.units <- row.names(match.matrix)[num.matches>0]
    c.units <- na.omit(as.vector(as.matrix(match.matrix)))
    matched <-c(t.units,unique(c.units))
    
    weights <- rep(0,length(treat))
    names(weights) <- labels
    weights[t.units] <- 1
    
    for (cont in clabels) {
      treats <- na.omit(row.names(match.matrix)[cont==match.matrix[,1]])
      if (dim(match.matrix)[2]>1) {
        for (j in 2:dim(match.matrix)[2]) 
          treats <- c(na.omit(row.names(match.matrix)[cont==match.matrix[,j]]),treats)	
      }
      weights[cont] <- length(unique(treats))
    }
    
    if (sum(weights[clabels])==0){
      weights[clabels] <- rep(0, length(weights[clabels]))
    } else {
      weights[clabels] <-length(weights[clabels][weights[clabels]>0])*weights[clabels]/sum(weights[clabels][weights[clabels]>0])
    }
    
  } else if(identical(exact,TRUE)) {
    
    weights <- rep(0,length(treat))
    names(weights) <- labels
    t.units <- names(psclass)[psclass>0 & treat==1]
    weights[t.units] <- 1
    for(j in 1:max(psclass)){
        qn0 <- sum(treat==0 & psclass==j)
        qn1 <- sum(treat==1 & psclass==j)
        weights[treat==0 & psclass==j] <- qn1/qn0
      }
    if (sum(weights[treat==0])==0) {
      weights[treat==0] <- rep(0, length(weights[clabels]))
    } else {
      weights[clabels] <- length(weights[treat==0][weights[treat==0]>0]) *
        weights[treat==0]/sum(weights[treat==0][weights[treat==0]>0])
    }
    matched <- names(treat)
    
  } else {
    matched <- names(treat)
    weights <- rep(1,length(treat))
    names(weights) <- names(treat)
  }
  weights[!in.sample] <- 0
  if (sum(weights)==0) {
    stop("Error: No units were matched")
  } else if (sum(weights[tlabels])==0){
    stop("No treated units were matched",call.=FALSE)
  } else if (sum(weights[clabels])==0){
    stop("No control units were matched",call.=FALSE)
  }
  cat("Done\n")
  list(weights=weights)
}

distance <- function(formula, model="logit", data,
                     discard=0,reestimate=FALSE,counter=TRUE,...){
  if(counter){cat("Calculating propensity score...")}
  mf<-match.call()
  if(!is.vector(discard) | !identical(round(length(discard)),1)){
    warning("discard=", mf$discard,
            " is invalid; used discard=0 instead", call.=FALSE)
    discard <- 0}
  if(!(identical(reestimate,TRUE) | identical(reestimate,FALSE))){
    warning("reestimate=",
            mf$reestimate," is invalid; used reestimate=FALSE instead",
            call.=FALSE)
    mf$reestimate=FALSE}
  mf$counter <- mf$discard<-mf$model<-mf$reestimate<-NULL
  dotsub <- eval(mf$subset, sys.frame(sys.parent()))
  a <- 0 
  while(a<2){
    if(a==1)
      if(is.null(dotsub))
        mf$subset<-in.sample
      else{
        dotsub[dotsub] <- in.sample
        mf$subset <- dotsub
      }
    if(model=="logit" || model=="probit"){
      mf[[1]]<-as.name("glm")
      if(model=="logit"){mf$family<-as.name("binomial")} else
      {mf$family <- binomial(link="probit")}  #is there a better way to define the link function here?
      res<-eval(as.call(mf), sys.frame(sys.parent()))
      if(a==0){
        pscore<-fitted(res)
        treat<-res$y
        covariates <- res$model[,-1,drop=F]
      }
      else{
        pscore[!in.sample] <- NA
        pscore[in.sample] <- fitted(res)
      }
      assign.model<-summary(res)
    }
    else if(model=="nnet"){
      require(nnet)
      mf[[1]] <- as.name("nnet.formula")
      mf$size <- 3  #is this the right default setting?
      res <- eval(as.call(mf),sys.frame(sys.parent()))
      if(a==0){
        pscore <- as.vector(fitted(res))
        nnames <- dimnames(fitted(res))[[1]]
        names(pscore) <- nnames
        if(is.null(mf$data)){
          treat <- model.frame(formula(res))[,1]
          covariates <- model.frame(delete.response(terms(formula(res))),drop=F)
        }else{
          treat <- model.frame(formula(res),data)[,1]
          covariates <- model.frame(delete.response(terms(formula(res))),data,drop=F)
        }
        if(!is.null(dotsub)){
          treat <- treat[dotsub]
          covariates <- covariates[dotsub,,drop=F]
        }
        names(treat) <- nnames
        assign.model <- summary(res)
      }
      else{
        pscore[!in.sample] <- NA
        pscore[in.sample] <- as.vector(fitted(res))
        assign.model <- res
        dan <- mf
      }
    }    
    else if(model=="GAM"){
      require(mgcv)
      mf[[1]] <- as.name("gam")
      mf$family<-as.name("binomial")
      res <- eval(as.call(mf),sys.frame(sys.parent()))
      if(a==0){
        treat <- res$y
        names(treat) <- row.names(data)  #assuming we have a data frame
        pscore <- fitted(res)
        names(pscore) <- names(treat)
        covariates <- res$model[,2:ncol(res$model),drop=F]
      }
      else{
        pscore[!in.sample] <- NA
        pscore[in.sample] <- fitted(res)
      }
      assign.model <- summary(res)
    }
    else if(model=="cart"){
      require(rpart)
      mf[[1]] <- as.name("rpart")
      res <- eval(as.call(mf),sys.frame(sys.parent()))
      if(a==0){
        pscore <- predict(res)
        treat <- res$y
        if(!is.null(dotsub)){
          covariates <- data[dotsub,attr(delete.response(terms(formula)),"term.labels"),drop=F]
        } else {covariates <- data[,attr(delete.response(terms(formula)),"term.labels"),drop=F]}
      }
      else{
        pscore[!in.sample] <- NA
        pscore[in.sample] <- predict(res)
      }
      assign.model <- print(res)
    }
    else {stop("Must specify valid model to estimate propensity score",call.=FALSE)}
    if(a==0){
      maxp0 <- max(pscore[treat==0])
      minp0 <- min(pscore[treat==0])
      maxp1 <- max(pscore[treat==1])
      minp1 <- min(pscore[treat==1])}
    if(discard==0){
      in.sample <- !is.na(treat) #just to make all TRUE
      a <- 2}
    else if(discard==1){
      if(a==0){in.sample <- (pscore>=max(minp1,minp0) & pscore<=min(maxp0,maxp1))}
      if(reestimate){(a <- a+1)}
      else {a <- 2}
    }
    else if(discard==2){
      if(a==0){in.sample <- (pscore>=minp1 & pscore<=maxp1)}
      if(reestimate){a <- a+1}
      else {a <- 2}
    }
    else if(discard==3){
      if(a==0){in.sample <- (pscore>=minp0 & pscore<=maxp0)}
      if(reestimate){a <- a+1}
      else {a <- 2}
    }
    else(stop("Invalid discard option",call.=FALSE))
  }
  ## return
  if(counter){cat("Done\n")}
  list(in.sample=in.sample, pscore=pscore, treat=treat, covariates=covariates,
                     assign.model=assign.model)
}
exactmatch <-  function(formula,data,counter=TRUE) {
  if(counter){cat("Exact matching...")}
  data <- eval(data,parent.frame())
  treata <- model.frame(formula,data)[,1,drop=FALSE]
  treat <- as.vector(treata[,1])
  names(treat) <- row.names(treata)
  covariates <- model.frame(delete.response(terms(formula)),data)[,,drop=FALSE]
  n <- length(treat)
  n1 <- length(treat[treat==1])
  n0 <- length(treat[treat==0])
  covariates <- as.data.frame(covariates)
  xx <- apply(covariates, 1, function(x) paste(x, collapse = "\r"))
  xx1 <- xx[treat==1]
  xx0 <- xx[treat==0]
  cc <- unique(xx1)
  cc <- cc[cc%in%xx0]
  ncc <- length(cc)
  psclass <- rep(0,n)
  names(psclass) <- names(treat)
  for(i in 1:ncc){
    psclass[xx==cc[i]] <- i
  }
  if(counter){cat("Done\n")}
  list(psclass = psclass)
}
help.matchit <- function (object=NULL) 
{
    under.unix <- !(version$os == "Microsoft Windows" || version$os == 
        "Win32" || version$os == "mingw32")
    sys <- function(command, text = NULL) {
        cmd <- if (length(text)) 
            paste(command, text)
        else command
        if (under.unix) 
            system(cmd)
        else shell(cmd, wait = TRUE)
    }
    browser <- .Options$help.browser
    if (!length(browser)) 
        browser <- .Options$browser
    if (!length(browser)) 
        browser <- getOption("browser")
    url <- NULL
 
    if (is.null(object))
      url <- c("http://gking.harvard.edu/matchit")
    if (!is.null(object)) {
    if (object == "matchit")
        url <- c("http://gking.harvard.edu/matchit/docs/Usage.html")
    if (object == "distance")
      url <- c("http://gking.harvard.edu/matchit/docs/_TT_distance_TT.html")
    if (object == "matchdef")
      url <- c("http://gking.harvard.edu/matchit/docs/_TT_matchdef_TT.html")
    if (object == "subclassify")
      url <- c("http://gking.harvard.edu/matchit/docs/_TT_subclassify_TT.html")
    if (object == "print.matchit")
      url <- c("http://gking.harvard.edu/matchit/docs/Print.html")
    if (object == "summary.matchit")
      url <- c("http://gking.harvard.edu/matchit/docs/Summary.html")
    if (object == "plot.matchit")
      url <- c("http://gking.harvard.edu/matchit/docs/Plot.html")
    }

    if (is.null(url)) {
        cat("Error:", object, "currently not documented in help.matchit. \n Please check http://gking.harvard.edu/matchit. \n", 
            sep = " ")
	url <- c("http://gking.harvard.edu/matchit")
    }

    if (under.unix) {
        sys(paste(browser, url, "&"))
        invisible()
    }
    if (!under.unix) {
        browseURL(url, browser = browser)
        invisible("")
    }
}
match.data <- function(object, group = "all") {
  data <- object$data
  treat <- eval(object$treat, data)
  if (group == "all")
    return(data[object$matched,])
  else if (group == "treat")
    return(data[object$matched & treat == 1,])
  else if (group == "control")
    return(data[object$matched & treat == 0,])
  else
    stop("error: invalid input for group.")
}
matchdef <-  function(formula,in.sample,pscore,nearest=TRUE,replace=FALSE,
                      m.order=2,ratio=1,caliper=0,calclosest=FALSE,mahvars=NULL,exact=FALSE,
                      data=NULL, counter=TRUE){
  treata <- model.frame(formula,data)[,1,drop=FALSE]
  treat <- as.vector(treata[,1])
  names(treat) <- row.names(treata)
  
  covariates <- model.frame(delete.response(terms(formula)),data)[,,drop=FALSE]

                                        # Total number of units
  n <- length(treat)
  n1 <- length(treat[treat==1])
  n0 <- length(treat[treat==0])
  
  p1 <- pscore[treat==1]
  p0 <- pscore[treat==0]
  
  if(nearest==TRUE)
    {
      
                                        # Generating separate indices for treated and control units
      if(is.null(names(treat))){names(treat) <- seq(1,n)}
      labels <- names(treat)
      clabels <- labels[treat==0]
      tlabels <- labels[treat==1]
      
                                        # Generating match matrix
      match.matrix <- matrix(0,n1,ratio)
      match.matrix <- as.data.frame(match.matrix)
      row.names(match.matrix) <- tlabels
      
                                        # Vectors of whether unit has been matched: = 0 if not matched (unit # of match if matched)
                                        # = -1 if can't be matched (if in.sample=0)
      matchedc <- rep(0,length(p0))
      names(matchedc) <- clabels
      
                                        # These are the units that are ineligible because of discard (in.sample==0)
      matchedc[in.sample[clabels]==0] <- -1
      match.matrix[in.sample[tlabels]==0,] <- -1
      matchedt <- match.matrix[,1] 
      
      names(matchedt) <- tlabels
                                        # total number of matches (including ratios) = ratio * n1
      tr <- length(match.matrix[match.matrix!=-1])
      r <- 1
      
                                        # Caliper for matching (=0 if caliper matching not done)
      sd.cal <- caliper*sqrt(var(pscore[in.sample==1]))
      
					# Var-covar matrix for Mahalanobis (currently set for full sample)      
      if (!is.null(mahvars)) {
        if(!sum(mahvars%in%names(data))==length(mahvars)){stop("Mahvars not contained in data",call.=FALSE)}
        mahvars <- as.matrix(data[,mahvars])
        row.names(mahvars) <- labels
        Sigma <- var(mahvars)
      }   

                                        # Now for exact matching
                                        # If exact=T, set exact to covariates from distance
                                        # Make new variable that is T/F for exact matching, then exact will have the vars if
                                        # exactmatch=T
      exactmatch <- exact
      if (!is.logical(exact)) {
        exactmatch <- TRUE} else if (identical(exact,TRUE)) {exact <- covariates}
      
      if (exactmatch==TRUE){
        if(!sum(exact%in%names(data))==length(exact)){stop("Exact variables not contained in data",call.=FALSE)}
        exact <- as.matrix(data[,exact])
        row.names(exact) <- labels
      }   
                                        # Looping through nearest neighbour matching for all treatment units
                                        # Only do matching for units with in.sample==1 (matched!=-1)
      if(counter==TRUE){
        trseq <- floor(seq(tr/10,tr,tr/10))
        cat("Matching Treated: ")
      }
      
      for(i in 1:tr){
					# Make new matchedc column to be used for exact matching
					# Will only be 0 (eligible for matching) if it's an exact match
        if(counter==TRUE) {if(i%in%trseq){cat(10*which(trseq==i),"%...",sep="")}}  # a counter
	matchedc2 <- matchedc
        if(!0%in%matchedc2){  #in cases there's no replacement and all controls have been used up
          match.matrix[match.matrix[,r]==0 & !is.na(match.matrix[,r]),r] <- NA
          if(r<ratio){match.matrix[,(r+1):ratio] <- NA}
          break}

                                        #in case there's replacement, but all units have been used in
                                        #previous ratios
        if(sum(!is.na(match.matrix[,r]))==0){
          if(r<ratio){match.matrix[,(r+1):ratio] <- NA}
          break}
        
                                        # Which ratio we're on
        if(r!=ceiling(i/(tr/ratio))) {r <- r+1; matchedt <- match.matrix[,r]}
        
        if(m.order==2) {iterp1 <- max(p1[matchedt==0],na.rm=T)}
        if(m.order==3) {iterp1 <- min(p1[matchedt==0],na.rm=T)}
        if(m.order==4) {iterp1 <- sample(p1[matchedt==0][!is.na(p1[matchedt==0])],1)}
        
                                        # The treatment unit for this iteration, again resolving ties randomly
        itert <- as.vector(na.omit(tlabels[iterp1==p1 & matchedt==0]))
        if(length(itert)>1){itert <- sample(itert,1)}
        
                                        # Calculating all the absolute deviations in propensity scores
                                        # Calculate only for those eligible to be matched (matchedc==0)
                                        # this first if statement only applies to replacement ratio
                                        # matching, so that each treatment unit is matched to a different
                                        # control unit than from the previous round
        
                                        # match number = NA if no units within caliper
        
					# Set things up for exact matching
					# Make matchedc2==-2 if it isn't an exact match
					# There might be a more efficient way to do this, but I couldn't figure
					# out another way to compare a vector with the matrix
	if (exactmatch==TRUE) {
          for (k in 1:dim(exact)[2]) matchedc2[exact[itert,k]!=exact[clabels,k]] <- -2
        }
        
	# Need to add a check in case there aren't any eligible matches left...
        if(replace==TRUE & r!=1) {
          if (sum(!clabels%in%match.matrix[itert,(1:r-1)] & matchedc2==0)==0) { 
            deviation <- NULL 
            mindev <- NA
          }
          else deviation <- abs(p0[!clabels%in%match.matrix[itert,(1:r-1)] & matchedc2==0]-iterp1)
	}	else { 
          if (sum(matchedc2==0)==0) { 
            deviation <- NULL
            mindev <- NA
          }
          else deviation <- abs(p0[matchedc2==0]-iterp1)
	}
        
        if (caliper!=0 & (!is.null(deviation))) {
          if(replace==TRUE & r!=1) pool <- clabels[!clabels%in%match.matrix[itert,(1:r-1)] & matchedc2==0][deviation 
                                                                                               <= sd.cal]
          else pool <- clabels[matchedc2==0][deviation <= sd.cal]
          
          if(length(pool)==0) { 
            if (calclosest==FALSE) mindev <- NA
            else { 
              if (replace==TRUE & r!= 1){ 
                mindev <- clabels[!clabels%in%match.matrix[itert,(1:r-1)]][min(deviation)==deviation]
              } else{mindev <- clabels[matchedc2==0][min(deviation)==deviation]}
	    }
          }
          else if (length(pool)==1) mindev <- pool[1]
          else if (is.null(mahvars)) mindev <- sample(pool, 1)
          else {
                                        # This has the important vars for the C's within the caliper
            poolvarsC <- mahvars[pool,,drop=F]
                                        # Sigma is the full group var/covar matrix of Mahalvars
            mahal <- mahalanobis(poolvarsC, mahvars[itert,],Sigma)
            mindev <- pool[mahal==min(mahal)]
          }
        }        else if(!is.null(deviation)) { 
          if (replace==TRUE & r!=1){ 
            mindev <- clabels[!clabels%in%match.matrix[itert,(1:r-1)] & matchedc2==0][min(deviation)==deviation]
          } else {mindev <- clabels[matchedc2==0][min(deviation)==deviation]}
        }
        
                                        # Resolving ties in minimum deviation by random draw
        if(length(mindev)>1){goodmatch <- sample(mindev,1)} else goodmatch <- mindev
        
                                        # Storing which treatment unit has been matched to control, and
                                        # vice versa
        matchedt[itert==tlabels] <- goodmatch
        matchedc[goodmatch==clabels] <- itert
        
                                        # instead of the in.sample, we now have an index with dimensions n1 by # of
                                        # matches (ratio)
        match.matrix[which(itert==tlabels),r] <- goodmatch
        
                                        # If matching with replacement, set matchedc back to 0 so it can be reused
        if (replace==TRUE) matchedc[goodmatch==clabels] <- 0
        
      }
    }
  if(counter==TRUE){cat("Done\n")}
  if(!nearest){match.matrix <- NULL}   else {
    x <- as.matrix(match.matrix)
    x[x==-1] <- NA
    match.matrix <- as.data.frame(x)}
  list(match.matrix = match.matrix, in.sample = in.sample)
}
matchit <- function(formula, model="logit", data, discard=0,
                    reestimate=FALSE, nearest=TRUE, replace=FALSE,
                    m.order=2, ratio=1, caliper=0, calclosest=FALSE,
                    subclass=0, sub.by="treat", mahvars=NULL, exact=FALSE,
                    counter=TRUE,...){ 
  cl <- match.call()
  #Checking input format
  #data input
  is.null(data)  #there must be a better way to get a warning for data input
  if(!is.data.frame(data)){
    stop("Data ", cl$data, " must be a dataframe",call.=FALSE)}
  #check pscore / psweights names in dataframe
  if("pscore"%in%names(data)){
    stop("Dataframe contains the variable 'pscore'.  Please change this name.",call.=FALSE)
  } 
  if("psweights"%in%names(data)){
    stop("Dataframe contains the variable 'psweights'.  Please change this name.",call.=FALSE)
  }
  if("psclass"%in%names(data)){
    stop("Dataframe contains the variable 'psclass'.  Please change this name.", call.=FALSE)
  }

  #sub.by
  if(!is.vector(sub.by)|length(sub.by)!=1){
    warning(cl$sub.by," is not a valid sub.by option; sub.by=\"treat\" used instead",call.=FALSE); sub.by <- "treat"}
  if(!sub.by%in%c("treat","control","all")){
    warning(cl$sub.by, " is not a valid sub.by option; sub.by=\"treat\" used instead",call.=FALSE); sub.by <- "treat"}
  #subclass
  if(length(subclass)==1){
    if(subclass<0 | subclass==1){
      warning(cl$subclass, " is not a valid subclass; subclass=0 used instead",call.=FALSE); subclass <- 0}
  } else {
    if(!is.vector(subclass)){
      warning(cl$subclass, " is not a valid subclass; subclass=0 used instead",call.=FALSE); subclass <- 0}
    if(sum(subclass<=1 & subclass>=0)!=length(subclass)){
      warning("Subclass ", cl$subclass, " is not bounded by 0 and 1; subclass=0 used instead",
              call.=FALSE); subclass <- 0}
  }
  #nearest
  if(!(identical(nearest,TRUE) | identical(nearest,FALSE))){
    warning("nearest=",cl$nearest," is invalid; used nearest=TRUE instead",call.=FALSE);nearest=TRUE}
  #replace
  if(!(identical(replace,TRUE) | identical(replace,FALSE))){
    warning("replace=",cl$replace," is invalid; used replace=FALSE instead",call.=FALSE);replace=0}
  #m.order
  if(!(identical(m.order,2) | identical(m.order,3) | identical(m.order,4))){
    warning("m.order=",cl$m.order," is invalid; used m.order=2 instead",call.=FALSE);m.order=2}
  #ratio
  ratio <- round(ratio)
  if(!is.numeric(ratio) | ratio[1]<1 | !identical(round(length(ratio)),1)){
    warning("ratio=",cl$ratio," is invalid; used ratio=1 instead",call.=FALSE);ratio=1}
  #caliper
  if(!is.vector(caliper) | !identical(round(length(caliper)),1)){
    warning("caliper=",cl$caliper," is invalid; used caliper=0 instead",call.=FALSE);caliper=0}
  if(caliper<0){
    warning("caliper=",cl$caliper," is less than 0; used caliper=0 instead",call.=FALSE);caliper=0}
  #calclosest
  if(!(identical(calclosest,TRUE)| identical(calclosest,FALSE))){
    warning("calclosest=",cl$calclosest," is invalid; used calclosest=FALSE instead",call.=FALSE)
    calclosest=FALSE}  
  #mahvars & caliper
  if (!is.null(mahvars) & caliper[1]==0){
    warning("Must specify caliper > 0 to use Mahalanobis matching. Mahalanobis matching not done",call. = FALSE)}
  #counter
  if(!(identical(counter,TRUE)| identical(counter,FALSE))){
    warning("counter=",cl$counter," is invalid; used counter=TRUE instead",call.=FALSE);counter=TRUE}    
  
  # Set up for output
  tt <- terms(formula(cl))

  if (identical(TRUE, exact)) {
    if(!is.null(c(cl$discard, cl$reestimate, cl$model, cl$nearest, cl$replace, cl$m.order, cl$ratio,
                  cl$caliper, cl$calclosest, cl$subclass, cl$sub.by, cl$mahvars)))
      warning("exact=TRUE chosen; all other matching options ignored")
    data <- eval(data,parent.frame())
    mf <- match.call()
    dotsub <- eval(mf$subset,data)
    if(!is.null(dotsub)) { 
      data <- data[dotsub,]
    }

    treata <- model.frame(formula,data)[,1,drop=FALSE]
    treat <- as.vector(treata[,1])
    names(treat) <- row.names(treata)
    covariates <- model.frame(delete.response(terms(formula)),data)[,,drop=FALSE]

    b <- exactmatch(formula, data, counter=counter)
    
    dummy.pscore <- rep(0.5, dim(covariates)[1])
    dummy.in.sample <- rep(TRUE, dim(covariates)[1])
    names(dummy.pscore) <- names(treat)
    names(dummy.in.sample) <- names(treat)
    c <- diagnose(formula, match.matrix=NULL, pscore=dummy.pscore,
                  in.sample = dummy.in.sample, data=data,exact=TRUE,
                  mahvars=NULL, subclass=max(b$psclass),
                  psclass=b$psclass, nearest=FALSE, counter=counter)
    
    psweights <- c$weights
    pscore <- dummy.pscore

    if (is.null(b$psclass)) psclass <- rep(1, length(pscore))
    else psclass <- b$psclass

    data <- cbind(data,psweights,pscore, psclass)
    
    z <- list(formula=formula, psweights=c$weights,
              treat=attr(terms(formula(cl)),"variables")[[2]],
              in.sample=NULL, match.matrix=b$match.matrix,
              covariates=attr(delete.response(tt),"term.labels"), psclass=b$psclass,
              q.cut=NULL, assign.model=NULL, call=cl, matched=(c$weights!=0),
              data= data)
    class(z) <- "matchit"
    return(z)
  }
  else {
    mf <- match.call()
    dotsub <- eval(mf$subset,data)
    mf[[1]] <- as.name("distance")
    mf$nearest <- mf$replace <- mf$m.order <- mf$ratio <- mf$caliper <- 
      mf$calclosest <- mf$subclass <- mf$sub.by <- mf$mahvars <-
        mf$exact <- NULL 
    a <- eval(as.call(mf), sys.frame(sys.parent()))
    if(!is.null(dotsub)) { 
      data <- data[dotsub,]
    }
    
    b <- matchdef(formula, a$in.sample, a$pscore,
                  nearest=nearest, replace=replace, m.order=m.order,
                  ratio=ratio, caliper=caliper, calclosest=calclosest,
                  mahvars=mahvars,data=data, exact=exact, counter=counter) 
    b1 <- subclassify(formula, data, a$in.sample, a$pscore, nearest=nearest,
                      b$match.matrix, subclass=subclass,
                      sub.by=sub.by, counter=counter)
    c <- diagnose(formula, b$match.matrix, a$pscore, a$in.sample,
                  data=data, exact=exact, mahvars=mahvars,
                  subclass=subclass, psclass=b1$psclass,nearest=nearest,
                  q.cut=b1$q.cut, counter=counter)
    #adding pscore and weights to dataframe
    pscore <- a$pscore
    psweights <- c$weights
    
    # Create psclass with subclasses
    if (is.null(b1$psclass)) psclass <- rep(1, length(pscore))
    else psclass <- b1$psclass
    
    data <- cbind(data,psweights,pscore, psclass)

    z <- list(formula=formula,psweights=psweights,
              treat=attr(terms(formula(cl)),"variables")[[2]],
              in.sample=a$in.sample, match.matrix=b$match.matrix,
              covariates=attr(delete.response(tt),"term.labels"),
              psclass=b1$psclass, q.cut=b1$q.cut, assign.model=a$assign.model,
              call=cl,matched=(c$weights!=0), data=data)
    class(z) <- "matchit"
    return(z)
  }
}
neyman<-function(Y, object, bootstrap=NULL, counter=TRUE){
  mf<-match.call()
  matchit.call <- object$call
  weighted.var<-function(x, w)
    sum(w*(x-weighted.mean(x, w))^2)/(sum(w)-1)
  est<-function(Y, object, data){
    treat<-eval(object$treat, data)
    Y<-eval(mf$Y, data)
    W<-object$psweights
    indx<-object$psclass
    if(is.null(indx)){
      m<-1
      indx<-rep(1, length(Y))
    } else {
      m<-max(na.omit(indx))
    }
    ate.est<-ate.var<-Ntrt<-Ncont<-matrix(NA, nrow=1, ncol=m)
    for(i in 1:m){
      Ntrt[i]<-length(Y[treat==1 & indx==i])
      Ncont[i]<-length(Y[treat==0 & indx==i])
      ate.est[i]<-weighted.mean(Y[treat==1 & indx==i],
                                W[treat==1 & indx==i]) -
                                weighted.mean(Y[treat==0 & indx==i],
                                              W[treat==0 & indx==i])
      ate.var[i]<-weighted.var(Y[treat==1 & indx==i],
                               W[treat==1 & indx==i])/sum(W[treat==1 & indx == i]) +
                                 weighted.var(Y[treat==0 & indx==i],
                                              W[treat==0 & indx==i])/sum(W[treat==0 & indx == i])
    }
    if(m==1) {
      Ntrt=length(Y[treat==1 & W!=0])
      Ncont=length(Y[treat==0 & W!=0])
    } else {
      colnames(ate.est)<-c(1:m)
      for(i in 1:m)
        colnames(ate.est)[i]<-paste("Subclass", i)
    }
    list(ate.est=ate.est, ate.var=ate.var, Ntrt=Ntrt, Ncont=Ncont, bootstrap=NULL)
  }
    
  if(class(object)!="matchit"){
    stop("the input object is not a matchit object")
  }

  data<-eval(object$data, sys.parent())

  if (identical(eval(matchit.call$exact), TRUE)) {
    Y <- eval(mf$Y, object$data)
    W <- object$psweights
    treat <- eval(object$treat, object$data)
    
    ate <- weighted.mean(Y[treat==1], W[treat==1]) -
      weighted.mean(Y[treat==0], W[treat==0])
    se <- sqrt(weighted.var(Y[treat==1], W[treat==1])/sum(W[treat==1])
               + weighted.var(Y[treat==0], W[treat==0])/sum(W[treat==0]))
    numtrt <- sum(W[treat==1])
    numcont <- sum(W[treat==0])
  }
  else {
    if(is.null(bootstrap)) {
      res<-est(Y, object, data)
    } else {
      ate.est<-Ntrt<-Ncont<-NULL
      cl<-object$call
      cl$counter<-FALSE
      cl$data<-as.name("dta")
      if(counter) {
        bseq<-floor(seq(bootstrap/10, bootstrap, bootstrap/10))
      }
      for(i in 1:bootstrap){
        dta<-data[sample(c(1:nrow(data)), size=nrow(data),
                         replace=T),!names(data) %in% c("pscore", "psweights", "psclass")]
        bobject<-eval(cl)
        tmp<-est(Y,bobject,dta)
        ate.est<-rbind(ate.est, tmp$ate.est)
        Ntrt<-rbind(Ntrt, tmp$Ntrt)
        Ncont<-rbind(Ncont, tmp$Ncont)
        if(counter & i%in%bseq)
          cat(10*which(bseq==i),"%...", sep="")
      }
      cat("Done\n")
      row.names(ate.est)<-c(1:bootstrap)
      res<-list(ate.est=ate.est, ate.var=NULL, Ntrt=Ntrt, Ncont=Ncont, bootstrap=bootstrap)
    }
    
  # Subclass estimates calculated above...now aggregate
    if(is.null(object$call$sub.by)){
      sub.by <- "treat"
    } else {
      sub.by <- object$call$sub.by
    }
    if(!is.null(bootstrap)){
      if(is.null(object$call$subclass) | identical(object$call$subclass,0)){
        ate<-apply(res$ate.est, 2, mean)
        numtrt <- apply(res$Ntrt, 2, mean)
        numcont <- apply(res$Ncont, 2, mean)
      } else{
        if(sub.by=="treat"){
          w<-res$Ntrt/apply(res$Ntrt,1,sum)
        } else if(sub.by=="control"){
          w<-res$Ncont/apply(res$Ncont,1,sum)
        } else if(sub.by=="all"){
          w <- (res$Ncont+res$Ntrt)/apply((res$Ncont+res$Ntrt),1,sum)
        }
        # Set w to 0 if block has fewer then 2 t or c units
        w[res$Ntrt<=1 | res$Ncont<=1] <- rep(0, sum(res$Ntrt<=1 | res$Ncont<=1)) 
        ate.boot <- NULL
        for(i in 1:bootstrap){
          ate.boot<-c(ate.boot,weighted.mean(res$ate.est[i,],w[i,]))
        }
        ate<-c(mean(ate.boot),apply(res$ate.est, 2, mean))
        numtrt <- c(mean(apply(res$Ntrt, 1, sum)), apply(res$Ntrt, 2, mean))
        numcont <- c(mean(apply(res$Ncont, 1, sum)), apply(res$Ncont, 2, mean))
      }
    }
    else {
      if(sub.by=="treat"){
        w<-res$Ntrt/sum(res$Ntrt)
      } else if(sub.by=="control"){
        w<-res$Ncont/sum(res$Ncont)
      } else if(sub.by=="all"){
        w <- (res$Ncont+res$Ntrt)/sum(res$Ncont+res$Ntrt)
      } 
      # Set w to 0 if block has fewer then 2 t or c units
      w[res$Ntrt<=1 | res$Ncont<=1] <- rep(0, sum(res$Ntrt<=1 | res$Ncont<=1)) 
      ate<-apply(res$ate.est, 2, weighted.mean, w)
      if (length(res$Ntrt)> 1) {
	numtrt <- c(sum(res$Ntrt), res$Ntrt)
        numcont <- c(sum(res$Ncont), res$Ncont)
      } else {
	numtrt <- res$Ntrt
	numcont <- res$Ncont
      }
    }
    if(is.null(res$ate.var)){
      if(is.null(object$call$subclass) | identical(object$call$subclass,0)){
        se<-sqrt(apply(res$ate.est, 2, var))
      } else {
        se <- c(sd(ate.boot),sqrt(apply(res$ate.est, 2, var)))
      }
    } else{
      se<-sqrt(res$ate.var)
    }
    if(is.null(bootstrap)){
      if(length(ate)>1){
        if(is.null(res$ate.var)){
          tmp<-apply(res$ate.est, 1, weighted.mean, w)
          ate<-c(mean(tmp), ate)
          se<-c(sqrt(var(tmp)), se)
        }    
        else{
	  # weight=0 if block has < 2 units in either group
          ate<-c(weighted.mean(ate, w, na.rm=T), ate)
	  # Only calculate se for blocks with > 1 unit (se!=nan)	
          se<-c(sqrt(sum((w[!is.nan(se)]*se[!is.nan(se)])^2)), se)
        }
      }
    }
  }
  
  ressum<-list(ate=ate, se=se, Ntrt=numtrt, Ncont=numcont)
  class(ressum) <- "neyman"
  ressum
}
plot.matchit <- function(x, ...)
  {
   if (identical(eval(x$call$exact),TRUE)){
    return(warning("Plot() not appropriate for exact matching.  No plots generated"))
    }
    match.matrix <- x$match.matrix
    xdata <- x$data
    treata <- model.frame(x$formula,xdata)[,1,drop=FALSE]
    pscore <- xdata[,"pscore"]
    in.sample <- x$in.sample
    covariates <- x$covariates
    treat <- as.vector(treata[,1])
    names(treat) <- row.names(treata)
   covariates <- model.frame(delete.response(terms(x$formula)),xdata)[,,drop=FALSE]
    mahvars <- eval(x$call$mahvars)
    exact <- eval(x$call$exact)

 if (!is.null(mahvars)){
      w <- mahvars%in%names(covariates)
      if(sum(w)!=length(mahvars)){
      md <- as.data.frame(as.matrix(xdata[,mahvars[!w]]))
      names(md) <- mahvars[!w]
      covariates <- cbind(covariates,md)
    }
  }
  if (!identical(TRUE, exact) & !identical(FALSE, exact)) {
    w <- exact%in%names(covariates)
    if(sum(w)!=length(exact)) {
      ed <- as.data.frame(as.matrix(xdata[,exact[!w]]))
      names(ed) <- exact[!w]
      covariates <- cbind(covariates,ed)
    }
  }
    psclass <- x$psclass
    if(is.null(psclass)){subclass <- 0}else{subclass <- qbins <- max(na.omit(x$psclass))}
    nearest <- !is.null(x$match.matrix)
    q.cut <- x$q.cut
    n <- length(treat)
    n0 <- length(treat[treat==0])
    n1 <- length(treat[treat==1])

    if(is.null(names(treat))){names(treat) <- seq(1,n)}
    labels <- names(treat)
    clabels <- labels[treat==0]
    tlabels <- labels[treat==1]
   if(is.null(in.sample)){in.sample <- rep(TRUE,n); names(in.sample) <- names(treat)}
    if(nearest==TRUE) {
      match.matrix <- match.matrix[in.sample[treat==1],,drop=F]
      num.matches <- dim(match.matrix)[2]-apply(as.matrix(match.matrix), 1, function(x) { sum(is.na(x)) })
      names(num.matches) <- tlabels[in.sample[treat==1]]
      t.units <- row.names(match.matrix)[num.matches>0]
      c.units <- na.omit(as.vector(as.matrix(match.matrix)))
      matched <-c(t.units,unique(c.units))
      weights <- rep(0,length(treat))
      names(weights) <- labels
      weights[t.units] <- 1
      
      for (cont in clabels) {
        treats <- na.omit(row.names(match.matrix)[cont==match.matrix[,1]])
        if (dim(match.matrix)[2]>1) {
          for (j in 2:dim(match.matrix)[2]) 
            treats <- c(na.omit(row.names(match.matrix)[cont==match.matrix[,j]]),treats)	
        }
	if (length(treats)==0) weights[cont] <- 0
	else for (k in 1:length(treats)) weights[cont] <- weights[cont]+1/num.matches[treats[k]] 
      }

    } else {
      	matched <- names(treat)
	weights <- rep(1,length(treat))
	names(weights) <- names(treat)
	}
    
    weights[!in.sample] <- 0 

    doverlay <- function(x,treat,weights=NULL,xlab="",main="",lines=FALSE){
      xmiss <- !is.na(x)
      x <- x[xmiss]
      treat <- treat[xmiss]
      weights <- weights[xmiss]
      minobs <- min(x)
      maxobs <- max(x)
      dx1 <- density(x[treat==1],from=minobs,to=maxobs)
      dx0 <- density(x[treat==0],from=minobs,to=maxobs)
      if(!is.null(weights))
        {
          x1 <- x[treat==1&weights!=0]
          x0 <- sample(x[treat==0], size=min(10000,(100*length(x[treat==0]))),
                       replace=TRUE,prob=weights[treat==0]/sum(weights[treat==0]))
          d1 <- density(x1,from=minobs,to=maxobs)
          bw <- d1$bw 
          d0 <- density(x0,from=minobs,to=maxobs,bw=bw)
          par(mfrow=c(2,1))
          matplot(dx0$x,cbind(dx1$y,dx0$y),type="l",ylim=range(c(dx0$y,dx1$y,d1$y,d0$y)),
                  ylab="Density",xlab=xlab,main=paste(main,": All Units",sep=""))
          legend(minobs,max(c(d1$y,d0$y,dx1$y,dx0$y)), lty=1:2, col=1:2,
                 legend=c("Treatment","Control"))
          matplot(dx0$x,cbind(d1$y,d0$y),type="l",ylim=range(c(dx0$y,dx1$y,d1$y,d0$y)),
                  ylab="Density",xlab=xlab,main=paste(main,": Matched Units",sep=""))
          legend(minobs,max(c(d1$y,d0$y,dx1$y,dx0$y)), lty=1:2, col=1:2,
                 legend=c("Treatment","Control"))
          if (lines==T) abline(v=q.cut,col="grey",lty=1)
          par(mfrow=c(1,1))
        } else{
          matplot(dx0$x,cbind(dx1$y,dx0$y),type="l",ylab="Density",xlab=xlab,main=main)
          legend(minobs,max(c(dx1$y,dx0$y)),lty=1:2,col=1:2,legend=c("Treatment","Control"))
	  if (lines==T) abline(v=q.cut, col="grey", lty=1)
        }
    }
    
    choice.menu <- function(choices,question)
      {
        k <- length(choices)-1
        Choices <- data.frame(choices)
        row.names(Choices) <- 0:k
        names(Choices) <- "Choices"
        print.data.frame(Choices,right=FALSE)
        ans <- readline(question)          
        while(!ans%in%c(0:k))
          {
            print("Not valid -- please pick one of the choices")
            print.data.frame(Choices,right=FALSE)
            ans <- readline(question)
          }
        return(ans)
      }
    if(!is.null(pscore))
      {
        densq <- ("Would you like to see density estimates of the propensity scores?")
        denschoice <- c("No","Yes")
        densplot <- choice.menu(denschoice,densq)
        if(densplot==1){
          if(nearest){
            doverlay(pscore,treat,weights,main="Propensity Score",lines=T)
            #if(length(subclass)!=1 | subclass!=0){abline(v=q.cut,col="grey",lty=1)}
	  } else {
              doverlay(pscore,treat,main="Propensity Score", lines=T) 
              #if(length(subclass)!=1 | subclass!=0){abline(v=q.cut,col="grey",lty=1)}
          }
        }
      
        jitq <- ("Would you like to see a jitterplot of the propensity scores?")
        jitchoice <- c("No","Yes")
        jitplot <- choice.menu(jitchoice,jitq)
        if(jitplot==1) {
          jitp <- jitter(rep(1,length(treat)),factor=6)-(treat==0)
          cwt <- sqrt(weights)
          plot(pscore,xlim=range(na.omit(pscore)),ylim=c(-1,2),type="n",ylab="",xlab="Propensity Score",axes=F,main="Distribution of Propensity Scores")
          if(length(subclass)!=1 | subclass!=0){abline(v=q.cut,col="grey",lty=1)}
          points(pscore[treat==1&weights!=0],jitp[treat==1&weights!=0],pch=18,cex=cwt[treat==1&weights!=0])
          points(pscore[treat==1&weights==0],jitp[treat==1&weights==0],pch=5,col="grey",cex=0.5)
          points(pscore[treat==0&weights!=0],jitp[treat==0&weights!=0],pch=18,cex=cwt[treat==0&weights!=0])
          points(pscore[treat==0&weights==0],jitp[treat==0&weights==0],pch=5,col="grey",cex=0.5)
          axis(1)
          text(sum(range(na.omit(pscore)))/2,1.5,"Treatment Units")
          text(sum(range(na.omit(pscore)))/2,-0.5,"Control Units")
          box()
          print("To identify the units, use first mouse button; to stop, use second.")
          identify(pscore,jitp,names(treat))
        }
      }

    if(nearest){
      choices <- c("No",paste("Yes : ",names(covariates)))
      question <- "Would you like to see density estimates of any other covariates?"
      ans <- -1
      while(ans!=0)
        {
          ans <- as.numeric(choice.menu(choices,question))
          if(ans!=0)
            if(nearest)
              {
                doverlay(covariates[,(ans)],treat,weights,main=names(covariates)[ans])
              } else {
                doverlay(covariates[,(ans)],treat,main=names(covariates)[ans])
              }
        }
    }
    if(length(subclass)!=1 | subclass!=0){
      choices <- c("No",paste("Yes : Subclass ", 1:qbins))
      question <- "Would you like to see density estimates of any subclass covariates?"
      ans <- -1
      while(ans!=0)
        {
          ans <- as.numeric(choice.menu(choices,question))
          if(ans!=0)
            {
              question2 <- "Which covariates?"
              choices2 <- c(paste(names(covariates)))
              ans2 <- as.numeric(choice.menu(choices2,question2))
              ans2 <- ans2+1
              if(sum(treat[psclass==ans]==1,na.rm=TRUE)<=2){
                print("Not enough treatment units in this subclass")
              } else if(sum(treat[psclass==ans]==0,na.rm=TRUE)<=2){
                print("Not enough control units in this subclass")} else{
                  doverlay(covariates[,(ans2)][psclass==ans],treat[psclass==ans],main=paste("Subclass ",ans," : ",names(covariates)[ans2]))}
            }
        }
      
    } 
    
  }

print.matchit <- function(x,  digits = max(3, getOption("digits") -
  3), ...)
  {
    t.test.wtd <- function(x, treat, weights) {
      trt <- cov.wt(as.matrix(x[treat==1]), wt=weights[treat==1])
      mean.t <- trt$center
      cov.t <- trt$cov    
      cont <- cov.wt(as.matrix(x[treat==0]), wt=weights[treat==0])
      mean.c <- cont$center
      cov.c <- cont$cov
      n.t <- sum(weights[treat==1])
      n.c <- sum(weights[treat==0])
      ratio.t <- cov.t/n.t
      ratio.c <- cov.c/n.c
      tt  <- (mean.t-mean.c)/sqrt(ratio.t+ratio.c)
    }
    weighted.var <- function(x, w) {
      sum(w * (x - weighted.mean(x,w))^2)/(sum(w) - 1)}
    qoi <- function(xx,tt,ww){
      xsum <- matrix(0,2,5)
      xsum <- as.data.frame(xsum)
      row.names(xsum) <- c("Full","Matched")
      names(xsum) <- c("Means Treated","Means Control","SD","T-stat","Bias")
      xn <- matrix(0,2,3)
      xn <- as.data.frame(xn)
      row.names(xn) <- c("Full","Matched")
      names(xn) <- c("Treated","Control","Total")
      xn[1,1] <- sum(tt==1)
      xn[1,2] <- sum(tt==0)
      xn[1,3] <- length(tt)
      x1 <- xx[tt==1]
      x0 <- xx[tt==0]
      xsum[1,1] <- mean(x1,na.rm=T)
      xsum[1,2] <- mean(x0,na.rm=T)
      X.t.m <- xx[tt==1][ww[tt==1]>0]
      X.c.m <- xx[tt==0][ww[tt==0]>0]
      xsum[2,1] <- weighted.mean(X.t.m, ww[tt==1][ww[tt==1]>0])
      xsum[2,2] <- weighted.mean(X.c.m, ww[tt==0][ww[tt==0]>0])
      if(sum(tt==1)<2|(sum(tt==0)<2)){
        xsum[c(1,2),c(3,4,5)] <- NA
      } else {
        xsum[1,3] <- sd(xx,na.rm=T) 
        xsum[1,4] <- -1*t.test(xx~tt)$sta
        xsum[1,5] <- (mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x1,na.rm=T)
        xsum[2,3] <- sqrt(weighted.var(xx[ww>0],ww[ww>0]))
        xsum[2,4] <- t.test.wtd(xx[ww>0],tt[ww>0],ww[ww>0])
        xsum[2,5] <- (xsum[2,1]-xsum[2,2])/sqrt(weighted.var(X.t.m,ww[tt==1][ww[tt==1]>0]))
      } 
      xn[2,1] <- sum(tt[ww>0]==1)
      xn[2,2] <- sum(tt[ww>0]==0)
      xn[2,3] <- length(tt[ww>0])
      list(xsum=xsum,xn=xn)
    }
    qoi.by.sub <- function(xx,tt,ww,qq){
      qbins <- max(qq,na.rm=TRUE)
      q.table <- matrix(0,5,qbins)
      qn <- matrix(0,3,qbins)
      matched <- ww!=0
      for (i in 1:qbins)
        {
          qi <- qq[matched]==i
          qx <- xx[matched][qi]
          qt <- tt[matched][qi]
          qw <- ww[matched][qi]
          if(sum(qt==1)<2|(sum(qt==0)<2)){
            if(sum(qt==1)<2){warning("Not enough treatment units in subclass ",i,call.=FALSE)}
            else if(sum(qt==0)<2){warning("Not enough control units in subclass ",i,call.=FALSE)}
          }
          qoi.i <- qoi(qx,qt,qw)
          q.table[,i] <- as.numeric(qoi.i$xsum[1,])
          qn[,i] <- as.numeric(qoi.i$xn[1,])
        }
      q.table <- as.data.frame(q.table)
      qn <- as.data.frame(qn)
      names(q.table) <- names(qn) <- paste("Subclass",1:qbins)
      row.names(q.table) <- names(qoi.i$xsum)
      row.names(qn) <- c("Treated","Control","Total")
      list(q.table=q.table,qn=qn)
    }
    
    pscore <- x$data[,"pscore"]
    treat <- eval(x$treat,x$data)
    names(treat) <- names(pscore)
    cat("\nAssignment model specification:\n", deparse(x$call),"\n\n",sep = "")
    if(identical(eval(x$call$exact),TRUE)) #for exact matching
      {
        cat("\nSample sizes for full and exactly matched data:\n\n")
        xqoi <- qoi(pscore, treat, x$psweights)
        print.data.frame(xqoi$xn, digits = digits)
        cat("\n\n")
      } else {        
        if(is.null(x$match.matrix))
          {
            cat("Summary of propensity score for full sample:\n\n")
            xqoi <- qoi(pscore, treat, x$psweights)
            print.data.frame(xqoi$xsum[1,], digits = digits)
            cat("\nSample sizes:\n\n")
            print.data.frame(xqoi$xn[1,], digits = digits)
          } else {
            cat("\nSummary of propensity score for full and matched samples:\n\n")
            xqoi <- qoi(pscore, treat, x$psweights)
            print.data.frame(xqoi$xsum, digits = digits)
            cat("\nSample sizes:\n\n")
            print.data.frame(xqoi$xn, digits = digits)
          }
        cat("\n\n")
        if(!is.null(x$psclass))
          {
            qqoi <- qoi.by.sub(pscore, treat, x$psweights, x$psclass)
            cat("\nSummary of propensity score by subclasses:\n\n")
            print.data.frame(qqoi$q.table, digits = digits)
            cat("\nSample sizes by subclasses:\n\n")
            print.data.frame(qqoi$qn, digits = digits)
            cat("\n\n")
            invisible(x)
          }
      } 
  }
print.neyman<-function(x, digits=max(3, getOption("digits") -3), ...){
  print.res<-cbind(x$ate, x$se)
  colnames(print.res)<-c("ATE", "SD")
  print.res <- print.res[1,]
                                                                                                                                                             
  cat("\n Average Treatment Effect:\n \n")
  print.matrix(print.res, digits=digits, ...)
  cat("\n")
}
print.summary.matchit <- function(x, digits = max(3, getOption("digits") - 3), ...){  
  verbose <- x$verbose
  sum.all <- x$sum.all
  sum.matched <- x$sum.matched
  q.table <- x$q.table
  xn <- x$xn
  sig <- x$sig
  cat("\nAssignment model specification:\n", deparse(x$call),"\n\n",sep = "")
  if(verbose){
    cat("Summary of covariates and interactions for all data:\n\n")
  } else{
    cat("Summary of covariates for all data:\n\n")
  }
  print(sum.all, digits=digits)
  cat("\n")
  xs1 <- sum.matched
  cc <- row.names(sum.all)
  if(!is.null(x$match.matrix))
    {
      if(verbose){
        cat("Summary of covariates and interactions for matched data:\n\n")
      } else {
        cat("Summary of covariates for matched data:\n\n")
      }
      print(xs1, digits=digits)
      cat("\nSample sizes:\n\n")
      print.data.frame(xn, digits=digits)
      cat("\n")
      cat("Problematic covariates:  ")
      xs2 <- xs1[1:(nrow(xs1)-1),]
      cat(row.names(xs2)[!is.na(xs2[,4]) & abs(xs2[,4])>sig])
      cat("\n")
    }
  if(!is.null(x$psclass))
    {
      if(identical(eval(x$call$exact),TRUE)){
        cat("\nSample sizes for full and exactly matched data:\n\n")
        print(xn,digits = digits)
        cat("\nSample sizes by covariates:\n\n")
        print(q.table, digits = digits)
        cat("\n")
      } else {
        print(q.table, digits = digits)
        cat("\nSample sizes by subclasses:\n\n")
        print.data.frame(x$qn, digits = digits)
        cat("\n")
        cat("Problematic covariates:\n")
        for(i in 1:dim(q.table)[3])
          {
            cat("Subclass ", i, ":  ")
            xs1 <- q.table[cc,,i]
            cat(row.names(xs1)[!is.na(xs1[,4]) & abs(xs1[,4])>sig])
            cat("\n")
          }
      }
    }
  
  cat("Number of units discarded:  ", sum(x$in.sample==FALSE))
  cat("\n\n")
  invisible(x)
}
print.summary.neyman<-function(x, digits=max(3, getOption("digits")-3),
                               signif.stars=getOption("show.signif.stars"), 
                               ...){
  cat("\n Average Treatment Effect:\n")
  printCoefmat(x$results, digits=digits, signif.starts = signif.stars, ...)
  cat("\n")

  cat("\n Sample Sizes:\n")
  print.matrix(x$sample.size)
  cat("\n")
}
subclassify <- function(formula,data,in.sample,pscore,nearest=TRUE,
                        match.matrix,subclass=0,sub.by="treat", counter=TRUE){
  data <- eval(data,parent.frame())
  treata <- model.frame(formula,data)[,1,drop=FALSE]
  treat <- as.vector(treata[,1])
  names(treat) <- row.names(treata)

  if(!identical(subclass,0))
    {
      if(counter){cat("Subclassifying...")}  
      n <- length(treat)
      if(nearest){
        match.matrix <- match.matrix[in.sample[treat==1],,drop=F]
        t.units <- row.names(match.matrix)[!is.na(match.matrix)]

        c.units <- na.omit(as.vector(as.matrix(match.matrix)))
        matched <-c(t.units,c.units)
        matched <- names(treat)%in%matched
       } else{
          matched <- rep(TRUE,n)}
      names(matched) <- names(treat)
      m1 <- matched[treat==1]
      m0 <- matched[treat==0]
      p1 <- pscore[treat==1][m1]
      p0 <- pscore[treat==0][m0]
      
      if(length(subclass)!=1 | (length(subclass)==1 & identical(subclass<1,TRUE)))
        {
	  subclass <- sort(subclass)
	  if (subclass[1]==0) subclass <- subclass[-1]
	  if (subclass[length(subclass)]==1) subclass <- subclass[-length(subclass)]
          if(sub.by=="treat")
            {
              q <- c(0,quantile(p1,probs=c(subclass)),1)
            }
          else if(sub.by=="control")
            {
            q <- c(0,quantile(p0,probs=c(subclass)),1)
            }
          else if(sub.by=="all")
            {
              q <- c(0,quantile(pscore,probs=c(subclass)),1)
	  }
#          else {stop("Must specify a valid sub.by",call.=FALSE)
#              }
        }
      else        {
#        if(subclass<=0){stop("Subclass must be a positive vector",call.=FALSE)}
        sprobs <- seq(0,1,length=(round(subclass)+1))
        sprobs <- sprobs[2:(length(sprobs)-1)]
        if(sub.by=="treat")            {
          q <- c(0,quantile(p1,probs=sprobs,na.rm=TRUE),1)
        }
        else if(sub.by=="control")
        {
          q <- c(0,quantile(p0,probs=sprobs,na.rm=TRUE),1)
        }
      else if(sub.by=="all")
        {
          q <- c(0,quantile(pscore,probs=sprobs,na.rm=TRUE),1)
        }
#      else {stop("Must specify a valid sub.by",call.=FALSE)}
      }

      #print(sprobs)
      #print(q)
      qbins <- length(q)-1
      psclass <- rep(0,n)
      names(psclass) <- names(treat)
      for (i in 1:qbins){
        q1 <- q[i]
        q2 <- q[i+1]
        psclass <- psclass+i*as.numeric(pscore<q2 & pscore>=q1)}
      
                                        # Make sure not to assign subclass to discarded units
      psclass[in.sample==0] <- 0
      psclass[!matched] <- 0
      if(counter){cat("Done\n")}
    }  
  
  else {psclass <- NULL; q=NULL}

  list(psclass = psclass, q.cut = q)
}
summary.matchit <- function(object, verbose=F, sig=2, ...) {
  t.test.wtd <- function(x, treat, weights) {
    trt <- cov.wt(as.matrix(x[treat==1]), wt=weights[treat==1])
    mean.t <- trt$center
    cov.t <- trt$cov    
    cont <- cov.wt(as.matrix(x[treat==0]), wt=weights[treat==0])
    mean.c <- cont$center
    cov.c <- cont$cov
    n.t <- sum(weights[treat==1])
    n.c <- sum(weights[treat==0])
    ratio.t <- cov.t/n.t
    ratio.c <- cov.c/n.c
    tt  <- (mean.t-mean.c)/sqrt(ratio.t+ratio.c)
  }
  weighted.var <- function(x, w) {
    sum(w * (x - weighted.mean(x,w))^2)/(sum(w) - 1)}
  qoi <- function(xx,tt,ww){
    xsum <- matrix(0,2,5)
    xsum <- as.data.frame(xsum)
    row.names(xsum) <- c("Full","Matched")
    names(xsum) <- c("Means Treated","Means Control","SD","T-stat","Bias")
    xn <- matrix(0,2,3)
    xn <- as.data.frame(xn)
    row.names(xn) <- c("Full","Matched")
    names(xn) <- c("Treated","Control","Total")
    xn[1,1] <- sum(tt==1)
    xn[1,2] <- sum(tt==0)
    xn[1,3] <- length(tt)
    x1 <- xx[tt==1]
    x0 <- xx[tt==0]
    xsum[1,1] <- mean(x1,na.rm=T)
    xsum[1,2] <- mean(x0,na.rm=T)
    X.t.m <- xx[tt==1][ww[tt==1]>0]
    X.c.m <- xx[tt==0][ww[tt==0]>0]
    xsum[2,1] <- weighted.mean(X.t.m, ww[tt==1][ww[tt==1]>0])
    xsum[2,2] <- weighted.mean(X.c.m, ww[tt==0][ww[tt==0]>0])
    if(sum(tt==1)<2|(sum(tt==0)<2)){
      xsum[c(1,2),c(3,4,5)] <- NA
    } else {
      xsum[1,3] <- sd(xx,na.rm=T) 
      xsum[1,4] <- -1*t.test(xx~tt)$sta
      xsum[1,5] <- (mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x1,na.rm=T)
      xsum[2,3] <- sqrt(weighted.var(xx[ww>0],ww[ww>0]))
      xsum[2,4] <- t.test.wtd(xx[ww>0],tt[ww>0],ww[ww>0])
      xsum[2,5] <- (xsum[2,1]-xsum[2,2])/sqrt(weighted.var(X.t.m,ww[tt==1][ww[tt==1]>0]))
    } 
    xn[2,1] <- sum(tt[ww>0]==1)
    xn[2,2] <- sum(tt[ww>0]==0)
    xn[2,3] <- length(tt[ww>0])
    list(xsum=xsum,xn=xn)
  }
  qoi.by.sub <- function(xx,tt,ww,qq){
    qbins <- max(qq,na.rm=TRUE)
    q.table <- matrix(0,5,qbins)
    qn <- matrix(0,3,qbins)
    matched <- ww!=0
    for (i in 1:qbins)
        {
          qi <- qq[matched]==i
          qx <- xx[matched][qi]
          qt <- tt[matched][qi]
          qw <- ww[matched][qi]
          if(sum(qt==1)<2|(sum(qt==0)<2)){
            if(sum(qt==1)<2){warning("Not enough treatment units in subclass ",i,call.=FALSE)}
            else if(sum(qt==0)<2){warning("Not enough control units in subclass ",i,call.=FALSE)}
          }
          qoi.i <- qoi(qx,qt,qw)
          q.table[,i] <- as.numeric(qoi.i$xsum[1,])
          qn[,i] <- as.numeric(qoi.i$xn[1,])
        }
    q.table <- as.data.frame(q.table)
    qn <- as.data.frame(qn)
    names(q.table) <- names(qn) <- paste("Subclass",1:qbins)
    row.names(q.table) <- names(qoi.i$xsum)
      row.names(qn) <- c("Treated","Control","Total")
    list(q.table=q.table,qn=qn)
  }  
  #setting up covariates
  pscore <- object$data[,"pscore"]
  treat <- eval(object$treat,object$data)
  names(treat) <- names(pscore)
  data <- object$data
  weights <- object$psweights
  covariates <- as.data.frame(model.matrix(terms(object$formula),data)[,,drop=FALSE])
  if ("(Intercept)"%in%colnames(covariates))
    covariates <- covariates[,-match("(Intercept)", colnames(covariates)),drop=FALSE]
  mahvars <- eval(object$call$mahvars)
  exact <- eval(object$call$exact)
  subclass <- eval(object$call$subclass)
  if(is.null(subclass)){subclass <- 0}
  psclass <- object$psclass
  if (!is.null(mahvars)){
    w <- mahvars%in%names(covariates)
    if(sum(w)!=length(mahvars)){
      md <- as.data.frame(as.matrix(data[,mahvars[!w]]))
      names(md) <- mahvars[!w]
      row.names(md) <- row.names(covariates)
      covariates <- cbind(covariates,md)
    }
  }
  if (!identical(TRUE, exact) & !identical(FALSE, exact) & !is.null(exact)) {
    w <- eval(exact)%in%names(covariates)
    if(sum(w)!=length(exact)) {
      ed <- as.data.frame(as.matrix(data[,exact[!w]]))
      names(ed) <- exact[!w]
      row.names(ed) <- row.names(covariates)
      covariates <- cbind(covariates,ed)
    }    
  }
  if(!identical(TRUE,exact)){ 
    covariates<- cbind(pscore,covariates)
  } 
  
  aa <- apply(covariates,2,qoi,tt=treat,ww=weights)
  k <- length(aa)
  sum.all <- as.data.frame(matrix(0,k,5))
  sum.matched <- as.data.frame(matrix(0,k,6))
  row.names(sum.all) <- row.names(sum.matched) <- names(aa)
  names(sum.all) <- names(sum.matched) <- names(aa[[1]]$xsum)
  names(sum.matched) <- c(names(aa[[1]]$xsum),"Reduction")
  sum.all.int <- sum.matched.int <- NULL
  for(i in 1:k){
    sum.all[i,] <- aa[[i]]$xsum[1,]
    sum.matched[i,1:5] <- aa[[i]]$xsum[2,]
    sum.matched[i,6] <- (abs(sum.all[i,5])-abs(sum.matched[i,5]))>0
    if(verbose){
      for(j in i:k){
        x2 <- covariates[,i]*as.matrix(covariates[,j])
        sum.all.int <- rbind(sum.all.int,qoi(x2,tt=treat,ww=weights)$xsum[1,])
        Reduction <- (abs(sum.all.int[nrow(sum.all.int),5])-abs(qoi(x2,tt=treat,ww=weights)$xsum[2,5]))>0
        sum.matched.int <- rbind(sum.matched.int,cbind(qoi(x2,tt=treat,ww=weights)$xsum[2,],Reduction))
        row.names(sum.all.int)[nrow(sum.all.int)] <-
          row.names(sum.matched.int)[nrow(sum.matched.int)] <-
            paste(names(covariates)[i],names(covariates)[j],sep="x")
      }
    }
  }
  xn <- aa[[1]]$xn
  sum.all <- rbind(sum.all,sum.all.int)
  sum.matched <- rbind(sum.matched,sum.matched.int)
  
  #now subclassification
  if(!identical(subclass,0))
    {
      qbins <- max(psclass,na.rm=TRUE)
      if(verbose){
        q.table <- array(0,dim=c(k+sum(1:k),5,qbins))
        ii <- 0
        nn <- NULL
      } else {
        q.table <- array(0,dim=c(k,5,qbins))
      }
      aa <- apply(covariates,2,qoi.by.sub,tt=treat,ww=weights,
                  qq=object$psclass)
      for(i in 1:k){
        if(!verbose){
          q.table[i,,] <- as.matrix(aa[[i]]$q.table)
          nn <- names(aa)
        } else {
          ii <- ii + 1 
          q.table[ii,,] <- as.matrix(aa[[i]]$q.table)
          nn <- c(nn,names(aa)[i])
          for(j in i:k){
            ii <- ii + 1 
            x2 <- covariates[,i]*as.matrix(covariates[,j])
            q.table[ii,,] <- as.matrix(qoi.by.sub(x2,tt=treat,ww=weights,qq=object$psclass)$q.table)
            nn <- c(nn,paste(names(covariates)[i],names(covariates)[j],sep="x"))
          }
        }
      }   
      qn <- aa[[1]]$qn
      dimnames(q.table) <- list(nn,row.names(aa[[i]]$q.table),paste("Subclass",1:qbins))
    } else {
      q.table <- NULL; q.bias <- NULL; qn <- NULL
    }
  if(identical(exact,TRUE)){
    qbins <- max(object$psclass,na.rm=TRUE)
    q.table <- as.data.frame(matrix(0,qbins,k+3))
    names(q.table) <- c(names(covariates),"Treated","Control","Total")
    for(i in 1:qbins){
      qi <- object$psclass==i
      q.table[i,] <- c(as.numeric(covariates[qi,,drop=F][1,]), sum(treat[qi]==1), sum(treat[qi]==0), length(treat[qi]))
    }
  }
  ans <- list(sum.all = sum.all, sum.matched = sum.matched,
              q.table=q.table, call=object$call, verbose=verbose,
              xn = xn, match.matrix = object$match.matrix, sig = sig,
              psclass=object$psclass, in.sample = object$in.sample,
              qn = qn)
  class(ans) <- "summary.matchit"
  ans
}  

summary.neyman<-function(object,...){
  tval<-object$ate/object$se
  res<-cbind(object$ate, object$se, tval, 2*pnorm(abs(tval), lower.tail=FALSE))
  colnames(res)<-c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
  rownames(res)[1]<-c("Overall")
  
  sampsize <- cbind(object$Ntrt, object$Ncont)
  colnames(sampsize) <- c("Ntrt", "Ncont")
  rownames(sampsize) <- rownames(res)

  output <- list(results=res, sample.size=sampsize)
  class(output)<-"summary.neyman"
  output
}
