[r-cran-matchit] 01/45: Import Upstream version 1.0-1
Andreas Tille
tille at debian.org
Fri Oct 20 06:17:18 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-matchit.
commit 01c6119114e55363eab75fca8687a4987aa01e51
Author: Andreas Tille <tille at debian.org>
Date: Fri Oct 20 07:40:48 2017 +0200
Import Upstream version 1.0-1
---
DESCRIPTION | 19 ++
NAMESPACE | 25 ++
R/diagnose.R | 82 +++++++
R/distance.R | 136 ++++++++++
R/exactmatch.R | 25 ++
R/help.matchit.R | 53 ++++
R/match.data.R | 12 +
R/matchdef.R | 215 ++++++++++++++++
R/matchit.R | 174 +++++++++++++
R/neyman.R | 164 +++++++++++++
R/plot.matchit.R | 203 +++++++++++++++
R/print.matchit.R | 146 +++++++++++
R/print.neyman.R | 9 +
R/print.summary.matchit.R | 61 +++++
R/print.summary.neyman.R | 11 +
R/subclassify.R | 124 ++++++++++
R/summary.matchit.R | 220 +++++++++++++++++
R/summary.neyman.R | 14 ++
data/incomedata.tab | 201 +++++++++++++++
data/lalonde.tab | 615 ++++++++++++++++++++++++++++++++++++++++++++++
demo/00Index | 3 +
demo/analysis.R | 82 +++++++
demo/diagnose.R | 35 +++
demo/lalonde.R | 112 +++++++++
man/diagnose.Rd | 85 +++++++
man/distance.Rd | 77 ++++++
man/exactmatch.Rd | 59 +++++
man/help.matchit.Rd | 42 ++++
man/incomedata.Rd | 23 ++
man/lalonde.Rd | 31 +++
man/match.data.Rd | 37 +++
man/matchdef.Rd | 85 +++++++
man/matchit.Rd | 146 +++++++++++
man/neyman.Rd | 47 ++++
man/plot.matchit.Rd | 39 +++
man/print.matchit.Rd | 36 +++
man/subclassify.Rd | 85 +++++++
man/summary.matchit.Rd | 65 +++++
man/summary.neyman.Rd | 44 ++++
39 files changed, 3642 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..325109c
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,19 @@
+Package: MatchIt
+Version: 1.0-1
+Date: 2005-1-3
+Title: MatchIt: Nonparametric Preprocessing for Parametric Casual Inference
+Author: Daniel Ho <daniel.ho at yale.edu>,
+ Kosuke Imai <kimai at Princeton.Edu>,
+ Gary King <king at harvard.edu>,
+ Elizabeth Stuart <stuart at stat.harvard.edu>
+Maintainer: Daniel Ho <daniel.ho at yale.edu>
+Depends: R (>= 2.0)
+Suggests: Zelig, optmatch
+Description: MatchIt selects matched samples of the
+ the original treated and control groups with similar
+ covariate distributions -- can be used to match exactly
+ on covariates, to match on propensity scores, or perform
+ a variety of other matching procedures.
+License: GPL version 2 or newer
+URL: http://gking.harvard.edu/matchit
+Packaged: Thu Jan 6 13:44:52 2005; estuart
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..7b19386
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,25 @@
+export(
+ matchdef,
+ match.data,
+ matchit,
+ diagnose,
+ distance,
+ subclassify,
+ exactmatch,
+ neyman,
+ help.matchit,
+ print.matchit,
+ print.summary.matchit,
+ print.summary.neyman,
+ plot.matchit,
+ summary.matchit,
+ summary.neyman
+ )
+
+S3method(print, matchit)
+S3method(print, summary.matchit)
+S3method(print, summary.neyman)
+S3method(plot, matchit)
+S3method(summary, matchit)
+S3method(summary, neyman)
+
diff --git a/R/diagnose.R b/R/diagnose.R
new file mode 100644
index 0000000..e268322
--- /dev/null
+++ b/R/diagnose.R
@@ -0,0 +1,82 @@
+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)
+}
+
diff --git a/R/distance.R b/R/distance.R
new file mode 100644
index 0000000..bac2983
--- /dev/null
+++ b/R/distance.R
@@ -0,0 +1,136 @@
+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)
+}
diff --git a/R/exactmatch.R b/R/exactmatch.R
new file mode 100644
index 0000000..75756dc
--- /dev/null
+++ b/R/exactmatch.R
@@ -0,0 +1,25 @@
+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)
+}
diff --git a/R/help.matchit.R b/R/help.matchit.R
new file mode 100644
index 0000000..014f9e1
--- /dev/null
+++ b/R/help.matchit.R
@@ -0,0 +1,53 @@
+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("")
+ }
+}
diff --git a/R/match.data.R b/R/match.data.R
new file mode 100644
index 0000000..e5a6689
--- /dev/null
+++ b/R/match.data.R
@@ -0,0 +1,12 @@
+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.")
+}
diff --git a/R/matchdef.R b/R/matchdef.R
new file mode 100644
index 0000000..8af675a
--- /dev/null
+++ b/R/matchdef.R
@@ -0,0 +1,215 @@
+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]
+
+ ## 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 <- as.data.frame(matrix(0,n1,ratio), row.names = 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
+
+ ## optimal ratio matching
+ if (m.order==1) {
+ cat("Optimal Matching Treated: ")
+ distance <- matrix(0, ncol=n0, nrow=n1)
+ rownames(distance) <- row.names(treata)[treat==1]
+ colnames(distance) <- row.names(treata)[treat==0]
+ for (i in 1:n1)
+ distance[i,] <- abs(p1[i]-p0)
+ full <- as.matrix(fullmatch(distance, min.controls = ratio,
+ max.controls = ratio,
+ omit.fraction = (n0-ratio*n1)/n0,...))
+ psclass <- full[pmatch(row.names(treata), row.names(full)),]
+ psclass <- as.numeric(as.factor(psclass))
+ names(psclass) <- row.names(treata)
+ for (i in 1:n1) {
+ match.matrix[i,] <-
+ names(which(psclass[tlabels[i]]==
+ psclass[-pmatch(tlabels[i],
+ names(psclass))]))
+ }
+ if(counter==TRUE){cat("Done\n")}
+
+ } else if(nearest) {
+ ## 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
+ ##in cases there's no replacement and all controls have been used up
+ if(!0%in%matchedc2){
+ 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)
+}
diff --git a/R/matchit.R b/R/matchit.R
new file mode 100644
index 0000000..faf5900
--- /dev/null
+++ b/R/matchit.R
@@ -0,0 +1,174 @@
+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, full = FALSE, full.options=list(), ...){
+ cl <- match.call()
+ if (m.order==1 | full)
+ require(optmatch)
+
+ #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) | identical(m.order,1))){
+ 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}
+ #full matching
+ if(!(identical(full,TRUE) | identical(full,FALSE))){
+ warning("full=",cl$full,") is invalid; used full=FALSE instead",call. = FALSE)
+ full <- FALSE
+ } else if(full & nearest){
+ warning("Full matching will ignore nearest neighbor inputs",call. = FALSE)
+ nearest <- F
+ }
+
+ # 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 <- mf$full <- mf$full.options <- 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, full = full, full.options=full.options)
+
+ 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)
+ }
+}
diff --git a/R/neyman.R b/R/neyman.R
new file mode 100644
index 0000000..f1f1bf9
--- /dev/null
+++ b/R/neyman.R
@@ -0,0 +1,164 @@
+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
+}
diff --git a/R/plot.matchit.R b/R/plot.matchit.R
new file mode 100644
index 0000000..b3807e5
--- /dev/null
+++ b/R/plot.matchit.R
@@ -0,0 +1,203 @@
+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]))}
+ }
+ }
+
+ }
+
+ }
+
diff --git a/R/print.matchit.R b/R/print.matchit.R
new file mode 100644
index 0000000..eb9fae4
--- /dev/null
+++ b/R/print.matchit.R
@@ -0,0 +1,146 @@
+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,psc = NULL, full = F){
+ 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)
+ if(!full){
+ 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)
+ if (is.null(x$call$exact))
+ xsum[1,4] <- -1*t.test(xx~tt)$statistic
+ else if(is.logical((eval(x$call$exact))))
+ if(eval(x$call$exact))
+ xsum[1,4] <- 0
+ else
+ xsum[1,4] <- -1*t.test(xx~tt)$statistic
+ else
+ xsum[1,4] <- -1*t.test(xx~tt)$statistic
+ 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])
+ } else {
+ xn[2,] <- xn[1,] #??
+ m0 <- lm(xx~ tt)
+ m1 <- lm(xx~ tt + as.factor(psc))
+ tdat <- m1$model
+ tdat$tt <- 1
+ xsum[2,1] <- mean(predict(m1,newdata=tdat))
+ tdat$tt <- 0
+ xsum[2,2] <- mean(predict(m1,newdata=tdat))
+ xsum[1,3] <- xsum[2,3] <- sd(xx)
+ xsum[1,4] <- summary(m0)$coef[2,3]
+ xsum[1,5] <- (xsum[1,1]-xsum[1,2])/sd(x1)
+ xsum[2,4] <- summary(m1)$coef[2,3]
+ xsum[2,5] <- (xsum[2,1]-xsum[2,2])/sd(x1)
+ }
+ 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) & !identical(eval(x$call$full),TRUE))
+ {
+ 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 if(identical(eval(x$call$full),TRUE)){
+ cat("\nSummary of propensity score for full and matched samples:\n\n")
+ xqoi <- qoi(pscore, treat, x$psweights,
+ psc = x$data$psclass, full = T)
+ print.data.frame(xqoi$xsum, digits = digits)
+ cat("\nSample sizes:\n\n")
+ print.data.frame(xqoi$xn, 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) & !identical(eval(x$call$full),TRUE)) {
+ 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)
+ }
+ }
+ }
diff --git a/R/print.neyman.R b/R/print.neyman.R
new file mode 100644
index 0000000..aa50268
--- /dev/null
+++ b/R/print.neyman.R
@@ -0,0 +1,9 @@
+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")
+}
diff --git a/R/print.summary.matchit.R b/R/print.summary.matchit.R
new file mode 100644
index 0000000..58ad985
--- /dev/null
+++ b/R/print.summary.matchit.R
@@ -0,0 +1,61 @@
+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) | identical(eval(x$call$full),TRUE))
+ {
+ 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) & !identical(eval(x$call$full),TRUE))
+ {
+ 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)
+}
diff --git a/R/print.summary.neyman.R b/R/print.summary.neyman.R
new file mode 100644
index 0000000..cffbc07
--- /dev/null
+++ b/R/print.summary.neyman.R
@@ -0,0 +1,11 @@
+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")
+}
diff --git a/R/subclassify.R b/R/subclassify.R
new file mode 100644
index 0000000..550fb08
--- /dev/null
+++ b/R/subclassify.R
@@ -0,0 +1,124 @@
+subclassify <- function(formula,data,in.sample,pscore,nearest=TRUE,
+ match.matrix,subclass=0,sub.by="treat",
+ counter=TRUE, full = FALSE, full.options=list()){
+
+ data <- eval(data,parent.frame())
+ treata <- model.frame(formula,data)[,1,drop=FALSE]
+ treat <- as.vector(treata[,1])
+ names(treat) <- row.names(treata)
+
+ if(full) { # full matching with propensity score
+ if(counter){cat("Full Matching...")}
+ full <- full.options
+ if(!is.list(full.options)){
+ warning("full.options must be a list; assuming defaults for full matching",call.=FALSE)
+ }
+ if(is.null(full$min.controls)){
+ full$min.controls <- 0
+ }
+ if(is.null(full$max.controls)){
+ full$max.controls <- Inf
+ }
+ if(is.null(full$omit.fraction)){
+ full$omit.fraction <- NULL
+ }
+ if(is.null(full$tol)){
+ full$tol <- 0.01
+ }
+ if(is.null(full$subclass.indices)){
+ full$subclass.indices <- NULL
+ }
+ notin <- names(full.options)[which(!names(full.options)%in%c("min.controls","max.controls",
+ "omit.fraction","omit.fraction",
+ "tol", "subclass.indices"))]
+ if(!is.null(notin)){
+ warning(paste(notin,collapse=" "), " in full.options invalid and ignored for full matching",call.=FALSE)
+ }
+ n1 <- length(treat[treat==1])
+ n0 <- length(treat[treat==0])
+ p1 <- pscore[treat==1]
+ p0 <- pscore[treat==0]
+ distance <- matrix(0, ncol=n0, nrow=n1)
+ rownames(distance) <- row.names(treata)[treat==1]
+ colnames(distance) <- row.names(treata)[treat==0]
+ for (i in 1:n1)
+ distance[i,] <- abs(p1[i]-p0)
+ full <- as.matrix(fullmatch(distance,subclass.indices = full$subclass.indices,
+ min.controls = full$min.controls,
+ max.controls = full$max.controls,
+ omit.fraction = full$omit.fraction,
+ tol = full$tol))
+ psclass <- full[pmatch(row.names(treata), row.names(full)),]
+ psclass <- as.numeric(as.factor(psclass))
+ names(psclass) <- row.names(treata)
+ q <- NULL
+ if(counter){cat("Done\n")}
+ }
+ else if(subclass) {
+ 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)}
+ }
+
+ 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}
+
+ return(list(psclass = psclass, q.cut = q))
+}
diff --git a/R/summary.matchit.R b/R/summary.matchit.R
new file mode 100644
index 0000000..4dbd586
--- /dev/null
+++ b/R/summary.matchit.R
@@ -0,0 +1,220 @@
+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, psc = NULL, full = F){
+ 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)
+ if(!full){
+ 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])
+ } else {
+ xn[2,] <- xn[1,] #??
+ m0 <- lm(xx~ tt)
+ m1 <- lm(xx~ tt + as.factor(psc))
+ tdat <- m1$model
+ tdat$tt <- 1
+ xsum[2,1] <- mean(predict(m1,newdata=tdat))
+ tdat$tt <- 0
+ xsum[2,2] <- mean(predict(m1,newdata=tdat))
+ xsum[1,3] <- xsum[2,3] <- sd(xx)
+ xsum[1,4] <- summary(m0)$coef[2,3]
+ xsum[1,5] <- (xsum[1,1]-xsum[1,2])/sd(x1)
+ xsum[2,4] <- summary(m1)$coef[2,3]
+ xsum[2,5] <- (xsum[2,1]-xsum[2,2])/sd(x1)
+ }
+ 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)
+ }
+
+ if(!identical(eval(object$call$full),TRUE)){
+ aa <- apply(covariates,2,qoi,tt=treat,ww=weights)
+ } else {
+ aa <- apply(covariates,2,qoi,tt=treat,ww=weights,
+ psc = object$data$psclass, full = T)
+ }
+ 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])
+ if(!identical(eval(object$call$full),TRUE)){
+ jqoi <- qoi(x2,tt=treat,ww=weights)
+ } else {
+ jqoi <- qoi(x2,tt=treat,ww=weights,
+ psc = object$data$psclass, full = T)
+ }
+ sum.all.int <- rbind(sum.all.int,jqoi$xsum[1,])
+ Reduction <- (abs(sum.all.int[nrow(sum.all.int),5])-abs(jqoi$xsum[2,5]))>0
+ sum.matched.int <- rbind(sum.matched.int,cbind(jqoi$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) & !identical(eval(object$call$full),TRUE))
+ {
+ 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
+}
+
diff --git a/R/summary.neyman.R b/R/summary.neyman.R
new file mode 100644
index 0000000..b8be27a
--- /dev/null
+++ b/R/summary.neyman.R
@@ -0,0 +1,14 @@
+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
+}
diff --git a/data/incomedata.tab b/data/incomedata.tab
new file mode 100644
index 0000000..64dd46c
--- /dev/null
+++ b/data/incomedata.tab
@@ -0,0 +1,201 @@
+"treat" "income74" "income75" "male" "educ" "income78"
+1 1441 2488 1 11 3561
+1 1869 1703 0 9 2579
+1 1516 2258 1 7 3714
+1 1690 2093 0 10 2806
+1 2336 3105 1 13 3423
+1 2154 2426 0 9 2471
+1 1653 2006 0 9 3777
+1 2252 3262 0 12 2686
+1 1589 1695 0 11 3605
+1 1763 3223 0 12 2433
+1 1950 2913 0 8 3075
+1 2017 2166 1 11 2914
+1 2135 2900 1 10 3574
+1 2000 1858 0 13 2700
+1 2228 2592 1 12 3439
+1 2052 1975 1 12 2593
+1 2174 1452 1 10 3734
+1 2922 3112 1 13 2915
+1 1493 2374 1 15 2975
+1 1543 2219 1 14 2926
+1 1502 3039 1 12 3141
+1 2503 3374 1 11 2501
+1 1776 2440 0 10 2632
+1 1635 1721 1 12 3007
+1 2140 2279 0 12 2261
+1 2574 2378 1 7 2717
+1 1285 2321 1 10 4377
+1 1689 2327 1 9 3197
+1 2041 2799 1 7 3696
+1 2854 1988 0 9 2570
+1 2695 1105 1 13 2196
+1 2490 2053 1 11 3556
+1 1861 1990 1 13 2372
+1 1419 1300 1 11 2911
+1 1237 812 1 12 3602
+1 1925 2585 0 8 3125
+1 2321 2985 1 12 3263
+1 1893 3209 1 14 2643
+1 2908 2462 1 11 2464
+1 2204 2186 0 9 2606
+1 2011 1768 0 9 3112
+1 1883 2510 1 14 2856
+1 2281 1955 1 11 3553
+1 1283 2837 0 10 2782
+1 1681 2097 0 9 2651
+1 2370 2767 1 12 3460
+1 2041 1939 0 10 2902
+1 1769 1670 1 12 3427
+1 1589 2605 0 12 2320
+1 2432 1902 1 10 3104
+1 1583 2005 1 14 3277
+1 2320 2719 1 14 3169
+1 2034 2151 0 13 2051
+1 2099 1958 0 6 3049
+1 2616 2801 0 14 2832
+1 2049 2546 0 13 2937
+1 2544 1694 0 7 3061
+1 672 2341 0 14 2505
+1 1284 2741 0 11 3410
+1 2104 2399 0 10 3655
+1 1479 2128 1 10 2773
+1 2360 2582 1 10 3180
+1 1809 2350 0 9 2654
+1 1378 2093 1 10 3356
+1 1851 2296 1 14 2126
+1 2953 1796 0 13 3351
+1 2485 1692 0 9 2497
+1 2420 2222 0 9 3512
+1 2197 1839 1 11 3726
+1 2779 2448 0 11 2524
+1 1132 2116 1 9 2804
+1 1641 3345 1 11 2811
+1 1436 2188 0 10 2766
+1 2422 2097 0 14 2781
+1 2202 2251 0 12 3580
+1 1140 2018 1 14 3126
+1 2038 2106 0 8 3440
+1 2076 1941 1 12 2670
+1 1839 1290 0 12 3067
+1 1853 1800 0 11 1986
+1 2721 1568 0 10 3238
+1 2190 2231 0 9 3047
+1 3164 2588 0 13 3423
+1 2136 2304 0 11 2546
+1 3143 2774 0 10 2152
+1 2187 2496 1 8 2847
+1 1945 2259 1 14 3251
+1 1089 2017 0 13 3593
+1 1964 2234 1 9 4722
+1 1817 1685 0 9 3707
+1 1765 2393 0 12 2598
+1 2103 2611 1 9 3456
+1 1827 2130 0 8 2921
+1 1521 2960 1 7 2674
+1 2026 2426 1 10 2634
+1 1594 3443 1 7 2739
+1 2720 1741 0 10 2559
+1 2146 2509 1 14 2578
+1 1444 2091 0 9 3868
+1 2551 2230 0 10 3477
+0 2721 2126 0 10 2072
+0 2305 1725 1 11 2587
+0 1442 2435 0 9 2631
+0 1836 2866 0 11 3295
+0 2471 1876 1 9 2453
+0 2176 2267 0 12 2770
+0 1735 2628 1 11 1559
+0 1856 1771 0 10 2131
+0 1860 2679 1 10 2019
+0 1810 1889 0 13 2600
+0 2299 1837 1 10 1429
+0 2391 2793 0 11 2504
+0 1878 1761 0 10 3778
+0 2800 2567 0 10 2204
+0 1729 2727 1 13 2545
+0 2388 2036 1 9 2920
+0 1509 2882 1 8 2572
+0 1878 1765 0 9 2514
+0 1997 2881 0 10 1954
+0 2327 2182 0 6 1726
+0 1385 2557 0 10 3083
+0 2172 2266 1 10 2879
+0 2435 2715 0 11 2989
+0 979 3556 1 9 1859
+0 1904 2371 1 9 2522
+0 2817 2793 1 13 2760
+0 2641 1462 0 12 2599
+0 1915 2634 1 7 2519
+0 2108 2338 0 8 3242
+0 2920 2026 1 10 1681
+0 1599 2709 1 7 1949
+0 2042 1808 0 12 2864
+0 2739 1805 1 10 2639
+0 2398 1906 0 16 2156
+0 1674 2573 0 13 1819
+0 1809 2233 0 8 3175
+0 1871 1874 0 11 1819
+0 2030 2033 1 11 2607
+0 1129 2427 1 10 2402
+0 3068 2391 1 6 3093
+0 1679 1418 0 11 2612
+0 1472 2432 1 9 2804
+0 1792 3312 0 11 2316
+0 2220 1786 1 11 2999
+0 2263 1701 0 11 3088
+0 2414 1845 1 10 2745
+0 2233 1954 1 10 2695
+0 1152 2250 0 8 2566
+0 2083 1785 0 10 2588
+0 575 1968 0 6 2419
+0 2161 1969 1 6 2533
+0 2036 2178 0 9 2508
+0 1701 1993 1 12 2520
+0 2263 1640 0 10 2993
+0 2391 2347 1 9 3449
+0 1254 1579 1 13 2639
+0 2415 2403 1 7 3172
+0 1848 2482 0 13 3709
+0 3216 2329 0 14 1494
+0 2276 1493 1 13 2635
+0 1767 2513 0 12 3317
+0 2395 2826 0 9 1685
+0 1773 2060 0 11 2708
+0 2723 2016 0 8 2394
+0 2523 2664 1 9 1988
+0 1587 1837 0 9 2407
+0 2253 2117 0 13 2160
+0 1504 2267 0 11 3059
+0 1746 3184 1 12 2491
+0 1941 2292 1 7 2275
+0 2411 2679 1 14 3262
+0 3127 2180 1 13 2383
+0 1801 3165 0 12 2278
+0 2120 2393 0 9 1754
+0 2074 2602 0 7 2025
+0 1793 1595 1 8 2139
+0 2399 2331 0 12 2674
+0 2043 2454 0 10 2398
+0 2477 2171 0 10 2543
+0 2061 3127 1 9 1783
+0 2182 2148 0 9 3201
+0 2355 2819 0 10 2712
+0 2063 2427 1 11 3145
+0 2250 1879 1 6 1759
+0 1332 1859 1 12 2037
+0 1198 2314 0 9 1635
+0 1504 2069 0 12 2763
+0 2482 2873 0 9 2728
+0 1257 2179 0 9 1303
+0 1553 2339 0 13 1741
+0 2368 1731 0 14 2578
+0 2432 2622 1 11 1864
+0 1866 2699 1 8 2979
+0 1562 2001 0 10 2471
+0 1856 2778 0 8 2579
+0 1572 1882 1 8 1515
+0 2636 2742 0 10 2160
+0 1453 2282 1 15 2825
+0 1558 2190 0 13 3293
+0 1941 2395 1 13 2832
diff --git a/data/lalonde.tab b/data/lalonde.tab
new file mode 100644
index 0000000..47d6b3c
--- /dev/null
+++ b/data/lalonde.tab
@@ -0,0 +1,615 @@
+"treat" "age" "educ" "black" "hispan" "married" "nodegree" "re74" "re75" "re78"
+"NSW1" 1 37 11 1 0 1 1 0 0 9930.046
+"NSW2" 1 22 9 0 1 0 1 0 0 3595.894
+"NSW3" 1 30 12 1 0 0 0 0 0 24909.45
+"NSW4" 1 27 11 1 0 0 1 0 0 7506.146
+"NSW5" 1 33 8 1 0 0 1 0 0 289.7899
+"NSW6" 1 22 9 1 0 0 1 0 0 4056.494
+"NSW7" 1 23 12 1 0 0 0 0 0 0
+"NSW8" 1 32 11 1 0 0 1 0 0 8472.158
+"NSW9" 1 22 16 1 0 0 0 0 0 2164.022
+"NSW10" 1 33 12 0 0 1 0 0 0 12418.07
+"NSW11" 1 19 9 1 0 0 1 0 0 8173.908
+"NSW12" 1 21 13 1 0 0 0 0 0 17094.64
+"NSW13" 1 18 8 1 0 0 1 0 0 0
+"NSW14" 1 27 10 1 0 1 1 0 0 18739.93
+"NSW15" 1 17 7 1 0 0 1 0 0 3023.879
+"NSW16" 1 19 10 1 0 0 1 0 0 3228.503
+"NSW17" 1 27 13 1 0 0 0 0 0 14581.86
+"NSW18" 1 23 10 1 0 0 1 0 0 7693.4
+"NSW19" 1 40 12 1 0 0 0 0 0 10804.32
+"NSW20" 1 26 12 1 0 0 0 0 0 10747.35
+"NSW21" 1 23 11 1 0 0 1 0 0 0
+"NSW22" 1 41 14 0 0 0 0 0 0 5149.501
+"NSW23" 1 38 9 0 0 0 1 0 0 6408.95
+"NSW24" 1 24 11 1 0 0 1 0 0 1991.4
+"NSW25" 1 18 10 1 0 0 1 0 0 11163.17
+"NSW26" 1 29 11 1 0 1 1 0 0 9642.999
+"NSW27" 1 25 11 1 0 0 1 0 0 9897.049
+"NSW28" 1 27 10 0 1 0 1 0 0 11142.87
+"NSW29" 1 17 10 1 0 0 1 0 0 16218.04
+"NSW30" 1 24 11 1 0 0 1 0 0 995.7002
+"NSW31" 1 17 10 1 0 0 1 0 0 0
+"NSW32" 1 48 4 1 0 0 1 0 0 6551.592
+"NSW33" 1 25 11 1 0 1 1 0 0 1574.424
+"NSW34" 1 20 12 1 0 0 0 0 0 0
+"NSW35" 1 25 12 1 0 0 0 0 0 3191.753
+"NSW36" 1 42 14 1 0 0 0 0 0 20505.93
+"NSW37" 1 25 5 1 0 0 1 0 0 6181.88
+"NSW38" 1 23 12 1 0 1 0 0 0 5911.551
+"NSW39" 1 46 8 1 0 1 1 0 0 3094.156
+"NSW40" 1 24 10 1 0 0 1 0 0 0
+"NSW41" 1 21 12 1 0 0 0 0 0 1254.582
+"NSW42" 1 19 9 0 0 0 1 0 0 13188.83
+"NSW43" 1 17 8 1 0 0 1 0 0 8061.485
+"NSW44" 1 18 8 0 1 1 1 0 0 2787.96
+"NSW45" 1 20 11 1 0 0 1 0 0 3972.54
+"NSW46" 1 25 11 1 0 1 1 0 0 0
+"NSW47" 1 17 8 1 0 0 1 0 0 0
+"NSW48" 1 17 9 1 0 0 1 0 0 0
+"NSW49" 1 25 5 1 0 0 1 0 0 12187.41
+"NSW50" 1 23 12 1 0 0 0 0 0 4843.176
+"NSW51" 1 28 8 1 0 0 1 0 0 0
+"NSW52" 1 31 11 1 0 1 1 0 0 8087.487
+"NSW53" 1 18 11 1 0 0 1 0 0 0
+"NSW54" 1 25 12 1 0 0 0 0 0 2348.973
+"NSW55" 1 30 11 1 0 1 1 0 0 590.7818
+"NSW56" 1 17 10 1 0 0 1 0 0 0
+"NSW57" 1 37 9 1 0 0 1 0 0 1067.506
+"NSW58" 1 41 4 1 0 1 1 0 0 7284.986
+"NSW59" 1 42 14 1 0 1 0 0 0 13167.52
+"NSW60" 1 22 11 0 0 0 1 0 0 1048.432
+"NSW61" 1 17 8 1 0 0 1 0 0 0
+"NSW62" 1 29 8 1 0 0 1 0 0 1923.938
+"NSW63" 1 35 10 1 0 0 1 0 0 4666.236
+"NSW64" 1 27 11 1 0 0 1 0 0 549.2984
+"NSW65" 1 29 4 1 0 0 1 0 0 762.9146
+"NSW66" 1 28 9 1 0 0 1 0 0 10694.29
+"NSW67" 1 27 11 1 0 0 1 0 0 0
+"NSW68" 1 23 7 0 0 0 1 0 0 0
+"NSW69" 1 45 5 1 0 1 1 0 0 8546.715
+"NSW70" 1 29 13 1 0 0 0 0 0 7479.656
+"NSW71" 1 27 9 1 0 0 1 0 0 0
+"NSW72" 1 46 13 1 0 0 0 0 0 647.2046
+"NSW73" 1 18 6 1 0 0 1 0 0 0
+"NSW74" 1 25 12 1 0 0 0 0 0 11965.81
+"NSW75" 1 28 15 1 0 0 0 0 0 9598.541
+"NSW76" 1 25 11 0 0 0 1 0 0 18783.35
+"NSW77" 1 22 12 1 0 0 0 0 0 18678.08
+"NSW78" 1 21 9 1 0 0 1 0 0 0
+"NSW79" 1 40 11 1 0 0 1 0 0 23005.6
+"NSW80" 1 22 11 1 0 0 1 0 0 6456.697
+"NSW81" 1 25 12 1 0 0 0 0 0 0
+"NSW82" 1 18 12 1 0 0 0 0 0 2321.107
+"NSW83" 1 38 12 0 0 0 0 0 0 4941.849
+"NSW84" 1 27 13 1 0 0 0 0 0 0
+"NSW85" 1 27 8 1 0 0 1 0 0 0
+"NSW86" 1 38 11 1 0 0 1 0 0 0
+"NSW87" 1 23 8 0 1 0 1 0 0 3881.284
+"NSW88" 1 26 11 1 0 0 1 0 0 17230.96
+"NSW89" 1 21 12 0 0 0 0 0 0 8048.603
+"NSW90" 1 25 8 1 0 0 1 0 0 0
+"NSW91" 1 31 11 1 0 1 1 0 0 14509.93
+"NSW92" 1 17 10 1 0 0 1 0 0 0
+"NSW93" 1 25 11 1 0 0 1 0 0 0
+"NSW94" 1 21 12 1 0 0 0 0 0 9983.784
+"NSW95" 1 44 11 1 0 0 1 0 0 0
+"NSW96" 1 25 12 0 0 0 0 0 0 5587.503
+"NSW97" 1 18 9 1 0 0 1 0 0 4482.845
+"NSW98" 1 42 12 1 0 0 0 0 0 2456.153
+"NSW99" 1 25 10 1 0 0 1 0 0 0
+"NSW100" 1 31 9 0 1 0 1 0 0 26817.6
+"NSW101" 1 24 10 1 0 0 1 0 0 0
+"NSW102" 1 26 10 1 0 0 1 0 0 9265.788
+"NSW103" 1 25 11 1 0 0 1 0 0 485.2298
+"NSW104" 1 18 11 1 0 0 1 0 0 4814.627
+"NSW105" 1 19 11 1 0 0 1 0 0 7458.105
+"NSW106" 1 43 9 1 0 0 1 0 0 0
+"NSW107" 1 27 13 1 0 0 0 0 0 34099.28
+"NSW108" 1 17 9 1 0 0 1 0 0 1953.268
+"NSW109" 1 30 11 1 0 0 1 0 0 0
+"NSW110" 1 26 10 1 0 1 1 2027.999 0 0
+"NSW111" 1 20 9 1 0 0 1 6083.994 0 8881.665
+"NSW112" 1 17 9 0 1 0 1 445.1704 74.34345 6210.67
+"NSW113" 1 20 12 1 0 0 0 989.2678 165.2077 0
+"NSW114" 1 18 11 1 0 0 1 858.2543 214.5636 929.8839
+"NSW115" 1 27 12 1 0 1 0 3670.872 334.0493 0
+"NSW116" 1 21 12 0 0 0 0 3670.872 334.0494 12558.02
+"NSW117" 1 27 12 1 0 0 0 2143.413 357.9499 22163.25
+"NSW118" 1 20 12 1 0 0 0 0 377.5686 1652.637
+"NSW119" 1 19 10 1 0 0 1 0 385.2741 8124.715
+"NSW120" 1 23 12 1 0 0 0 5506.308 501.0741 671.3318
+"NSW121" 1 29 14 1 0 0 0 0 679.6734 17814.98
+"NSW122" 1 18 10 1 0 0 1 0 798.9079 9737.154
+"NSW123" 1 19 9 1 0 0 1 0 798.9079 17685.18
+"NSW124" 1 27 13 0 0 1 0 9381.566 853.7225 0
+"NSW125" 1 18 11 0 0 0 1 3678.231 919.5579 4321.705
+"NSW126" 1 27 9 1 0 1 1 0 934.4454 1773.423
+"NSW127" 1 22 12 1 0 0 0 5605.852 936.1773 0
+"NSW128" 1 23 10 1 0 1 1 0 936.4386 11233.26
+"NSW129" 1 23 12 0 1 0 0 9385.74 1117.439 559.4432
+"NSW130" 1 20 11 1 0 0 1 3637.498 1220.836 1085.44
+"NSW131" 1 17 9 1 0 0 1 1716.509 1253.439 5445.2
+"NSW132" 1 28 11 1 0 0 1 0 1284.079 60307.93
+"NSW133" 1 26 11 1 0 1 1 0 1392.853 1460.36
+"NSW134" 1 20 11 1 0 0 1 16318.62 1484.994 6943.342
+"NSW135" 1 24 11 1 0 1 1 824.3886 1666.113 4032.708
+"NSW136" 1 31 9 1 0 0 1 0 1698.607 10363.27
+"NSW137" 1 23 8 0 0 1 1 0 1713.15 4232.309
+"NSW138" 1 18 10 1 0 0 1 2143.411 1784.274 11141.39
+"NSW139" 1 29 12 1 0 0 0 10881.94 1817.284 0
+"NSW140" 1 26 11 0 0 0 1 0 2226.266 13385.86
+"NSW141" 1 24 9 1 0 0 1 9154.7 2288.675 4849.559
+"NSW142" 1 25 12 1 0 0 0 14426.79 2409.274 0
+"NSW143" 1 24 10 1 0 0 1 4250.402 2421.947 1660.508
+"NSW144" 1 46 8 1 0 0 1 3165.658 2594.723 0
+"NSW145" 1 31 12 0 0 0 0 0 2611.218 2484.549
+"NSW146" 1 19 11 1 0 0 1 2305.026 2615.276 4146.603
+"NSW147" 1 19 8 1 0 0 1 0 2657.057 9970.681
+"NSW148" 1 27 11 1 0 0 1 2206.94 2666.274 0
+"NSW149" 1 26 11 1 0 1 1 0 2754.646 26372.28
+"NSW150" 1 20 10 1 0 0 1 5005.731 2777.355 5615.189
+"NSW151" 1 28 10 1 0 0 1 0 2836.506 3196.571
+"NSW152" 1 24 12 1 0 0 0 13765.75 2842.764 6167.681
+"NSW153" 1 19 8 1 0 0 1 2636.353 2937.264 7535.942
+"NSW154" 1 23 12 1 0 0 0 6269.341 3039.96 8484.239
+"NSW155" 1 42 9 1 0 1 1 0 3058.531 1294.409
+"NSW156" 1 25 13 1 0 0 0 12362.93 3090.732 0
+"NSW157" 1 18 9 1 0 0 1 0 3287.375 5010.342
+"NSW158" 1 21 12 1 0 0 0 6473.683 3332.409 9371.037
+"NSW159" 1 27 10 1 0 0 1 1001.146 3550.075 0
+"NSW160" 1 21 8 1 0 0 1 989.2678 3695.897 4279.613
+"NSW161" 1 22 9 1 0 0 1 2192.877 3836.986 3462.564
+"NSW162" 1 31 4 1 0 0 1 8517.589 4023.211 7382.549
+"NSW163" 1 24 10 1 0 1 1 11703.2 4078.152 0
+"NSW164" 1 29 10 1 0 0 1 0 4398.95 0
+"NSW165" 1 29 12 1 0 0 0 9748.387 4878.937 10976.51
+"NSW166" 1 19 10 0 0 0 1 0 5324.109 13829.62
+"NSW167" 1 19 11 0 1 1 1 5424.485 5463.803 6788.463
+"NSW168" 1 31 9 1 0 0 1 10717.03 5517.841 9558.501
+"NSW169" 1 22 10 1 0 1 1 1468.348 5588.664 13228.28
+"NSW170" 1 21 9 1 0 0 1 6416.47 5749.331 743.6666
+"NSW171" 1 17 10 1 0 0 1 1291.468 5793.852 5522.788
+"NSW172" 1 26 12 1 0 1 0 8408.762 5794.831 1424.944
+"NSW173" 1 20 9 0 1 0 1 12260.78 5875.049 1358.643
+"NSW174" 1 19 10 1 0 0 1 4121.949 6056.754 0
+"NSW175" 1 26 10 1 0 0 1 25929.68 6788.958 672.8773
+"NSW176" 1 28 11 1 0 0 1 1929.029 6871.856 0
+"NSW177" 1 22 12 0 1 1 0 492.2305 7055.702 10092.83
+"NSW178" 1 33 11 1 0 0 1 0 7867.916 6281.433
+"NSW179" 1 22 12 0 0 0 0 6759.994 8455.504 12590.71
+"NSW180" 1 29 10 0 1 0 1 0 8853.674 5112.014
+"NSW181" 1 33 12 1 0 1 0 20279.95 10941.35 15952.6
+"NSW182" 1 25 14 1 0 1 0 35040.07 11536.57 36646.95
+"NSW183" 1 35 9 1 0 1 1 13602.43 13830.64 12803.97
+"NSW184" 1 35 8 1 0 1 1 13732.07 17976.15 3786.628
+"NSW185" 1 33 11 1 0 1 1 14660.71 25142.24 4181.942
+"PSID1" 0 30 12 0 0 1 0 20166.73 18347.23 25564.67
+"PSID2" 0 26 12 0 0 1 0 25862.32 17806.55 25564.67
+"PSID3" 0 25 16 0 0 1 0 25862.32 15316.21 25564.67
+"PSID4" 0 42 11 0 0 1 1 21787.05 14265.29 15491.01
+"PSID5" 0 25 9 1 0 1 1 14829.69 13776.53 0
+"PSID6" 0 37 9 1 0 1 1 13685.48 12756.05 17833.2
+"PSID7" 0 32 12 0 0 1 0 19067.58 12625.35 14146.28
+"PSID8" 0 20 12 1 0 0 0 7392.314 12396.19 17765.23
+"PSID9" 0 38 9 0 1 1 1 16826.18 12029.18 0
+"PSID10" 0 39 10 0 0 1 1 16767.41 12022.02 4433.18
+"PSID11" 0 41 5 0 0 1 1 10785.76 11991.58 19451.31
+"PSID12" 0 31 14 0 0 1 0 17831.29 11563.69 22094.97
+"PSID13" 0 34 8 0 0 1 1 8038.872 11404.35 5486.799
+"PSID14" 0 29 12 0 0 1 0 14768.95 11146.55 6420.722
+"PSID15" 0 22 14 1 0 1 0 748.4399 11105.37 18208.55
+"PSID16" 0 42 0 0 1 1 1 2797.833 10929.92 9922.934
+"PSID17" 0 25 9 0 1 0 1 5460.477 10589.76 7539.361
+"PSID18" 0 28 9 0 0 1 1 11091.41 10357.02 15406.78
+"PSID19" 0 40 13 0 0 1 0 3577.621 10301.52 11911.95
+"PSID20" 0 35 9 0 0 1 1 11475.43 9397.403 11087.38
+"PSID21" 0 27 10 0 1 1 1 15711.36 9098.419 17023.41
+"PSID22" 0 27 6 0 1 1 1 7831.189 9071.565 5661.171
+"PSID23" 0 36 12 0 0 1 0 25535.12 8695.597 21905.82
+"PSID24" 0 47 8 1 0 1 1 9275.169 8543.419 0
+"PSID25" 0 40 11 0 0 1 1 20666.35 8502.242 25564.67
+"PSID26" 0 27 7 0 0 1 1 3064.293 8461.065 11149.45
+"PSID27" 0 36 9 1 0 1 1 13256.4 8457.484 0
+"PSID28" 0 39 6 0 1 1 1 13279.91 8441.371 25048.94
+"PSID29" 0 21 9 0 0 1 1 11156.07 8441.371 1213.214
+"PSID30" 0 29 12 0 0 1 0 11199.17 8081.516 0
+"PSID31" 0 22 13 0 1 0 0 6404.843 7882.79 9453.017
+"PSID32" 0 25 10 0 0 1 1 13634.54 7793.274 11688.82
+"PSID33" 0 27 12 0 0 1 0 12270.89 7709.129 7806.829
+"PSID34" 0 45 8 0 0 1 1 22415.97 7635.726 15931.37
+"PSID35" 0 26 12 0 0 1 0 2345.242 7565.903 2838.713
+"PSID36" 0 27 12 0 0 1 0 9788.497 7496.081 14038.4
+"PSID37" 0 33 8 0 0 1 1 12312.03 7474.597 25514.43
+"PSID38" 0 25 12 0 0 1 0 11381.38 7467.435 4162.756
+"PSID39" 0 49 8 0 0 1 1 6459.703 7431.629 7503.896
+"PSID40" 0 40 3 0 1 1 1 7576.485 7426.258 12104.06
+"PSID41" 0 22 12 1 0 1 0 9729.719 7372.548 2231.367
+"PSID42" 0 25 5 0 0 1 1 7891.927 7293.774 14617.67
+"PSID43" 0 25 12 0 0 1 0 11516.57 7263.339 19588.74
+"PSID44" 0 21 12 0 0 1 0 13601.23 7202.468 10746.03
+"PSID45" 0 33 9 0 1 1 1 11959.36 7087.887 25564.67
+"PSID46" 0 20 12 1 0 1 0 9555.344 7055.661 0
+"PSID47" 0 19 11 0 0 1 1 4306.468 6978.677 837.871
+"PSID48" 0 25 12 1 0 1 0 295.8493 6942.871 461.0507
+"PSID49" 0 29 12 0 0 1 0 15303.83 6932.129 24290.87
+"PSID50" 0 20 12 0 0 1 0 3558.029 6797.855 6680.802
+"PSID51" 0 29 6 0 1 1 1 8542.403 6701.177 7196.528
+"PSID52" 0 25 13 0 0 1 0 19259.59 6652.839 13015.82
+"PSID53" 0 41 15 0 0 1 0 25862.32 6563.323 24647
+"PSID54" 0 39 10 0 0 1 1 22745.13 6493.5 25564.67
+"PSID55" 0 33 12 0 0 1 0 10819.07 6369.968 2936.243
+"PSID56" 0 29 8 0 0 1 1 9169.369 6352.065 20575.86
+"PSID57" 0 21 11 0 0 1 1 10679.96 6276.871 10923.35
+"PSID58" 0 31 12 0 0 1 0 23652.27 6228.532 22403.81
+"PSID59" 0 36 12 1 0 1 0 11040.47 6221.371 7215.739
+"PSID60" 0 25 7 0 0 1 1 5597.625 6099.629 122.6513
+"PSID61" 0 35 7 0 0 1 1 10715.23 6087.097 15177.73
+"PSID62" 0 22 9 0 0 1 1 5683.833 6038.758 4742.025
+"PSID63" 0 31 2 0 1 1 1 3262.179 5965.355 9732.307
+"PSID64" 0 40 15 0 0 1 0 10907.24 5922.387 6238.962
+"PSID65" 0 47 3 0 0 1 1 9047.894 5911.645 6145.865
+"PSID66" 0 26 8 0 1 0 1 3168.134 5872.258 11136.15
+"PSID67" 0 42 7 0 0 1 1 10971.89 5806.016 9241.702
+"PSID68" 0 53 12 0 0 0 0 17104.4 5775.581 19965.56
+"PSID69" 0 30 17 1 0 0 0 17827.37 5546.419 14421.13
+"PSID70" 0 28 10 0 0 1 1 10415.46 5544.629 10289.41
+"PSID71" 0 46 11 0 0 1 1 14753.28 5299.355 0
+"PSID72" 0 28 12 0 0 0 0 8256.35 5279.661 21602.88
+"PSID73" 0 27 12 0 1 1 0 17604.01 5222.371 25564.67
+"PSID74" 0 25 10 0 0 1 1 4335.857 5181.194 12418.81
+"PSID75" 0 38 8 0 0 1 1 11242.27 5174.032 0
+"PSID76" 0 26 12 0 1 0 0 7968.338 5109.581 4181.966
+"PSID77" 0 54 12 0 0 0 0 7165.039 5012.903 0
+"PSID78" 0 38 8 0 1 1 1 22606.02 4978.887 8720.065
+"PSID79" 0 23 17 0 0 0 0 0 4876.839 16747.08
+"PSID80" 0 23 8 0 0 1 1 3595.255 4866.097 2782.559
+"PSID81" 0 23 12 0 0 1 0 11690.95 4764.048 14065
+"PSID82" 0 25 12 0 1 1 0 8746.167 4762.258 379.7757
+"PSID83" 0 25 15 0 0 1 0 7386.436 4738.984 12705.49
+"PSID84" 0 37 11 0 1 0 1 615.2098 4713.919 0
+"PSID85" 0 40 12 0 0 1 0 18389.68 4688.855 21857.05
+"PSID86" 0 19 10 0 0 0 1 5777.878 4672.742 135.9508
+"PSID87" 0 48 7 0 0 1 1 13326.93 4636.935 0
+"PSID88" 0 19 12 0 0 0 0 8530.648 4620.823 0
+"PSID89" 0 16 9 0 0 0 1 2539.21 4579.645 0
+"PSID90" 0 29 10 0 0 1 1 713.1731 4542.048 7781.708
+"PSID91" 0 30 16 0 0 0 0 3093.682 4468.645 15538.29
+"PSID92" 0 22 11 0 0 1 1 8761.841 4463.274 10642.59
+"PSID93" 0 22 10 0 0 0 1 17268.98 4400.613 2453.026
+"PSID94" 0 47 10 1 0 1 1 13311.26 4397.032 19330.14
+"PSID95" 0 25 12 0 1 1 0 2266.872 4361.226 3020.473
+"PSID96" 0 47 10 1 0 0 1 21918.32 4323.629 19438.02
+"PSID97" 0 24 12 1 0 1 0 8573.752 4293.194 0
+"PSID98" 0 20 12 1 0 1 0 2648.929 4273.5 0
+"PSID99" 0 28 12 1 0 0 0 16722.34 4253.806 7314.747
+"PSID100" 0 47 11 0 0 0 1 8060.424 4232.323 3358.873
+"PSID101" 0 50 0 0 0 1 1 10162.72 4218 220.1813
+"PSID102" 0 18 12 0 0 0 0 2217.89 4191.145 8957.978
+"PSID103" 0 21 12 0 0 0 0 9665.063 4110.581 1687.564
+"PSID104" 0 47 11 0 0 1 1 23924.61 4096.258 17358.85
+"PSID105" 0 21 12 0 0 0 0 2827.222 4056.871 5937.505
+"PSID106" 0 34 11 0 0 1 1 0 4010.323 18133.18
+"PSID107" 0 19 12 0 0 1 0 5817.063 3919.016 1066.919
+"PSID108" 0 44 13 0 0 1 0 8032.994 3881.419 3104.704
+"PSID109" 0 21 15 0 0 1 0 6951.479 3879.629 0
+"PSID110" 0 20 12 1 0 0 0 5099.971 3842.032 12718.79
+"PSID111" 0 51 11 0 0 0 1 48.98167 3813.387 1525.014
+"PSID112" 0 28 13 0 0 0 0 5260.631 3790.113 9253.524
+"PSID113" 0 24 15 0 0 0 0 12746.99 3743.565 0
+"PSID114" 0 28 8 0 1 1 1 8305.332 3718.5 0
+"PSID115" 0 20 11 0 0 1 1 5822.941 3532.306 11075.56
+"PSID116" 0 29 12 0 0 1 0 14288.93 3503.661 8133.407
+"PSID117" 0 23 12 0 0 1 0 14347.71 3482.177 3818.445
+"PSID118" 0 20 11 1 0 0 1 0 3480.387 5495.665
+"PSID119" 0 42 7 0 0 1 1 4324.102 3457.113 9856.436
+"PSID120" 0 43 12 0 0 1 0 14328.12 3453.532 18781.9
+"PSID121" 0 27 13 0 0 0 0 16406.9 3426.677 5344.937
+"PSID122" 0 27 4 0 1 1 1 626.9654 3410.565 3367.739
+"PSID123" 0 25 12 0 0 1 0 21469.65 3405.194 7981.201
+"PSID124" 0 18 12 0 0 0 0 4729.67 3328.21 12602.05
+"PSID125" 0 31 16 0 0 1 0 25862.32 3254.806 25564.67
+"PSID126" 0 27 12 0 0 1 0 4043.927 3231.532 7240.86
+"PSID127" 0 18 11 0 0 0 1 0 3226.161 15814.63
+"PSID128" 0 24 7 0 0 1 1 7860.578 3213.629 0
+"PSID129" 0 23 12 0 0 1 0 7856.66 3213.629 5535.564
+"PSID130" 0 50 12 0 0 1 0 19929.66 3190.355 18597.19
+"PSID131" 0 19 12 0 0 0 0 99.92261 3172.452 15436.33
+"PSID132" 0 23 10 0 0 1 1 15811.28 3145.597 6398.556
+"PSID133" 0 51 12 0 0 1 0 21001.38 3140.226 16015.6
+"PSID134" 0 19 11 1 0 0 1 5607.422 3054.29 94.5745
+"PSID135" 0 20 10 0 0 1 1 3099.56 2970.145 21141.83
+"PSID136" 0 20 11 0 1 0 1 2868.367 2968.355 7403.41
+"PSID137" 0 21 12 0 0 0 0 8128.998 2939.71 0
+"PSID138" 0 39 10 0 0 1 1 0 2886 18761.22
+"PSID139" 0 36 5 0 0 0 1 3814.692 2873.468 2751.527
+"PSID140" 0 19 9 1 0 0 1 1079.556 2873.468 14344.29
+"PSID141" 0 42 6 0 1 1 1 2425.572 2832.29 1907.745
+"PSID142" 0 20 7 0 0 0 1 1902.448 2792.903 6098.578
+"PSID143" 0 23 12 0 0 1 0 4954.986 2771.419 0
+"PSID144" 0 35 12 0 0 1 0 1469.45 2719.5 0
+"PSID145" 0 18 12 0 0 0 0 881.6701 2696.226 12120.31
+"PSID146" 0 43 8 0 0 1 1 18338.74 2674.742 6395.601
+"PSID147" 0 37 14 0 0 1 0 18501.36 2638.935 13429.58
+"PSID148" 0 24 10 0 0 1 1 4719.874 2565.532 2173.736
+"PSID149" 0 51 12 0 0 0 0 20742.76 2538.677 1019.631
+"PSID150" 0 22 11 0 1 0 1 7341.373 2535.097 14187.65
+"PSID151" 0 19 12 0 0 0 0 336.9939 2518.984 7118.209
+"PSID152" 0 52 0 0 1 1 1 773.9104 2506.452 0
+"PSID153" 0 21 12 0 0 0 0 2903.633 2456.323 4787.834
+"PSID154" 0 24 12 0 0 0 0 9784.578 2413.355 0
+"PSID155" 0 35 8 0 0 1 1 2241.401 2399.032 9460.406
+"PSID156" 0 20 13 0 0 0 0 0 2352.484 0
+"PSID157" 0 17 7 1 0 0 1 1054.086 2286.242 1613.677
+"PSID158" 0 18 10 1 0 0 1 311.5234 2284.452 8154.095
+"PSID159" 0 28 12 1 0 0 0 6285.328 2255.806 7310.313
+"PSID160" 0 25 14 0 1 1 0 1622.273 2239.694 1892.968
+"PSID161" 0 40 12 0 1 0 0 13616.9 2228.952 876.2919
+"PSID162" 0 50 3 0 0 1 1 3136.786 2203.887 13976.34
+"PSID163" 0 48 8 0 0 1 1 16050.31 2116.161 11600.15
+"PSID164" 0 17 7 0 1 0 1 0 2082.145 6460.621
+"PSID165" 0 30 12 0 0 1 0 7347.251 2080.355 14475.81
+"PSID166" 0 30 7 0 0 1 1 574.0652 2010.532 366.4762
+"PSID167" 0 22 11 0 0 1 1 3030.986 1976.516 0
+"PSID168" 0 27 12 0 0 1 0 11493.06 1906.694 13419.24
+"PSID169" 0 25 9 0 0 1 1 23377.97 1901.323 1898.879
+"PSID170" 0 21 14 0 0 0 0 80.32994 1890.581 6389.69
+"PSID171" 0 17 10 0 0 0 1 0 1888.79 19993.64
+"PSID172" 0 39 7 0 0 0 1 7786.126 1844.032 9206.237
+"PSID173" 0 18 9 1 0 0 1 1183.397 1822.548 803.8833
+"PSID174" 0 25 12 0 0 1 0 2721.422 1754.516 1037.364
+"PSID175" 0 20 8 0 0 1 1 2360.916 1741.984 0
+"PSID176" 0 19 13 0 0 0 0 2366.794 1709.758 0
+"PSID177" 0 19 11 0 0 0 1 0 1693.645 9853.481
+"PSID178" 0 22 12 0 0 0 0 10137.25 1679.323 25564.67
+"PSID179" 0 18 11 1 0 0 1 2068.986 1623.823 20243.38
+"PSID180" 0 21 10 0 0 0 1 1767.259 1555.79 7675.312
+"PSID181" 0 24 12 0 0 1 0 7643.1 1546.839 3262.82
+"PSID182" 0 18 11 0 0 0 1 1273.523 1532.516 12489.75
+"PSID183" 0 17 10 0 0 0 1 568.1874 1525.355 6231.573
+"PSID184" 0 17 10 0 0 0 1 0 1503.871 7843.773
+"PSID185" 0 18 10 0 0 0 1 0 1491.339 237.914
+"PSID186" 0 53 10 0 1 0 1 7878.212 1489.548 13170.98
+"PSID187" 0 18 11 1 0 0 1 1191.234 1478.806 3683.972
+"PSID188" 0 17 10 0 1 0 1 0 1453.742 6918.716
+"PSID189" 0 26 12 1 0 0 0 0 1448.371 0
+"PSID190" 0 39 5 0 0 1 1 13082.02 1434.048 18323.81
+"PSID191" 0 18 12 1 0 0 0 1579.169 1408.984 3057.416
+"PSID192" 0 23 13 0 0 0 0 601.4949 1394.661 4975.505
+"PSID193" 0 18 8 0 0 0 1 5023.56 1391.081 6756.166
+"PSID194" 0 28 10 0 0 1 1 7578.444 1383.919 2404.261
+"PSID195" 0 32 4 0 0 1 1 0 1378.548 0
+"PSID196" 0 18 11 1 0 0 1 0 1367.806 33.98771
+"PSID197" 0 40 10 0 0 1 1 1543.902 1342.742 0
+"PSID198" 0 21 14 0 0 0 0 8456.196 1330.21 16967.26
+"PSID199" 0 29 10 0 1 0 1 3732.403 1323.048 6694.101
+"PSID200" 0 31 6 0 0 0 1 2666.562 1321.258 0
+"PSID201" 0 46 7 0 0 1 1 19171.43 1317.677 0
+"PSID202" 0 20 9 0 1 1 1 0 1283.661 0
+"PSID203" 0 36 18 0 0 1 0 3273.935 1269.339 18227.76
+"PSID204" 0 45 12 0 0 1 0 16559.72 1265.758 7987.112
+"PSID205" 0 16 10 0 0 0 1 1026.656 1224.581 6847.785
+"PSID206" 0 18 12 0 0 0 0 818.9735 1208.468 2232.845
+"PSID207" 0 40 12 0 1 0 0 11867.28 1195.935 3873.121
+"PSID208" 0 16 9 0 0 0 1 0 1188.774 2451.548
+"PSID209" 0 16 10 0 0 0 1 574.0652 1181.613 5578.418
+"PSID210" 0 28 5 0 1 1 1 10967.98 1178.032 239.3917
+"PSID211" 0 20 12 0 0 0 0 0 1147.597 15554.55
+"PSID212" 0 19 8 0 0 1 1 39.18534 1136.855 5327.204
+"PSID213" 0 16 8 0 0 0 1 0 1113.581 542.3257
+"PSID214" 0 20 11 0 0 1 1 2547.047 1099.258 0
+"PSID215" 0 35 10 0 0 1 1 4964.782 1086.726 1745.195
+"PSID216" 0 32 6 0 1 1 1 979.6334 1036.597 0
+"PSID217" 0 32 16 1 0 0 0 17135.75 1031.226 0
+"PSID218" 0 17 9 1 0 0 1 0 981.0968 8900.347
+"PSID219" 0 16 7 0 0 0 1 0 975.7258 4728.725
+"PSID220" 0 32 15 0 0 0 0 489.8167 968.5645 7684.178
+"PSID221" 0 19 12 0 0 0 0 815.055 964.9839 12059.73
+"PSID222" 0 40 12 0 0 1 0 16851.65 961.4032 17717.94
+"PSID223" 0 50 7 0 0 1 1 11473.47 956.0323 0
+"PSID224" 0 39 11 0 0 0 1 0 930.9677 0
+"PSID225" 0 18 8 0 1 0 1 0 902.3226 1306.31
+"PSID226" 0 39 10 1 0 0 1 844.444 889.7903 701.9201
+"PSID227" 0 17 11 0 1 0 1 0 873.6774 7759.542
+"PSID228" 0 17 5 1 0 0 1 96.00407 868.3065 0
+"PSID229" 0 19 12 0 0 0 0 2425.572 861.1452 2587.499
+"PSID230" 0 27 15 0 0 0 0 0 857.5645 3392.86
+"PSID231" 0 18 11 1 0 0 1 587.78 841.4516 7933.914
+"PSID232" 0 20 14 0 0 1 0 0 805.6452 1454.083
+"PSID233" 0 20 12 0 0 1 0 12145.49 791.3226 13683.75
+"PSID234" 0 19 13 1 0 0 0 1714.358 785.9516 9067.33
+"PSID235" 0 24 8 0 0 1 1 213.5601 760.8871 2340.719
+"PSID236" 0 27 12 0 0 1 0 4222.22 751.9355 0
+"PSID237" 0 19 9 0 0 0 1 773.9104 676.7419 5647.871
+"PSID238" 0 52 8 1 0 1 1 5454.599 666 0
+"PSID239" 0 18 11 0 1 0 1 0 630.1935 0
+"PSID240" 0 16 10 0 1 0 1 0 630.1935 3892.332
+"PSID241" 0 18 12 0 1 0 0 0 630.1935 4843.988
+"PSID242" 0 45 12 0 0 0 0 4473.006 608.7097 0
+"PSID243" 0 21 14 0 0 0 0 9708.167 594.3871 2256.488
+"PSID244" 0 36 8 0 0 1 1 2715.544 585.4355 0
+"PSID245" 0 21 13 0 0 0 0 513.3279 578.2742 0
+"PSID246" 0 41 7 0 0 1 1 19573.08 565.7419 0
+"PSID247" 0 18 7 0 0 0 1 491.776 558.5806 642.8111
+"PSID248" 0 39 9 0 0 0 1 11230.52 537.0968 5752.79
+"PSID249" 0 19 3 0 0 1 1 0 537.0968 0
+"PSID250" 0 32 13 0 0 1 0 12553.02 524.5645 15353.58
+"PSID251" 0 16 9 0 0 0 1 0 485.1774 4112.513
+"PSID252" 0 16 7 0 0 0 1 658.3136 479.8065 6210.885
+"PSID253" 0 21 9 1 0 0 1 1030.574 470.8548 1223.558
+"PSID254" 0 22 12 0 0 1 0 12096.51 469.0645 14289.62
+"PSID255" 0 23 11 0 1 1 1 8946.012 469.0645 4776.012
+"PSID256" 0 17 8 1 0 0 1 0 451.1613 0
+"PSID257" 0 21 8 0 0 1 1 5699.507 388.5 8844.194
+"PSID258" 0 18 10 0 0 0 1 0 386.7097 0
+"PSID259" 0 24 12 0 0 1 0 9051.813 327.629 8547.171
+"PSID260" 0 24 12 1 0 1 0 4232.016 320.4677 1273.8
+"PSID261" 0 16 9 0 0 0 1 0 320.4677 3707.616
+"PSID262" 0 20 8 0 0 1 1 621.0876 306.1452 5551.819
+"PSID263" 0 42 8 0 0 0 1 17925.33 300.7742 14116.72
+"PSID264" 0 17 8 0 1 0 1 391.8534 300.7742 18891.26
+"PSID265" 0 19 8 0 1 0 1 368.3422 300.7742 18510
+"PSID266" 0 17 9 1 0 0 1 0 297.1935 54.67588
+"PSID267" 0 21 14 0 0 0 0 107.7597 293.6129 7698.955
+"PSID268" 0 16 9 1 0 0 1 0 277.5 3983.951
+"PSID269" 0 23 13 1 0 0 0 172.4155 272.129 582.2243
+"PSID270" 0 16 9 0 0 0 1 411.446 254.2258 1725.985
+"PSID271" 0 17 11 0 1 0 1 803.2994 248.8548 5173.521
+"PSID272" 0 46 7 0 0 0 1 1081.515 245.2742 0
+"PSID273" 0 32 10 0 0 1 1 4145.809 238.1129 8245.714
+"PSID274" 0 18 11 0 0 0 1 131.2709 218.4194 7503.896
+"PSID275" 0 23 12 0 1 1 0 0 216.629 0
+"PSID276" 0 18 10 0 0 1 1 0 211.2581 14053.18
+"PSID277" 0 19 10 1 0 0 1 1056.045 205.8871 0
+"PSID278" 0 16 7 1 0 0 1 133.2301 205.8871 6145.865
+"PSID279" 0 26 7 0 0 1 1 1538.024 189.7742 650.1997
+"PSID280" 0 16 10 0 0 0 1 0 189.7742 2136.793
+"PSID281" 0 17 10 0 0 0 1 0 182.6129 6423.677
+"PSID282" 0 17 10 0 0 0 1 0 171.871 1483.637
+"PSID283" 0 23 8 0 0 1 1 33.30754 166.5 0
+"PSID284" 0 29 12 0 0 1 0 14641.6 162.9194 9473.705
+"PSID285" 0 17 10 0 0 0 1 0 152.1774 10301.23
+"PSID286" 0 49 8 0 0 1 1 14684.7 136.0645 14963.46
+"PSID287" 0 20 10 0 0 1 1 6563.544 134.2742 15363.92
+"PSID288" 0 40 16 0 0 1 0 0 114.5806 0
+"PSID289" 0 19 10 0 0 0 1 1933.796 112.7903 675.321
+"PSID290" 0 18 11 0 0 0 1 1481.206 57.29032 1421.573
+"PSID291" 0 16 6 1 0 0 1 0 44.75806 0
+"PSID292" 0 22 8 0 0 1 1 105.8004 42.96774 209.8372
+"PSID293" 0 31 12 1 0 1 0 0 42.96774 11023.84
+"PSID294" 0 20 11 0 0 1 1 4478.884 39.3871 6280.338
+"PSID295" 0 17 11 0 1 0 1 601.4949 10.74194 1913.656
+"PSID296" 0 50 12 0 0 1 0 25862.32 0 25564.67
+"PSID297" 0 49 14 0 0 1 0 25862.32 0 25564.67
+"PSID298" 0 47 9 0 0 1 1 25862.32 0 25564.67
+"PSID299" 0 34 11 0 1 1 1 22198.49 0 0
+"PSID300" 0 22 8 1 0 1 1 16961.37 0 959.0445
+"PSID301" 0 27 12 0 0 1 0 15509.56 0 12593.19
+"PSID302" 0 30 10 0 0 1 1 14913.94 0 11563.21
+"PSID303" 0 52 12 0 0 1 0 14780.71 0 25564.67
+"PSID304" 0 43 12 0 0 1 0 13321.05 0 16860.86
+"PSID305" 0 27 9 0 1 1 1 12829.28 0 0
+"PSID306" 0 35 13 0 0 0 0 9537.711 0 11269.14
+"PSID307" 0 45 12 0 0 1 0 9277.128 0 12108.49
+"PSID308" 0 22 11 1 0 1 1 9049.853 0 9088.018
+"PSID309" 0 22 12 0 0 1 0 9022.424 0 3342.618
+"PSID310" 0 23 11 0 0 1 1 8910.745 0 4183.444
+"PSID311" 0 55 7 0 0 1 1 8832.375 0 0
+"PSID312" 0 26 14 0 0 0 0 8411.132 0 0
+"PSID313" 0 34 12 0 0 0 0 8125.079 0 6032.08
+"PSID314" 0 22 11 0 0 0 1 8013.401 0 5748.356
+"PSID315" 0 31 12 0 0 0 0 6156.016 0 4094.78
+"PSID316" 0 19 12 0 0 0 0 5797.47 0 2160.436
+"PSID317" 0 24 10 0 0 1 1 5523.173 0 5040.525
+"PSID318" 0 36 12 0 0 1 0 5374.269 0 0
+"PSID319" 0 20 9 0 0 1 1 5229.283 0 15892.95
+"PSID320" 0 23 8 0 0 1 1 4610.155 0 0
+"PSID321" 0 35 11 0 0 1 1 3975.352 0 21963.45
+"PSID322" 0 23 12 0 0 0 0 3893.063 0 16324.45
+"PSID323" 0 29 10 0 0 0 1 3751.996 0 251.2135
+"PSID324" 0 24 9 0 0 1 1 3438.513 0 818.6605
+"PSID325" 0 18 10 0 0 0 1 3360.143 0 0
+"PSID326" 0 45 8 1 0 0 1 3299.405 0 31.03226
+"PSID327" 0 21 13 0 1 0 0 3015.312 0 17627.8
+"PSID328" 0 29 13 0 0 1 0 2780.2 0 14339.86
+"PSID329" 0 21 15 0 0 0 0 2629.336 0 1717.118
+"PSID330" 0 22 16 1 0 0 0 2564.68 0 116.7404
+"PSID331" 0 24 12 1 0 1 0 2355.039 0 2448.593
+"PSID332" 0 20 14 0 0 0 0 2210.053 0 2813.591
+"PSID333" 0 19 6 1 0 0 1 1955.348 0 14998.92
+"PSID334" 0 19 9 0 1 0 1 1822.118 0 3372.172
+"PSID335" 0 19 12 1 0 0 0 1681.051 0 0
+"PSID336" 0 20 13 0 0 0 0 1657.54 0 913.235
+"PSID337" 0 19 12 1 0 0 0 1655.58 0 0
+"PSID338" 0 26 5 0 0 1 1 1573.291 0 3700.227
+"PSID339" 0 26 9 0 1 0 1 1563.495 0 2862.356
+"PSID340" 0 23 12 0 0 0 0 1504.717 0 0
+"PSID341" 0 20 9 0 1 0 1 1500.798 0 12618.31
+"PSID342" 0 20 10 0 0 0 1 1412.631 0 6290.682
+"PSID343" 0 36 11 0 0 1 1 1404.794 0 0
+"PSID344" 0 39 12 0 0 1 0 1289.198 0 1202.869
+"PSID345" 0 17 9 1 0 0 1 1222.582 0 422.6298
+"PSID346" 0 55 3 0 0 0 1 1208.868 0 0
+"PSID347" 0 28 8 0 0 1 1 1202.99 0 19516.33
+"PSID348" 0 19 12 0 1 0 0 1058.004 0 8923.991
+"PSID349" 0 37 7 0 0 1 1 963.9593 0 0
+"PSID350" 0 16 9 0 0 1 1 920.8554 0 15997.87
+"PSID351" 0 17 10 0 0 0 1 646.558 0 9438.24
+"PSID352" 0 24 12 1 0 0 0 566.2281 0 2284.565
+"PSID353" 0 19 11 0 0 0 1 540.7576 0 3406.16
+"PSID354" 0 50 5 1 0 1 1 411.446 0 9166.338
+"PSID355" 0 19 9 1 0 0 1 384.0163 0 0
+"PSID356" 0 36 1 1 0 0 1 348.7495 0 0
+"PSID357" 0 18 11 0 0 0 1 321.3198 0 7722.599
+"PSID358" 0 16 7 0 1 0 1 289.9715 0 7515.717
+"PSID359" 0 21 11 0 0 1 1 246.8676 0 6708.879
+"PSID360" 0 55 6 0 0 1 1 111.6782 0 0
+"PSID361" 0 37 12 0 0 0 0 48.98167 0 877.7696
+"PSID362" 0 26 12 0 1 1 0 47.0224 0 0
+"PSID363" 0 54 12 0 0 1 0 0 0 0
+"PSID364" 0 50 12 0 0 1 0 0 0 0
+"PSID365" 0 16 8 0 0 0 1 0 0 2559.422
+"PSID366" 0 16 9 0 1 0 1 0 0 0
+"PSID367" 0 18 10 1 0 0 1 0 0 2281.61
+"PSID368" 0 40 11 1 0 1 1 0 0 0
+"PSID369" 0 16 8 0 0 0 1 0 0 0
+"PSID370" 0 16 9 1 0 0 1 0 0 2158.959
+"PSID371" 0 26 14 0 0 0 0 0 0 6717.745
+"PSID372" 0 20 9 1 0 0 1 0 0 6083.8
+"PSID373" 0 20 12 1 0 0 0 0 0 0
+"PSID374" 0 18 11 1 0 0 1 0 0 0
+"PSID375" 0 46 11 1 0 1 1 0 0 2820.98
+"PSID376" 0 17 8 1 0 0 1 0 0 12760.17
+"PSID377" 0 16 9 0 0 0 1 0 0 4974.028
+"PSID378" 0 30 10 0 0 1 1 0 0 3151.991
+"PSID379" 0 33 12 0 1 1 0 0 0 5841.453
+"PSID380" 0 34 12 1 0 1 0 0 0 18716.88
+"PSID381" 0 21 13 1 0 0 0 0 0 17941.08
+"PSID382" 0 29 11 0 0 1 1 0 0 0
+"PSID383" 0 19 12 0 0 0 0 0 0 0
+"PSID384" 0 31 4 0 1 0 1 0 0 1161.493
+"PSID385" 0 19 12 0 1 0 0 0 0 18573.55
+"PSID386" 0 20 12 1 0 0 0 0 0 11594.24
+"PSID387" 0 55 4 1 0 0 1 0 0 0
+"PSID388" 0 19 11 1 0 0 1 0 0 16485.52
+"PSID389" 0 18 11 1 0 0 1 0 0 7146.286
+"PSID390" 0 48 13 0 0 1 0 0 0 0
+"PSID391" 0 16 9 0 1 1 1 0 0 6821.186
+"PSID392" 0 17 10 1 0 0 1 0 0 0
+"PSID393" 0 38 12 0 0 1 0 0 0 18756.78
+"PSID394" 0 34 8 0 0 1 1 0 0 2664.341
+"PSID395" 0 53 12 0 0 0 0 0 0 0
+"PSID396" 0 48 14 0 0 1 0 0 0 7236.427
+"PSID397" 0 16 9 0 0 0 1 0 0 6494.608
+"PSID398" 0 17 8 1 0 0 1 0 0 4520.366
+"PSID399" 0 27 14 1 0 0 0 0 0 10122.43
+"PSID400" 0 37 8 1 0 0 1 0 0 648.722
+"PSID401" 0 17 10 1 0 0 1 0 0 1053.619
+"PSID402" 0 16 8 0 0 0 1 0 0 0
+"PSID403" 0 48 12 0 0 1 0 0 0 1491.026
+"PSID404" 0 55 7 0 0 0 1 0 0 0
+"PSID405" 0 21 15 0 0 0 0 0 0 0
+"PSID406" 0 16 10 1 0 0 1 0 0 1730.418
+"PSID407" 0 23 12 0 0 0 0 0 0 3902.676
+"PSID408" 0 46 11 1 0 1 1 0 0 0
+"PSID409" 0 17 10 0 0 0 1 0 0 14942.77
+"PSID410" 0 42 16 0 0 0 0 0 0 23764.8
+"PSID411" 0 18 10 1 0 0 1 0 0 5306.516
+"PSID412" 0 53 12 1 0 0 0 0 0 0
+"PSID413" 0 17 10 0 0 1 1 0 0 3859.822
+"PSID414" 0 17 6 0 0 0 1 0 0 0
+"PSID415" 0 43 6 0 0 1 1 0 0 0
+"PSID416" 0 34 12 1 0 0 0 0 0 0
+"PSID417" 0 16 8 0 1 0 1 0 0 12242.96
+"PSID418" 0 27 12 0 0 1 0 0 0 1533.88
+"PSID419" 0 51 4 1 0 0 1 0 0 0
+"PSID420" 0 39 2 1 0 1 1 0 0 964.9555
+"PSID421" 0 55 8 0 0 1 1 0 0 0
+"PSID422" 0 16 9 0 0 0 1 0 0 5551.819
+"PSID423" 0 27 10 1 0 0 1 0 0 7543.794
+"PSID424" 0 25 14 0 0 0 0 0 0 0
+"PSID425" 0 18 11 0 0 0 1 0 0 10150.5
+"PSID426" 0 24 1 0 1 1 1 0 0 19464.61
+"PSID427" 0 21 18 0 0 0 0 0 0 0
+"PSID428" 0 32 5 1 0 1 1 0 0 187.6713
+"PSID429" 0 16 9 0 0 0 1 0 0 1495.459
diff --git a/demo/00Index b/demo/00Index
new file mode 100644
index 0000000..8a58bec
--- /dev/null
+++ b/demo/00Index
@@ -0,0 +1,3 @@
+lalonde General examples of MatchIt methods, using Lalonde dataset
+analysis Demo of analysis procedures to estimate treatment effects
+diagnose Demo of how to diagnose a propensity score model
diff --git a/demo/analysis.R b/demo/analysis.R
new file mode 100644
index 0000000..9351ca5
--- /dev/null
+++ b/demo/analysis.R
@@ -0,0 +1,82 @@
+# Examples of analyses, with and without Zelig
+# Without zelig
+
+data(lalonde)
+
+# First just neyman on matched sample
+print("Neyman estimates using a matched sample")
+match.out1 <- matchit(treat ~ age + educ + black + hispan + nodegree +
+ married + re74 + re75, data = lalonde)
+neyman.out <- neyman(re78, match.out1)
+print(neyman.out)
+print(summary(neyman.out))
+silent <- readline("\nPress <return> to continue: ")
+
+# Neyman using propensity score subclasses
+print("Neyman estimates with propensity score subclassification")
+match.out1s <- matchit(treat ~ age + educ + black + hispan + nodegree
+ + married + re74 + re75, data = lalonde, subclass=4)
+neyman.outs <- neyman(re78, match.out1s)
+print(neyman.outs)
+print(summary(neyman.outs))
+silent <- readline("\nPress <return> to continue: ")
+
+# Neyman using exact matching
+print("Neyman estimates using exact matching")
+match.out1e <- matchit(treat ~ black + hispan, data = lalonde, exact=TRUE)
+neyman.oute <- neyman(re78, match.out1e)
+print(neyman.oute)
+print(summary(neyman.oute))
+silent <- readline("\nPress <return> to continue: ")
+
+# Run linear regression on matched data: use match.out1
+print("Linear regression on matched data set")
+lm.1 <- lm(re78 ~ treat + pscore, data=match.data(match.out1))
+
+print(summary(lm.1)$coef[2,])
+silent <- readline("\nPress <return> to continue: ")
+
+# Zelig part from zelig documentation (match.R in demo directory)
+if ("Zelig"%in%.packages(all=T)) {
+ library(Zelig)
+
+ ## an example for propensity score matching
+ print("Using Zelig for analysis")
+ z.out1 <- zelig(re78 ~ pscore, data = match.data(match.out1,
+ "control"), model = "ls")
+ x.out1 <- setx(z.out1, data = match.data(match.out1, "treat"),
+ fn = NULL, cond = TRUE)
+ s.out1 <- sim(z.out1, x = x.out1)
+ print(summary(s.out1))
+ user.prompt()
+
+ ## an example for subclassification
+ print("Using Zelig for analysis, with propensity score subclassification")
+ z.out2 <- zelig(re78 ~ pscore, data = match.data(match.out1s,
+ "control"), model="ls", by="psclass")
+ x.out2 <- setx(z.out2, data = match.data(match.out1s, "treat"),
+ fn = NULL, cond = TRUE)
+ s.out2 <- sim(z.out2, x = x.out2, num = 100)
+ print(summary(s.out2)) # overall results
+ print(summary(s.out2, subset = 1)) # subclass 1
+ print(summary(s.out2, subset = 2)) # subclass 2
+ print(summary(s.out2, subset = 3)) # subclass 3
+ user.prompt()
+}
+if(!("Zelig"%in%.packages(all=T))){
+ cat("Zelig is not installed. \nFor further illustrations of analyses,\nyou may install Zelig from \nhttp://gking.harvard.edu/zelig/\n")
+}
+
+# Full matching: if package installed
+if ("optmatch"%in%.packages(all=T)) {
+print("Analysis after full matching: using fixed effects")
+foo1 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, full=T)
+m1 <- lm(re78~ treat + age + educ + black + hispan + married +
+ nodegree + re74 + re75 + as.factor(psclass),
+ data=foo1$data)
+print(summary(m1))
+}
+if(!("optmatch"%in%.packages(all=T)))
+ cat("Optmatch is not installed. \nMore information is available from \nhttp://www.stat.lsa.umich.edu/~bbh/optmatch.html \n")
+
diff --git a/demo/diagnose.R b/demo/diagnose.R
new file mode 100644
index 0000000..4f63257
--- /dev/null
+++ b/demo/diagnose.R
@@ -0,0 +1,35 @@
+# Demo of propensity score specification diagnostics
+# Initial model
+data(lalonde)
+match.out1 <- matchit(treat ~ age + married + black + hispan +
+ nodegree + educ + re74 + re75, data=lalonde,
+ subclass=c(0, .5, 1))
+
+summary(match.out1, verbose=F)
+
+# Create more blocks if t-stat of pscore > 2.5
+match.out1 <- matchit(treat ~ age + educ + married + nodegree + hispan
+ + black + re74 + re75, data=lalonde,
+ subclass=c(0,.25, .5,.75,1))
+
+summary(match.out1, verbose=T)
+
+# Add terms that are significant in > 1 block
+match.out2 <- matchit(treat ~ age + educ + married + nodegree + hispan
+ + black + re74 + re75 + I(educ*black) +
+ I(black*re74) + I(age*married) + I(educ*married)
+ + I(educ*nodegree) + I(married*nodegree) +
+ I(married*black) + I(age*nodegree),
+ data=lalonde, subclass=c(0,.5,1))
+
+summary(match.out2, verbose=F)
+
+match.out2 <- matchit(treat ~ age + educ + married + nodegree + hispan
+ + black + re74 + re75 + I(educ*black) +
+ I(black*re74) + I(age*married) + I(educ*married)
+ + I(educ*nodegree) + I(married*nodegree) +
+ I(married*black) + I(age*nodegree),
+ data=lalonde, subclass=c(0,.25,.5,1))
+
+summary(match.out2, verbose=F)
+
diff --git a/demo/lalonde.R b/demo/lalonde.R
new file mode 100644
index 0000000..3b143a0
--- /dev/null
+++ b/demo/lalonde.R
@@ -0,0 +1,112 @@
+# MatchIt Demo script with Lalonde data
+
+data(lalonde)
+
+#Exact matching
+foo1 <- matchit(treat ~ black + hispan, exact=TRUE, data=lalonde)
+print(foo1)
+summary(foo1)
+summary(foo1,verbose=T)
+row.names(foo1$data)[foo1$psclass==1][1:20]
+foo1$data[foo1$psclass==2,c("black", "hispan")][1:20,]
+silent <- readline("\nPress <return> to continue: ")
+
+#Propensity score matching
+foo2 <- matchit(treat ~ re74 + re75, data=lalonde)
+print(foo2)
+summary(foo2)
+summary(foo2, verbose=T)
+plot(foo2)
+
+# Greedy vs. optimal matching
+if ("optmatch"%in%.packages(all=T)) {
+foo <- matchit(treat ~ re74+re75+age+educ, data=lalonde)
+print(summary(foo))
+foo <- matchit(treat ~ re74+re75+age+educ, data=lalonde, m.order=1)
+print(summary(foo))
+}
+
+#Propensity score with exact restriction
+foo3 <- matchit(treat ~ re74 + re75, exact=c("black","hispan"), data=lalonde)
+print(foo3)
+sum.foo3 <- summary(foo3)
+sum.foo3$sum.all[c("hispan","black"),]
+sum.foo3$sum.matched[c("hispan","black"),]
+dta.matched <- match.data(foo3)
+dta.treated <- match.data(foo3, group = "treat")
+dta.control <- match.data(foo3, group = "control")
+
+#Replication of Dehejia & Wahba
+foo <- matchit(treat ~ age + I(age^2) + educ + I(educ^2) + black +
+ hispan + married + nodegree + re74 + I(re74^2) +
+ re75 + I(re75^2) + I(as.numeric(re74==0)) +
+ I(as.numeric(re75==0)), data=lalonde, replace=TRUE,
+ discard=1)
+
+#Non-matched and discarded units
+foo$in.sample[!foo$in.sample]
+foo$psweights[foo$in.sample & foo$psweights==0]
+sdisc <- apply(lalonde[!foo$in.sample,],2,mean) #taking column-wise means
+streat <- apply(lalonde[lalonde$treat==1,],2,mean)
+round(rbind(sdisc,streat),2) #rounding
+
+#Look at weights
+print(foo$psweights[foo$data$treat==1][1:10])
+print(foo$psweights[foo$data$treat==0][1:20])
+print((foo$psweights[foo$data$treat==0]*sum(foo$psweights[foo$data$treat==1])/sum(foo$psweights[foo$data$treat==0]))[1:20])
+
+#Subclassification
+foo1 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, replace=TRUE,
+ subclass=6)
+foo2 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, nearest=FALSE,
+ replace=TRUE, subclass=6)
+print(foo2)
+plot(foo2)
+
+#Exact matching on all covariates
+foo <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, exact=TRUE)
+summary(foo)
+
+# Full matching
+if ("optmatch"%in%.packages(all=T)) {
+foo1 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, full=T)
+print(summary(foo1))
+}
+
+# Caliper matching
+foo <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data=lalonde, caliper=0.25,
+ replace=TRUE)
+
+#Mahalanobis matching
+foo <- matchit(treat ~ age + educ, data=lalonde, mahvars = c("age", "educ"),
+ caliper=Inf, replace=TRUE)
+
+#Assignment Model specification
+foo1 <- matchit(treat~educ+re74+re75,data=lalonde,model="probit")
+foo2 <- matchit(treat~educ+re74+re75,data=lalonde,model="nnet")
+foo3 <- matchit(treat~educ+re74+re75,data=lalonde,model="GAM")
+foo4 <- matchit(treat~educ+re74+re75,data=lalonde,model="cart")
+
+#Using Observation Names
+test <- lalonde[1:4,] #taking a lalonde subset
+row.names(test) <- c("Dan","Kosuke","Liz","Gary") #assigning row names
+print(test)
+
+#Matching on One Covariate
+index <- rep(TRUE,nrow(lalonde)) #keeping all units
+names(index) <- row.names(lalonde) #assigning obs names
+treat <- lalonde$treat
+names(treat) <- row.names(lalonde)
+re75 <- lalonde$re75
+names(re75) <- row.names(lalonde)
+foo <- matchdef(treat~re75,index,pscore=re75,replace=TRUE,data=lalonde)
+cbind(lalonde[row.names(foo$match.matrix),]$re75,
+ lalonde[as.vector(foo$match.matrix[,1]),]$re75) #matched re75
+plot(lalonde[row.names(foo$match.matrix),]$re75,
+ lalonde[as.vector(foo$match.matrix[,1]),]$re75, pch=16,
+ xlab="Treated", ylab="Controls") #plotting matched data
diff --git a/man/diagnose.Rd b/man/diagnose.Rd
new file mode 100644
index 0000000..8121ac1
--- /dev/null
+++ b/man/diagnose.Rd
@@ -0,0 +1,85 @@
+\name{diagnose}
+
+\alias{diagnose}
+
+\title{Diagnostics for matching procedure}
+
+\description{\code{diagnose} is a sub-function of \code{matchit} which calculates summary
+statistics for the matching, including the weight for each unit.}
+
+
+\details{This is a sub-function of \code{matchit} which calculates summary statistics for the
+matching procedure in \code{matchit}, calculating t-statistics for covariates as well as
+weights for each unit. This function is called directly by \code{matchit} and does not
+generally need to be called directly by users; these details are included for advanced users. }
+
+\usage{ diagnose <- diagnose(formula, match.matrix, pscore, in.sample, data, exact=FALSE,
+ mahvars=NULL, subclass=0, psclass=NULL, nearest=TRUE, q.cut=NULL, counter=TRUE)
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. The \code{+} symbol means
+"inclusion" not "addition." You may also include interaction terms in the form if
+\code{I(X1*X2)} or squared terms in the form of \code{I(X1^2)}.}
+
+ \item{match.matrix}{(required). n1 by ratio data frame where the rows correspond to treated
+units and the columns store the names of the control units matched to each treated unit. NA
+indicates that treated unit was not matched. Generally created in \code{nearest}.}
+
+ \item{pscore}{(required). Vector of estimated propensity scores. Generally calculated in
+\code{distance}.}
+
+ \item{in.sample}{(required). Vector of length n showing whether each unit was eligible for
+matching due to common support restrictions with \code{discard}. Generally calculated in
+\code{distance}.}
+
+ \item{data}{(required). Data frame containing the variables called in the \code{formula}.
+The dataframe should not include variables with the names \code{psclass}, \code{psweights}, or
+\code{pscore}, as these are expressly reserved in the output dataframe for MatchIt.}
+
+ \item{exact}{"FALSE" (default)=no exact matching. "TRUE"=exact matching on all
+variables in \code{formula}. A vector of variable names (that are in \code{data} to indicate
+separate variables on which to exact match, in combination with matching on the propensity
+score.}
+
+ \item{mahvars}{Variables on which to perform Mahalanobis matching within each caliper
+(default=NULL). Should be entered as a vector of names of variables in \code{data}.}
+
+ \item{subclass}{Either a scaler specifying the number of subclasses (default=0) or a
+vector of probabilities to create quantiles based on \code{sub.by}.}
+
+ \item{psclass}{Subclass index in an ordinal scale from 1 to the number of subclasses.
+Unmatched units have subclass=0. Generally computed in \code{subclassify}.}
+
+ \item{nearest}{Whether to perform nearest-neighbor matching (default=TRUE). }
+ \item{q.cut}{Subclass cut points. Generally calculated in \code{subclassify}.}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).} }
+
+\value{
+ \item{psweights}{Vector of length n giving the weight assigned to each unit in the matching
+process. Each weight is proportional to the number of times that unit was matched.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
+
+
+
+
+
+
+
+
diff --git a/man/distance.Rd b/man/distance.Rd
new file mode 100644
index 0000000..4364841
--- /dev/null
+++ b/man/distance.Rd
@@ -0,0 +1,77 @@
+\name{distance}
+
+\alias{distance}
+
+\title{Distance function: estimating propensity scores}
+
+\description{ The distance function calculates the distance measure to be used in the matching,
+usually the propensity score. It is a sub-function of \code{matchit}.}
+
+\details{This is a sub-function of the \code{matchit} command, which calculates the distance
+measure used in the matching, usually the propensity score. This function is called directly by
+\code{matchit} and does not generally need to be called directly by users; these details are
+included for advanced users. }
+
+
+\usage{
+distance <- distance(formula, model="logit", data, discard=0, reestimate=FALSE, counter=TRUE, ...)
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. The \code{+} symbol means
+"inclusion" not "addition." You may also include interaction terms in the form if
+\code{I(X1*X2)} or squared terms in the form of \code{I(X1^2)}.}
+
+ \item{data}{(required). Data frame containing the variables called in the \code{formula}.
+The dataframe should not include variables with the names \code{psclass}, \code{psweights}, or
+\code{pscore}, as these are expressly reserved in the output dataframe for MatchIt.}
+
+ \item{model}{Method used to estimate the propensity score. May be "logit" (default),
+"probit", "nnet", "GAM", or "cart".}
+
+ \item{discard}{Whether to discard units that fall outside some measure of support of the
+distance score. 0 (default)=keep all units. 1=keep all units with common support. 2=discard
+only control units outside the support of the distance measure of the treated units. 3=discard
+only treated units outside the support of the distance measure of the control units.}
+
+ \item{reestimate}{Specifies whether to reestimate the propensity score model after
+discarding units (default=FALSE).}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).}
+
+ \item{...}{Additional arguments to be passed to \code{distance}, depending on the model
+to be used.} }
+
+\value{
+ \item{in.sample}{Vector of length n showing whether each unit was eligible for matching due to
+common support restrictions with \code{discard}.}
+ \item{pscore}{Vector of estimated propensity scores.}
+ \item{treat}{The treatment indicator from \code{data}.}
+ \item{covariates}{Covariates used in the right-hand side of the assignment model.}
+ \item{assign.model}{Output of the assignment model.}
+}
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{internal}
+
+
+
+
+
+
+
+
diff --git a/man/exactmatch.Rd b/man/exactmatch.Rd
new file mode 100644
index 0000000..c09d4ea
--- /dev/null
+++ b/man/exactmatch.Rd
@@ -0,0 +1,59 @@
+\name{exactmatch}
+
+\alias{exactmatch}
+
+\title{Exact matching}
+
+\description{\code{exactmatch} is a sub-function of \code{matchit} that does exact matching
+on the set of covariates.}
+
+
+\details{This sub-function is called automatically by the \code{matchit} command when
+\code{exact=TRUE} is chosen, meaning exact matching on all covariates in the assignment model.
+It creates subclasses defined by the values of the covariates, where each subclass must contain
+some treated and some control units. Units with no exact matches in the other treatment group
+are discarded. This function is called directly by \code{matchit} and does not generally need
+to be called directly by users; these details are included for advanced users. }
+
+\usage{
+exactmatch <- exactmatch(formula, data, counter=TRUE)
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. Exact matching is done on the
+variables on the right-hand side of \code{formula}.}
+
+ \item{data}{(required). Data frame containing the variables called in the \code{formula}.
+The dataframe should not include variables with the names \code{psclass}, \code{psweights}, or
+\code{pscore}, as these are expressly reserved in the output dataframe for MatchIt.}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).} }
+
+\value{
+ \item{psclass}{Subclass index in an ordinal scale from 1 to the number of subclasses.
+Unmatched units have subclass=0. With exact matching, these subclasses are defined by the
+categories of the variables in \code{formula}.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
+
+
+
+
+
+
+
+
diff --git a/man/help.matchit.Rd b/man/help.matchit.Rd
new file mode 100644
index 0000000..957fdcf
--- /dev/null
+++ b/man/help.matchit.Rd
@@ -0,0 +1,42 @@
+\name{help.matchit}
+
+\alias{help.matchit}
+
+\title{HTML Help for Matchit Commands and Models}
+
+\description{
+ The \code{help.matchit} command launches html help for Matchit commands
+ and supported models. The full manual is available online at
+ \url{http://gking.harvard.edu/matchit}.
+ }
+
+\usage{
+help.matchit(object)
+}
+
+\arguments{
+ \item{object}{a character string representing a Matchit command or model.
+ \code{help.matchit("command")} will take you to an index of Matchit
+ commands and \code{help.matchit("model")} will take you to a list of
+ models. }
+ }
+
+\seealso{The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{documentation}
+
+
+
+
+
+
+
+
diff --git a/man/incomedata.Rd b/man/incomedata.Rd
new file mode 100644
index 0000000..f0b288d
--- /dev/null
+++ b/man/incomedata.Rd
@@ -0,0 +1,23 @@
+\name{incomedata}
+\docType{data}
+\alias{incomedata}
+\title{Simulated data to illustrate matching procedure for Matchit}
+
+\description{
+ This data set is a fake data set to be used to illustrate matching methods in the
+Matchit library. It has data on 200 individuals, 100 treated and 100 control, with treatment
+assignemt, four covariates (income in 1974, income in 1975, male, and education), and one
+outcome variable (income in 1978) observed. }
+
+\usage{data(incomedata)}
+
+\format{
+ A data frame with 200 observations (rows) and six variables (columns). The first
+column, "treat", is the treatment assignment (1=treatment, 0=control). The second column,
+"income74", is income in 1974 (in U.S. dollars). The third column, "income75", is income in
+1975, also in U.S. dollars. The fourth column, "male", is an indicator for male (1=male,
+0=female). The fifth column, "educ", is the education level in years of schooling. The sixth
+column is the outcome variable, income in 1978 ("income78"). }
+
+\source{Simulated data.}
+\keyword{datasets}
diff --git a/man/lalonde.Rd b/man/lalonde.Rd
new file mode 100644
index 0000000..c086551
--- /dev/null
+++ b/man/lalonde.Rd
@@ -0,0 +1,31 @@
+\name{lalonde}
+\docType{data}
+\alias{lalonde}
+\title{Data from National Supported Work Demonstration and PSID, as analyzed by Dehejia and Wahba (1999).}
+
+\description{ This is a subsample of the data from the treated group in the National Supported
+Work Demonstration (NSW) and the comparison sample from the Current Population Survey (CPS).
+This data was previously analyzed extensively by Lalonde (1986) and Dehejia and Wahba (1999).
+The full dataset is available at \url{http://www.columbia.edu/~rd247/nswdata.html}. }
+
+\usage{data(incomedata)}
+
+\format{ A data frame with 313 observations (185 treated, 429 control). There are 10 variables
+measured for each individual. "treat" is the treatment assignment (1=treated, 0=control).
+"age" is age in years. "educ" is education in number of years of schooling. "black" is an
+indicator for African-American (1=African-American, 0=not). "hispan" is an indicator for being
+of Hispanic origin (1=Hispanic, 0=not). "married" is an indicator for married (1=married, 0=not
+married). "nodegree" is an indicator for whether the individual has a high school degree (1=no
+degree, 0=degree). "re74" is income in 1974, in U.S. dollars. "re75" is income in 1975, in U.S.
+dollars. "re78" is income in 1978, in U.S. dollars. }
+
+\references{
+ Lalonde, R. (1986). Evaluating the econometric evaluations of training programs
+ with experimental data. {\it American Economic Review} 76: 604-620. \\
+
+ Dehejia, R.H. and Wahba, S. (1999). Causal Effects in Nonexperimental Studies:
+ Re-Evaluating the Evaluation of Training Programs. {\it Journal of the American
+ Statistical Association} 94: 1053-1062.
+}
+\source{\url{http://www.columbia.edu/~rd247/nswdata.html}}
+\keyword{datasets}
diff --git a/man/match.data.Rd b/man/match.data.Rd
new file mode 100644
index 0000000..c7eda3b
--- /dev/null
+++ b/man/match.data.Rd
@@ -0,0 +1,37 @@
+\name{match.data}
+
+\alias{match.data}
+
+\title{Output matched data sets}
+
+\description{The code \code{match.data} creates output data sets from a \code{matchit}
+matching algorithm.}
+
+\usage{
+match.data <- match.data(object, group="all")
+}
+
+\arguments{
+\item{object}{(required). Stored output from \code{matchit}.}
+
+\item{group}{Which units to output. "all" (default) gives all matched units (treated and
+control), "treat" gives just the matched treated units, and "control" gives just the matched
+control units.} }
+
+\value{
+\item{Returns a subset of the original data set sent to \code{matchit}, with just the matched units. The data set also contains the additional
+variables \code{psclass}, \code{pscore}, and \code{psweights}. The variable \code{psclass} gives the subclass index for each unit (if applicable).
+The variable \code{pscore} gives the propensity scores, and \code{psweights} gives the weights for each unit, generated in the matching procedure.
+See the \code{matchit} documentation for more details.}
+}
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
diff --git a/man/matchdef.Rd b/man/matchdef.Rd
new file mode 100644
index 0000000..f5c78e2
--- /dev/null
+++ b/man/matchdef.Rd
@@ -0,0 +1,85 @@
+\name{matchdef}
+
+\alias{matchdef}
+
+\title{Nearest neighbor matching algorithm}
+
+\description{ This is a sub-function of \code{matchit} which performs the nearest neighbor
+matching; \code{matchdef} is called automatically by \code{matchit}. }
+
+
+\details{The matching is done using nearest neighbor matching on the propensity score, estimated
+in the sub-function \ref{distance}. This function is called directly by \code{matchit} and
+does not generally need to be called directly by users; these details are included for advanced
+users. }
+
+\usage{
+matchdef <- matchdef(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)
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. The \code{+} symbol means
+"inclusion" not "addition." You may also include interaction terms in the form if
+\code{I(X1*X2)} or squared terms in the form of \code{I(X1^2)}.}
+
+ \item{data}{Data frame containing the variables in \code{formula}.}
+ \item{in.sample}{Vector of length n showing whether each unit was eligible for matching due to
+common support restrictions with \code{discard}. Computed in \code{distance}.}
+
+ \item{pscore}{Vector of propensity scores, calculated in \code{distance}.}
+
+ \item{nearest}{Whether to perform nearest-neighbor matching (default=TRUE). }
+ \item{replace}{Whether to match with replacement (default=FALSE). }
+ \item{m.order}{Order in which to match treated units with control units. 1=optimal (requires ``optmatch" package),
+2 (default)=from high to low, 3=from low to high, 4=random order.}
+ \item{ratio}{The number of control units to be matched to each treated unit (default=1).}
+
+ \item{caliper}{Standard deviations of the propensity score within which to draw control
+units (default=0).}
+
+ \item{calclosest}{If \code{caliper!=0}, whether to take the nearest available match if
+no matches are available within \code{caliper} (default=FALSE).}
+
+ \item{mahvars}{Variables on which to perform Mahalanobis matching within each caliper
+(default=NULL). Should be entered as a vector of names of variables in \code{data}.}
+
+ \item{exact}{"FALSE" (default)=no exact matching. "TRUE"=exact matching on all
+variables in \code{formula}. A vector of variable names (that are in \code{data}) to indicate
+separate variables on which to exact match, in combination with matching on the propensity
+score.}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).} }
+
+\value{
+ \item{match.matrix}{n1 by ratio data frame where the rows correspond to treated units and the
+columns store the names of the control units matched to each treated unit. NA indicates that
+treated unit was not matched.}
+
+ \item{in.sample}{Vector of length n showing whether each unit was eligible for matching due to
+common support restrictions with \code{discard}.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{internal}
+
+
+
+
+
+
+
+
diff --git a/man/matchit.Rd b/man/matchit.Rd
new file mode 100644
index 0000000..7df07a6
--- /dev/null
+++ b/man/matchit.Rd
@@ -0,0 +1,146 @@
+\name{matchit}
+
+\alias{matchit}
+\alias{MatchIt}
+\alias{Matchit}
+\alias{i}
+
+\title{Matchit: Matching Software for Causal Inference}
+
+\description{
+\emph{Matchit} enables parametric models for causal inference to work better by selecting
+well-matched subsets of the original treated and control groups. MatchIt implements the
+suggestions of Ho, Imai, King, and Stuart (2004) for improving parametric statistical models by
+preprocessing data with nonparametric matching methods. MatchIt implements a wide range of
+sophisticated matching methods, making it possible to greatly reduce the dependence of causal
+inferences on hard-to-justify, but commonly made, statistical modeling assumptions. The
+software also easily fits into existing research practices since, after preprocessing with
+MatchIt, researchers can use whatever parametric model they would have used without MatchIt, but
+produce inferences with substantially more robustness and less sensitivity to modeling
+assumptions. Matched data sets created by MatchIt can be entered easily in Zelig
+(\url{http://gking.harvard.edu/zelig}) for subsequent parametric analyses. Full documentation is
+available online at \url{http://gking.harvard.edu/matchit}, and help for specific commands is
+available through \code{help.matchit}.}
+
+\details{The matching is done using the \code{matchit(treat ~ X, ...)} command, where
+\code{treat} is the vector of treatment assignments and \code{X} are the covariates to be used
+in the matching. There are a number of matching options, detailed below. The full syntax is
+\code{matchit(formula, data=NULL, discard=0, exact=FALSE, replace=FALSE, ratio=1, model="logit",
+reestimate=FALSE, nearest=TRUE, m.order=2, caliper=0, calclosest=FALSE, mahvars=NULL,
+subclass=0, sub.by="treat", counter=TRUE, full=FALSE, full.options=list(), \dots)} A summary of the results can be seen
+graphically
+using \code{plot(matchitobject)}, or numerically using \code{summary(matchitobject)}.
+\code{print(matchitobject)} also prints out the output. }
+
+\usage{
+matchit <- matchit(formula, data, model="logit", 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, full=FALSE, full.options=list(),...)
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. The \code{+} symbol means
+"inclusion" not "addition." You may also include interaction terms in the form if
+\code{I(X1*X2)} or squared terms in the form of \code{I(X1^2)}.}
+
+ \item{data}{(required). Data frame containing the variables called in the \code{formula}.
+The dataframe should not include variables with the names \code{psclass}, \code{psweights}, or
+\code{pscore}, as these are expressly reserved in the output dataframe for MatchIt.}
+
+ \item{model}{Method used to estimate the propensity score. May be "logit" (default),
+"probit", "nnet", "GAM", or "cart".}
+
+ \item{discard}{Whether to discard units that fall outside some measure of support of the
+distance score. 0 (default)=keep all units. 1=keep all units with common support. 2=discard
+only control units outside the support of the distance measure of the treated units. 3=discard
+only treated units outside the support of the distance measure of the control units.}
+
+ \item{reestimate}{Specifies whether to reestimate the propensity score model after
+discarding units (default=FALSE).}
+ \item{nearest}{Whether to perform nearest-neighbor matching (default=TRUE). }
+ \item{replace}{Whether to match with replacement (default=FALSE). }
+ \item{m.order}{Order in which to match treated units with control units. 1=optimal (requires ``optmatch"
+package), 2 (default)=from high to low, 3=from low to high, 4=random order.}
+
+ \item{ratio}{The number of control units to be matched to each treated unit (default=1).}
+ \item{caliper}{Standard deviations of the propensity score within which to draw control
+units (default=0).}
+
+ \item{calclosest}{If \code{caliper!=0}, whether to take the nearest available match if
+no matches are available within \code{caliper} (default=FALSE).}
+
+ \item{subclass}{Either a scaler specifying the number of subclasses (default=0) or a
+vector of probabilities to create quantiles based on \code{sub.by}.}
+
+ \item{sub.by}{If \code{subclass!=0}, by what criteria to subclassify. "treat" (default)
+=by the number of treated units, "control"=by the number of control units, "all"=by the total
+number of units.}
+
+ \item{mahvars}{Variables on which to perform Mahalanobis matching within each caliper
+(default=NULL). Should be entered as a vector of names of variables in \code{data}.}
+
+ \item{exact}{"FALSE" (default)=no exact matching. "TRUE"=exact matching on all
+variables in \code{formula}. A vector of variable names (that are in \code{data} to indicate
+separate variables on which to exact match, in combination with matching on the propensity
+score.}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).}
+
+ \item{full}{Whether to do full matching (default=FALSE). Requires ``optmatch" package.}
+ \item{full.options}{Additional options for full matching.}
+
+ \item{...}{Additional arguments to be passed to \code{matchit}, depending on the model
+to be used.}
+}
+
+\value{
+ \item{call}{The original \code{matchit} call.}
+ \item{formula}{Formula used to specify the propensity score.}
+ \item{match.matrix}{n1 by ratio data frame where the rows correspond to treated units and the
+columns store the names of the control units matched to each treated unit. NA indicates that
+treated unit was not matched.}
+
+ \item{in.sample}{Vector of length n showing whether each unit was eligible for matching due to
+common support restrictions with \code{discard}.}
+
+ \item{matched}{Vector of length n showing whether each unit was matched.}
+
+ \item{psweights}{Vector of length n giving the weight assigned to each unit in the matching
+process. Each weight is proportional to the number of times that unit was matched.}
+
+ \item{psclass}{Subclass index in an ordinal scale from 1 to the number of subclasses.
+Unmatched units have subclass=0.}
+ \item{q.cut}{Subclass cut points.}
+ \item{assign.model}{Output of the assignment model.}
+ \item{data}{The original data set, with \code{psclass}, \code{psweights}, and \code{pscore}
+(propensity scores) added.}
+ \item{treat}{The treatment indicator from \code{data}.}
+ \item{covariates}{Covariates used in the right-hand side of the assignment model.}
+}
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\references{Daniel Ho, Kosuke Imai, Gary King, and Elizabeth Stuart (2004)
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{environment}
+
+
+
+
+
+
+
+
diff --git a/man/neyman.Rd b/man/neyman.Rd
new file mode 100644
index 0000000..d6f5c04
--- /dev/null
+++ b/man/neyman.Rd
@@ -0,0 +1,47 @@
+\name{neyman}
+
+\alias{neyman}
+\alias{print.neyman}
+
+\title{Estimating treatment effects: Neyman's method}
+
+\description{Calculates basic estimate of treatment effect, using simple difference in means and
+Neyman's estimate of the variance.}
+
+\details{Calculates overall treatment effect estimate using matched samples generated by
+\code{matchit}. If subclasses were generated in the matching procedure, the estimate is a
+weighted average over the subclasses, with subclass weights defined by the \code{sub.by} option
+in the \code{matchit} call. The standard deviation can also be estimated using the bootstrap
+procedure, which is often used when matching done with replacement. For more details on the
+calculations, see the complete documentation (link below). The \code{summary} command on a
+\code{neyman} object will provide a test of significance of the outcome as well as sample sizes
+of the matched treated and control groups, and subclass estimates if applicable. }
+
+\usage{
+neyman <- neyman(Y, object, bootstrap=NULL, counter=TRUE)
+}
+
+\arguments{
+\item{Y}{(required). Outcome variable of interest, included in \code{data} from
+\code{matchit} object.}
+\item{object}{(required). Stored output from \code{matchit}.}
+\item{bootstrap}{Use the bootstrap to estimate standard errors (default=FALSE). Often used when
+matching done with replacement.}
+\item{counter}{Counter indicating progress of the bootstrap
+(default=TRUE, which shows the counter).}
+}
+
+\value{
+\item{Returns estimate of the effect of the treatment, as well as the standard deviation
+of that estimate.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{htest}
diff --git a/man/plot.matchit.Rd b/man/plot.matchit.Rd
new file mode 100644
index 0000000..0cf8529
--- /dev/null
+++ b/man/plot.matchit.Rd
@@ -0,0 +1,39 @@
+\name{plot.matchit}
+
+\alias{plot.matchit}
+
+\title{Graphing Quantities of Interest}
+
+\description{ The \code{matchit} method for the generic \code{plot}
+ command generates default plots.}
+
+\usage{
+\method{plot}{matchit}(x, ...)
+}
+
+\arguments{
+\item{x}{stored output from \code{\link{matchit}}. }
+\item{\dots}{Additional parameters passed to \code{plot.default}.
+ Because \code{plot.matchit} primarily produces diagnostic plots, many
+ of these parameters are hard-coded for convenience and
+ presentation. }
+}
+
+\value{
+\item{\code{plot.matchit} allows the user to check the distribution of all covariates
+in the assignment model, squares, and interactions, within all subclasses. The graphs present
+density estimate graphs of the propensity score of treated and control units in the full and
+matched samples, jitter plots of the propensity score for treated and control units, density
+estimate graphs of any covariates, and density estimate graphs of any covariates by subclass.}}
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{hplot}
diff --git a/man/print.matchit.Rd b/man/print.matchit.Rd
new file mode 100644
index 0000000..591c692
--- /dev/null
+++ b/man/print.matchit.Rd
@@ -0,0 +1,36 @@
+\name{print.matchit}
+
+\alias{print.matchit}
+
+\title{Printing Quantities of Interest}
+
+\description{ The \code{matchit} method for the generic \code{print}
+ command generates default output.
+}
+
+\usage{
+\method{print}{matchit}(x, digits = max(3, getOption("digits") - 3), ...)
+}
+
+\arguments{
+\item{x}{stored output from \code{matchit}. }
+\item{digits}{Optionally specify the significant digits to be displayed.}
+\item{\dots}{Additional parameters passed to \code{print.default}.}
+}
+
+\value{
+\item{Returns the original assignment model call, summary statistics of the propensity
+score for the full and matched samples, summary statistics of the propensity score for each
+subclass (if applicable), and sample sizes for the full and matched samples (full and by
+subclass).} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
diff --git a/man/subclassify.Rd b/man/subclassify.Rd
new file mode 100644
index 0000000..264d86d
--- /dev/null
+++ b/man/subclassify.Rd
@@ -0,0 +1,85 @@
+\name{subclassify}
+
+\alias{subclassify}
+
+\title{Subclassify units}
+
+\description{This is a sub-function of \code{matchit} which generates subclasses of units
+given a propensity score.}
+
+\details{This function creates subclasses of units based on their covariate values, usually
+defined by the propensity score. This function is called directly by \code{matchit} and does
+not generally need to be called directly by users; these details are included for advanced
+users. }
+
+\usage{
+subclassify <- subclassify(formula, data, in.sample, pscore, nearest=TRUE, match.matrix,
+ subclass=0, sub.by="treat", counter=TRUE, full=FALSE, full.options=list())
+}
+
+
+\arguments{
+ \item{formula}{(required). Takes the form of \code{T ~ X1 + X2}, where \code{T} is a binary
+treatment indicator and \code{X1} and \code{X2} are the pre-treatment covariates, and \code{T},
+\code{X1}, and \code{X2} are contained in the same data frame. The \code{+} symbol means
+"inclusion" not "addition." You may also include interaction terms in the form if
+\code{I(X1*X2)} or squared terms in the form of \code{I(X1^2)}.}
+
+ \item{data}{(required). Data frame containing the variables called in the \code{formula}.
+The dataframe should not include variables with the names \code{psclass}, \code{psweights}, or
+\code{pscore}, as these are expressly reserved in the output dataframe for MatchIt.}
+
+ \item{in.sample}{(required). Vector of length n showing whether each unit was eligible for
+matching due to common support restrictions with \code{discard}. Generally created by
+\code{distance}.}
+
+ \item{pscore}{(required). Vector of estimated propensity scores. Generally created by
+\code{distance}.}
+
+ \item{nearest}{Whether to perform nearest-neighbor matching (default=TRUE). }
+
+ \item{match.matrix}{(required) n1 by ratio data frame where the rows correspond to treated
+units and the columns store the names of the control units matched to each treated unit. NA
+indicates that treated unit was not matched. Generally created by \code{matchdef}.}
+
+ \item{subclass}{Either a scaler specifying the number of subclasses (default=0) or a
+vector of probabilities to create quantiles based on \code{sub.by}.}
+
+ \item{sub.by}{If \code{subclass!=0}, by what criteria to subclassify. "treat" (default)
+=by the number of treated units, "control"=by the number of control units, "all"=by the total
+number of units.}
+
+ \item{counter}{Whether to display counter indicating the progress of the matching
+(default=TRUE).}
+
+ \item{full}{Whether to do full matching (default=FALSE). Requires ``optmatch" package.}
+ \item{full.options}{Additional options for full matching.}
+
+
+}
+
+\value{
+ \item{psclass}{Subclass index in an ordinal scale from 1 to the number of subclasses.
+Unmatched units have subclass=0.}
+ \item{q.cut}{Subclass cut points.}
+}
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.
+}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{internal}
+
+
+
+
+
+
+
+
diff --git a/man/summary.matchit.Rd b/man/summary.matchit.Rd
new file mode 100644
index 0000000..08c7098
--- /dev/null
+++ b/man/summary.matchit.Rd
@@ -0,0 +1,65 @@
+\name{summary.matchit}
+
+\alias{summary.matchit}
+\alias{print.summary.matchit}
+
+\title{Summarizing Quantities of Interest}
+
+\description{ The \code{matchit} method for the generic \code{summary}
+ command generates default output.
+}
+
+\usage{
+\method{summary}{matchit}(object, verbose=F, sig=2, ...)
+}
+
+\arguments{
+\item{object}{stored output from \code{matchit}.}
+
+\item{sig}{Optionally specifies the size of the t-statistic for a covariate to be labeled
+"problematic" (default=2).}
+
+\item{verbose}{Option to show summary statistics in \code{sum.all} and \code{sum.matched} for
+all covariates, their squares, and two-way interactions when \code{verbose=TRUE} and only the
+covariates themselves when \code{verbose=FALSE} (default).}
+
+\item{\dots}{Additional parameters
+passed to \code{summary.default}.} }
+
+\value{
+\item{call}{The original assignment model call.}
+
+\item{sum.all}{Data frame containing summary statistics on all observations in the treated and
+control groups (means, standard deviations, t-statistics and bias statistics for all
+covariates).}
+
+\item{sum.matched}{Data frame containing summary statistics on the matched observations in the
+treated and control groups (means, standard deviations, t-statistics, bias statistics, and
+whether bias was reduced by the matching).}
+
+\item{Problematic}{Problematic covariates that remain imbalanced in the matched samples.
+Variables listed when the absolute value of the t-statistic in the matched sample exceeds
+\code{sig}.}
+
+\item{xn}{Sample sizes in the full and matched samples.}
+
+\item{q.table}{Array that contains the same information as \code{sum.all} and \code{sum.matched}
+for each subclass.}
+
+\item{qn}{Sample sizes in the full and matched samples, by subclass.}
+
+\item{match.matrix}{From the \code{matchit} output.}
+
+\item{psclass}{From the \code{matchit} output.}
+\item{in.sample}{From the \code{matchit} output.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
diff --git a/man/summary.neyman.Rd b/man/summary.neyman.Rd
new file mode 100644
index 0000000..18134f9
--- /dev/null
+++ b/man/summary.neyman.Rd
@@ -0,0 +1,44 @@
+\name{summary.neyman}
+
+\alias{summary.neyman}
+\alias{print.summary.neyman}
+
+\title{Summarizing the estimated treatment effects}
+
+\description{Calculates basic estimate of treatment effect, using simple difference in means and
+Neyman's estimate of the variance. Also provides test of significance and subclass estimates if
+applicable.}
+
+\details{Calculates overall treatment effect estimate using matched samples generated by
+\code{matchit}. Also gives test of significance of the estimated treatment effect, and subclass
+estimates, if applicable. If subclasses were generated in the matching procedure, an overall
+estimate is also provided, which is a weighted average over the subclasses, with subclass
+weights defined by the \code{sub.by} option in the \code{matchit} call. The standard deviation
+can also be estimated using the bootstrap procedure, which is often used when matching done with
+replacement. For more details on the calculations, see the complete documentation (link below).
+}
+
+\usage{
+\method{summary}{neyman}(object, ...)
+}
+
+\arguments{
+\item{object}{(required). Stored output from \code{neyman}. }
+\item{...}{Further arguments passed to or from other methods.}
+}
+
+\value{ \item{Returns estimate of the effect of the treatment, as well as the standard deviation
+of that estimate, a test of significance of the estimate, and subclass estimates if applicable.
+Also shows the sample sizes of the matched treated and control groups used to estimate the
+treatment effect.} }
+
+\seealso{Please use \code{help.matchit} to access the matchit reference
+ manual. The complete document is available online at
+ \url{http://gking.harvard.edu/matchit}.}
+
+\author{
+ Daniel Ho <\email{deho at fas.harvard.edu}>; Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Elizabeth Stuart<\email{stuart at stat.harvard.edu}>
+}
+
+\keyword{methods}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-matchit.git
More information about the debian-science-commits
mailing list