[Gmm-commits] r181 - in pkg: causalGel causalGel/R causalGel/man causalGel/src gmm gmm/man gmm/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 9 23:33:59 CET 2021
Author: chaussep
Date: 2021-02-09 23:33:58 +0100 (Tue, 09 Feb 2021)
New Revision: 181
Added:
pkg/causalGel/R/otherCausal.R
pkg/causalGel/man/causalfit-class.Rd
pkg/causalGel/man/otherCausal.Rd
pkg/causalGel/src/
pkg/causalGel/src/Makevars
pkg/causalGel/src/causal.f
pkg/causalGel/src/causalGel.h
pkg/causalGel/src/src.c
Modified:
pkg/causalGel/DESCRIPTION
pkg/causalGel/NAMESPACE
pkg/causalGel/R/allClasses.R
pkg/causalGel/R/causalGel.R
pkg/causalGel/man/causalGEL.Rd
pkg/gmm/DESCRIPTION
pkg/gmm/man/charStable.Rd
pkg/gmm/vignettes/empir.bib
pkg/gmm/vignettes/gmm_with_R.pdf
pkg/gmm/vignettes/gmm_with_R.rnw
Log:
added other causal methods to causaGel and fixed vignette of gmm
Modified: pkg/causalGel/DESCRIPTION
===================================================================
--- pkg/causalGel/DESCRIPTION 2020-11-09 13:42:09 UTC (rev 180)
+++ pkg/causalGel/DESCRIPTION 2021-02-09 22:33:58 UTC (rev 181)
@@ -1,6 +1,6 @@
Package: causalGel
-Version: 0.1-1
-Date: 2020-11-09
+Version: 0.2
+Date: 2020-11-23
Title: Causal Inference using Generalized Empirical
Likelihood Methods
Author: Pierre Chausse <pchausse at uwaterloo.ca>
@@ -7,12 +7,12 @@
Maintainer: Pierre Chausse <pchausse at uwaterloo.ca>
Description: Methods for causal inference in which covariates are balanced using generalized empirical likelihod methods.
Depends: R (>= 3.0.0), momentfit (>= 0.2)
-Imports: stats, methods
+Imports: stats, methods, quadprog
Suggests: lmtest, knitr, texreg
Collate: 'allClasses.R' 'causalMethods.R' 'rcausalMethods.R' 'causalGel.R'
- 'causalfitMethods.R'
+ 'causalfitMethods.R' 'otherCausal.R'
License: GPL (>= 2)
-NeedsCompilation: no
+NeedsCompilation: yes
VignetteBuilder: knitr
Modified: pkg/causalGel/NAMESPACE
===================================================================
--- pkg/causalGel/NAMESPACE 2020-11-09 13:42:09 UTC (rev 180)
+++ pkg/causalGel/NAMESPACE 2021-02-09 22:33:58 UTC (rev 181)
@@ -1,8 +1,11 @@
+useDynLib(causalGel, .registration = TRUE, .fixes="F_")
import("momentfit")
-importFrom("stats", "lm", "model.response", "terms", "model.frame", "reformulate", "as.formula")
+importFrom("stats", "lm", "model.response", "terms", "model.frame",
+ "reformulate", "as.formula", "binomial", "fitted", "glm",
+ "lm.wfit", "optim", "plogis")
importFrom("utils", "head", "tail")
-
+importFrom("quadprog", "solve.QP")
importFrom("methods", is, new, show, "slot<-", "slotNames", "validObject",
"getClassDef", "selectMethod", "callNextMethod", "as", "setAs",
"getMethod")
@@ -10,11 +13,12 @@
### S4 Methods and Classes
exportClasses()
-exportClasses("causalData", "causalModel", "causalGelfit", "rcausalModel")
+exportClasses("causalData", "causalModel", "causalGelfit", "rcausalModel",
+ "causalfit")
exportMethods("causalMomFct", "checkConv")
-export("causalModel", "causalGEL")
+export("causalModel", "causalGEL", "matching", "LLmatching", "ipw")
Modified: pkg/causalGel/R/allClasses.R
===================================================================
--- pkg/causalGel/R/allClasses.R 2020-11-09 13:42:09 UTC (rev 180)
+++ pkg/causalGel/R/allClasses.R 2021-02-09 22:33:58 UTC (rev 181)
@@ -17,6 +17,14 @@
setClass("causalGelfit", contains="gelfit")
+setClass("causalfit", representation(estim="numeric",
+ type="character",
+ method="character",
+ form="list",
+ details="list",
+ info="list",
+ data="data.frame",
+ call="callORNULL"))
## converters
Modified: pkg/causalGel/R/causalGel.R
===================================================================
--- pkg/causalGel/R/causalGel.R 2020-11-09 13:42:09 UTC (rev 180)
+++ pkg/causalGel/R/causalGel.R 2021-02-09 22:33:58 UTC (rev 181)
@@ -16,6 +16,53 @@
U[, 1:r, drop = FALSE]
}
+.getLamEEL <- function (gmat, k=1, restrictedLam = integer(), ...)
+{
+ gmat <- as.matrix(gmat)
+ if (length(restrictedLam))
+ {
+ if (length(restrictedLam) > ncol(gmat))
+ stop("The number of restricted Lambda exceeds the number of moments")
+ if (!all(restrictedLam %in% (1:ncol(gmat))))
+ stop(paste("restrictedLam must be a vector of integers between 1 and ",
+ ncol(gmat), sep = ""))
+ gmat <- gmat[, -restrictedLam, drop = FALSE]
+ } else {
+ restrictedLam <- integer()
+ }
+ res <- .EEL_quad(gmat, k = 1)
+ if (res$convergence == 2)
+ {
+ res <- EEL_lam(gmat, k = 1)
+ res$convergence <- 2
+ }
+ if (length(restrictedLam)) {
+ lambda <- numeric(ncol(gmat) + length(restrictedLam))
+ lambda[-restrictedLam] <- res$lambda
+ res$lambda <- lambda
+ }
+ res
+}
+
+.EEL_quad <- function (gmat, k = 1)
+{
+ Dmat <- crossprod(gmat)/nrow(gmat)
+ dvec <- -colMeans(gmat)
+ Amat <- t(gmat)
+ bvec <- rep(-1, nrow(gmat))
+ res <- try(solve.QP(Dmat, dvec, Amat, bvec), silent=TRUE)
+ if (any(class(res) == "try-error"))
+ {
+ conv <- list(convergence = 2)
+ lambda0 <- rep(0, ncol(gmat))
+ } else {
+ conv <- list(convergence = 0)
+ lambda0 <- res$solution
+ }
+ list(lambda = lambda0, convergence = conv,
+ obj = mean(rhoEEL(gmat, lambda0, 0, k)))
+}
+
causalModel <- function(g, balm, data,theta0=NULL,
momType=c("ACE","ACT","ACC", "uncondBal"),
popMom = NULL, ACTmom=1L, orthoBases=FALSE)
@@ -101,9 +148,10 @@
initTheta <- match.arg(initTheta)
coefSlv <- match.arg(coefSlv)
gelType <- match.arg(gelType)
+ if (gelType == "REEL")
+ lamSlv <- .getLamEEL
if (initTheta=="theta0" & is.null(theta0))
stop("theta0 is required when initTheta='theta0'")
-
model <- causalModel(g, balm, data, theta0, momType, popMom, ACTmom,
orthoBases)
Added: pkg/causalGel/R/otherCausal.R
===================================================================
--- pkg/causalGel/R/otherCausal.R (rev 0)
+++ pkg/causalGel/R/otherCausal.R 2021-02-09 22:33:58 UTC (rev 181)
@@ -0,0 +1,531 @@
+###################################
+######## Matching methods
+##################################
+
+## units in x2 are matched to each element of x1
+## minn is the minimum number of matches
+## tol: if other individuals are within tol of
+## the highest distance in the first minn matches, then they are also included.
+## y2 is y from group 2
+## it returns:
+## list in indices for the match
+## E(y2 | z=1, x) (n1 x 1 vector of estimated missing Y(2))
+
+getNN <- function(x,y,z,fromTo, minn, tol=1e-7)
+{
+ if (length(fromTo)!=2)
+ stop("fromTo must be a vector of 2 with 0s or 1s")
+ if (!all(fromTo %in% c(0,1)))
+ stop("fromTo must be a vector of 2 with 0s or 1s")
+ x <- as.matrix(x)
+ self <- z==fromTo[1]
+ selt <- z==fromTo[2]
+ x1 <- x[self,,drop=FALSE]
+ x2 <- x[selt,,drop=FALSE]
+ y2 <- y[selt]
+ n1 <- nrow(x1)
+ n2 <- nrow(x2)
+ k <- ncol(x1)
+ type <- paste("Nearest neighbour of the ",
+ ifelse(fromTo[1]==1, "treated", "control"), " among the ",
+ ifelse(fromTo[2]==1, "treated", "control"), sep="")
+ res <- .Fortran(F_findnn, as.double(x1), as.double(x2),as.double(y2),
+ as.double(tol),
+ as.integer(n1), as.integer(n2), as.integer(k),
+ as.integer(minn), K=double(n2), nn=integer(n1),
+ ind=integer(n1*n2), y2hat=double(n1))
+ nn <- res$nn
+ ind <- matrix(res$ind,nrow=n1)
+ list(matches=lapply(1:n1, function(i) ind[i,1:nn[i]]), K=res$K,
+ y2hat = res$y2hat, type=type)
+}
+
+## Function to get the missing counterfactual in y
+## if which="y0", y[z==1] are replace by an estimate Y(0)|z=1
+## if which="y1", y[z==0] are replace by an estimate Y(1)|z=0
+## if type="match", only match method is used
+## if type="reg", only regression is used
+## if type="bc_match" the bias-corrected method is used.
+## x is the matrix used for matching
+## regMat is the regression matrix for E(Y|X,Z=1) or E(Y|X,Z=0)
+
+getMissingY <- function(y, z, x, which=c("y0", "y1"),
+ type=c("match","bc_match","reg","all"), minn=4, tol=1e-7,
+ regMat=cbind(1,x))
+{
+ which <- match.arg(which)
+ type <- match.arg(type)
+ ## regMat should have an intercept included
+ x <- as.matrix(x)
+ if (which == "y0")
+ {
+ sel <- z==1
+ fromTo <- c(1,0)
+ } else {
+ sel <- z==0
+ fromTo <- c(0,1)
+ }
+ res <- getNN(x,y, z, fromTo, minn, tol)
+ if (type != "match" | type == "all")
+ b <- lm.wfit(regMat[!sel,,drop=FALSE], y[!sel], res$K)$coefficients
+ if (type == "match")
+ {
+ y[sel] <- res$y2hat
+ y <- as.matrix(y)
+ colnames(y) <- "match"
+ } else if (type == "reg") {
+ y[sel] <- c(regMat[sel,,drop=FALSE]%*%b)
+ y <- as.matrix(y)
+ colnames(y) <- "reg"
+ } else if (type=="bc_match") {
+ fity <- c(regMat%*%b)
+ bc <- fity[sel] - sapply(res$matches, function(l) mean(fity[!sel][l]))
+ y[sel] <- res$y2hat + bc
+ y <- as.matrix(y)
+ colnames(y) <- "bc"
+ } else {
+ ym <- y
+ ym[sel] <- res$y2hat
+ yreg <- y
+ yreg[sel] <- c(regMat[sel,,drop=FALSE]%*%b)
+ ybc <- y
+ fity <- c(regMat%*%b)
+ bc <- fity[sel] - sapply(res$matches, function(l) mean(fity[!sel][l]))
+ ybc[sel] <- res$y2hat + bc
+ y <- cbind(ym, ybc, yreg)
+ colnames(y) <- c("match","bc","reg")
+ }
+ attr(y, "which") <- which
+ y
+}
+
+## PSForm is the formula used to compute the propensity score.
+## BCorForm is the bias correction regression model
+
+matching <- function(form, balm, data, type=c("ACE","ACT","ACC"), M=4,
+ adjust=FALSE, matchPS=FALSE, psForm=NULL,
+ bcForm=NULL)
+{
+ ## If bcForm does not have an intercept it will be added.
+ type <- match.arg(type)
+ Call <- match.call()
+ mf <- model.frame(form, data)
+ Y <- c(mf[[1]])
+ T <- c(mf[[2]])
+ if (is.null(bcForm) & adjust)
+ bcForm <- balm
+ if (is.null(psForm) & matchPS)
+ {
+ psForm <- attr(terms(balm), "term.labels")
+ psForm <- reformulate(psForm, response=colnames(mf)[2])
+ }
+ BC <- NA
+ if (!matchPS)
+ {
+ balm <- attr(terms(balm), "term.labels")
+ balm <- reformulate(balm, intercept=FALSE)
+ X <- model.matrix(balm, data)
+ } else {
+ res <- glm(psForm, data=data, family=binomial())
+ X <- as.matrix(fitted(res))
+ data$PScore <- X
+ balm <- ~PScore-1
+ }
+ X <- scale(X)
+ if (adjust)
+ {
+ bcForm <- reformulate(attr(terms(bcForm), "term.labels"))
+ Z <- model.matrix(bcForm, data)
+ } else {
+ Z <- NULL
+ }
+ typeCor <- ifelse(adjust, "all", "match")
+ if (type == "ACT")
+ {
+ Y0 <- getMissingY(Y, T, X, which="y0", type=typeCor, minn=M,
+ regMat=Z)
+ NN <- mean(Y[T==1]-Y0[T==1,"match"])
+ if (adjust)
+ BC <- mean(Y[T==1]-Y0[T==1,"bc"])
+ } else if (type == "ACC") {
+ Y1 <- getMissingY(Y, T, X, which="y1", type=typeCor, minn=M,
+ regMat=Z)
+ NN <- mean(Y1[T==0,"match"] - Y[T==0])
+ if (adjust)
+ BC <- mean(Y1[T==0,"bc"] - Y[T==0])
+ } else {
+ Y0 <- getMissingY(Y, T, X, which="y0", type=typeCor, minn=M,
+ regMat=Z)
+ Y1 <- getMissingY(Y, T, X, which="y1", type=typeCor, minn=M,
+ regMat=Z)
+ NN <- mean(Y1[,"match"]-Y0[,"match"])
+ if (adjust)
+ BC <- mean(Y1[,"bc"]-Y0[,"bc"])
+ }
+ estim <- c(NN=NN, NN.BiasCor=BC)
+ form <- list(balForm=balm, psForm=psForm, bcForm=bcForm)
+ info <- list()
+ method <- ifelse(matchPS, "PS Matching Method", "Covariate Matching Method")
+ details <- list(ifelse(matchPS,
+ paste("PS formula: ", deparse(psForm), sep=""),
+ paste("Covariate formula: ",
+ deparse(balm), sep="")),
+ paste("Min. number of matches: ", M, sep=""))
+ if (adjust)
+ details <- c(details, list(paste("Bias-cor. formula: ",
+ deparse(bcForm),sep="")))
+ shortm <- paste(c("","BiasCor_"),
+ paste(ifelse(matchPS, "ps","cov"),
+ "Matching", M, "_", type, sep=""),sep="")
+ names(shortm) <- c("NN","NN.BiasCor")
+ shortMet <- shortm
+ new("causalfit", estim=estim, type=type, method=method,
+ form=form, details=details, info=info, data=data, call=Call)
+}
+
+
+
+### ACE using local linear regression matching
+### if h is set, it is used, if not, the optimal h is computed
+### by minimizing the CV. A grid in seq(from, to, length=nh) is
+### first obtain and the brent method is used to complete.
+### In the output, if info=0, it converged, if info=1, the minimum
+### is at from or to, so the brent merhod was not launch, and for info = 2
+### the brent reached the minimum iteration.
+### the method requires propensity score so PSForm must be provided
+
+LLmatching <- function(form, psForm, data, type=c("ACE","ACT","ACC"),
+ kern=c("Gaussian","Epanechnikov"),tol=1e-4,
+ h=NULL, from=.00001, to=5, ngrid=10, maxit=100,
+ hMethod=c("Brent","Grid"))
+{
+ type <- match.arg(type)
+ Call <- match.call()
+ kern <- match.arg(kern)
+ hMethod <- match.arg(hMethod)
+ mf <- model.frame(form, data)
+ Y <- c(mf[[1]])
+ T <- c(mf[[2]])
+ res <- glm(psForm, data=data, family=binomial())
+ X <- as.matrix(fitted(res))
+ if (type == "ACT")
+ {
+ Y0 <- getMissingY.LL(Y, T, X, "y0", h, kern,
+ from, to, ngrid, tol, maxit, hMethod)
+ ACE <- mean(Y[T==1]-Y0[T==1], na.rm=TRUE)
+ info <- attr(Y0, "h")
+ } else if (type == "ACC") {
+ Y1 <- getMissingY.LL(Y, T, X, "y1", h, kern,
+ from, to, ngrid, tol, maxit, hMethod)
+ info <- attr(Y1, "h")
+ ACE <- mean(Y1[T==0] - Y[T==0], na.rm=TRUE)
+ } else {
+ Y0 <- getMissingY.LL(Y, T, X, "y0", h, kern,
+ from, to, ngrid, tol, maxit, hMethod)
+ Y1 <- getMissingY.LL(Y, T, X, "y1", h, kern,
+ from, to, ngrid, tol, maxit, hMethod)
+ info <- list(Y0=attr(Y0, "h"), Y1=attr(Y1, "h"))
+ ACE <- mean(Y1-Y0, na.rm=TRUE)
+ }
+ estim <- c(LL=ACE)
+ info2 <- unlist(info)
+ method <- "Local-Linear Regression Method"
+ details <- list(paste("Kernel: ", kern, sep=""),
+ paste("PS formula: ", deparse(psForm), sep=""),
+ ifelse(is.null(h), paste("Bandwidth select: ",hMethod,sep=""),
+ "Bandwidth select: Fixed"))
+ if (type == "ACE")
+ {
+ tmp <- paste("h (Y(0)) = ", round(info[[1]][[1]], 6), sep="")
+ tmp2 <- paste("h (Y(1)) = ", round(info[[2]][[1]], 6), sep="")
+ details <- c(details, tmp, tmp2)
+ } else {
+ tmp <- paste("h = ", round(info[[1]], 6), sep="")
+ details <- c(details, tmp)
+ }
+ if (is.null(h))
+ {
+ if (type == "ACE")
+ {
+ tmp <- paste("Conv. code for bandwidth (Y(0)): ",
+ info[[1]][[2]], sep="")
+ tmp2 <- paste("Conv. code for bandwidth (Y(1)): ",
+ info[[2]][[2]], sep="")
+ details <- c(details, tmp, tmp2)
+ } else {
+ details <- c(details,
+ list(paste("Conv. code for bandwidth: ",
+ info[[2]],sep="")))
+ }
+ }
+ if (type == "ACE")
+ info2 <- list(meanh = mean(c(info[[1]][[1]], info[[2]][[1]])),
+ convergence = max(info[[1]][[2]], info[[2]][[2]]))
+ shortm <- paste("LL", "_", strtrim(kern,3),"_",
+ ifelse(is.null(h), "autoh", "fixedh"),"_", type, sep="")
+ shortMet <- shortm
+ new("causalfit", estim=estim, type=type, method=method,
+ form=list(psForm=psForm),data=data,
+ details=details, info=info2, call=Call)
+}
+
+## p is the fitted propensity score
+## The function estimate the counterfactual Y(0) or Y(1)
+
+getMissingY.LL <- function(y, z, p, which=c("y0", "y1"), h=NULL,
+ kern=c("Gaussian","Epanechnikov"),
+ from, to, ngrid, tol, maxit,
+ hMethod=c("Brent","Grid"))
+{
+ which <- match.arg(which)
+ hMethod <- match.arg(hMethod)
+ kern <- match.arg(kern)
+ if (which == "y0")
+ sel <- z==1
+ else
+ sel <- z==0
+ missY <- llrF(p[sel], p[!sel], y[!sel], h,
+ kern, from, to, ngrid, tol, maxit, hMethod)
+ y[sel] <- missY
+ attr(y, "which") <- which
+ attr(y, "h") <- attr(missY, "h")
+ y
+}
+
+cvF <- function(p, y, h, kern=c("Gaussian","Epanechnikov"))
+{
+ kern <- match.arg(kern)
+ nh <- length(h)
+ kernF <- tolower(strtrim(kern,1))
+ n <- length(p)
+ res <- .Fortran(F_cvfct, as.double(p), as.double(y),
+ as.integer(n), as.double(h), as.character(kernF),
+ as.integer(nh), cv=double(nh))
+ res$cv
+}
+
+llrF <- function(p1, p0, y0, h=NULL, kern=c("Gaussian","Epanechnikov"),
+ from, to, ngrid, tol, maxit, hMethod=c("Brent","Grid"))
+{
+ hMethod <- match.arg(hMethod)
+ kern <- match.arg(kern)
+ kernF <- ifelse(kern=="Gaussian", 1L, 2L)
+ n1 <- length(p1)
+ n2 <- length(p0)
+ convergence <- NA
+ if (is.null(h))
+ {
+ resh <- optCVF(p0, y0, kern, from, to, ngrid,
+ tol, maxit, hMethod)
+ h <- resh$h
+ convergence <- resh$info
+ }
+ res <- .Fortran(F_llr, as.double(p1), as.double(p0), as.double(y0),
+ as.integer(n1), as.integer(n2), as.double(h),
+ as.integer(kernF), info = integer(n1),
+ y0hat = double(n1))
+ y0hat <- res$y0hat
+ y0hat[res$info == 1] <- NA
+ attr(y0hat, "h") <- list(h=h,convergence=convergence)
+ y0hat
+}
+
+
+gridcvF <- function(p, y, kern=c("Gaussian","Epanechnikov"), from, to, nh)
+{
+ kern <- match.arg(kern)
+ kernF <- ifelse(kern=="Gaussian", 1L, 2L)
+ n <- length(p)
+ res <- .Fortran(F_gridcv, as.double(p), as.double(y),
+ as.integer(n), as.integer(kernF),
+ as.double(from), as.double(to), as.integer(nh),
+ info=integer(1), h=double(3), cv=double(3))
+ cbind(res$h,res$cv)
+}
+
+### The gricv fortran subroutine returns, if info=0, a 3x1 vector of h and cv,
+### such that cv[1]>cv[2] and cv[3]>cv[2]. The minimum is therefore the middle one.
+
+optCVF <- function(p, y, kern=c("Gaussian","Epanechnikov"), from, to, nh,
+ tol=1e-4, maxit=100, method=c("Brent","Grid"))
+{
+ kern <- match.arg(kern)
+ method <- match.arg(method)
+ kernF <- ifelse(kern=="Gaussian", 1L, 2L)
+ n <- length(p)
+ res <- .Fortran(F_gridcv, as.double(p), as.double(y),
+ as.integer(n), as.integer(kernF),
+ as.double(from), as.double(to), as.integer(nh),
+ info=integer(1), h=double(3), cv=double(3))
+ if (res$info == 1)
+ {
+ return(list(cv=res$cv[1], h=res$h[1], info=res$info))
+ }
+ if (method == "Grid")
+ {
+ return(list(cv=res$cv[2], h=res$h[2], info=res$info))
+ }
+ res2 <- .Fortran(F_brentcv, as.double(p), as.double(y),
+ as.integer(n), as.integer(kernF),
+ as.double(tol), as.integer(maxit),
+ as.double(res$h), as.double(res$cv),
+ info=integer(1), h=double(1), cv=double(1))
+ info <- ifelse(res2$info == 0, 0, 2)
+ list(cv=res2$cv, h=res2$h, info=info)
+}
+
+
+### Weighting Methods
+########################
+
+ipw <- function(form, psForm, data, type=c("ACE","ACT","ACC"),
+ normalized=FALSE, ...)
+{
+ type <- match.arg(type)
+ Call <- match.call()
+ mf <- model.frame(form, data)
+ Y <- c(mf[[1]])
+ T <- c(mf[[2]])
+ mnorm <- normalized
+ info <- list()
+ if (normalized == "GPE")
+ {
+ res <- getGPE.PS(form, psForm, data, ...)
+ X <- res$phat
+ info <- list(res$info)
+ normalized <- TRUE
+ } else {
+ res <- glm(psForm, data=data, family=binomial())
+ X <- as.matrix(fitted(res))
+ }
+ if (type == "ACT")
+ {
+ s1 <- sum(T)
+ w1 <- rep(1/s1, s1)
+ w0 <- X[T==0]/(1-X[T==0])/s1
+ if (normalized)
+ w0 <- w0/sum(w0)
+ } else if (type == "ACC") {
+ s0 <- sum(1-T)
+ w1 <- (1-X[T==1])/X[T==1]/s0
+ w0 <- rep(1/s0,s0)
+ if (normalized)
+ w1 <- w1/sum(w1)
+ } else {
+ n <- length(X)
+ w1 <- 1/(n*X[T==1])
+ w0 <- 1/(n*(1-X[T==0]))
+ if (normalized)
+ {
+ w1 <- w1/sum(w1)
+ w0 <- w0/sum(w0)
+ }
+ }
+ ACE <- sum(w1*Y[T==1]) - sum(w0*Y[T==0])
+ estim <- c(IPW=ACE)
+ form <- list(psForm=psForm)
+ method <- "Inverse Probability Weight Method"
+ details <- list(paste("PS formula: ", deparse(psForm), sep=""))
+ if (mnorm == "GPE")
+ details <- c(details,
+ list(paste("Method: normalized with GPE method")))
+ else
+ details <- c(details,
+ list(paste("Method: ", ifelse(mnorm, "normalized",
+ "non-normalized"), sep="")))
+ if (mnorm == "GPE")
+ details <- c(details,
+ paste("Convergence of optim for PS: ",
+ info[[1]][2], sep=""))
+ shortm <- paste("IPW", "_",
+ ifelse(mnorm=="GPE", "GPE",
+ ifelse(mnorm, "normalized", "unnormalized")),
+ "_", type, sep="")
+ shortMet <- shortm
+ new("causalfit", estim=estim, type=type, method=method,
+ form=list(psForm=psForm),
+ details=details, info=info, call=Call)
+}
+
+
+getGPE.PS <- function(form, PSForm=z~kX, data, ...)
+ {
+ mf <- model.frame(form, data)
+ Y <- c(mf[[1]])
+ T <- c(mf[[2]])
+ X <- model.matrix(PSForm, data)
+ g <- function(theta, T, X)
+ {
+ kX <- c(X%*%theta)
+ G <- plogis(kX)
+ e <- (T-G)/(1-G)
+ m <- as.matrix(e*X)
+ mean(colMeans(m, na.rm=TRUE)^2)
+ }
+ res <- glm(PSForm, data=data, family=binomial())
+ res1 <- optim(coef(res), g, method="BFGS", T=T, X=X, ...)
+ b <- res1$par
+ fit <- c(X%*%b)
+ phat <- plogis(fit)
+ phat[phat>.999] <- .999
+ list(phat=phat, info=c(obj=res1$value, convergence=res1$convergence))
+ }
+
+
+## print
+
+setMethod("print", "gmmfit",
+ function(x, model=TRUE, ...) {
+ theta <- coef(x)
+ if (model)
+ print(x at model)
+ ntype <- matrix(c("Two-Step GMM", "Iterated GMM", "CUE",
+ "One-Step GMM with fixed weights","Two-Stage Least Squares",
+ "Evaluated at a fixed Theta; No estimation",
+ "One-Step Efficient M.D.E.",
+ "twostep","iter","cue","onestep","tsls", "eval","mde"),
+ ncol=2)
+ type <- ntype[match(x at type, ntype[,2]),1]
+ spec <- modelDims(x at model)
+ if (spec$q==spec$k && x at type != "eval")
+ type <- "One-Step, Just-Identified"
+ cat("\nEstimation: ", type,"\n")
+ if (!is.null(x at convergence))
+ cat("Convergence Optim: ", x at convergence, "\n")
+ if (!is.null(x at convIter))
+ cat("Convergence Iteration: ", x at convIter, "\n")
+ cat("coefficients:\n")
+ print.default(format(theta, ...), print.gap=2L, quote=FALSE)
+ })
+
+## show
+
+setMethod("show","causalfit", function(object) print(object))
+
+## print
+
+setMethod("print", "causalfit",
+ function(x, digits=5, ...) {
+ cat(x at method, "\n")
+ cat("**************************\n")
+ cat("Type of causal effect: ", x at type, "\n", sep="")
+ if (length(x at estim)==1)
+ {
+ cat("Causal estimate: ",
+ format(x at estim, digits=digits, ...), "\n")
+ } else {
+ cat("Causal estimates: \n")
+ for (i in 1:length(x at estim))
+ cat("\t", names(x at estim)[i], " = ",
+ format(x at estim[i], digits=digits, ...),
+ "\n", sep="")
+ }
+ if (length(x at details))
+ {
+ cat("Details:\n")
+ for (i in 1:length(x at details))
+ cat("\t", x at details[[i]], "\n")
+ }
+ invisible(x)
+ })
+
Modified: pkg/causalGel/man/causalGEL.Rd
===================================================================
--- pkg/causalGel/man/causalGEL.Rd 2020-11-09 13:42:09 UTC (rev 180)
+++ pkg/causalGel/man/causalGEL.Rd 2021-02-09 22:33:58 UTC (rev 181)
@@ -118,5 +118,11 @@
fit2 <- causalGEL(g, balm, nsw, gelType="EL")
fit2
+## This is REEL, which is EEL with restriction in pt to be positive.
+## It is based on the quadprog package
+
+fit3 <- causalGEL(g, balm, nsw, gelType="REEL")
+fit3
+
}
Added: pkg/causalGel/man/causalfit-class.Rd
===================================================================
--- pkg/causalGel/man/causalfit-class.Rd (rev 0)
+++ pkg/causalGel/man/causalfit-class.Rd 2021-02-09 22:33:58 UTC (rev 181)
@@ -0,0 +1,33 @@
+\name{causalfit-class}
+\docType{class}
+\alias{causalfit-class}
+
+\title{Class \code{"causalfit"}}
+\description{
+This is the class for causal effect estimates other than GEL.
+}
+\section{Objects from the Class}{
+Objects can be created by calls of the form \code{new("causalfit", ...)}.
+But, it is created by \code{\link{matching}}, \code{\link{LLmatching}},
+or \code{\link{ipw}}.
+}
+\section{Slots}{
+ \describe{
+ \item{\code{estim}:}{Object of class \code{"numeric"} ~~ }
+ \item{\code{type}:}{Object of class \code{"character"} ~~ }
+ \item{\code{method}:}{Object of class \code{"character"} ~~ }
+ \item{\code{form}:}{Object of class \code{"list"} ~~ }
+ \item{\code{details}:}{Object of class \code{"list"} ~~ }
+ \item{\code{info}:}{Object of class \code{"list"} ~~ }
+ \item{\code{data}:}{Object of class \code{"data.frame"} ~~ }
+ \item{\code{call}:}{Object of class \code{"callORNULL"} ~~ }
+ }
+}
+\section{Methods}{
+No methods defined with class "causalfit" in the signature.
+}
+
+\examples{
+showClass("causalfit")
+}
+\keyword{classes}
Added: pkg/causalGel/man/otherCausal.Rd
===================================================================
--- pkg/causalGel/man/otherCausal.Rd (rev 0)
+++ pkg/causalGel/man/otherCausal.Rd 2021-02-09 22:33:58 UTC (rev 181)
@@ -0,0 +1,106 @@
+\name{otherCausal}
+\alias{otherCausal}
+\alias{matching}
+\alias{LLmatching}
+\alias{ipw}
+
+\title{Popular Causal Effect estimation methods}
+
+\description{ This documentation file presents a collection of popular
+methods used to estimate the average causal effect, the causal effect on
+the treated and the causal effect on the non-treated.
+}
+
+\usage{
+
+matching(form, balm, data, type=c("ACE","ACT","ACC"), M=4,
+ adjust=FALSE, matchPS=FALSE, psForm=NULL,
+ bcForm=NULL)
+
+LLmatching(form, psForm, data, type=c("ACE","ACT","ACC"),
+ kern=c("Gaussian","Epanechnikov"),tol=1e-4,
+ h=NULL, from=.00001, to=5, ngrid=10, maxit=100,
+ hMethod=c("Brent","Grid"))
+
+ipw(form, psForm, data, type=c("ACE","ACT","ACC"),
+ normalized=FALSE, ...)
+}
+\arguments{
+
+ \item{form}{A formula that links the outcome variable to the treatment
+ indicator. For the moment, only one treatment group is allowed.}
+
+ \item{balm}{A formula or a matrix with balancing covariates to be
+ matched.}
+
+ \item{data}{A data.frame or a matrix with column names.}
+
+ \item{type}{The type of causal effect to compute. \code{ACE} stands
+ for average causal effect, \code{ACT} for causal effect on the treated
+ and \code{ACC} for causal effect on the control or non-treated.}
+
+ \item{M}{The minimum number of matches}
+
+ \item{adjust}{Should we apply a bias correction to the estimate.}
+
+ \item{matchPS}{Should we match the propensity scores instead of the
+ covariates.}
+
+ \item{psForm}{It is the \code{formula} argument passed to
+ \code{\link{glm}} to estimate the propensity scores by a logistic
+ regression.}
+
+ \item{bcForm}{A formula that represents the right hand side in the
+ regression used for the bias correction.}
+
+ \item{kern}{The type of kernel to use in the local linear regression
+ method.}
+
+ \item{tol}{The tolerance level for the stopping rule used to compute
+ the optimal bandwidth.}
+
+ \item{h}{A fixed bandwidth. By default, the optimal bandwidth is found
+ by minimizing a cross-validation.}
+
+ \item{from}{The lower bound for the search of the optimal bandwidth.}
+
+ \item{to}{The upper bound for the search of the optimal bandwidth.}
+
+ \item{ngrid}{The number of grid points if the optimal bandwidth is
+ obtained by grid search.}
+
+ \item{maxit}{The maximum number of iterations for the minimization of
+ the cross-validation.}
+
+ \item{hMethod}{The method used to find the optimal bandwidth.}
+
+ \item{normalized}{Should the weights be normalized. If set to
+ \code{GPE}, the GPE method is used.}
+
+ \item{...}{Additional arguments to pass to \code{\link{optim}} when
+ the \code{GPE} method is used.}
+
+}
+
+\value{
+ All methods return an object of classes \code{"causalfit"}.
+ }
+
+\examples{
+data(nsw)
+
+balm <- ~age+ed+black+hisp:married+nodeg+re75+I(re75^2)
+g <- re78~treat
+ps <- treat~age+ed+black+hisp:married+nodeg+re75+I(re75^2)
+
+fit1 <- matching(g, balm, nsw, "ACE", adjust=TRUE)
+fit1
+
+fit2 <- LLmatching(g, ps, nsw, "ACE")
+fit2
+
+fit3 <- ipw(g, ps, nsw, "ACE", normalized=TRUE)
+fit3
+
+}
+
Added: pkg/causalGel/src/Makevars
===================================================================
--- pkg/causalGel/src/Makevars (rev 0)
+++ pkg/causalGel/src/Makevars 2021-02-09 22:33:58 UTC (rev 181)
@@ -0,0 +1 @@
+PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lgomp
\ No newline at end of file
Added: pkg/causalGel/src/causal.f
===================================================================
--- pkg/causalGel/src/causal.f (rev 0)
+++ pkg/causalGel/src/causal.f 2021-02-09 22:33:58 UTC (rev 181)
@@ -0,0 +1,322 @@
+ subroutine edist(x1,x2,n,k,dist)
+ integer n, i
+ double precision x1(k), x2(n,k), dist(n)
+ do i=1,n
+ call exprsum((x1-x2(i,:))**2, k, dist(i))
+ end do
+ end
+
+ subroutine allmin(x,n,nm,loc)
+ integer n, nm, loc(n), indi
+ double precision x(n), tmp1, tmp2
+ logical sel(n)
+
+ nm = 0
+ sel = .true.
+ do
+ indi = minloc(x,dim=1,mask=sel)
+ sel(indi) = .false.
+ if (nm == 0) then
+ tmp1 = x(indi)
+ nm = 1
+ loc(1) = indi
+ cycle
+ end if
+ tmp2 = x(indi)
+ if (tmp1 == tmp2) then
+ nm = nm+1
+ loc(nm) = indi
+ tmp1 = tmp2
+ else
+ exit
+ end if
+ end do
+ end
+
+c x1 is for group 1 and x2 fro group 2 (either could be treated)
+c Here units in x2 are matched to THE unit in x1
+
+ subroutine findnni(x1, x2, tol, n, k, minn, nn, ind)
+ integer n, k, minn, nn, ind(n), indi
+ double precision x1(k), x2(n,k), dist(n), tmp1, tmp2, tol
+ logical sel(n)
+
+ nn = 0
+ sel = .true.
+ ind = 0
+
+ call edist(x1,x2,n,k,dist)
+ do
+ indi = minloc(dist,dim=1,mask=sel)
+ if (indi == 0) then
+ exit
+ end if
+ sel(indi) = .false.
+ if (nn==0) then
+ tmp1 = sqrt(dist(indi))
+ ind(1) = indi
+ nn = 1
+ cycle
+ end if
+ tmp2 = sqrt(dist(indi))
+ if (nn < minn) then
+ nn = nn+1
+ ind(nn) = indi
+ else
+ if (abs(tmp2-tmp1)<tol) then
+ nn = nn+1
+ ind(nn) = indi
+ else
+ exit
+ end if
+ end if
+ tmp1 = tmp2
+ end do
+ end
+
+c Here units in x2 are matched to every unit in x1
+c Here we also estimate the missing Y(2) [ Y(2)|Z=1 ] by averaging the Y's
+c from group 2 matched to group 1.
+c Y(2) is therefore n1 x 1 using Yi from group 2.
+c Also, rk is the weights attached to y2. it sums to n1
+
+ subroutine findnn(x1, x2, y2, tol, n1, n2, k, minn, rk, nn, ind,
+ * y2h)
+ integer n1, n2, k, minn, nn(n1), ind(n1,n2), i
+ double precision x1(n1,k), x2(n2,k), tol, rk(n2), y2(n2), y2h(n1)
+
+ rk = 0.0d0
+ do i=1,n1
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 181
More information about the Gmm-commits
mailing list