causalGel.R
In causalGel: Causal Inference using Generalized Empirical Likelihood Methods

## ----extract, message=FALSE, warning=FALSE, echo=FALSE------------------------
library(causalGel)    
library(texreg)
setMethod("extract", "causalGelfit", 
          function(model, includeSpecTest=FALSE, 
                   specTest=c("LR","LM","J"), include.nobs=TRUE, 
                   include.obj.fcn=TRUE, ...)
              {
                  specTest <- match.arg(specTest)
                  s <- summary(model, ...)
                  wspecTest <- grep(specTest, rownames(s@specTest@test))
                  spec <- modelDims(model@model)
                  coefs <- s@coef
                  names <- rownames(coefs)
                  coef <- coefs[, 1]
                  se <- coefs[, 2]
                  pval <- coefs[, 4]
                  n <- model@model@n
                  gof <- numeric()
                  gof.names <- character()
                  gof.decimal <- logical()
                  if (includeSpecTest) {
                      if (spec$k == spec$q)
                          {
                              obj.fcn <- NA
                              obj.pv <- NA
                          } else {
                              obj.fcn <- s@specTest@test[wspecTest,1]
                              obj.pv <- s@specTest@test[wspecTest,3]
                          }
                      gof <- c(gof, obj.fcn, obj.pv)                      
                      gof.names <- c(gof.names, 
                                     paste(specTest,"-test Statistics", sep=""),
                                     paste(specTest,"-test p-value", sep=""))
                      gof.decimal <- c(gof.decimal, TRUE, TRUE)
                  }
                  if (include.nobs == TRUE) {
                      gof <- c(gof, n)
                      gof.names <- c(gof.names, "Num.\\ obs.")
                      gof.decimal <- c(gof.decimal, FALSE)
                  }
                  nbal <- length(model@model@X@balCov)
                  gof.names <- c(gof.names, "Num. Bal. Cov.")
                  gof <- c(gof, nbal)
                  gof.decimal <- c(gof.decimal, FALSE)
                  tr <- createTexreg(coef.names = names, coef = coef, 
                                     se = se, pvalues = pval, 
                                     gof.names = gof.names, gof = gof, 
                                     gof.decimal = gof.decimal)
                  return(tr)
              })

## ----echo=FALSE, message=FALSE, warning=FALSE---------------------------------
library(knitr)
opts_chunk$set(size='footnotesize')
library(causalGel)

## -----------------------------------------------------------------------------
library(causalGel)
data(nsw)
## We express income in thousands for better stability
nsw$re78 <- nsw$re78/1000
nsw$re75 <- nsw$re75/1000

## -----------------------------------------------------------------------------
g <- re78~treat

## -----------------------------------------------------------------------------
balm <- ~age+ed+black+hisp+re75

## -----------------------------------------------------------------------------
ace <- causalModel(g, balm, nsw, momType="ACE")
act <- causalModel(g, balm, nsw, momType="ACT")
aceRT <- causalModel(g, balm, nsw, momType="uncondBal")

## -----------------------------------------------------------------------------
ace

## -----------------------------------------------------------------------------
print(act, printBalCov=TRUE)

## -----------------------------------------------------------------------------
balm2 <- ~age*ed+black+hisp+re75+I(re75^2)
ace2 <- causalModel(g, balm2, nsw, momType="ACE")
print(ace2, printBalCov=TRUE)

## -----------------------------------------------------------------------------
fit1 <- gelFit(ace, gelType="EL") ## EL is the default
print(fit1, model=FALSE, lambda=TRUE)

## -----------------------------------------------------------------------------
vcov(fit1)$vcov_par

## -----------------------------------------------------------------------------
confint(fit1)

## -----------------------------------------------------------------------------
confint(fit1, 2, type="invLR")

## -----------------------------------------------------------------------------
cr <- confint(fit1, 1:2, area=TRUE)
cr

## ----fig.height=4.5-----------------------------------------------------------
plot(cr, col="lightblue", density=20)

## -----------------------------------------------------------------------------
summary(fit1)

## -----------------------------------------------------------------------------
ace <- causalModel(re78~treat, 
                   ~(age+black+ed)*(age+black+ed) + I(age^2) + I(ed^2),
                   data=nsw)
print(ace, TRUE)

## -----------------------------------------------------------------------------
print(ace[-c(3,4)], TRUE)

## -----------------------------------------------------------------------------
print(ace[1:5], TRUE)

## -----------------------------------------------------------------------------
print(ace[,1:100], TRUE)
print(ace[1:3,nsw$re75>0], TRUE)

## -----------------------------------------------------------------------------
fit <- gelFit(ace)
fit2 <- fit[,nsw$age<48]
fit3 <- fit[1:3,1:500]

## ----echo=FALSE, results='asis'-----------------------------------------------
texreg(list(fit,fit2,fit3), digits=4, label="t1")

## -----------------------------------------------------------------------------
checkConv(fit)

