1# Jasjeet S. Sekhon <sekhon@berkeley.edu> 2# HTTP://sekhon.berkeley.edu/ 3# UC Berkeley 4 5# Match(): function to estimate treatments using a matching estimator. 6# Currently only the ability to estimate average treatment effects 7# using the approach of Abadie and Imbens is implemented. In the 8# future, quantile treatment effects will be implemented along with 9# the ability to use robust estimation when estimating the propensity 10# score. MatchBalance(), and balanceUV() test for balance. 11 12Match <- function(Y=NULL,Tr,X,Z=X,V=rep(1,length(Y)), estimand="ATT", M=1, 13 BiasAdjust=FALSE,exact=NULL,caliper=NULL, replace=TRUE, ties=TRUE, 14 CommonSupport=FALSE,Weight=1,Weight.matrix=NULL, weights=NULL, 15 Var.calc=0, sample=FALSE, restrict=NULL, match.out=NULL, 16 distance.tolerance=0.00001, tolerance=sqrt(.Machine$double.eps), 17 version="standard") 18 { 19 20 BiasAdj <- as.double(BiasAdjust) 21 sample <- as.double(sample) 22 23 if ( (BiasAdj != 0) & (BiasAdj != 1) ) 24 { 25 warning("User set 'BiasAdjust' to a non-logical value. Resetting to the default which is FALSE.") 26 BiasAdj <- 0 27 } 28 29 #we don't need to use a Y 30 if (is.null(Y)) 31 { 32 Y = rep(0, length(Tr)) 33 version <- "fast" 34 35 if(BiasAdj) 36 { 37 warning("'BiasAdjust' set to FALSE because Y is NULL") 38 BiasAdj <- FALSE 39 } 40 } 41 42 Y <- as.double(Y) 43 Tr <- as.double(Tr) 44 X <- as.matrix(X) 45 Z <- as.matrix(Z) 46 V <- as.matrix(V) 47 48 orig.nobs <- length(Y) 49 nobs <- orig.nobs 50 xvars <- ncol(X) 51 orig.tr.nobs <- length(Tr) 52 if (orig.tr.nobs != orig.nobs) 53 { 54 stop("length(Y) != length(Tr)") 55 } 56 if( orig.tr.nobs != nrow(X)) 57 { 58 stop("length(Tr) != nrow(X)") 59 } 60 61 if( orig.nobs != nrow(X)) 62 { 63 stop("length(Y) != nrow(X)") 64 } 65 if( orig.nobs != nrow(V)) 66 { 67 stop("length(Y) != nrow(V)") 68 } 69 if( orig.nobs != nrow(Z)) 70 { 71 stop("length(Y) != nrow(Z)") 72 } 73 74 if (is.null(weights)) 75 { 76 weights <- rep(1,length(Y)) 77 weights.flag <- FALSE 78 } else { 79 weights.flag <- TRUE 80 weights <- as.double(weights) 81 if( orig.tr.nobs != length(weights)) 82 { 83 stop("length(Tr) != length(weights)") 84 } 85 } 86 87 isna <- sum(is.na(Y)) + sum(is.na(Tr)) + sum(is.na(X)) + sum(is.na(weights)) + sum(is.na(Z)) + sum(is.na(V)) 88 if (isna!=0) 89 { 90 stop("Match(): input includes NAs") 91 return(invisible(NULL)) 92 } 93 94 if (sum(Tr !=1 & Tr !=0) > 0) { 95 stop("Treatment indicator ('Tr') must be a logical variable---i.e., TRUE (1) or FALSE (0)") 96 } 97 if (var(Tr)==0) { 98 stop("Treatment indicator ('Tr') must contain both treatment and control observations") 99 } 100 if (distance.tolerance < 0) 101 { 102 warning("User set 'distance.tolerance' to less than 0. Resetting to the default which is 0.00001.") 103 distance.tolerance <- 0.00001 104 } 105 106 #CommonSupport 107 if (CommonSupport !=1 & CommonSupport !=0) { 108 stop("'CommonSupport' must be a logical variable---i.e., TRUE (1) or FALSE (0)") 109 } 110 if(CommonSupport==TRUE) 111 { 112 tr.min <- min(X[Tr==1,1]) 113 tr.max <- max(X[Tr==1,1]) 114 115 co.min <- min(X[Tr==0,1]) 116 co.max <- max(X[Tr==0,1]) 117 118 if(tr.min >= co.min) 119 { 120 indx1 <- X[,1] < (tr.min-distance.tolerance) 121 } else { 122 indx1 <- X[,1] < (co.min-distance.tolerance) 123 } 124 125 if(co.max <= tr.max) 126 { 127 indx2 <- X[,1] > (co.max+distance.tolerance) 128 } else { 129 indx2 <- X[,1] > (tr.max+distance.tolerance) 130 } 131 132 indx3 <- indx1==0 & indx2==0 133 Y <- as.double(Y[indx3]) 134 Tr <- as.double(Tr[indx3]) 135 X <- as.matrix(X[indx3,]) 136 Z <- as.matrix(Z[indx3,]) 137 V <- as.matrix(V[indx3,]) 138 weights <- as.double(weights[indx3]) 139 140 #let's recalculate these for common support 141 orig.nobs <- length(Y) 142 nobs <- orig.nobs 143 }#end of CommonSupport 144 145 #check additional inputs 146 if (tolerance < 0) 147 { 148 warning("User set 'tolerance' to less than 0. Resetting to the default which is 0.00001.") 149 tolerance <- 0.00001 150 } 151 if (M < 1) 152 { 153 warning("User set 'M' to less than 1. Resetting to the default which is 1.") 154 M <- 1 155 } 156 if ( M != round(M) ) 157 { 158 warning("User set 'M' to an illegal value. Resetting to the default which is 1.") 159 M <- 1 160 } 161 if (Var.calc < 0) 162 { 163 warning("User set 'Var.calc' to less than 0. Resetting to the default which is 0.") 164 Var.calc <- 0 165 } 166 if ( (sample != 0) & (sample != 1) ) 167 { 168 warning("User set 'sample' to a non-logical value. Resetting to the default which is FALSE.") 169 sample <- 0 170 } 171 if (Weight != 1 & Weight != 2 & Weight != 3) 172 { 173 warning("User set 'Weight' to an illegal value. Resetting to the default which is 1.") 174 Weight <- 1 175 } 176 if (version!="fast" & version != "standard" & version != "legacy" & version != "Matchby" & version != "MatchbyAI") 177 { 178 warning("User set 'version' to an illegal value. Resetting to the default which is 'standard'.") 179 version <- "standard" 180 } 181 182 if(version=="Matchby") 183 { 184 version <- "fast" 185 Matchby.call <- TRUE 186 MatchbyAI <- FALSE 187 } else if(version=="MatchbyAI") 188 { 189 version <- "standard" 190 Matchby.call <- TRUE 191 MatchbyAI <- TRUE 192 } else { 193 Matchby.call <- FALSE 194 MatchbyAI <- FALSE 195 } 196 197 if (Var.calc !=0 & version=="fast") 198 { 199 warning("Var.calc cannot be estimate when version=='fast'") 200 Var.calc=0 201 } 202 if (BiasAdj!=FALSE & version=="fast") 203 { 204 warning("Bias Adjustment cannot be estimated when version=='fast'") 205 BiasAdj=0 206 } 207 if (replace!=FALSE & replace!=TRUE) 208 { 209 warning("'replace' must be TRUE or FALSE. Setting to TRUE") 210 replace <- TRUE 211 } 212 if(replace==FALSE) 213 { 214 ties <- FALSE 215 version="fast" 216 if (version=="legacy") 217 warning("'version' is set to 'fast' because replace==FALSE") 218 } 219 if (ties!=FALSE & ties!=TRUE) 220 { 221 warning("'ties' must be TRUE or FALSE. Setting to TRUE") 222 ties <- TRUE 223 } 224 if(ties==FALSE) 225 { 226 version="fast" 227 if (version=="legacy") 228 warning("'version' is set to 'fast' because ties==FALSE") 229 230 if(BiasAdjust==TRUE) { 231 warning("Bias Adjustment can only be estimated when ties==TRUE and replace=TRUE. Setting BiasAdjust=FALSE") 232 BiasAdjust <- FALSE 233 BiasAdj <- 0 234 } 235 } 236 237 if (!is.null(match.out) & class(match.out) != "Match") { 238 warning("match.out object not of class 'Match'") 239 return(invisible(NULL)) 240 } 241 242 ccc <- tolerance 243 cdd <- distance.tolerance 244 245 orig.treated.nobs <- sum(Tr==1) 246 orig.control.nobs <- sum(Tr==0) 247 orig.wnobs <- sum(weights) 248 orig.weighted.treated.nobs <- sum( weights[Tr==1] ) 249 orig.weighted.control.nobs <- sum( weights[Tr==0] ) 250 weights.orig <- as.matrix(weights) 251 zvars <- ncol(Z); 252 253 estimand.orig <- estimand 254 if (estimand=="ATT") 255 { 256 estimand <- 0 257 258 if(BiasAdj==1 & orig.treated.nobs<zvars) 259 { 260 warning("Fewer treated obs than variables in 'Z': BiasAdjust set to FALSE") 261 BiasAdj=0 262 } 263 264 } else if(estimand=="ATE") { 265 estimand <- 1 266 267 if(BiasAdj==1 & orig.nobs<zvars) 268 { 269 warning("Fewer obs than variables in 'Z': BiasAdjust set to FALSE") 270 BiasAdj=0 271 } 272 } else if(estimand=="ATC") { 273 estimand <- 2 274 275 if(BiasAdj==1 & orig.control.nobs<zvars) 276 { 277 warning("Fewer control obs than variables in 'Z': BiasAdjust set to FALSE") 278 BiasAdj=0 279 } 280 281 } else { 282 estimand <- 0 283 warning("User set 'estimand' to an illegal value. Resetting to the default which is 'ATT'") 284 } 285 286 if (!is.null(Weight.matrix)) 287 { 288 289 if(inherits(Weight.matrix, "GenMatch")) 290 { 291 Weight.matrix = Weight.matrix$Weight.matrix 292 } 293 294 if (Weight==2) 295 { 296 warning("User supplied 'Weight.matrix' is being used even though 'Weight' is not set equal to 3") 297 } 298 Weight <- 3 299 } else { 300 Weight.matrix <- dim(X)[2] 301 } 302 303 if(Var.calc > orig.weighted.treated.nobs) 304 { 305 warning("'Var.calc' > the number of treated obs: 'Var.calc' reset to ", 306 orig.weighted.treated.nobs, immediate.=Matchby.call) 307 Var.calc <- orig.weighted.treated.nobs 308 } 309 if(Var.calc > orig.weighted.control.nobs) 310 { 311 warning("'Var.calc' > the number of control obs: 'Var.calc' reset to ", 312 orig.weighted.control.nobs, immediate.=Matchby.call) 313 Var.calc <- orig.weighted.control.nobs 314 } 315 316 if(orig.nobs > 20000 & version!="fast" & !Matchby.call) 317 { 318 warning("The version='fast' option is recommended for large datasets if speed is desired. For additional speed, you may also consider using the ties=FALSE option.", immediate.=TRUE) 319 } 320 321 #check the restrict matrix input 322 if(!is.null(restrict)) 323 { 324 if(!is.matrix(restrict)) 325 stop("'restrict' must be a matrix of restricted observations rows and three columns: c(i,j restriction)") 326 327 if(ncol(restrict)!=3 ) 328 stop("'restrict' must be a matrix of restricted observations rows and three columns: c(i,j restriction)") 329 } 330 331 if (!is.null(exact)) 332 { 333 exact = as.vector(exact) 334 nexacts = length(exact) 335 if ( (nexacts > 1) & (nexacts != xvars) ) 336 { 337 warning("length of exact != ncol(X). Ignoring exact option") 338 exact <- NULL 339 } else if (nexacts==1 & (xvars > 1) ){ 340 exact <- rep(exact, xvars) 341 } 342 } 343 344 if (!is.null(caliper)) 345 { 346 caliper = as.vector(caliper) 347 ncalipers = length(caliper) 348 if ( (ncalipers > 1) & (ncalipers != xvars) ) 349 { 350 warning("length of caliper != ncol(X). Ignoring caliper option") 351 caliper <- NULL 352 } else if (ncalipers==1 & (xvars > 1) ){ 353 caliper <- rep(caliper, xvars) 354 } 355 } 356 357 if (!is.null(caliper)) 358 { 359 ecaliper <- vector(mode="numeric", length=xvars) 360 sweights <- sum(weights.orig) 361 for (i in 1:xvars) 362 { 363 meanX <- sum( X[,i]*weights.orig )/sweights 364 sdX <- sqrt(sum( (X[,i]-meanX)^2 )/sweights) 365 ecaliper[i] <- caliper[i]*sdX 366 } 367 } else { 368 ecaliper <- NULL 369 } 370 371 if (!is.null(exact)) 372 { 373 if(is.null(caliper)) 374 { 375 max.diff <- abs(max(X)-min(X) + tolerance * 100) 376 ecaliper <- matrix(max.diff, nrow=xvars, ncol=1) 377 } 378 379 for (i in 1:xvars) 380 { 381 if (exact[i]) 382 ecaliper[i] <- tolerance; 383 } 384 } 385 386 if(replace==FALSE) 387 { 388 #replace==FALE, needs enough observation 389 #ATT 390 orig.weighted.control.nobs <- sum(weights[Tr!=1]) 391 if(estimand==0) 392 { 393 if (orig.weighted.treated.nobs > orig.weighted.control.nobs) 394 { 395 warning("replace==FALSE, but there are more (weighted) treated obs than control obs. Some treated obs will not be matched. You may want to estimate ATC instead.") 396 } 397 } else if(estimand==1) 398 { 399 #ATE 400 if (orig.weighted.treated.nobs > orig.weighted.control.nobs) 401 { 402 warning("replace==FALSE, but there are more (weighted) treated obs than control obs. Some treated obs will not be matched. You may want to estimate ATC instead.") 403 } 404 if (orig.weighted.treated.nobs < orig.weighted.control.nobs) 405 { 406 warning("replace==FALSE, but there are more (weighted) control obs than treated obs. Some control obs will not be matched. You may want to estimate ATT instead.") 407 } 408 } else 409 { 410 #ATC 411 if (orig.weighted.treated.nobs < orig.weighted.control.nobs) 412 { 413 warning("replace==FALSE, but there are more (weighted) control obs than treated obs. Some obs will be dropped. You may want to estimate ATC instead") 414 } 415 } 416 417 #we need a restrict matrix if we are going to not do replacement 418 if(is.null(restrict)) 419 { 420 restrict <- t(as.matrix(c(0,0,0))) 421 } 422 if(version!="fast" & version!="standard") 423 { 424 warning("reverting to 'standard' version because replace=FALSE") 425 version="standard" 426 } 427 }#end of replace==FALSE 428 429# if(version=="fast" & is.null(ecaliper) & sum(weights==1)==orig.nobs) 430 if(version=="fast" | version=="standard") 431 { 432 if(!is.null(match.out)) 433 { 434 ret <- RmatchLoop(Y=Y, Tr=Tr, X=X, Z=Z, V=V, All=estimand, M=M, BiasAdj=BiasAdj, 435 Weight=Weight, Weight.matrix=Weight.matrix, Var.calc=Var.calc, 436 weight=weights, SAMPLE=sample, ccc=ccc, cdd=cdd, 437 ecaliper=ecaliper, exact=exact, caliper=caliper, 438 restrict=restrict, MatchLoopC.indx=match.out$MatchLoopC, 439 weights.flag=weights.flag, replace=replace, ties=ties, 440 version=version, MatchbyAI=MatchbyAI) 441 } else { 442 ret <- RmatchLoop(Y=Y, Tr=Tr, X=X, Z=Z, V=V, All=estimand, M=M, BiasAdj=BiasAdj, 443 Weight=Weight, Weight.matrix=Weight.matrix, Var.calc=Var.calc, 444 weight=weights, SAMPLE=sample, ccc=ccc, cdd=cdd, 445 ecaliper=ecaliper, exact=exact, caliper=caliper, 446 restrict=restrict, weights.flag=weights.flag, replace=replace, ties=ties, 447 version=version, MatchbyAI=MatchbyAI) 448 } 449 } else { 450 ret <- Rmatch(Y=Y, Tr=Tr, X=X, Z=Z, V=V, All=estimand, M=M, BiasAdj=BiasAdj, 451 Weight=Weight, Weight.matrix=Weight.matrix, Var.calc=Var.calc, 452 weight=weights, SAMPLE=sample, ccc=ccc, cdd=cdd, 453 ecaliper=ecaliper, restrict=restrict) 454 } 455 456 if(is.null(ret$est)) 457 { 458 if(!Matchby.call) 459 { 460 if(ret$valid < 1) 461 { 462 if (ret$sum.caliper.drops > 0) { 463 warning("'Match' object contains no valid matches (probably because of the caliper or the exact option).") 464 } else { 465 warning("'Match' object contains no valid matches") 466 } 467 } else { 468 if (ret$sum.caliper.drops > 0) { 469 warning("'Match' object contains only 1 valid match (probably because of the caliper or the exact option).") 470 } else { 471 warning("'Match' object contains only one valid match") 472 } 473 } 474 } #endof if(!Matchby.call) 475 476 z <- NA 477 class(z) <- "Match" 478 return(z) 479 } 480 481 indx <- cbind(ret$art.data[,1], ret$art.data[,2], ret$W) 482 483 index.treated <- indx[,1] 484 index.control <- indx[,2] 485 weights <- indx[,3] 486 sum.caliper.drops <- ret$sum.caliper.drops 487 488 #RESET INDEX.TREATED 489 indx <- as.matrix(cbind(index.treated,index.control)) 490 if (estimand==0) { 491 #"ATT" 492 index.treated <- indx[,1] 493 index.control <- indx[,2] 494 } else if(estimand==1) { 495 #"ATE" 496 tmp.index.treated <- indx[,1] 497 tmp.index.control <- indx[,2] 498 499 tl <- length(tmp.index.treated) 500 index.treated <- vector(length=tl, mode="numeric") 501 index.control <- vector(length=tl, mode="numeric") 502 trt <- Tr[tmp.index.treated]==1 503 for (i in 1:tl) 504 { 505 if (trt[i]) { 506 index.treated[i] <- tmp.index.treated[i] 507 index.control[i] <- tmp.index.control[i] 508 } else { 509 index.treated[i] <- tmp.index.control[i] 510 index.control[i] <- tmp.index.treated[i] 511 } 512 } 513 } else if(estimand==2) { 514 #"ATC" 515 index.treated <- indx[,2] 516 index.control <- indx[,1] 517 } 518 519 520 mdata <- list() 521 mdata$Y <- c(Y[index.treated],Y[index.control]) 522 mdata$Tr <- c(Tr[index.treated],Tr[index.control]) 523 mdata$X <- rbind(X[index.treated,],X[index.control,]) 524 mdata$orig.weighted.treated.nobs <- orig.weighted.treated.nobs 525 526 #naive standard errors 527 mest <- sum((Y[index.treated]-Y[index.control])*weights)/sum(weights) 528 v1 <- Y[index.treated] - Y[index.control] 529 varest <- sum( ((v1-mest)^2)*weights)/(sum(weights)*sum(weights)) 530 se.standard <- sqrt(varest) 531 532 wnobs <- sum(weights) 533 534 if(estimand==0) 535 { 536 #ATT 537 actual.drops <- orig.weighted.treated.nobs-wnobs 538 } else if (estimand==1) 539 { 540 #ATE 541 actual.drops <- orig.wnobs-wnobs 542 } else { 543 #ATC 544 actual.drops <- (orig.wnobs-orig.weighted.treated.nobs)-wnobs 545 } 546 547 #What obs were dropped? 548 index.dropped <- NULL #nothing was dropped 549 if (sum.caliper.drops > 0 ) 550 { 551 if(estimand.orig=="ATT") 552 { 553 matched.index <- which(Tr==1) 554 matched <- !(matched.index %in% index.treated) 555 556 } else if(estimand.orig=="ATC") 557 { 558 matched.index <- which(Tr==0) 559 matched <- !(matched.index %in% index.control) 560 } else if(estimand.orig=="ATE") 561 { 562 matched.index <- 1:length(Tr) 563 matched <- !(matched.index %in% c(index.treated,index.control)) 564 } 565 index.dropped <- matched.index[matched] #obs not matched 566 } #end of sum.caliper.drops > 0 567 568 z <- list(est=ret$est, se=ret$se, est.noadj=mest, se.standard=se.standard, 569 se.cond=ret$se.cond, 570 mdata=mdata, 571 index.treated=index.treated, index.control=index.control, index.dropped=index.dropped, 572 weights=weights, orig.nobs=orig.nobs, orig.wnobs=orig.wnobs, 573 orig.treated.nobs=orig.treated.nobs, 574 nobs=nobs, wnobs=wnobs, 575 caliper=caliper, ecaliper=ecaliper, exact=exact, 576 ndrops=actual.drops, ndrops.matches=sum.caliper.drops, 577 MatchLoopC=ret$MatchLoopC, version=version, 578 estimand=estimand.orig) 579 580 if(MatchbyAI) 581 { 582 z$YCAUS <- ret$YCAUS 583 z$ZCAUS <- ret$ZCAUS 584 z$Kcount <- ret$Kcount 585 z$KKcount <- ret$KKcount 586 z$Sigs <- ret$Sigs 587 } 588 589 class(z) <- "Match" 590 return(z) 591 } #end of Match 592 593summary.Match <- function(object, ..., full=FALSE, digits=5) 594 { 595 if(!is.list(object)) { 596 warning("'Match' object contains less than two valid matches. Cannot proceed.") 597 return(invisible(NULL)) 598 } 599 600 if (class(object) != "Match") { 601 warning("Object not of class 'Match'") 602 return(invisible(NULL)) 603 } 604 605 if(object$version!="fast") 606 { 607 cat("\n") 608 cat("Estimate... ",format(object$est,digits=digits),"\n") 609 cat("AI SE...... ",format(object$se,digits=digits),"\n") 610 cat("T-stat..... ",format(object$est/object$se,digits=digits),"\n") 611 cat("p.val...... ",format.pval((1-pnorm(abs(object$est/object$se)))*2,digits=digits),"\n") 612 cat("\n") 613 } else { 614 cat("\n") 615 cat("Estimate... ",format(object$est,digits=digits),"\n") 616 cat("SE......... ",format(object$se.standard,digits=digits),"\n") 617 cat("T-stat..... ",format(object$est/object$se.standard,digits=digits),"\n") 618 cat("p.val...... ",format.pval((1-pnorm(abs(object$est/object$se.standard)))*2,digits=digits),"\n") 619 cat("\n") 620 } 621 622 if(full) 623 { 624 cat("Est noAdj.. ",format(object$est.noadj,digits=digits),"\n") 625 cat("SE......... ",format(object$se.standard,digits=digits),"\n") 626 cat("T-stat..... ",format(object$est.noadj/object$se.standard,digits=digits),"\n") 627 cat("p.val...... ",format.pval((1-pnorm(abs(object$est.noadj/object$se.standard)))*2,digits=digits),"\n") 628 cat("\n") 629 } 630 631 632 if(object$orig.wnobs!=object$orig.nobs) 633 cat("Original number of observations (weighted)... ", round(object$orig.wnobs, 3),"\n") 634 cat("Original number of observations.............. ", object$orig.nobs,"\n") 635 if(object$mdata$orig.weighted.treated.nobs!=object$orig.treated.nobs) 636 cat("Original number of treated obs (weighted).... ", round(object$mdata$orig.weighted.treated.nobs, 3),"\n") 637 if(object$estimand!="ATC") 638 { 639 cat("Original number of treated obs............... ", object$orig.treated.nobs,"\n") 640 } else { 641 cat("Original number of control obs............... ", 642 object$orig.nobs-object$orig.treated.nobs,"\n") 643 } 644 cat("Matched number of observations............... ", round(object$wnobs, 3),"\n") 645 cat("Matched number of observations (unweighted). ", length(object$index.treated),"\n") 646 647 cat("\n") 648 if(!is.null(object$exact)) 649 { 650 cat("Number of obs dropped by 'exact' or 'caliper' ", object$ndrops.matches, "\n") 651 if (object$ndrops.matches!=round(object$ndrops)) 652 cat("Weighted #obs dropped by 'exact' or 'caliper' ", round(object$ndrops, 3),"\n") 653 cat("\n") 654 }else if(!is.null(object$caliper)) 655 { 656 cat("Caliper (SDs)........................................ ",object$caliper,"\n") 657 cat("Number of obs dropped by 'exact' or 'caliper' ", object$ndrops.matches, "\n") 658 if (object$ndrops.matches!=round(object$ndrops)) 659 cat("Weighted #obs dropped by 'exact' or 'caliper' ", round(object$ndrops, 3),"\n") 660 cat("\n") 661 } 662 663 z <- list() 664 class(z) <- "summary.Match" 665 return(invisible(z)) 666 } #end of summary.Match 667 668print.summary.Match <- function(x, ...) 669 { 670 invisible(NULL) 671 } 672 673 674Rmatch <- function(Y, Tr, X, Z, V, All, M, BiasAdj, Weight, Weight.matrix, Var.calc, weight, 675 SAMPLE, ccc, cdd, ecaliper=NULL, restrict=NULL) 676 { 677 sum.caliper.drops <- 0 678 X.orig <- X 679 680#are we using the restriction matrix? 681 if(is.matrix(restrict)) { 682 restrict.trigger <- TRUE 683 } else { 684 restrict.trigger <- FALSE 685 } 686 687# if SATC is to be estimated the treatment indicator is reversed 688 if (All==2) { 689 Tr <- 1-Tr 690 } 691 692 693 694# check on the number of matches, to make sure the number is within the limits 695# feasible given the number of observations in both groups. 696 if (All==1) 697 { 698 M <- min(M,min(sum(Tr),sum(1-Tr))); 699 } else { 700 M <- min(M,sum(1-Tr)); 701 } 702 703# two slippage parameters that are used to determine whether distances are equal 704# distances less than ccc or cdd are interpeted as zero. 705# these are passed in. ccc, cdd 706 707 708# I. set up 709# I.a. vector for which observations we want the average effect 710# iot_t is the vector with weights in the average treatment effects 711# iot_c is the vector of indicators for potential controls 712 713 if (All==1) 714 { 715 iot.t <- weight; 716 iot.c <- as.matrix(rep(1,length(Tr))) 717 } else { 718 iot.t <- Tr*weight; 719 iot.c <- 1-Tr 720 } 721 722# I.b. determine sample and covariate vector sizes 723 N <- nrow(X) 724 Kx <- ncol(X) 725 Kz <- ncol(Z) 726 727# K covariates 728# N observations 729# Nt <- sum(Tr) 730# Nc <- sum(1-Tr) 731# on <- as.matrix(rep(1,N)) 732 733# I.c. normalize regressors to have mean zero and unit variance. 734# If the standard deviation of a variable is zero, its normalization 735# leads to a variable with all zeros. 736# The matrix AA enables one to transform the user supplied weight matrix 737# to take account of this transformation. BUT THIS IS NOT USED!! 738# Mu_X and Sig_X keep track of the original mean and variances 739# AA <- diag(Kx) 740 Mu.X <- matrix(0, Kx, 1) 741 Sig.X <- matrix(0, Kx, 1) 742 743 for (k in 1:Kx) 744 { 745 Mu.X[k,1] <- sum(X[,k]*weight)/sum(weight) 746 eps <- X[,k]-Mu.X[k,1] 747 Sig.X[k,1] <- sqrt(max(ccc, sum(X[,k]*X[,k]*weight))/sum(weight)-Mu.X[k,1]^2) 748 Sig.X[k,1] <- Sig.X[k,1]*sqrt(N/(N-1)) 749 if(Sig.X[k,1] < ccc) 750 Sig.X[k,1] <- ccc 751 X[,k]=eps/Sig.X[k,1] 752# AA[k,k]=Sig.X[k,1] 753 } #end of k loop 754 755# Nv <- nrow(V) 756 Mv <- ncol(V) 757 Mu.V <- matrix(0, Mv, 1) 758 Sig.V <- matrix(0, Mv, 1) 759 760 for (j in 1:Mv) 761 { 762 Mu.V[j,1]= ( t(V[,j])%*%weight ) /sum(weight) 763# dv <- V[,j]-Mu.V[j,1] 764 sv <- sum(V[,j]*V[,j]*weight)/sum(weight)-Mu.V[j,1]^2 765 if (sv > 0) 766 { 767 sv <- sqrt(sv) 768 } else { 769 sv <- 0 770 } 771 sv <- sv * sqrt(N/(N-1)) 772 Sig.V[j,1] <- sv 773 } #end of j loop 774 775# I.d. define weight matrix for metric, taking into account normalization of 776# regressors. 777# If the eigenvalues of the regressors are too close to zero the Mahalanobis metric 778# is not used and we revert back to the default inverse of variances. 779 if (Weight==1) 780 { 781 Weight.matrix=diag(Kx) 782 } else if (Weight==2) { 783 if (min (eigen( t(X)%*%X/N, only.values=TRUE)$values) > ccc) 784 { 785 Weight.matrix= solve(t(X)%*%X/N) 786 } else { 787 Weight.matrix <- diag(Kx) 788 } 789 } 790 # DO NOT RESCALE THE Weight.matrix!! 791 #else if (Weight==3) 792 # { 793 # Weight.matrix <- AA %*% Weight.matrix %*% AA 794 # } 795 796# if (exact==1) 797# { 798# Weight.matrix <- cbind(Weight.matrix, matrix(0,nrow=Kx,ncol=Mv)) 799# Weight.matrix <- rbind(Weight.matrix, cbind(matrix(0,nrow=Mv,ncol=Kx), 800# 1000*solve(diag(as.vector(Sig.V*Sig.V), nrow=length(Sig.V))))) 801# Weight.matrix <- as.matrix(Weight.matrix) 802# X <- cbind(X,V) 803# Mu.X <- rbind(Mu.X, matrix(0, nrow(Mu.V), 1)) 804# Sig.X <- rbind(Sig.X, matrix(1, nrow(Sig.V), 1)) 805# } #end of exact 806 807 Nx <- nrow(X) 808 Kx <- ncol(X) 809 810 if ( min(eigen(Weight.matrix, only.values=TRUE)$values) < ccc ) 811 Weight.matrix <- Weight.matrix + diag(Kx)*ccc 812 813 # I.fg. initialize matrices before looping through sample 814 YCAUS <- as.matrix(rep(0, N)) 815 SCAUS <- as.matrix(rep(0, N)) 816 XCAUS <- matrix(0, nrow=N, ncol=Kx) 817 ZCAUS <- matrix(0, nrow=N, ncol=Kz) 818 Kcount <- as.matrix(rep(0, N)) 819 KKcount <- as.matrix(rep(0, N)) 820 MMi <- as.matrix(rep(0, N)) 821 822# II. Loop through all observations that need to be matched. 823 INN <- as.matrix(1:N) 824 ww <- chol(Weight.matrix) # so that ww*ww=w.m 825# TT <- as.matrix(1:N) 826 827 # initialize some data objects 828 DD <- NULL 829 I <- NULL 830 IM <- NULL 831 IT <- NULL 832 IX <- NULL 833 IZ <- NULL 834 IY <- NULL 835 W <- NULL 836 ADist <- NULL 837 WWi <- NULL 838 Xt <- NULL 839 Zt <- NULL 840 Yt <- NULL 841 Xc <- NULL 842 Zc <- NULL 843 Yc <- NULL 844 845 for (i in 1:N) 846 { 847 #treatment indicator for observation to be matched 848 TREATi <- Tr[i] 849 850 # proceed with all observations if All==1 851 # but only with treated observations if All=0 852 if ( (TREATi==1 & All!=1) | All==1 ) 853 { 854 # covariate value for observation to be matched 855 xx <- t(as.matrix(X[i,])) 856 857 # covariate value for observation to be matched 858 zz <- t(as.matrix(Z[i,])) 859 860 # outcome value for observation to be matched 861 yy <- Y[i] 862 863 #JSS: check * 864 foo <- as.matrix(rep(1, Nx)) 865 DX <- (X - foo %*% xx) %*% t(ww) 866 867 if (Kx>1) 868 { 869 #JSS 870 foo <- t(DX*DX) 871 Dist <- as.matrix(apply(foo, 2, sum)) 872 } else { 873 Dist <- as.matrix(DX*DX) 874 } #end of Kx 875 876 # Dist distance to observation to be matched 877 # is N by 1 vector 878 879 #use the restriction matrix 880 if (restrict.trigger) 881 { 882 for (j in 1:nrow(restrict)) 883 { 884 if (restrict[j,1]==i) 885 { 886 if (restrict[j,3] < 0) { 887 Dist[restrict[j,2]] = .Machine$double.xmax 888 } else { 889 Dist[restrict[j,2]] = restrict[j,3] 890 } 891 } else if (restrict[j,2]==i) { 892 if (restrict[j,3] < 0) { 893 Dist[restrict[j,1]] = .Machine$double.xmax 894 } else { 895 Dist[restrict[j,1]] = restrict[j,3] 896 } 897 } 898 } 899 } #end if restrict.trigger 900 901 # set of potential matches (all observations with other treatment) 902 # JSS, note:logical vector 903 POTMAT <- Tr == (1-TREATi) 904 905 # X's for potential matches 906# XPOT <- X[POTMAT,] 907 DistPot <- Dist[POTMAT,1] 908 909# TTPotMat <- TT[POTMAT,1] 910 weightPot <- as.matrix(weight[POTMAT,1]) 911 912 # sorted distance of potential matches 913 S <- sort(DistPot) 914 L <- order(DistPot) 915 weightPot.sort <- weightPot[L,1] 916 weightPot.sum <- cumsum(weightPot.sort) 917 tt <- 1:(length(weightPot.sum)) 918 MMM <- min(tt[weightPot.sum>=M]) 919 920 # distance at last match 921 Distmax <- S[MMM] 922 923 # selection of actual matches 924 #logical index 925 if (restrict.trigger) 926 { 927 ACTMAT <- POTMAT & ( (Dist <= (Distmax+cdd)) & (Dist < .Machine$double.xmax) ) 928 929 if (sum(ACTMAT) < 1) 930 next; 931 } else { 932 ACTMAT <- POTMAT & ( Dist <= (Distmax+cdd) ) 933 } 934 935 Ii <- i * matrix(1, nrow=sum(ACTMAT), ncol=1) 936 IMi <- as.matrix(INN[ACTMAT,1]) 937 938 if(!is.null(ecaliper)) 939 { 940 for (j in 1:length(Ii)) 941 { 942 for( x in 1:Kx) 943 { 944 diff <- abs(X.orig[i,x] - X.orig[IMi[j], x]) 945 if (diff > ecaliper[x]) 946 { 947# print(diff) 948 949 ACTMAT[IMi[j]] <- FALSE 950 sum.caliper.drops <- sum.caliper.drops+1 951 952 break 953 } 954 } #x loop 955 } #j loop 956 957 if (sum(ACTMAT) < 1) 958 next; 959 } #ecaliper check 960 961 # distance to actual matches 962 ACTDIST <- as.matrix(Dist[ACTMAT,1]) 963 964 # counts how many times each observation is matched. 965 Kcount <- Kcount + weight[i] * weight*ACTMAT/sum(ACTMAT*weight) 966 967 KKcount <- KKcount+weight[i,1] * weight*weight*ACTMAT / 968 (sum(ACTMAT*weight)*sum(ACTMAT*weight)) 969 970 # counts how many times each observation is matched. 971 # counts number of matches for observation i 972 # Unless there are ties this should equal M 973 974 Mi <- sum(weight*ACTMAT) 975 MMi[i,1] <- Mi 976 Wi <- as.matrix(weight[ACTMAT,1]) 977 978 # mean of Y's for actual matches 979 ymat <- t(Y[ACTMAT,1]) %*% Wi/Mi 980 981 # mean of X's for actual matches 982 # mean of Z's for actual matches 983 if (length(Wi)>1) 984 { 985 xmat <- t(t(X[ACTMAT,]) %*% Wi/Mi) 986 zmat <- t(t(Z[ACTMAT,]) %*% Wi/Mi) 987 } else { 988 xmat <- t(X[ACTMAT,]) * as.double(Wi)/Mi 989 zmat <- t(Z[ACTMAT,]) * as.double(Wi)/Mi 990 } 991 992 # estimate causal effect on y for observation i 993 YCAUS[i,1] <- TREATi %*% (yy-ymat)+(1-TREATi) %*% (ymat-yy) 994 995 # difference between x and actual matches for observation i 996 XCAUS[i,] <- TREATi %*% (xx-xmat)+(1-TREATi) %*% (xmat-xx) 997 998 ZCAUS[i,] <- TREATi %*% (zz-zmat)+(1-TREATi) %*% (zmat-zz) 999 1000 # collect results 1001 I <- rbind(I, i * matrix(1, nrow=sum(ACTMAT), ncol=1)) 1002 DD <- rbind(DD, ACTDIST) 1003 IM <- rbind(IM, as.matrix(INN[ACTMAT,1])) 1004 IT <- rbind(IT, TREATi * as.matrix(rep(1, sum(ACTMAT)))) 1005 IX <- rbind(IX, as.matrix(rep(1, sum(ACTMAT))) %*% xx) 1006 IZ <- rbind(IZ, as.matrix(rep(1, sum(ACTMAT))) %*% zz) 1007 IY <- rbind(IY, as.matrix(rep(1, sum(ACTMAT))) * yy) 1008 1009 # weight for matches 1010 W <- as.matrix(c(W, weight[i,1] * Wi/Mi)) 1011 ADist <- as.matrix(c(ADist, ACTDIST)) 1012 WWi <- as.matrix(c(WWi, Wi)) 1013 1014 if (TREATi==1) 1015 { 1016 1017 if (ncol(X) > 1) 1018 { 1019 # covariates for treated 1020 Xt <- rbind(Xt, as.matrix(rep(1, sum(ACTMAT))) %*% xx) 1021 # covariate for matches 1022 Xc <- rbind(Xc, X[ACTMAT,]) 1023 1024 } else { 1025 # covariates for treated 1026 Xt <- as.matrix(c(Xt, as.matrix(rep(1, sum(ACTMAT))) %*% xx)) 1027 # covariate for matches 1028 Xc <- as.matrix(c(Xc, X[ACTMAT,])) 1029 } 1030 1031 if (ncol(Z) > 1) 1032 { 1033 # covariates for treated 1034 Zt <- rbind(Zt, as.matrix(rep(1, sum(ACTMAT))) %*% zz) 1035 # covariate for matches 1036 Zc <- rbind(Zc, Z[ACTMAT,]) 1037 } else { 1038 # covariates for treated 1039 Zt <- as.matrix(c(Zt, as.matrix(rep(1, sum(ACTMAT))) %*% zz)) 1040 # covariate for matches 1041 Zc <- as.matrix(c(Zc, Z[ACTMAT,])) 1042 } 1043 1044 # outcome for treated 1045 Yt <- as.matrix(c(Yt, yy * as.matrix(rep(1, sum(ACTMAT))) )) 1046 # outcome for matches 1047 Yc <- as.matrix(c(Yc, Y[ACTMAT,1])) 1048 } else { 1049 1050 if (ncol(X) > 1) 1051 { 1052 # covariates for controls 1053 Xc <- rbind(Xc, as.matrix(rep(1, sum(ACTMAT))) %*% xx) 1054 # covariate for matches 1055 Xt <- rbind(Xt, X[ACTMAT,]) 1056 } else { 1057 # covariates for controls 1058 Xc <- as.matrix(c(Xc, as.matrix(rep(1, sum(ACTMAT))) %*% xx)) 1059 # covariate for matches 1060 Xt <- as.matrix(c(Xt, X[ACTMAT,])) 1061 } 1062 1063 if (ncol(Z) > 1) 1064 { 1065 # covariates for controls 1066 Zc <- rbind(Zc, as.matrix(rep(1, sum(ACTMAT))) %*% zz) 1067 # covariate for matches 1068 Zt <- rbind(Zt, Z[ACTMAT,]) 1069 } else { 1070 # covariates for controls 1071 Zc <- as.matrix(c(Zc, as.matrix(rep(1, sum(ACTMAT))) %*% zz)) 1072 # covariate for matches 1073 Zt <- as.matrix(c(Zt, Z[ACTMAT,])) 1074 } 1075 # outcome for controls 1076 Yc <- as.matrix(c(Yc, as.matrix(rep(1, sum(ACTMAT))) * yy)) 1077 # outcome for matches 1078 Yt <- as.matrix(c(Yt, Y[ACTMAT,1])) 1079 1080 } #end of TREATi 1081 } #end of if 1082 } #i loop 1083 1084 # transform matched covariates back for artificial data set 1085 Xt.u <- Xt 1086 Xc.u <- Xc 1087 IX.u <- IX 1088 for (k in 1:Kx) 1089 { 1090 Xt.u[,k] <- Mu.X[k,1]+Sig.X[k,1] * Xt.u[,k] 1091 Xc.u[,k] <- Mu.X[k,1]+Sig.X[k,1] * Xc.u[,k] 1092 IX.u[,k] <- Mu.X[k,1]+Sig.X[k,1] * IX.u[,k] 1093 } 1094 if (All!=1) 1095 { 1096 I <- as.matrix(I[IT==1,]) 1097 IM <- as.matrix(IM[IT==1,]) 1098 IT <- as.matrix(IT[IT==1,]) 1099 IY <- as.matrix(IY[IT==1,]) 1100 Yc <- as.matrix(Yc[IT==1,]) 1101 Yt <- as.matrix(Yt[IT==1,]) 1102 W <- as.matrix(W[IT==1,]) 1103 ADist <- as.matrix(ADist[IT==1,]) 1104 WWi <- as.matrix(WWi[IT==1,]) 1105 IX.u <- as.matrix(IX.u[IT==1,]) 1106 Xc.u <- as.matrix(Xc.u[IT==1,]) 1107 Xt.u <- as.matrix(Xt.u[IT==1,]) 1108 Xc <- as.matrix(Xc[IT==1,]) 1109 Xt <- as.matrix(Xt[IT==1,]) 1110 Zc <- as.matrix(Zc[IT==1,]) 1111 Zt <- as.matrix(Zt[IT==1,]) 1112 IZ <- as.matrix(IZ[IT==1,]) 1113 } #end of if 1114 1115 if (length(I) < 1) 1116 { 1117 return(list(sum.caliper.drops=sum.caliper.drops,valid=0)) 1118 } else if(length(I) < 2) 1119 { 1120 return(list(sum.caliper.drops=sum.caliper.drops,valid=1)) 1121 } 1122 1123 if (BiasAdj==1) 1124 { 1125 # III. Regression of outcome on covariates for matches 1126 if (All==1) 1127 { 1128 # number of observations 1129 NNt <- nrow(Z) 1130 # add intercept 1131 ZZt <- cbind(rep(1, NNt), Z) 1132 # number of covariates 1133 Nx <- nrow(ZZt) 1134 Kx <- ncol(ZZt) 1135 xw <- ZZt*(sqrt(Tr*Kcount) %*% t(as.matrix(rep(1,Kx)))) 1136 1137 foo <- min(eigen(t(xw)%*%xw, only.values=TRUE)$values) 1138 foo <- as.double(foo<=ccc) 1139 foo2 <- apply(xw, 2, sd) 1140 1141 options(show.error.messages = FALSE) 1142 wout <- NULL 1143 try(wout <- solve( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 1144 (t(xw) %*% (Y*sqrt(Tr*Kcount)))) 1145 if(is.null(wout)) 1146 { 1147 wout2 <- NULL 1148 try(wout2 <- ginv( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 1149 (t(xw) %*% (Y*sqrt(Tr*Kcount)))) 1150 if(!is.null(wout2)) 1151 { 1152 wout <-wout2 1153 warning("using generalized inverse to calculate Bias Adjustment probably because of singular 'Z'") 1154 } 1155 } 1156 options(show.error.messages = TRUE) 1157 if(is.null(wout)) 1158 { 1159 warning("unable to calculate Bias Adjustment probably because of singular 'Z'") 1160 BiasAdj <- 0 1161 } else { 1162 NW <- nrow(wout) 1163# KW <- ncol(wout) 1164 Alphat <- wout[2:NW,1] 1165 } 1166 } else { 1167 Alphat <- matrix(0, nrow=Kz, ncol=1) 1168 } #end if ALL 1169 } 1170 1171 if(BiasAdj==1) 1172 { 1173 # III.b. Controls 1174 NNc <- nrow(Z) 1175 ZZc <- cbind(matrix(1, nrow=NNc, ncol=1),Z) 1176 Nx <- nrow(ZZc) 1177 Kx <- ncol(ZZc) 1178 1179 xw <- ZZc*(sqrt((1-Tr)*Kcount) %*% matrix(1, nrow=1, ncol=Kx)) 1180 1181 foo <- min(eigen(t(xw)%*%xw, only.values=TRUE)$values) 1182 foo <- as.double(foo<=ccc) 1183 foo2 <- apply(xw, 2, sd) 1184 1185 options(show.error.messages = FALSE) 1186 wout <- NULL 1187 try(wout <- solve( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 1188 (t(xw) %*% (Y*sqrt((1-Tr)*Kcount)))) 1189 if(is.null(wout)) 1190 { 1191 wout2 <- NULL 1192 try(wout2 <- ginv( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 1193 (t(xw) %*% (Y*sqrt((1-Tr)*Kcount)))) 1194 if(!is.null(wout2)) 1195 { 1196 wout <-wout2 1197 warning("using generalized inverse to calculate Bias Adjustment probably because of singular 'Z'") 1198 } 1199 } 1200 options(show.error.messages = TRUE) 1201 if(is.null(wout)) 1202 { 1203 warning("unable to calculate Bias Adjustment probably because of singular 'Z'") 1204 BiasAdj <- 0 1205 } else { 1206 NW <- nrow(wout) 1207# KW <- ncol(wout) 1208 Alphac <- as.matrix(wout[2:NW,1]) 1209 } 1210 } 1211 1212 1213 if(BiasAdj==1) 1214 { 1215 # III.c. adjust matched outcomes using regression adjustment for bias adjusted matching estimator 1216 1217 SCAUS <- YCAUS-Tr*(ZCAUS %*% Alphac)-(1-Tr)*(ZCAUS %*% Alphat) 1218 # adjusted control outcome 1219 Yc.adj <- Yc+BiasAdj * (IZ-Zc) %*% Alphac 1220 # adjusted treated outcome 1221 Yt.adj <- Yt+BiasAdj*(IZ-Zt) %*% Alphat 1222 Tau.i <- Yt.adj - Yc.adj 1223 } else { 1224 Yc.adj <- Yc 1225 Yt.adj <- Yt 1226 Yt.adj <- Yt 1227 Tau.i <- Yt.adj - Yc.adj 1228 } 1229 1230 art.data <- cbind(I,IM,IT,DD,IY,Yc,Yt,W,WWi,ADist,IX.u,Xc.u,Xt.u, 1231 Yc.adj,Yt.adj,Tau.i) 1232 1233 1234 # III. If conditional variance is needed, initialize variance vector 1235 # and loop through all observations 1236 1237 Nx <- nrow(X) 1238 Kx <- ncol(X) 1239 1240# ww <- chol(Weight.matrix) 1241# NN <- as.matrix(1:N) 1242 if (Var.calc>0) 1243 { 1244 Sig <- matrix(0, nrow=N, ncol=1) 1245 # overall standard deviation of outcome 1246 # std <- sd(Y) 1247 for (i in 1:N) 1248 { 1249 # treatment indicator observation to be matched 1250 TREATi <- Tr[i,1] 1251 # covariate value for observation to be matched 1252 xx <- X[i,] 1253 # potential matches are all observations with the same treatment value 1254 POTMAT <- (Tr==TREATi) 1255 POTMAT[i,1] <- 0 1256 weightPOT <- as.matrix(weight[POTMAT==1,1]) 1257 DX <- (X - matrix(1, Nx,1) %*% xx) %*% t(ww) 1258 if (Kx>1) 1259 { 1260 foo <- apply(t(DX*DX), 2, sum) 1261 Dist <- as.matrix(foo) 1262 } else { 1263 Dist <- DX*DX 1264 } 1265 1266 # distance to observation to be matched 1267 1268 # Distance vector only for potential matches 1269 DistPot <- Dist[POTMAT==1,1] 1270 # sorted distance of potential matches 1271 S <- sort(DistPot) 1272 L <- order(DistPot) 1273 weightPOT.sort <- weightPOT[L,1] 1274 weightPOT.sum <- cumsum(weightPOT.sort) 1275 tt <- 1:(length(weightPOT.sum)) 1276 MMM <- min(tt[weightPOT.sum >= Var.calc]) 1277 MMM <- min(MMM,length(S)) 1278 Distmax=S[MMM] 1279 1280 # distance of Var_calc-th closest match 1281 ACTMAT <- (POTMAT==1) & (Dist<= (Distmax+ccc)) 1282 1283 # indicator for actual matches, that is all potential 1284 # matches closer than, or as close as the Var_calc-th 1285 # closest 1286 1287 Yactmat <- as.matrix(c(Y[i,1], Y[ACTMAT,1])) 1288 weightactmat <- as.matrix(c(weight[i,1], weight[ACTMAT,1])) 1289 fm <- t(Yactmat) %*% weightactmat/sum(weightactmat) 1290 sm <- sum(Yactmat*Yactmat*weightactmat)/sum(weightactmat) 1291 sigsig <- (sm-fm %*% fm)*sum(weightactmat)/(sum(weightactmat)-1) 1292 1293 # standard deviation of actual matches 1294 Sig[i,1] <- sqrt(sigsig) 1295 }# end of i loop 1296 #variance estimate 1297 Sigs <- Sig*Sig 1298 } #end of var.calc > 0 1299 1300 # matching estimator 1301 est <- t(W) %*% Tau.i/sum(W) 1302# est.t <- sum((iot.t*Tr+iot.c*Kcount*Tr)*Y)/sum(iot.t*Tr+iot.c*Kcount*Tr) 1303# est.c <- sum((iot.t*(1-Tr)+iot.c*Kcount*(1-Tr))*Y)/sum(iot.t*(1-Tr)+iot.c*Kcount*(1-Tr)) 1304 1305 if (Var.calc==0) 1306 { 1307 eps <- Tau.i - as.double(est) 1308 eps.sq <- eps*eps 1309 Sigs <- 0.5 * matrix(1, N, 1) %*% (t(eps.sq) %*% W)/sum(W) 1310# sss <- sqrt(Sigs[1,1]) 1311 } #end of Var.calc==0 1312 1313 SN <- sum(iot.t) 1314 var.sample <- sum((Sigs*(iot.t+iot.c*Kcount)*(iot.t+iot.c*Kcount))/(SN*SN)) 1315 1316 if (All==1) 1317 { 1318 var.pop <- sum((Sigs*(iot.c*Kcount*Kcount+2*iot.c*Kcount-iot.c*KKcount))/(SN*SN)) 1319 } else { 1320 var.pop=sum((Sigs*(iot.c*Kcount*Kcount-iot.c*KKcount))/(SN*SN)) 1321 } 1322 1323 if (BiasAdj==1) 1324 { 1325 dvar.pop <- sum(iot.t*(SCAUS-as.double(est))*(SCAUS-as.double(est)))/(SN*SN) 1326 } else { 1327 dvar.pop <- sum(iot.t*(YCAUS-as.double(est))*(YCAUS-as.double(est)))/(SN*SN) 1328 } 1329 1330 var.pop <- var.pop + dvar.pop 1331 1332 if (SAMPLE==1) 1333 { 1334 var <- var.sample 1335 } else { 1336 var <- max(var.sample, var.pop) 1337 var <- var.pop 1338 } 1339 1340 var.cond <- max(var.sample,var.pop)-var.sample 1341 se <- sqrt(var) 1342 se.cond <- sqrt(var.cond) 1343 1344# Sig <- sqrt(Sigs) 1345# aug.data <- cbind(Y,Tr,X,Z,Kcount,Sig,weight) 1346 1347 if (All==2) 1348 est <- -1*est 1349 1350# if (exact==1) 1351# { 1352# Vt.u <- Xt.u[,(Kx-Mv+1):Kx] 1353# Vc.u <- Xc.u[,(Kx-Mv+1):Kx] 1354# Vdif <- abs(Vt.u-Vc.u) 1355# 1356# if (Mv>1) 1357# Vdif <- as.matrix(apply(t(Vdif), 2, sum)) 1358# 1359# em[1,1] <- length(Vdif) 1360# em[2,1] <- sum(Vdif>0.000001) 1361# }#end of exact==1 1362 1363 return(list(est=est, se=se, se.cond=se.cond, W=W, 1364 sum.caliper.drops=sum.caliper.drops, 1365 art.data=art.data)) 1366 }# end of Rmatch 1367 1368 1369# 1370# See Rosenbaum and Rubin (1985) and Smith and Todd Rejoinder (JofEconometrics) p.9 1371# 1372sdiff.pooled <- function(Tr, Co, weights=rep(1,length(Co)), 1373 weights.Tr=rep(1,length(Tr)), 1374 weights.Co=rep(1,length(Co)), 1375 match=FALSE) 1376 { 1377 if (!match) 1378 { 1379 obs.Tr <- sum(weights.Tr) 1380 obs.Co <- sum(weights.Co) 1381# obs.total <- obs.Tr+obs.Co 1382 1383 mean.Tr <- sum(Tr*weights.Tr)/obs.Tr 1384 mean.Co <- sum(Co*weights.Co)/obs.Co 1385 diff <- mean.Tr - mean.Co 1386 1387#match Rubin 1388# mean.total <- sum(Tr*weights.Tr)/obs.total + sum(Co*weights.Co)/obs.total 1389# var.total <- sum( ( (Tr - mean.total)^2 )*weights.Tr)/(obs.total-1) + 1390# sum( ( (Co - mean.total)^2 )*weights.Co)/(obs.total-1) 1391 var.pooled <- ( sum( ( (Tr - mean.Tr)^2)*weights.Tr)/(obs.Tr-1) + 1392 sum( ( (Co - mean.Co)^2 )*weights.Co)/(obs.Co-1) )/2 1393 1394 1395 if(var.pooled==0 & diff==0) 1396 { 1397 sdiff <- 0 1398 } else { 1399 sdiff <- 100*diff/sqrt( var.pooled ) 1400 } 1401 1402 } else{ 1403 diff <- sum( (Tr-Co)*weights ) /sum(weights) 1404 mean.Tr <- sum(Tr*weights)/sum(weights) 1405 mean.Co <- sum(Co*weights)/sum(weights) 1406 1407#match Rubin 1408 obs <- sum(weights) 1409 1410# obs for total 1411# obs = sum(weights)*2 1412# mean.total <- (mean.Tr + mean.Co)/2 1413# var.total <- sum( ( (Tr - mean.total)^2 )*weights)/(obs-1) + 1414# sum( ( (Co - mean.total)^2 )*weights)/(obs-1) 1415 1416 var.pooled <- ( sum( ( (Tr - mean.Tr)^2 )*weights)/(obs-1) + 1417 sum( ( (Co - mean.Co)^2 )*weights)/(obs-1) )/2 1418 1419 if(var.pooled==0 & diff==0) 1420 { 1421 sdiff <- 0 1422 } else { 1423 sdiff <- 100*diff/sqrt(var.pooled) 1424 } 1425 } 1426 1427 return(sdiff) 1428 } 1429 1430 1431# 1432# STANDARDIZED BY THE VARIANCE OF THE TREATMENT GROUP 1433# See sdiff.rubin for Rosenbaum and Rubin (1985) and Smith and Todd Rejoinder (JofEconometrics) p.9 1434# 1435sdiff <- function(Tr, Co, weights=rep(1,length(Co)), 1436 weights.Tr=rep(1,length(Tr)), 1437 weights.Co=rep(1,length(Co)), 1438 match=FALSE, 1439 estimand="ATT") 1440 1441 { 1442 if (!match) 1443 { 1444 obs.Tr <- sum(weights.Tr) 1445 obs.Co <- sum(weights.Co) 1446 1447 mean.Tr <- sum(Tr*weights.Tr)/obs.Tr 1448 mean.Co <- sum(Co*weights.Co)/obs.Co 1449 diff <- mean.Tr - mean.Co 1450 1451 if(estimand=="ATC") 1452 { 1453 var <- sum( ( (Co - mean.Co)^2 )*weights.Co)/(obs.Co-1) 1454 } else { 1455 var <- sum( ( (Tr - mean.Tr)^2 )*weights.Tr)/(obs.Tr-1) 1456 } 1457 1458 if(var==0 & diff==0) 1459 { 1460 sdiff=0 1461 } else { 1462 sdiff <- 100*diff/sqrt( (var) ) 1463 } 1464 } else{ 1465 mean.Tr <- sum(Tr*weights)/sum(weights) 1466 mean.Co <- sum(Co*weights)/sum(weights) 1467 diff <- mean.Tr - mean.Co 1468 1469 if(estimand=="ATC") 1470 { 1471 var <- sum( ( (Co - mean.Co)^2 )*weights)/(sum(weights)-1) 1472 } else { 1473 var <- sum( ( (Tr - mean.Tr)^2 )*weights)/(sum(weights)-1) 1474 } 1475 1476 if(var==0 & diff==0) 1477 { 1478 sdiff <- 0 1479 } else { 1480 sdiff <- 100*diff/sqrt( (var) ) 1481 } 1482 } 1483 1484 return(sdiff) 1485 } 1486 1487# function runs sdiff and t.test 1488# 1489balanceUV <- function(Tr, Co, weights=rep(1,length(Co)), exact=FALSE, ks=FALSE, nboots=1000, 1490 paired=TRUE, match=FALSE, 1491 weights.Tr=rep(1,length(Tr)), weights.Co=rep(1,length(Co)), 1492 estimand="ATT") 1493 { 1494 ks.out <- NULL 1495 1496 if (!match) 1497 { 1498 sbalance.pooled <- sdiff.pooled(Tr=Tr, Co=Co, 1499 weights.Tr=weights.Tr, 1500 weights.Co=weights.Co, 1501 match=FALSE) 1502 1503 sbalance.constvar <- sdiff(Tr=Tr, Co=Co, 1504 weights.Tr=weights.Tr, 1505 weights.Co=weights.Co, 1506 match=FALSE, 1507 estimand=estimand) 1508 1509 obs.Tr <- sum(weights.Tr) 1510 obs.Co <- sum(weights.Co) 1511 1512 mean.Tr <- sum(Tr*weights.Tr)/obs.Tr 1513 mean.Co <- sum(Co*weights.Co)/obs.Co 1514 var.Tr <- sum( ( (Tr - mean.Tr)^2 )*weights.Tr)/(obs.Tr-1) 1515 var.Co <- sum( ( (Co - mean.Co)^2 )*weights.Co)/(obs.Co-1) 1516 var.ratio <- var.Tr/var.Co 1517 1518 qqsummary <- qqstats(x=Tr, y=Co, standardize=TRUE) 1519 qqsummary.raw <- qqstats(x=Tr, y=Co, standardize=FALSE) 1520 1521 tt <- Mt.test.unpaired(Tr, Co, 1522 weights.Tr=weights.Tr, 1523 weights.Co=weights.Co) 1524 1525 if (ks) 1526 ks.out <- ks.boot(Tr, Co,nboots=nboots) 1527 } else { 1528 sbalance.pooled <- sdiff(Tr=Tr, Co=Co, weights=weights, match=TRUE) 1529 sbalance.constvar <- sdiff(Tr=Tr, Co=Co, weights=weights, match=TRUE, estimand=estimand) 1530 1531 mean.Tr <- sum(Tr*weights)/sum(weights); 1532 mean.Co <- sum(Co*weights)/sum(weights); 1533 var.Tr <- sum( ( (Tr - mean.Tr)^2 )*weights)/sum(weights); 1534 var.Co <- sum( ( (Co - mean.Co)^2 )*weights)/sum(weights); 1535 var.ratio <- var.Tr/var.Co 1536 1537 qqsummary <- qqstats(x=Tr, y=Co, standardize=TRUE) 1538 qqsummary.raw <- qqstats(x=Tr, y=Co, standardize=FALSE) 1539 1540 if(paired) 1541 { 1542 tt <- Mt.test(Tr, Co, weights) 1543 } else { 1544 tt <- Mt.test.unpaired(Tr, Co, weights.Tr=weights, weights.Co=weights) 1545 } 1546 1547 if (ks) 1548 ks.out <- ks.boot(Tr, Co, nboots=nboots) 1549 } 1550 1551 ret <- list(sdiff=sbalance.constvar, 1552 sdiff.pooled=sbalance.pooled, 1553 mean.Tr=mean.Tr,mean.Co=mean.Co, 1554 var.Tr=var.Tr,var.Co=var.Co, p.value=tt$p.value, 1555 var.ratio=var.ratio, 1556 ks=ks.out, tt=tt, 1557 qqsummary=qqsummary, 1558 qqsummary.raw=qqsummary.raw) 1559 1560 class(ret) <- "balanceUV" 1561 return(ret) 1562 } #balanceUV 1563 1564summary.balanceUV <- function(object, ..., digits=5) 1565 { 1566 1567 if (!inherits(object, "balanceUV")) { 1568 warning("Object not of class 'balanceUV'") 1569 return(NULL) 1570 } 1571 1572 cat("mean treatment........", format(object$mean.Tr, digits=digits),"\n") 1573 cat("mean control..........", format(object$mean.Co, digits=digits),"\n") 1574 cat("std mean diff.........", format(object$sdiff, digits=digits),"\n\n") 1575 1576 cat("mean raw eQQ diff.....", format(object$qqsummary.raw$meandiff, digits=digits),"\n") 1577 cat("med raw eQQ diff.....", format(object$qqsummary.raw$mediandiff, digits=digits),"\n") 1578 cat("max raw eQQ diff.....", format(object$qqsummary.raw$maxdiff, digits=digits),"\n\n") 1579 1580 cat("mean eCDF diff........", format(object$qqsummary$meandiff, digits=digits),"\n") 1581 cat("med eCDF diff........", format(object$qqsummary$mediandiff, digits=digits),"\n") 1582 cat("max eCDF diff........", format(object$qqsummary$maxdiff, digits=digits),"\n\n") 1583 1584 cat("var ratio (Tr/Co).....", format(object$var.ratio, digits=digits),"\n") 1585 cat("T-test p-value........", format.pval(object$tt$p.value,digits=digits), "\n") 1586 if (!is.null(object$ks)) 1587 { 1588 if(!is.na(object$ks$ks.boot.pvalue)) 1589 { 1590 cat("KS Bootstrap p-value..", format.pval(object$ks$ks.boot.pvalue, digits=digits), "\n") 1591 } 1592 cat("KS Naive p-value......", format(object$ks$ks$p.value, digits=digits), "\n") 1593 cat("KS Statistic..........", format(object$ks$ks$statistic, digits=digits), "\n") 1594 } 1595 cat("\n") 1596 } #end of summary.balanceUV 1597 1598 1599#print Before and After balanceUV together 1600PrintBalanceUV <- function(BeforeBalance, AfterBalance, ..., digits=5) 1601 { 1602 1603 if (!inherits(BeforeBalance, "balanceUV")) { 1604 warning("BeforeBalance not of class 'balanceUV'") 1605 return(NULL) 1606 } 1607 1608 if (!inherits(AfterBalance, "balanceUV")) { 1609 warning("AfterBalance not of class 'balanceUV'") 1610 return(NULL) 1611 } 1612 1613 space.size <- digits*2 1614# space <- rep(" ",space.size) 1615 1616 cat(" ", "Before Matching", "\t \t After Matching\n") 1617 cat("mean treatment........", format(BeforeBalance$mean.Tr, digits=digits, width=space.size), "\t \t", 1618 format(AfterBalance$mean.Tr, digits=digits, width=space.size), 1619 "\n") 1620 cat("mean control..........", format(BeforeBalance$mean.Co, digits=digits, width=space.size), "\t \t", 1621 format(AfterBalance$mean.Co, digits=digits, width=space.size), 1622 "\n") 1623 cat("std mean diff.........", format(BeforeBalance$sdiff, digits=digits, width=space.size), "\t \t", 1624 format(AfterBalance$sdiff, digits=digits, width=space.size), 1625 "\n\n") 1626 1627 cat("mean raw eQQ diff.....", format(BeforeBalance$qqsummary.raw$meandiff, digits=digits, width=space.size), "\t \t", 1628 format(AfterBalance$qqsummary.raw$meandiff, digits=digits, width=space.size), 1629 "\n") 1630 cat("med raw eQQ diff.....", format(BeforeBalance$qqsummary.raw$mediandiff, digits=digits, width=space.size), "\t \t", 1631 format(AfterBalance$qqsummary.raw$mediandiff, digits=digits, width=space.size), 1632 "\n") 1633 cat("max raw eQQ diff.....", format(BeforeBalance$qqsummary.raw$maxdiff, digits=digits, width=space.size), "\t \t", 1634 format(AfterBalance$qqsummary.raw$maxdiff, digits=digits, width=space.size), 1635 "\n\n") 1636 1637 cat("mean eCDF diff........", format(BeforeBalance$qqsummary$meandiff, digits=digits, width=space.size), "\t \t", 1638 format(AfterBalance$qqsummary$meandiff, digits=digits, width=space.size), 1639 "\n") 1640 cat("med eCDF diff........", format(BeforeBalance$qqsummary$mediandiff, digits=digits, width=space.size), "\t \t", 1641 format(AfterBalance$qqsummary$mediandiff, digits=digits, width=space.size), 1642 "\n") 1643 cat("max eCDF diff........", format(BeforeBalance$qqsummary$maxdiff, digits=digits, width=space.size), "\t \t", 1644 format(AfterBalance$qqsummary$maxdiff, digits=digits, width=space.size), 1645 "\n\n") 1646 1647 cat("var ratio (Tr/Co).....", format(BeforeBalance$var.ratio, digits=digits, width=space.size), "\t \t", 1648 format(AfterBalance$var.ratio, digits=digits, width=space.size), 1649 "\n") 1650 cat("T-test p-value........", format(format.pval(BeforeBalance$tt$p.value,digits=digits), justify="right", width=space.size), "\t \t", 1651 format(format.pval(AfterBalance$tt$p.value,digits=digits), justify="right", width=space.size), 1652 "\n") 1653 if (!is.null(BeforeBalance$ks)) 1654 { 1655 if(!is.na(BeforeBalance$ks$ks.boot.pvalue)) 1656 { 1657 cat("KS Bootstrap p-value..", format(format.pval(BeforeBalance$ks$ks.boot.pvalue, digits=digits), justify="right",width=space.size), "\t \t", 1658 format(format.pval(AfterBalance$ks$ks.boot.pvalue, digits=digits), justify="right", width=space.size), 1659 "\n") 1660 } 1661 cat("KS Naive p-value......", format(format.pval(BeforeBalance$ks$ks$p.value, digits=digits), justify="right",width=space.size), "\t \t", 1662 format(format.pval(AfterBalance$ks$ks$p.value, digits=digits), justify="right",width=space.size), 1663 "\n") 1664 cat("KS Statistic..........", format(BeforeBalance$ks$ks$statistic, digits=digits, width=space.size), "\t \t", 1665 format(AfterBalance$ks$ks$statistic, digits=digits, width=space.size), 1666 "\n") 1667 } 1668 cat("\n") 1669 } #end of PrintBalanceUV 1670 1671 1672#removed as of 0.99-7 (codetools) 1673#McNemar <- function(Tr, Co, weights=rep(1,length(Tr))) 1674 1675McNemar2 <- function (Tr, Co, correct = TRUE, weights=rep(1,length(Tr))) 1676{ 1677 x <- Tr 1678 y <- Co 1679 if (is.matrix(x)) { 1680 stop("this version of McNemar cannot handle x being a matrix") 1681 } else { 1682 if (is.null(y)) 1683 stop("if x is not a matrix, y must be given") 1684 if (length(x) != length(y)) 1685 stop("x and y must have the same length") 1686 DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y))) 1687 OK <- complete.cases(x, y) 1688 x <- factor(x[OK]) 1689 y <- factor(y[OK]) 1690 r <- nlevels(x) 1691 if ((r < 2) || (nlevels(y) != r)) 1692 { 1693 stop("x and y must have the same number of levels (minimum 2)") 1694 } 1695 } 1696 tx <- table(x, y) 1697 facs <- levels(x) 1698 txw <- tx 1699 for(i in 1:r) 1700 { 1701 for(j in 1:r) 1702 { 1703 indx <- x==facs[i] & y==facs[j] 1704 txw[i,j] <- sum(weights[indx]); 1705 } 1706 } 1707 pdiscordant <- sum( ( (x!=y)*weights )/sum(weights) ) 1708 x <- txw 1709 1710 PARAMETER <- r * (r - 1)/2 1711 METHOD <- "McNemar's Chi-squared test" 1712 if (correct && (r == 2) && any(x - t(x))) { 1713 y <- (abs(x - t(x)) - 1) 1714 METHOD <- paste(METHOD, "with continuity correction") 1715 } else y <- x - t(x) 1716 x <- x + t(x) 1717 STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)]) 1718 PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE) 1719 names(STATISTIC) <- "McNemar's chi-squared" 1720 names(PARAMETER) <- "df" 1721 1722 RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, 1723 p.value = PVAL, method = METHOD, data.name = DNAME, 1724 pdiscordant=pdiscordant) 1725 class(RVAL) <- "htest" 1726 return(RVAL) 1727} 1728 1729ks<-function(x,y,w=F,sig=T){ 1730# Compute the Kolmogorov-Smirnov test statistic 1731# 1732# Code for the Kolmogorov-Smirnov test is adopted from the Splus code 1733# written by Rand R. Wilcox for his book titled "Introduction to 1734# Robust Estimation and Hypothesis Testing." Academic Press, 1997. 1735# 1736# Also see 1737#@book( knuth1998, 1738# author= {Knuth, Donald E.}, 1739# title= {The Art of Computer Programming, Vol. 2: Seminumerical Algorithms}, 1740# edition= "3rd", 1741# publisher= "Addison-Wesley", address= "Reading: MA", 1742# year= 1998 1743#) 1744# and Wilcox 1997 1745# 1746# w=T computes the weighted version instead. 1747# sig=T indicates that the exact significance level is to be computed. 1748# If there are ties, the reported significance level is exact when 1749# using the unweighted test, but for the weighted test the reported 1750# level is too high. 1751# 1752# This function uses the functions ecdf, kstiesig, kssig and kswsig 1753# 1754# This function returns the value of the test statistic, the approximate .05 1755# critical value, and the exact significance level if sig=T. 1756# 1757# Missing values are automatically removed 1758# 1759 1760 ecdf<-function(x,val){ 1761 # compute empirical cdf for data in x evaluated at val 1762 # That is, estimate P(X <= val) 1763 # 1764 ecdf<-length(x[x<=val])/length(x) 1765 ecdf 1766 }#end ecdf 1767 1768 x<-x[!is.na(x)] 1769 y<-y[!is.na(y)] 1770 w<-as.logical(w) 1771 sig<-as.logical(sig) 1772 tie<-logical(1) 1773 tie<-F 1774 siglevel<-NA 1775 z<-sort(c(x,y)) # Pool and sort the observations 1776 for (i in 2:length(z))if(z[i-1]==z[i])tie<-T #check for ties 1777 v<-1 # Initializes v 1778 for (i in 1:length(z))v[i]<-abs(ecdf(x,z[i])-ecdf(y,z[i])) 1779 ks<-max(v) 1780 crit<-1.36*sqrt((length(x)+length(y))/(length(x)*length(y))) # Approximate 1781# .05 critical value 1782 if(!w && sig && !tie)siglevel<-kssig(length(x),length(y),ks) 1783 if(!w && sig && tie)siglevel<-kstiesig(x,y,ks) 1784 if(w){ 1785 crit<-(max(length(x),length(y))-5)*.48/95+2.58+abs(length(x)-length(y))*.44/95 1786 if(length(x)>100 || length(y)>100){ 1787 print("When either sample size is greater than 100,") 1788 print("the approximate critical value can be inaccurate.") 1789 print("It is recommended that the exact significance level be computed.") 1790 } 1791 for (i in 1:length(z)){ 1792 temp<-(length(x)*ecdf(x,z[i])+length(y)*ecdf(y,z[i]))/length(z) 1793 temp<-temp*(1.-temp) 1794 v[i]<-v[i]/sqrt(temp) 1795 } 1796 v<-v[!is.na(v)] 1797 ks<-max(v)*sqrt(length(x)*length(y)/length(z)) 1798 if(sig)siglevel<-kswsig(length(x),length(y),ks) 1799 if(tie && sig){ 1800 print("Ties were detected. The reported significance level") 1801 print("of the weighted Kolmogorov-Smirnov test statistic is not exact.") 1802 }} 1803 #round off siglevel in a nicer way 1804 if(is.double(ks) & is.double(crit) & !is.na(ks) & !is.na(crit)) 1805 { 1806 if (is.na(siglevel) & ks < crit) 1807 { 1808 siglevel <- 0.99999837212332 1809 } 1810 if (is.double(siglevel) & !is.na(siglevel)) 1811 { 1812 if (siglevel < 0) 1813 siglevel <- 0 1814 } 1815 } 1816 list(test=ks,critval=crit,pval=siglevel) 1817} 1818 1819 1820kssig<-function(m,n,val){ 1821# 1822# Compute significance level of the Kolmogorov-Smirnov test statistic 1823# m=sample size of first group 1824# n=sample size of second group 1825# val=observed value of test statistic 1826# 1827 cmat<-matrix(0,m+1,n+1) 1828 umat<-matrix(0,m+1,n+1) 1829 for (i in 0:m){ 1830 for (j in 0:n)cmat[i+1,j+1]<-abs(i/m-j/n) 1831 } 1832 cmat<-ifelse(cmat<=val,1e0,0e0) 1833 for (i in 0:m){ 1834 for (j in 0:n)if(i*j==0)umat[i+1,j+1]<-cmat[i+1,j+1] 1835 else umat[i+1,j+1]<-cmat[i+1,j+1]*(umat[i+1,j]+umat[i,j+1]) 1836 } 1837 term<-lgamma(m+n+1)-lgamma(m+1)-lgamma(n+1) 1838 kssig<-1.-umat[m+1,n+1]/exp(term) 1839 return(kssig) 1840} 1841 1842kstiesig<-function(x,y,val){ 1843# 1844# Compute significance level of the Kolmogorov-Smirnov test statistic 1845# for the data in x and y. 1846# This function allows ties among the values. 1847# val=observed value of test statistic 1848# 1849m<-length(x) 1850n<-length(y) 1851z<-c(x,y) 1852z<-sort(z) 1853cmat<-matrix(0,m+1,n+1) 1854umat<-matrix(0,m+1,n+1) 1855for (i in 0:m){ 1856for (j in 0:n){ 1857if(abs(i/m-j/n)<=val)cmat[i+1,j+1]<-1e0 1858k<-i+j 1859if(k > 0 && k<length(z) && z[k]==z[k+1])cmat[i+1,j+1]<-1 1860} 1861} 1862for (i in 0:m){ 1863for (j in 0:n)if(i*j==0)umat[i+1,j+1]<-cmat[i+1,j+1] 1864else umat[i+1,j+1]<-cmat[i+1,j+1]*(umat[i+1,j]+umat[i,j+1]) 1865} 1866term<-lgamma(m+n+1)-lgamma(m+1)-lgamma(n+1) 1867kstiesig<-1.-umat[m+1,n+1]/exp(term) 1868kstiesig 1869} 1870 1871 1872 1873kswsig<-function(m,n,val){ 1874# 1875# Compute significance level of the weighted 1876# Kolmogorov-Smirnov test statistic 1877# 1878# m=sample size of first group 1879# n=sample size of second group 1880# val=observed value of test statistic 1881# 1882 mpn<-m+n 1883 cmat<-matrix(0,m+1,n+1) 1884 umat<-matrix(0,m+1,n+1) 1885 for (i in 1:m-1){ 1886 for (j in 1:n-1)cmat[i+1,j+1]<-abs(i/m-j/n)*sqrt(m*n/((i+j)*(1-(i+j)/mpn))) 1887 } 1888 cmat<-ifelse(cmat<=val,1,0) 1889 for (i in 0:m){ 1890 for (j in 0:n)if(i*j==0)umat[i+1,j+1]<-cmat[i+1,j+1] 1891 else umat[i+1,j+1]<-cmat[i+1,j+1]*(umat[i+1,j]+umat[i,j+1]) 1892 } 1893 term<-lgamma(m+n+1)-lgamma(m+1)-lgamma(n+1) 1894 kswsig<-1.-umat[m+1,n+1]/exp(term) 1895 kswsig 1896} 1897 1898Mks.test.handler <- function(w) 1899 { 1900 #suppress the following warning: 1901 #In ks.test() : 1902 # p-values will be approximate in the presence of ties 1903 if( any( grepl( "ties", w) ) ) 1904 invokeRestart( "muffleWarning" ) 1905 1906 #invoke as: 1907 #withCallingHandlers( ks.test(round(x1), round(x2)), warning = Mks.test.handler ) 1908 } 1909 1910Mks.test <- function(...) 1911 { 1912 RVAL <- withCallingHandlers( ks.test(...), warning = Mks.test.handler ) 1913 return(RVAL) 1914 } 1915 1916Mt.test <- function(Tr, Co, weights) 1917 { 1918 v1 <- Tr-Co 1919 estimate <- sum(v1*weights)/sum(weights) 1920 var1 <- sum( ((v1-estimate)^2)*weights )/( sum(weights)*sum(weights) ) 1921 1922 parameter <- Inf 1923 1924 #get rid of NA for t.test! 1925 if (estimate==0 & var1==0) 1926 { 1927 statistic = 0 1928 p.value = 1 1929 } else { 1930 statistic <- estimate/sqrt(var1) 1931 1932 p.value <- (1-MATCHpt(abs(statistic), df=sum(weights)-1))*2 1933 } 1934 1935 z <- list(statistic=statistic, parameter=parameter, p.value=p.value, 1936 estimate=estimate) 1937 return(z) 1938 } #end of Mt.test 1939 1940Mt.test.unpaired <- function(Tr, Co, 1941 weights.Tr=rep(1,length(Tr)), 1942 weights.Co=rep(1,length(Co))) 1943 { 1944 obs.Tr <- sum(weights.Tr) 1945 obs.Co <- sum(weights.Co) 1946 1947 mean.Tr <- sum(Tr*weights.Tr)/obs.Tr 1948 mean.Co <- sum(Co*weights.Co)/obs.Co 1949 estimate <- mean.Tr-mean.Co 1950 var.Tr <- sum( ( (Tr - mean.Tr)^2 )*weights.Tr)/(obs.Tr-1) 1951 var.Co <- sum( ( (Co - mean.Co)^2 )*weights.Co)/(obs.Co-1) 1952 dim <- sqrt(var.Tr/obs.Tr + var.Co/obs.Co) 1953 1954 parameter <- Inf 1955 1956 #get rid of NA for t.test! 1957 if (estimate==0 & dim==0) 1958 { 1959 statistic = 0 1960 p.value = 1 1961 } else { 1962 statistic <- estimate/dim 1963 1964 a1 <- var.Tr/obs.Tr 1965 a2 <- var.Co/obs.Co 1966 dof <- ((a1 + a2)^2)/( (a1^2)/(obs.Tr - 1) + (a2^2)/(obs.Co - 1) ) 1967 1968 p.value <- (1-MATCHpt(abs(statistic), df=dof))*2 1969 } 1970 1971 z <- list(statistic=statistic, parameter=parameter, p.value=p.value, 1972 estimate=estimate) 1973 return(z) 1974 } #end of Mt.test.unpaired 1975 1976MatchBalance <- function(formul, data=NULL, match.out=NULL, ks=TRUE, 1977 nboots=500, weights=NULL, 1978 digits=5, paired=TRUE, print.level=1) 1979 { 1980 if(!is.list(match.out) & !is.null(match.out)) { 1981 warning("'Match' object contains no valid matches") 1982 match.out <- NULL 1983 } 1984 1985 if ((!inherits(match.out, "Match")) & (!inherits(match.out, "Matchby")) & (!is.null(match.out)) ) { 1986 warning("Object not of class 'Match'") 1987 match.out <- NULL 1988 } 1989 1990 orig.na.action <- as.character(options("na.action")) 1991 options("na.action"=na.pass) 1992 if (is.null(data)) 1993 { 1994 xdata <- as.data.frame(get.xdata(formul,datafr=environment(formul))) 1995 Tr <- as.double(get.ydata(formul,datafr=environment(formul))) 1996 1997 } else { 1998 data <- as.data.frame(data) 1999 2000 xdata <- as.data.frame(get.xdata(formul, data)) 2001 Tr <- as.double(get.ydata(formul, data)) 2002 } 2003 options("na.action"=orig.na.action) 2004 2005 if (is.null(weights)) 2006 weights <- rep(1,length(Tr)) 2007 2008 if(!is.numeric(weights)) 2009 stop("'weights' must be a numeric vector") 2010 2011 if( sum(is.na(xdata))!=0 | sum(is.na(Tr))!=0) 2012 { 2013 2014 if(orig.na.action!="na.omit" & orig.na.action!="na.exclude" & orig.na.action!="na.fail") 2015 warning("'na.action' should be set to 'na.omit', 'na.exclude' or 'na.fail' see 'help(na.action)'") 2016 2017 if (orig.na.action=="na.fail") 2018 { 2019 stop("NA's found in data input.") 2020 } else { 2021 warning("NA's found in data input. IT IS HIGHLY RECOMMENDED THAT YOU TEST IF THE MISSING VALUES ARE BALANCED INSTEAD OF JUST DELETING THEM.") 2022 } 2023 2024 indx1 <- apply(is.na(xdata),1,sum)==0 & is.na(Tr)==0 2025 Tr <- Tr[indx1] 2026 xdata = xdata[indx1,] 2027 weights <- weights[indx1] 2028 } #end of na 2029 2030 if (sum(Tr !=1 & Tr !=0) > 0) { 2031 stop("Treatment indicator must be a logical variable---i.e., TRUE (1) or FALSE (0)") 2032 } 2033 2034 nvars <- ncol(xdata) 2035 names.xdata <- names(xdata) 2036 2037 findx <- 1 2038 if (sum(xdata[,1]==rep(1,nrow(xdata)))==nrow(xdata)) 2039 { 2040 findx <- 2 2041 } 2042 2043 if(nboots < 10 & nboots > 0) 2044 nboots <- 10 2045 2046 if (ks) 2047 { 2048 ks.bm <- KSbootBalanceSummary(index.treated=(Tr==0), 2049 index.control=(Tr==1), 2050 X=xdata[,findx:nvars], 2051 nboots=nboots) 2052 2053 if (!is.null(match.out)) 2054 { 2055 ks.am <- KSbootBalanceSummary(index.treated=match.out$index.treated, 2056 index.control=match.out$index.control, 2057 X=xdata[,findx:nvars], 2058 nboots=nboots) 2059 } 2060 } 2061 2062 BeforeMatchingBalance <- list() 2063 AfterMatchingBalance <- list() 2064 2065 BMsmallest.p.value <- 1 2066 BMsmallest.number <- 1 2067 BMsmallest.name <- names.xdata[findx] 2068 2069 AMsmallest.p.value <- NULL 2070 AMsmallest.number <- NULL 2071 AMsmallest.name <- NULL 2072 2073 if (!is.null(match.out)) 2074 { 2075 AMsmallest.p.value <- 1 2076 AMsmallest.number <- 1 2077 AMsmallest.name <- names.xdata[findx] 2078 } 2079 2080 for( i in findx:nvars) 2081 { 2082 count <- i-findx+1 2083 if(print.level > 0) 2084 cat("\n***** (V",count,") ", names.xdata[i]," *****\n",sep="") 2085 2086 ks.do <- FALSE 2087 is.dummy <- length(unique( xdata[,i] )) < 3 2088 if (ks & !is.dummy) 2089 ks.do <- TRUE 2090 2091 BeforeMatchingBalance[[count]] <- balanceUV(xdata[,i][Tr==1], xdata[,i][Tr==0], nboots=0, 2092 weights.Tr=weights[Tr==1], weights.Co=weights[Tr==0], 2093 match=FALSE) 2094 2095 if (BeforeMatchingBalance[[count]]$tt$p.value < BMsmallest.p.value) 2096 { 2097 BMsmallest.p.value <- BeforeMatchingBalance[[count]]$tt$p.value 2098 BMsmallest.number <- count 2099 BMsmallest.name <- names.xdata[i] 2100 } else if (BeforeMatchingBalance[[count]]$tt$p.value == BMsmallest.p.value) 2101 { 2102 BMsmallest.number <- c(BMsmallest.number,count) 2103 BMsmallest.name <- c(BMsmallest.name,names.xdata[i]) 2104 } 2105 2106 if (ks.do) 2107 { 2108 BeforeMatchingBalance[[count]]$ks <- list() 2109 BeforeMatchingBalance[[count]]$ks$ks <- list() 2110 BeforeMatchingBalance[[count]]$ks$ks$p.value <- ks.bm$ks.naive.pval[count] 2111 BeforeMatchingBalance[[count]]$ks$ks$statistic <- ks.bm$ks.stat[count] 2112 if (nboots > 0) 2113 { 2114 BeforeMatchingBalance[[count]]$ks$ks.boot.pvalue <- ks.bm$ks.boot.pval[count] 2115 2116 if (ks.bm$ks.boot.pval[count] < BMsmallest.p.value) 2117 { 2118 BMsmallest.p.value <- ks.bm$ks.boot.pval[count] 2119 BMsmallest.number <- count 2120 BMsmallest.name <- names.xdata[i] 2121 } else if ( (ks.bm$ks.boot.pval[count] == BMsmallest.p.value) & (sum(BMsmallest.number==count)==0) ) 2122 { 2123 BMsmallest.number <- c(BMsmallest.number,count) 2124 BMsmallest.name <- c(BMsmallest.name,names.xdata[i]) 2125 } 2126 } else { 2127 BeforeMatchingBalance[[count]]$ks$ks.boot.pvalue <- NA 2128 2129 if (ks.bm$ks.naive.pval[count] < BMsmallest.p.value) 2130 { 2131 BMsmallest.p.value <- ks.bm$ks.naive.pval[count] 2132 BMsmallest.number <- count 2133 BMsmallest.name <- names.xdata[i] 2134 } else if ( (ks.bm$ks.naive.pval[count] == BMsmallest.p.value) & (sum(BMsmallest.number==count)==0) ) 2135 { 2136 BMsmallest.number <- c(BMsmallest.number,count) 2137 BMsmallest.name <- c(BMsmallest.name,names.xdata[i]) 2138 } 2139 } 2140 2141 } else { 2142 BeforeMatchingBalance[[count]]$ks <- NULL 2143 } 2144 2145 if (!is.null(match.out)) 2146 { 2147 AfterMatchingBalance[[count]] <- balanceUV(xdata[,i][match.out$index.treated], 2148 xdata[,i][match.out$index.control], 2149 weights=match.out$weights, nboots=0, 2150 paired=paired, match=TRUE) 2151 2152 if (AfterMatchingBalance[[count]]$tt$p.value < AMsmallest.p.value) 2153 { 2154 AMsmallest.p.value <- AfterMatchingBalance[[count]]$tt$p.value 2155 AMsmallest.number <- count 2156 AMsmallest.name <- names.xdata[i] 2157 } else if ( (AfterMatchingBalance[[count]]$tt$p.value == AMsmallest.p.value) & (sum(AMsmallest.number==count)==0) ) 2158 { 2159 AMsmallest.number <- c(AMsmallest.number,count) 2160 AMsmallest.name <- c(AMsmallest.name,names.xdata[i]) 2161 } 2162 2163 if (ks.do) 2164 { 2165 AfterMatchingBalance[[count]]$ks <- list() 2166 AfterMatchingBalance[[count]]$ks$ks <- list() 2167 AfterMatchingBalance[[count]]$ks$ks$p.value <- ks.am$ks.naive.pval[count] 2168 AfterMatchingBalance[[count]]$ks$ks$statistic <- ks.am$ks.stat[count] 2169 if (nboots > 0) 2170 { 2171 AfterMatchingBalance[[count]]$ks$ks.boot.pvalue <- ks.am$ks.boot.pval[count] 2172 2173 if (ks.am$ks.boot.pval[count] < AMsmallest.p.value) 2174 { 2175 AMsmallest.p.value <- ks.am$ks.boot.pval[count] 2176 AMsmallest.number <- count 2177 AMsmallest.name <- names.xdata[i] 2178 } else if ( (ks.am$ks.boot.pval[count] == AMsmallest.p.value) & (sum(AMsmallest.number==count)==0) ) 2179 { 2180 AMsmallest.number <- c(AMsmallest.number,count) 2181 AMsmallest.name <- c(AMsmallest.name,names.xdata[i]) 2182 } 2183 } else { 2184 AfterMatchingBalance[[count]]$ks$ks.boot.pvalue <- NA 2185 2186 if (ks.am$ks.naive.pval[count] < AMsmallest.p.value) 2187 { 2188 AMsmallest.p.value <- ks.am$ks.naive.pval[count] 2189 AMsmallest.number <- count 2190 AMsmallest.name <- names.xdata[i] 2191 } else if ( (ks.am$ks.naive.pval[count] == AMsmallest.p.value) & (sum(AMsmallest.number==count)==0) ) 2192 { 2193 AMsmallest.number <- c(AMsmallest.number,count) 2194 AMsmallest.name <- c(AMsmallest.name,names.xdata[i]) 2195 } 2196 } 2197 } else { 2198 AfterMatchingBalance[[count]]$ks <- NULL 2199 } 2200 if(print.level > 0) 2201 PrintBalanceUV(BeforeMatchingBalance[[count]], AfterMatchingBalance[[count]], digits=digits) 2202 } else { 2203 if(print.level > 0) 2204 { 2205 cat("before matching:\n") 2206 summary(BeforeMatchingBalance[[count]], digits=digits) 2207 } 2208 } #end of if match.out 2209 } #end of for loop 2210 2211 if(print.level & ( (nvars-findx+1) > 1)) 2212 { 2213 cat("\n") 2214 2215 if (BMsmallest.p.value < 1) 2216 { 2217 cat("Before Matching Minimum p.value:", format.pval(BMsmallest.p.value, digits=digits),"\n") 2218 cat("Variable Name(s):",BMsmallest.name, " Number(s):",BMsmallest.number,"\n\n") 2219 } else { 2220 cat("Before Matching Minimum p.value: 1\n\n") 2221 } 2222 2223 if (!is.null(match.out)) 2224 { 2225 if(AMsmallest.p.value < 1) 2226 { 2227 cat("After Matching Minimum p.value:", format.pval(AMsmallest.p.value, digits=digits),"\n") 2228 cat("Variable Name(s):",AMsmallest.name, " Number(s):",AMsmallest.number,"\n\n") 2229 } else { 2230 cat("After Matching Minimum p.value: 1\n\n") 2231 } 2232 } #end of !is.null(match.out) 2233 }#end of print.level & (nvars > 1) 2234 2235 return(invisible(list(BeforeMatching=BeforeMatchingBalance, 2236 AfterMatching=AfterMatchingBalance, 2237 BMsmallest.p.value=BMsmallest.p.value, 2238 BMsmallestVarName=BMsmallest.name, 2239 BMsmallestVarNumber=BMsmallest.number, 2240 AMsmallest.p.value=AMsmallest.p.value, 2241 AMsmallestVarName=AMsmallest.name, 2242 AMsmallestVarNumber=AMsmallest.number))) 2243 } #end of MatchBalance 2244 2245 2246get.xdata <- function(formul, datafr) { 2247 t1 <- terms(formul, data=datafr); 2248 if (length(attr(t1, "term.labels"))==0 & attr(t1, "intercept")==0) { 2249 m <- NULL; # no regressors specified for the model matrix 2250 } else { 2251 m <- model.matrix(formul, data=datafr, drop.unused.levels = TRUE) 2252 } 2253 return(m); 2254} 2255 2256 2257# get.ydata: 2258# Return response vector corresponding to the formula in formul 2259# 2260get.ydata <- function(formul, datafr) { 2261 t1 <- terms(formul, data=datafr); 2262 if (length(attr(t1, "response"))==0) { 2263 m <- NULL; # no response variable specified 2264 } else { 2265 m <- model.response(model.frame(formul, data=datafr)) 2266 } 2267 return(m); 2268} 2269 2270# 2271# bootstrap ks test implemented. Fast version 2272# 2273 2274ks.boot <- function(Tr, Co, nboots=1000, alternative = c("two.sided", "less", "greater"), print.level=0) 2275 { 2276 alternative <- match.arg(alternative) 2277 tol <- sqrt(.Machine$double.eps) 2278 Tr <- Tr[!is.na(Tr)] 2279 Co <- Co[!is.na(Co)] 2280 2281 w <- c(Tr, Co) 2282 obs <- length(w) 2283 n.x <- length(Tr) 2284 n.y <- length(Co) 2285 cutp <- n.x 2286 ks.boot.pval <- NULL 2287 bbcount <- 0 2288 2289 if (nboots < 10) 2290 { 2291 nboots <- 10 2292 warning("At least 10 'nboots' must be run; seting 'nboots' to 10") 2293 } 2294 2295 if (nboots < 500) 2296 warning("For publication quality p-values it is recommended that 'nboots'\n be set equal to at least 500 (preferably 1000)") 2297 2298 fs.ks <- Mks.test(Tr, Co, alternative=alternative) 2299 2300 if (alternative=="two.sided") 2301 { 2302 if (print.level > 0) 2303 cat("ks.boot: two.sided test\n") 2304 for (bb in 1:nboots) 2305 { 2306 if (print.level > 1) 2307 cat("s:", bb, "\n") 2308 2309 sindx <- sample(1:obs, obs, replace=TRUE) 2310 2311 X1tmp <- w[sindx[1:cutp]] 2312 X2tmp <- w[sindx[(cutp+1):obs]] 2313 2314 s.ks <- ks.fast(X1tmp, X2tmp, n.x=n.x, n.y=n.y, n=obs) 2315 2316 if (s.ks >= (fs.ks$statistic - tol) ) 2317 bbcount <- bbcount + 1 2318 } 2319 } else { 2320 if (print.level > 0) 2321 cat("ks.boot:",alternative,"test\n") 2322 for (bb in 1:nboots) 2323 { 2324 if (print.level > 1) 2325 cat("s:", bb, "\n") 2326 2327 sindx <- sample(1:obs, obs, replace=TRUE) 2328 2329 X1tmp <- w[sindx[1:cutp]] 2330 X2tmp <- w[sindx[(cutp+1):obs]] 2331 2332 s.ks <- Mks.test(X1tmp, X2tmp, alternative=alternative)$statistic 2333 2334 if (s.ks >= (fs.ks$statistic - tol) ) 2335 bbcount <- bbcount + 1 2336 } 2337 } 2338 ks.boot.pval <- bbcount/nboots 2339 2340 ret <- list(ks.boot.pvalue=ks.boot.pval, ks=fs.ks, nboots=nboots) 2341 class(ret) <- "ks.boot" 2342 2343 return(ret) 2344 } #end of ks.boot 2345 2346summary.ks.boot <- function(object, ..., digits=5) 2347 { 2348 if(!is.list(object)) { 2349 warning("object not a valid 'ks.boot' object") 2350 return() 2351 } 2352 2353 if (!inherits(object, "ks.boot")) { 2354 warning("Object not of class 'ks.boot'") 2355 return() 2356 } 2357 2358 cat("\n") 2359 cat("Bootstrap p-value: ", format.pval(object$ks.boot.pvalue, digits=digits), "\n") 2360 cat("Naive p-value: ", format(object$ks$p.value, digits=digits), "\n") 2361 cat("Full Sample Statistic:", format(object$ks$statistic, digits=digits), "\n") 2362# cat("nboots completed ", object$nboots, "\n") 2363 cat("\n") 2364 2365 z <- list() 2366 class(z) <- "summary.ks.boot" 2367 return(invisible(z)) 2368 } #end of summary.ks.boot 2369 2370print.summary.ks.boot <- function(x, ...) 2371 { 2372 invisible(NULL) 2373 } 2374 2375 2376RmatchLoop <- function(Y, Tr, X, Z, V, All, M, BiasAdj, Weight, Weight.matrix, Var.calc, weight, 2377 SAMPLE, ccc, cdd, ecaliper=NULL, exact=NULL, caliper=NULL, restrict=NULL, 2378 MatchLoopC.indx=NULL, weights.flag, replace=TRUE, ties=TRUE, 2379 version="standard", MatchbyAI=FALSE) 2380 { 2381 s1 <- MatchGenoudStage1caliper(Tr=Tr, X=X, All=All, M=M, weights=weight, 2382 exact=exact, caliper=caliper, 2383 distance.tolerance=cdd, 2384 tolerance=ccc) 2385 2386 sum.caliper.drops <- 0 2387 X.orig <- X 2388 2389#are we using the restriction matrix? 2390 if(is.matrix(restrict)) { 2391 restrict.trigger <- TRUE 2392 } else { 2393 restrict.trigger <- FALSE 2394 } 2395 2396# if SATC is to be estimated the treatment indicator is reversed 2397 if (All==2) 2398 Tr <- 1-Tr 2399 2400# check on the number of matches, to make sure the number is within the limits 2401# feasible given the number of observations in both groups. 2402 if (All==1) 2403 { 2404 M <- min(M,min(sum(Tr),sum(1-Tr))); 2405 } else { 2406 M <- min(M,sum(1-Tr)); 2407 } 2408 2409# two slippage parameters that are used to determine whether distances are equal 2410# distances less than ccc or cdd are interpeted as zero. 2411# these are passed in. ccc, cdd 2412 2413 2414# I. set up 2415# I.a. vector for which observations we want the average effect 2416# iot_t is the vector with weights in the average treatment effects 2417# iot_c is the vector of indicators for potential controls 2418 2419 if (All==1) 2420 { 2421 iot.t <- weight; 2422 iot.c <- as.matrix(rep(1,length(Tr))) 2423 } else { 2424 iot.t <- Tr*weight; 2425 iot.c <- 1-Tr 2426 } 2427 2428# I.b. determine sample and covariate vector sizes 2429 N <- nrow(X) 2430 Kx <- ncol(X) 2431 Kz <- ncol(Z) 2432 2433# K covariates 2434# N observations 2435# Nt <- sum(Tr) 2436# Nc <- sum(1-Tr) 2437# on <- as.matrix(rep(1,N)) 2438 2439# I.c. normalize regressors to have mean zero and unit variance. 2440# If the standard deviation of a variable is zero, its normalization 2441# leads to a variable with all zeros. 2442# The matrix AA enables one to transform the user supplied weight matrix 2443# to take account of this transformation. BUT THIS IS NOT USED!! 2444# Mu_X and Sig_X keep track of the original mean and variances 2445# AA <- diag(Kx) 2446 Mu.X <- matrix(0, Kx, 1) 2447 Sig.X <- matrix(0, Kx, 1) 2448 2449 for (k in 1:Kx) 2450 { 2451 Mu.X[k,1] <- sum(X[,k]*weight)/sum(weight) 2452 eps <- X[,k]-Mu.X[k,1] 2453 Sig.X[k,1] <- sqrt(max(ccc, sum(X[,k]*X[,k]*weight))/sum(weight)-Mu.X[k,1]^2) 2454 Sig.X[k,1] <- Sig.X[k,1]*sqrt(N/(N-1)) 2455 if(Sig.X[k,1] < ccc) 2456 Sig.X[k,1] <- ccc 2457 X[,k]=eps/Sig.X[k,1] 2458# AA[k,k]=Sig.X[k,1] 2459 } #end of k loop 2460 2461# Nv <- nrow(V) 2462 Mv <- ncol(V) 2463 Mu.V <- matrix(0, Mv, 1) 2464 Sig.V <- matrix(0, Mv, 1) 2465 2466 for (j in 1:Mv) 2467 { 2468 Mu.V[j,1]= ( t(V[,j])%*%weight ) /sum(weight) 2469# dv <- V[,j]-Mu.V[j,1] 2470 sv <- sum(V[,j]*V[,j]*weight)/sum(weight)-Mu.V[j,1]^2 2471 if (sv > 0) 2472 { 2473 sv <- sqrt(sv) 2474 } else { 2475 sv <- 0 2476 } 2477 sv <- sv * sqrt(N/(N-1)) 2478 Sig.V[j,1] <- sv 2479 } #end of j loop 2480 2481# I.d. define weight matrix for metric, taking into account normalization of 2482# regressors. 2483# If the eigenvalues of the regressors are too close to zero the Mahalanobis metric 2484# is not used and we revert back to the default inverse of variances. 2485 if (Weight==1) 2486 { 2487 Weight.matrix=diag(Kx) 2488 } else if (Weight==2) { 2489 if (min (eigen( t(X)%*%X/N, only.values=TRUE)$values) > 0.0000001) 2490 { 2491 Weight.matrix= solve(t(X)%*%X/N) 2492 } else { 2493 Weight.matrix <- diag(Kx) 2494 } 2495 } 2496 # DO NOT RESCALE THE Weight.matrix!! 2497 #else if (Weight==3) 2498 # { 2499 # Weight.matrix <- AA %*% Weight.matrix %*% AA 2500 # } 2501 2502# if (exact==1) 2503# { 2504# Weight.matrix <- cbind(Weight.matrix, matrix(0,nrow=Kx,ncol=Mv)) 2505# Weight.matrix <- rbind(Weight.matrix, cbind(matrix(0,nrow=Mv,ncol=Kx), 2506# 1000*solve(diag(as.vector(Sig.V*Sig.V), nrow=length(Sig.V))))) 2507# Weight.matrix <- as.matrix(Weight.matrix) 2508# X <- cbind(X,V) 2509# Mu.X <- rbind(Mu.X, matrix(0, nrow(Mu.V), 1)) 2510# Sig.X <- rbind(Sig.X, matrix(1, nrow(Sig.V), 1)) 2511# } #end of exact 2512 2513 if ( min(eigen(Weight.matrix, only.values=TRUE)$values) < ccc ) 2514 Weight.matrix <- Weight.matrix + diag(Kx)*ccc 2515 2516 ww <- chol(Weight.matrix) # so that ww*ww=w.m 2517 2518 if(is.null(s1$ecaliper)) 2519 { 2520 caliperflag <- 0 2521 use.ecaliper <- 0 2522 } else { 2523 caliperflag <- 1 2524 use.ecaliper <- s1$ecaliper 2525 } 2526 2527 #if we have a diagonal matrix we can void cblas_dgemm 2528 if (Kx > 1) 2529 { 2530 DiagWeightMatrixFlag <- as.double(sum( (Weight.matrix!=diag(diag(Weight.matrix))) )==0) 2531 } else { 2532 DiagWeightMatrixFlag <- 1 2533 } 2534 2535 if(is.null(MatchLoopC.indx)) 2536 { 2537 #indx: 2538 # 1] I (unadjusted); 2] IM (unadjusted); 3] weight; 4] I (adjusted); 5] IM (adjusted) 2539 if(weights.flag==TRUE) 2540 { 2541 MatchLoopC.indx <- MatchLoopC(N=s1$N, xvars=Kx, All=s1$All, M=s1$M, 2542 cdd=cdd, caliperflag=caliperflag, replace=replace, ties=ties, 2543 ww=ww, Tr=s1$Tr, Xmod=s1$X, 2544 weights=weight, 2545 CaliperVec=use.ecaliper, Xorig=X.orig, 2546 restrict.trigger=restrict.trigger, restrict=restrict, 2547 DiagWeightMatrixFlag=DiagWeightMatrixFlag) 2548 } else { 2549 MatchLoopC.indx <- MatchLoopCfast(N=s1$N, xvars=Kx, All=s1$All, M=s1$M, 2550 cdd=cdd, caliperflag=caliperflag, replace=replace, ties=ties, 2551 ww=ww, Tr=s1$Tr, Xmod=s1$X, 2552 CaliperVec=use.ecaliper, Xorig=X.orig, 2553 restrict.trigger=restrict.trigger, restrict=restrict, 2554 DiagWeightMatrixFlag=DiagWeightMatrixFlag) 2555 } 2556 } 2557 2558 indx <- MatchLoopC.indx 2559 2560 if(indx[1,1]==0) 2561 { 2562 ret <- list() 2563 ret$valid <- 0 2564 if (caliperflag) 2565 { 2566 ret$sum.caliper.drops <- indx[1,6] 2567 } else { 2568 ret$sum.caliper.drops <- 0 2569 } 2570 return(ret) 2571 } 2572 #we now keep going if we only have 1 valid match 2573 #else if (nrow(indx)< 2) 2574 # { 2575 # ret <- list() 2576 # ret$valid <- 1 2577 # if (caliperflag) 2578 # { 2579 # ret$sum.caliper.drops <- indx[1,6] 2580 # } else { 2581 # ret$sum.caliper.drops <- 0 2582 # } 2583 # return(ret) 2584 # } 2585 2586 if (All==2) 2587 { 2588 foo <- indx[,5] 2589 indx[,5] <- indx[,4] 2590 indx[,4] <- foo 2591 } 2592 2593 if (caliperflag) 2594 { 2595 sum.caliper.drops <- indx[1,6] 2596 } else { 2597 sum.caliper.drops <- 0 2598 } 2599 2600 # 2601 # Generate variables which we need later on 2602 # 2603 2604 I <- indx[,1] 2605 IT <- Tr[indx[,1]] 2606 IM <- indx[,2] 2607 2608# IX <- X[indx[,1],] 2609# Xt <- X[indx[,4],] 2610# Xc <- X[indx[,5],] 2611 2612# IY <- Y[indx[,1]] 2613 Yt <- Y[indx[,4]] 2614 Yc <- Y[indx[,5]] 2615 2616 W <- indx[,3] 2617 2618 if(BiasAdj==1 & sum(W) < ncol(Z)) 2619 { 2620 warning("Fewer (weighted) matches than variables in 'Z': BiasAdjust set to FALSE") 2621 BiasAdj=0 2622 } 2623 2624 if(BiasAdj==1) 2625 { 2626 if(sum(W) < ncol(Z)) 2627 { 2628 warning("Fewer matches than variables for Bias Adjustment") 2629 } 2630 2631 IZ <- Z[indx[,1],] 2632 Zt <- Z[indx[,4],] 2633 Zc <- Z[indx[,5],] 2634 } 2635 2636 est.func <- function(N, All, Tr, indx, weight, BiasAdj, Kz) 2637 { 2638 Kcount <- as.matrix(rep(0, N)) 2639 KKcount <- as.matrix(rep(0, N)) 2640 YCAUS <- matrix(0, nrow=N, ncol=1) 2641 2642 if (BiasAdj==1) 2643 { 2644 ZCAUS <- matrix(0, nrow=N, ncol=Kz) 2645 } else { 2646 ZCAUS <- NULL 2647 } 2648 2649 for (i in 1:N) 2650 { 2651 if ( ( Tr[i]==1 & All!=1) | All==1 ) 2652 { 2653 2654 foo.indx <- indx[,1]==i 2655 foo.indx2 <- foo.indx 2656 2657 sum.foo <- sum(foo.indx) 2658 2659 if (sum.foo < 1) 2660 next; 2661 2662 foo <- rep(FALSE, N) 2663 foo.indx <- indx[foo.indx,2] 2664 foo[foo.indx] <- rep(TRUE,sum.foo) 2665 2666 # inner.func <- function(N, weight, indx, foo.indx2, Y, Tr, foo) 2667 2668 Kcount <- Kcount + weight[i] * weight*foo/sum(foo*weight) 2669 2670 KKcount <- KKcount + weight[i]*weight*weight*foo/ 2671 (sum(foo*weight)*sum(foo*weight)) 2672 2673 foo.indx2.2 <- indx[foo.indx2,2]; 2674 foo.indx2.3 <- indx[foo.indx2,3]; 2675 sum.foo.indx2.3 <- sum(foo.indx2.3) 2676 2677 if(Tr[i]==1) 2678 { 2679 YCAUS[i] <- Y[i] - sum((Y[foo.indx2.2]*foo.indx2.3))/sum.foo.indx2.3 2680 } else { 2681 YCAUS[i] <- sum((Y[foo.indx2.2]*foo.indx2.3))/sum.foo.indx2.3 - Y[i] 2682 } 2683 2684 if (BiasAdj==1) 2685 { 2686 2687 if(Tr[i]==1) 2688 { 2689 if (sum.foo > 1) 2690 { 2691 ZCAUS[i,] <- 2692 Z[i,] - t(Z[foo.indx2.2,]) %*% foo.indx2.3/sum.foo.indx2.3 2693 } else { 2694 ZCAUS[i,] <- 2695 Z[i,] - Z[foo.indx2.2,]*foo.indx2.3/sum.foo.indx2.3 2696 } 2697 } else { 2698 if (sum.foo > 1) 2699 { 2700 ZCAUS[i,] <- t(Z[foo.indx2.2,]) %*% foo.indx2.3/sum.foo.indx2.3 - Z[i,] 2701 } else { 2702 ZCAUS[i,] <- Z[foo.indx2.2,]*foo.indx2.3/sum.foo.indx2.3 - Z[i,] 2703 } 2704 } 2705 } #endof BiasAdj 2706 } 2707 } #end of if 2708 return(list(YCAUS=YCAUS,ZCAUS=ZCAUS,Kcount=Kcount,KKcount=KKcount)) 2709 } #end of est.func 2710 2711 if(version=="standard" & BiasAdj==0) 2712 { 2713 ret <- .Call("EstFuncC", as.integer(N), as.integer(All), as.integer(nrow(indx)), 2714 as.double(Y), as.double(Tr), 2715 as.double(weight), as.double(indx), 2716 PACKAGE="Matching") 2717 YCAUS <- ret[,1]; 2718 Kcount <- ret[,2]; 2719 KKcount <- ret[,3]; 2720 } else if (version=="standard") { 2721 ret.est <- est.func(N=N, All=All, Tr=Tr, indx=indx, weight=weight, BiasAdj=BiasAdj, Kz=Kz) 2722 YCAUS <- ret.est$YCAUS 2723 ZCAUS <- ret.est$ZCAUS 2724 Kcount <- ret.est$Kcount 2725 KKcount <- ret.est$KKcount 2726 } 2727 2728 if (All!=1) 2729 { 2730 I <- as.matrix(I[IT==1]) 2731 IT <- as.matrix(IT[IT==1]) 2732 Yc <- as.matrix(Yc[IT==1]) 2733 Yt <- as.matrix(Yt[IT==1]) 2734 W <- as.matrix(W[IT==1]) 2735 if (BiasAdj==1) 2736 { 2737 if (Kz > 1) 2738 { 2739 Zc <- as.matrix(Zc[IT==1,]) 2740 Zt <- as.matrix(Zt[IT==1,]) 2741 IZ <- as.matrix(IZ[IT==1,]) 2742 } else{ 2743 Zc <- as.matrix(Zc[IT==1]) 2744 Zt <- as.matrix(Zt[IT==1]) 2745 IZ <- as.matrix(IZ[IT==1]) 2746 } 2747 } 2748# IM <- as.matrix(IM[IT==1,]) 2749# IY <- as.matrix(IY[IT==1]) 2750# IX.u <- as.matrix(IX.u[IT==1,]) 2751# Xc.u <- as.matrix(Xc.u[IT==1,]) 2752# Xt.u <- as.matrix(Xt.u[IT==1,]) 2753# Xc <- as.matrix(Xc[IT==1,]) 2754# Xt <- as.matrix(Xt[IT==1,]) 2755 } #end of if 2756 2757 if (length(I) < 1) 2758 { 2759 return(list(sum.caliper.drops=N)) 2760 } 2761 2762 if (BiasAdj==1) 2763 { 2764 # III. Regression of outcome on covariates for matches 2765 if (All==1) 2766 { 2767 # number of observations 2768 NNt <- nrow(Z) 2769 # add intercept 2770 ZZt <- cbind(rep(1, NNt), Z) 2771 # number of covariates 2772 Kx <- ncol(ZZt) 2773 xw <- ZZt*(sqrt(Tr*Kcount) %*% t(as.matrix(rep(1,Kx)))) 2774 2775 foo <- min(eigen(t(xw)%*%xw, only.values=TRUE)$values) 2776 foo <- as.double(foo<=ccc) 2777 foo2 <- apply(xw, 2, sd) 2778 2779 options(show.error.messages = FALSE) 2780 wout <- NULL 2781 try(wout <- solve( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 2782 (t(xw) %*% (Y*sqrt(Tr*Kcount)))) 2783 if(is.null(wout)) 2784 { 2785 wout2 <- NULL 2786 try(wout2 <- ginv( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 2787 (t(xw) %*% (Y*sqrt(Tr*Kcount)))) 2788 if(!is.null(wout2)) 2789 { 2790 wout <-wout2 2791 warning("using generalized inverse to calculate Bias Adjustment probably because of singular 'Z'") 2792 } 2793 } 2794 options(show.error.messages = TRUE) 2795 if(is.null(wout)) 2796 { 2797 warning("unable to calculate Bias Adjustment probably because of singular 'Z'") 2798 BiasAdj <- 0 2799 } else { 2800 NW <- nrow(wout) 2801# KW <- ncol(wout) 2802 Alphat <- wout[2:NW,1] 2803 } 2804 } else { 2805 Alphat <- matrix(0, nrow=Kz, ncol=1) 2806 } #end if ALL 2807 } 2808 2809 if(BiasAdj==1) 2810 { 2811 # III.b. Controls 2812 NNc <- nrow(Z) 2813 ZZc <- cbind(matrix(1, nrow=NNc, ncol=1),Z) 2814 Kx <- ncol(ZZc) 2815 2816 xw <- ZZc*(sqrt((1-Tr)*Kcount) %*% matrix(1, nrow=1, ncol=Kx)) 2817 2818 foo <- min(eigen(t(xw)%*%xw, only.values=TRUE)$values) 2819 foo <- as.double(foo<=ccc) 2820 foo2 <- apply(xw, 2, sd) 2821 2822 options(show.error.messages = FALSE) 2823 wout <- NULL 2824 try(wout <- solve( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 2825 (t(xw) %*% (Y*sqrt((1-Tr)*Kcount)))) 2826 if(is.null(wout)) 2827 { 2828 wout2 <- NULL 2829 try(wout2 <- ginv( t(xw) %*% xw + diag(Kx) * ccc * (foo) * max(foo2)) %*% 2830 (t(xw) %*% (Y*sqrt((1-Tr)*Kcount)))) 2831 if(!is.null(wout2)) 2832 { 2833 wout <-wout2 2834 warning("using generalized inverse to calculate Bias Adjustment probably because of singular 'Z'") 2835 } 2836 } 2837 options(show.error.messages = TRUE) 2838 if(is.null(wout)) 2839 { 2840 warning("unable to calculate Bias Adjustment probably because of singular 'Z'") 2841 BiasAdj <- 0 2842 } else { 2843 NW <- nrow(wout) 2844# KW <- ncol(wout) 2845 Alphac <- as.matrix(wout[2:NW,1]) 2846 } 2847 } 2848 2849 if(BiasAdj==1) 2850 { 2851 # III.c. adjust matched outcomes using regression adjustment for bias adjusted matching estimator 2852 IZ <- as.matrix(IZ) 2853 Zc <- as.matrix(Zc) 2854 Zt <- as.matrix(Zt) 2855 Alphat <- as.matrix(Alphat) 2856 2857 SCAUS <- YCAUS-Tr*(ZCAUS %*% Alphac)-(1-Tr)*(ZCAUS %*% Alphat) 2858 # adjusted control outcome 2859 Yc.adj <- Yc+BiasAdj * (IZ-Zc) %*% Alphac 2860 # adjusted treated outcome 2861 Yt.adj <- Yt+BiasAdj*(IZ-Zt) %*% Alphat 2862 Tau.i <- Yt.adj - Yc.adj 2863 } else { 2864 Yc.adj <- Yc 2865 Yt.adj <- Yt 2866 Tau.i <- Yt.adj - Yc.adj 2867 } 2868 2869 art.data <- cbind(I,IM) 2870 2871 # III. If conditional variance is needed, initialize variance vector 2872 # and loop through all observations 2873 2874 if (Var.calc>0) 2875 { 2876 #For R version of this function see Matching version < 4.5-0008 2877 Sigs <- VarCalcMatchC(N=N, xvars=ncol(X), Var.calc=Var.calc, 2878 cdd=cdd, caliperflag=caliperflag, 2879 ww=ww, Tr=Tr, Xmod=s1$X, 2880 CaliperVec=use.ecaliper, Xorig=X.orig, 2881 restrict.trigger=restrict.trigger, restrict=restrict, 2882 DiagWeightMatrixFlag=DiagWeightMatrixFlag, 2883 Y=Y, weightFlag=weights.flag, weight=weight) 2884 } #end of var.calc > 0 2885 2886 est <- t(W) %*% Tau.i/sum(W) # matching estimator 2887 2888 if(version=="standard") 2889 { 2890 if (Var.calc==0) 2891 { 2892 eps <- Tau.i - as.double(est) 2893 eps.sq <- eps*eps 2894 Sigs <- 0.5 * matrix(1, N, 1) %*% (t(eps.sq) %*% W)/sum(W) #sss <- sqrt(Sigs[1,1]) 2895 } #end of Var.calc==0 2896 2897 SN <- sum(iot.t) 2898 var.sample <- sum((Sigs*(iot.t+iot.c*Kcount)*(iot.t+iot.c*Kcount))/(SN*SN)) 2899 2900 if (All==1) 2901 { 2902 var.pop <- sum((Sigs*(iot.c*Kcount*Kcount+2*iot.c*Kcount-iot.c*KKcount))/(SN*SN)) 2903 } else { 2904 var.pop=sum((Sigs*(iot.c*Kcount*Kcount-iot.c*KKcount))/(SN*SN)) 2905 } 2906 2907 if (BiasAdj==1) 2908 { 2909 dvar.pop <- sum(iot.t*(SCAUS-as.double(est))*(SCAUS-as.double(est)))/(SN*SN) 2910 } else { 2911 dvar.pop <- sum(iot.t*(YCAUS-as.double(est))*(YCAUS-as.double(est)))/(SN*SN) 2912 } 2913 2914 var.pop <- var.pop + dvar.pop 2915 2916 if (SAMPLE==1) 2917 { 2918 var <- var.sample 2919 } else { 2920 #var <- max(var.sample, var.pop) 2921 var <- var.pop 2922 } 2923 2924 var.cond <- max(var.sample,var.pop)-var.sample 2925 se <- sqrt(var) 2926 se.cond <- sqrt(var.cond) 2927 #Sig <- sqrt(Sigs) 2928 } else { 2929 se=NULL 2930 se.cond=NULL 2931 } 2932 2933 if (All==2) 2934 est <- -1*est 2935 2936# if (exact==1) 2937# { 2938# Vt.u <- Xt.u[,(Kx-Mv+1):Kx] 2939# Vc.u <- Xc.u[,(Kx-Mv+1):Kx] 2940# Vdif <- abs(Vt.u-Vc.u) 2941# 2942# if (Mv>1) 2943# Vdif <- as.matrix(apply(t(Vdif), 2, sum)) 2944# 2945# em[1,1] <- length(Vdif) 2946# em[2,1] <- sum(Vdif>0.000001) 2947# }#end of exact==1 2948 2949 if(!MatchbyAI) 2950 { 2951 return(list(est=est, se=se, se.cond=se.cond, W=W, 2952 sum.caliper.drops=sum.caliper.drops, 2953 art.data=art.data, 2954 MatchLoopC=MatchLoopC.indx)) 2955 } else { 2956 if(Var.calc==0) 2957 Sigs <- NULL 2958 return(list(est=est, se=se, se.cond=se.cond, W=W, 2959 sum.caliper.drops=sum.caliper.drops, 2960 art.data=art.data, 2961 MatchLoopC=MatchLoopC.indx, 2962 YCAUS=YCAUS, Kcount=Kcount, KKcount=KKcount, 2963 Sigs=Sigs)) 2964 } 2965 }# end of RmatchLoop 2966 2967MatchLoopC <- function(N, xvars, All, M, cdd, caliperflag, replace, ties, ww, Tr, Xmod, weights, CaliperVec, Xorig, 2968 restrict.trigger, restrict, DiagWeightMatrixFlag) 2969 { 2970 2971 if(restrict.trigger) 2972 { 2973 restrict.nrow <- nrow(restrict) 2974 } else { 2975 restrict.nrow <- 0 2976 } 2977 2978 ret <- .Call("MatchLoopC", as.integer(N), as.integer(xvars), as.integer(All), as.integer(M), 2979 as.double(cdd), as.integer(caliperflag), as.integer(replace), as.integer(ties), 2980 as.double(ww), as.double(Tr), 2981 as.double(Xmod), as.double(weights), as.double(CaliperVec), as.double(Xorig), 2982 as.integer(restrict.trigger), as.integer(restrict.nrow), as.double(restrict), 2983 as.double(DiagWeightMatrixFlag), 2984 PACKAGE="Matching") 2985 return(ret) 2986 } #end of MatchLoopC 2987 2988MatchLoopCfast <- function(N, xvars, All, M, cdd, caliperflag, replace, ties, ww, Tr, Xmod, CaliperVec, 2989 Xorig, restrict.trigger, restrict, DiagWeightMatrixFlag) 2990 { 2991 2992 if(restrict.trigger) 2993 { 2994 restrict.nrow <- nrow(restrict) 2995 } else { 2996 restrict.nrow <- 0 2997 } 2998 2999 ret <- .Call("MatchLoopCfast", as.integer(N), as.integer(xvars), as.integer(All), as.integer(M), 3000 as.double(cdd), as.integer(caliperflag), as.integer(replace), as.integer(ties), 3001 as.double(ww), as.double(Tr), 3002 as.double(Xmod), as.double(CaliperVec), as.double(Xorig), 3003 as.integer(restrict.trigger), as.integer(restrict.nrow), as.double(restrict), 3004 as.double(DiagWeightMatrixFlag), 3005 PACKAGE="Matching") 3006 return(ret) 3007 } #end of MatchLoopCfast 3008 3009VarCalcMatchC <- function(N, xvars, Var.calc, cdd, caliperflag, ww, Tr, Xmod, CaliperVec, 3010 Xorig, restrict.trigger, restrict, DiagWeightMatrixFlag, 3011 Y, weightFlag, weight) 3012 { 3013 3014 if(restrict.trigger) 3015 { 3016 restrict.nrow <- nrow(restrict) 3017 } else { 3018 restrict.nrow <- 0 3019 } 3020 3021 ret <- .Call("VarCalcMatchC", as.integer(N), as.integer(xvars), as.integer(Var.calc), 3022 as.double(cdd), as.integer(caliperflag), 3023 as.double(ww), as.double(Tr), 3024 as.double(Xmod), as.double(CaliperVec), as.double(Xorig), 3025 as.integer(restrict.trigger), as.integer(restrict.nrow), as.double(restrict), 3026 as.double(DiagWeightMatrixFlag), 3027 as.double(Y), 3028 as.integer(weightFlag), as.double(weight), 3029 PACKAGE="Matching") 3030 return(ret) 3031 } #end of VarCalcMatchC 3032 3033 3034 3035.onAttach <- function( ... ) 3036{ 3037 MatchLib <- dirname(system.file(package = "Matching")) 3038 version <- packageDescription("Matching", lib.loc = MatchLib)$Version 3039 BuildDate <- packageDescription("Matching", lib.loc = MatchLib)$Date 3040 3041 foo <- paste("## \n## Matching (Version ", version, ", Build Date: ", BuildDate, ")\n", 3042 "## See http://sekhon.berkeley.edu/matching for additional documentation.\n", 3043 "## Please cite software as:\n", 3044 "## Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching\n", 3045 "## Software with Automated Balance Optimization: The Matching package for R.''\n", 3046 "## Journal of Statistical Software, 42(7): 1-52. \n##\n", 3047 sep = "") 3048 packageStartupMessage(foo) 3049} 3050 3051 3052