## -----------------------------------------------------------------------------
data(nsw)
nsw$re78 <- nsw$re78/1000
nsw$re75 <- nsw$re75/1000
fit1 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                 momType="uncondBal")

## -----------------------------------------------------------------------------
fit2 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                  momType="ACE")
fit3 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                  momType="ACT")
fit4 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                  momType="ACC")

## ----echo=FALSE, results='asis'-----------------------------------------------
texreg(list(fit1,fit2,fit3,fit4), digits=4, label="t2", 
       custom.model.names=c("ACE(rand.)","ACE(non-random)", "ACT", "ACC"),
       caption="Causal Effect for a Training Program", ci.force=TRUE,
       fontsize='footnotesize')

## -----------------------------------------------------------------------------
causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
          momType="uncondBal", cstLHS="causalEffect=1")

## -----------------------------------------------------------------------------
causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
          momType="uncondBal", cstLHS=c("causalEffect=1", "probTreatment=0.5"),
          tControl=list(method="Brent", lower=0, upper=10))

## -----------------------------------------------------------------------------
causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
          momType="uncondBal", cstLHS=2, cstRHS=1)@theta
causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
          momType="uncondBal", cstLHS=2:3, cstRHS=c(1,.5),
          tControl=list(method="Brent", lower=0, upper=10))@theta

## -----------------------------------------------------------------------------
Un_model <- causalModel(re78~treat, ~age+ed+black+hisp+re75, nsw,
                    momType="uncondBal")

## -----------------------------------------------------------------------------
restModel(Un_model, causalEffect~1)
# or restModel(Un_model, "causalEffect=1")
restModel(Un_model, list(causalEffect~1, probTreatment~0.5))
# or restModel(Un_model, c("causalEffect=1", "probTreatment=0.5"))

## -----------------------------------------------------------------------------
fit1 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                 momType="uncondBal")
fit1@lambda[1:2]

## -----------------------------------------------------------------------------
fit2 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                 momType="uncondBal", restrictLam=TRUE)
rbind(coef(fit1), coef(fit2))
rbind(fit1@lambda, fit2@lambda)

## -----------------------------------------------------------------------------
fit1 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                 momType="ACT")
fit2 <- causalGEL(re78~treat, ~age+ed+black+hisp+re75, nsw, gelType="EL",
                 momType="ACT", orthoBases=TRUE)

## ----echo=FALSE, results='asis'-----------------------------------------------
texreg(list(fit1, fit2), label='ortho', digits=6,
       custom.model.names=c("Original","Orthogonal Bases"),
       caption="Comparing estimates with and without the orthogonal bases")

## -----------------------------------------------------------------------------
checkConv(fit2)

## ----extract, echo=TRUE, eval=FALSE-------------------------------------------
#  library(causalGel)
#  library(texreg)
#  setMethod("extract", "causalGelfit",
#            function(model, includeSpecTest=FALSE,
#                     specTest=c("LR","LM","J"), include.nobs=TRUE,
#                     include.obj.fcn=TRUE, ...)
#                {
#                    specTest <- match.arg(specTest)
#                    s <- summary(model, ...)
#                    wspecTest <- grep(specTest, rownames(s@specTest@test))
#                    spec <- modelDims(model@model)
#                    coefs <- s@coef
#                    names <- rownames(coefs)
#                    coef <- coefs[, 1]
#                    se <- coefs[, 2]
#                    pval <- coefs[, 4]
#                    n <- model@model@n
#                    gof <- numeric()
#                    gof.names <- character()
#                    gof.decimal <- logical()
#                    if (includeSpecTest) {
#                        if (spec$k == spec$q)
#                            {
#                                obj.fcn <- NA
#                                obj.pv <- NA
#                            } else {
#                                obj.fcn <- s@specTest@test[wspecTest,1]
#                                obj.pv <- s@specTest@test[wspecTest,3]
#                            }
#                        gof <- c(gof, obj.fcn, obj.pv)
#                        gof.names <- c(gof.names,
#                                       paste(specTest,"-test Statistics", sep=""),
#                                       paste(specTest,"-test p-value", sep=""))
#                        gof.decimal <- c(gof.decimal, TRUE, TRUE)
#                    }
#                    if (include.nobs == TRUE) {
#                        gof <- c(gof, n)
#                        gof.names <- c(gof.names, "Num.\\ obs.")
#                        gof.decimal <- c(gof.decimal, FALSE)
#                    }
#                    nbal <- length(model@model@X@balCov)
#                    gof.names <- c(gof.names, "Num. Bal. Cov.")
#                    gof <- c(gof, nbal)
#                    gof.decimal <- c(gof.decimal, FALSE)
#                    tr <- createTexreg(coef.names = names, coef = coef,
#                                       se = se, pvalues = pval,
#                                       gof.names = gof.names, gof = gof,
#                                       gof.decimal = gof.decimal)
#                    return(tr)
#                })