---
title: "Semiparametric Least Squares Inference for Causal Effects with R"
author: "Pierre Chausse^[University of Waterloo, pchausse@uwaterloo.ca], Mihai Giurcanu^[University of Chicago, giurcanu@uchicago.edu], Marinela Capanu^[Memorial Sloan Kettering Cancer Center, capanum@mskcc.org], George Luta^[Georgetown University, George.luta@georgetown.edu]"
date: ""
bibliography: slse.bib 
output:
 pdf_document:
  number_sections: true
abstract: "This vignette explains how to use the causalSLSE package to estimate causal effects using the semiparametric least squares methods developed by @cslse23. We describe the classes and methods implemented in the package as well as how they can be used to analyze synthetic and real data."
vignette: >
  %\VignetteIndexEntry{Semiparametric Least Squares Inference for Causal Effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteDepends{causalSLSE}
  %\VignettePackage{causalSLSE}
  %\VignetteEncoding{UTF-8}
header-includes:
- \DeclareMathOperator{\Ex}{E}
- \DeclareMathOperator{\var}{var}
- \DeclareMathOperator{\cov}{cov}
- \DeclareMathOperator{\pr}{Pr}
- \newcommand{\tp}[1]{#1^T}
- \newcommand{\reals}{\mathbb{R}}
---


# Introduction \label{sec:intro}

This document presents the `causalSLSE` package describing the
functions implemented in the package. It is intended for users
interested in the details about the methods presented in @cslse23 and
how they are implemented. We first present the theory and then present
the package in the following sections.

```{r, include=FALSE}
def.chunk.hook  <- knitr::knit_hooks$get("chunk")
knitr::knit_hooks$set(chunk = function(x, options) {
  x <- def.chunk.hook(x, options)
  paste0("\n \\", "footnotesize","\n\n", x, "\n\n \\normalsize")
})
knitr::opts_chunk$set(comment="")
```

The general semiparametric additive regression model is 

\begin{equation} \label{cslseTrue}
\begin{split}
Y & =  \beta_0 (1-Z) + \beta_1 Z + 
\sum_{l=1}^q f_{l,0}(X_l)(1-Z) + \sum_{l=1}^q f_{l,1}(X_l)Z + \xi\\
 & \equiv \beta_0 (1-Z) + \beta_1 Z +  f_0(X)(1-Z) + f_1(X)Z + \xi\,, 
 \end{split}
\end{equation} 

where $Y\in \reals$ is the response variable, $Z$ is the treatment
indicator defined as $Z=1$ for the treated and $Z=0$ for the
nontreated, and $X\in \reals^q$ is a $q$-vector of confounders. We 
approximate this model by the following regression model:

\begin{equation} \label{cslseApp}
\begin{split}
Y & =  \beta_0 (1-Z) + \beta_1 Z + 
\sum_{l=1}^q \tp{\psi_{l,0}}U_{l,0}(1-Z) + 
\sum_{l=1}^q \tp{\psi_{l,1}}U_{l,1}Z +  \zeta\\
& \equiv \beta_0 (1-Z) + \beta_1 Z + \tp{\psi_0}U_0(1-Z) + 
\tp{\psi_1}U_1Z + \zeta\,, 
 \end{split}
\end{equation} 

where $U_{l,k}=u_{l,k}(X_l)=(u_{j,l,k}(X_l)~:~ 1 \leq j \leq
p_{l,k})\in\reals^{p_{l,k}}$ is a vector of basis functions
corresponding to the $l^\mathrm{th}$ nonparametric component of the
$k^\mathrm{th}$ group $f_{l,k}(X_l)$, $\psi_{l,k}\in\reals^{p_{l,k}}$
is an unknown vector of regression coefficients,
$U_k=u_k(X)=(u_{l,k}(X_l)~:~ 1 \leq l \leq q)\in\reals^{p_k}$ and
$\psi_k=(\psi_{l,k}~:~ 1 \leq l \leq q)\in\reals^{p_k}$, with
$p_k=\sum_{l=1}^q p_{l,k}$. In this paper, we propose a data-driven
method for selecting the vectors of basis functions $u_0(X)$ and
$u_1(X)$. Note that we allow the number of basis functions ($p_{l,k}$)
to differ across confounders and groups.

Let the following be the regression model estimated by least squares:

\begin{equation}\label{cslseReg}
Y_i = \beta_0 (1-Z_i) + \beta_1 Z_i + \tp{\psi_0}U_{i,0}(1-Z_i) + 
\tp{\psi_1}U_{i,1}Z_i + \zeta_i\mbox{ for } i=1,...,n\,,
\end{equation} 

and $\hat\beta_0$, $\hat\beta_1$, $\hat\psi_0$ and $\hat\psi_1$ be the
least squares estimators of the regression parameters. Then, the
semiparametric least squares estimators (SLSE) of the average causal
effect (ACE), causal effect on the treated (ACT) and causal effect on
the nontreated (ACN) are defined respectively as follows:

\begin{equation} \label{causalEst}
\begin{split}
\mathrm{ACE} & = \hat\beta_1-\hat\beta_0 + \tp{\hat{\psi}_1}\bar{U}_1 -
\tp{\hat{\psi}_0}\bar{U}_0 \\
\mathrm{ACT} & = \hat\beta_1-\hat\beta_0 + \tp{\hat{\psi}_1}\bar{U}_{1,1} -
\tp{\hat{\psi}_0}\bar{U}_{0,1} \\
\mathrm{ACN} & = \hat\beta_1-\hat\beta_0 + \tp{\hat{\psi}_1}\bar{U}_{1,0} -
\tp{\hat{\psi}_0}\bar{U}_{0,0} \,,
\end{split}
\end{equation} 


where $\bar{U}_k = \frac{1}{n}\sum_{i=1}^n U_{i,k}$, $\bar{U}_{k, 1} =
\frac{1}{n_1}\sum_{i=1}^n U_{i, k}Z_i$, $\bar{U}_{k, 0} =
\frac{1}{n_0}\sum_{i=1}^n U_{i, k}(1-Z_i)$, for $k=0,1$, and $n_0$ and
$n_1$ are the sample sizes of the nontreated and treated groups
respectively. As shown by @cslse23, under some regularity conditions
these estimators are consistent and asymptotically normal.

To derive the variance of these causal effect estimators, note that
they can be expressed as a linear combination of the vector of least
squares estimators. Let $\hat{\theta}=\tp{\{\hat\beta_0, \hat\beta_1,
\tp{\hat\psi_0}, \tp{\hat\psi_1}\}}$. Then, the causal effect
estimators can be written as $\tp{\hat{D}_{\mathrm{c}}}\hat{\theta}$
for c=ACE, ACT or ACN, with $\hat{D}_{\mathrm{ACE}}=\tp{\{-1, 1,
\tp{-\bar{U}_0}, \tp{\bar{U}_1}\}}$, $\hat{D}_{\mathrm{ACT}}=\tp{\{-1,
1, \tp{-\bar{U}_{0,1}}, \tp{\bar{U}_{1,1}}\}}$ and
$\hat{D}_{\mathrm{ACN}}=\tp{\{-1, 1, \tp{-\bar{U}_{0,0}},
\tp{\bar{U}_{1,0}}\}}$. Since $\hat{D}_\mathrm{c}$ is random, we need
a first order Taylor expansion to derive the variance of the
estimators. Assuming that the data set is iid and using the asymptotic
properties of least squares estimators, we can show that the variance of
ACE=$\tp{\hat{D}_{\mathrm{ACE}}}\hat{\theta}$ can be consistently
estimated as follows (we can derive a similar expression for the ACT
and ACN):

\begin{equation} \label{causalVar}
\hat{V}_\mathrm{ACE} = \begin{pmatrix} -\hat\beta_0 & \hat\beta_1 & 
\tp{\hat{D}_\mathrm{ACE}} \end{pmatrix}
\begin{pmatrix}
\hat\Sigma_0 & \hat\Sigma_{0,1} & \hat\Sigma_{0,\hat\theta} \\
\hat\Sigma_{1,0} & \hat\Sigma_1 & \hat\Sigma_{1,\hat\theta} \\
\hat\Sigma_{\hat\theta,0} & \hat\Sigma_{\hat\theta,1} & \hat\Sigma_{\hat\theta} 
\end{pmatrix}
\begin{pmatrix} -\hat\beta_0 \\ \hat\beta_1 \\ \hat{D}_\mathrm{ACE} \end{pmatrix}\,,
\end{equation} 

where $\hat\Sigma_k=\widehat{\var(\bar{U}_k)}$,
$\hat\Sigma_{k,l}=\widehat{\cov(\bar{U}_k, \bar{U}_l)}$,
$\hat\Sigma_{k, \hat\theta}=\tp{\hat\Sigma_{\hat\theta,
k}}=\widehat{\cov(\bar{U}_k, \hat\theta)}$, for $k,l=0,1$, and
$\hat\Sigma_{\hat\theta}$ is a consistent estimator of the variance of
$\hat\theta$. We will discuss the choice of the covariance matrix
estimator $\hat\Sigma_{\hat\theta}$ in the next section.

To understand the package, it is important to know how the
$u_{l,k}(X_l)$'s are defined. For clarity, let's write
$U_{l,k}=u_{l,k}(X_l)$ as $U=u(X) = (u_j(X)~:~ 1 \leq j \leq
p)\in\reals^{p}$. We just need to keep in mind that it is different
for the treated and nontreated groups and also for different
confounders. We describe here how to construct the local linear
splines for a given confounder $X$ in a given group. To this end,
let $\{\kappa_{1}, \ldots, \kappa_{p-1}\}$ be a set of $p-1$ knots
strictly inside the support of $X$ satisfying
$\kappa_{1}<\kappa_2<\ldots < \kappa_{p-1}$. In the case of local
linear splines described in the paper, we have:

\begin{equation}\label{basisFct}
\begin{split}
u_1(x) &= xI(x\leq \kappa_{1}) + \kappa_{1}I(x> \kappa_{1}) \\
u_j(x) &= (x -  \kappa_{j-1})I(\kappa_{j-1} \leq x \leq \kappa_{j}) +
(\kappa_{j} - \kappa_{j-1})I(x> \kappa_{j})\,,~~2\leq j \leq p-1\\
u_p(x) &= (x - \kappa_{p-1})I(x > \kappa_{p-1})
\end{split}
\end{equation}

Therefore, if the number of knots is equal to 1, we only have two
local linear splines. Since the knots must be strictly inside the
support of $X$, for any categorical variable with two levels, the
number of knots must be equal to zero. In this case, $u(x)=x$. For
general ordinal variables, the number of knots cannot exceed the
number of levels minus two. The following illustrates local linear
spline functions when the number of knots is equal to 3:

```{r, echo=FALSE, fig.align='center', out.width='50%', fig.width=10, fig.height=6}
knots <- c(2,5,7.5)
maxX <- 11.5
minX <- -2
plot(NULL, xlim=c(minX,maxX+1), ylim=c(minX,knots[2]), axes=FALSE, ylab="u(x)", xlab="x")
axis(1, c("min(X)", expression(k[1]), expression(k[2]), expression(k[3]),
          "max(X)"), at=c(minX, knots, maxX))
axis(2, c("min(X)", 0), at=c(minX, 0))
curve(x*(x<=knots[1])+knots[1]*(x>knots[1]), minX, maxX, add=TRUE,lwd=2)
curve((x-knots[3])*(x>knots[3]), minX, maxX, add=TRUE, col=4, lty=4,lwd=2)
curve((x-knots[1])*(knots[1]<=x)*(x<=knots[2])+
      (knots[2]-knots[1])*(x>knots[2]), minX,
      maxX, add=TRUE, col=2, lty=2,lwd=2)
curve((x-knots[2])*(knots[2]<=x)*(x<=knots[3])+
      (knots[3]-knots[2])*(x>knots[3]), minX,
      maxX, add=TRUE, col=3, lty=3,lwd=2)
legend("topleft", c(expression(u[1](x)), expression(u[2](x)),
                    expression(u[3](x)), expression(u[4](x))),
       col=1:4, lty=1:4, bty='n', lwd=2)
box()
```

Note that for the sample regression, the knots of $X_l$ for group $k$,
$l=1,...,q$, must be strictly inside the sample range of
$(X_{i,l}~:~1\leq i\leq n,~ Z_i=k)\in\reals^{n_k}$, where $n_k$ is
the sample size in group $k$, instead of inside the support of $X_l$.

The following section explains in details how to use the package to
estimate the causal effects using this method, and the last section
summarizes the package by providing a list of all objects and methods.


# The `causalSLSE` package \label{sec:causalSLSE}

## The Semiparametric LSE model

Note that the regression model presented by Equation \eqref{cslseReg}
can expressed as:

\begin{equation}\label{slseTwo} 
\begin{split}
Y_i & = \beta_0 + \tp{\psi_0}U_{i,0}+  \zeta_{i,0} \mbox{ for } i  \mbox{ s.t. } Z_i=0\\
Y_i & = \beta_1 + \tp{\psi_1}U_{i,1}+  \zeta_{i,1} \mbox{ for } i  \mbox{ s.t. } Z_i=1\,.
\end{split}
\end{equation} 

Estimating Equation \eqref{cslseReg} is identical to estimating the
previous two models separately. The latter may even be numerically
more accurate since it avoids many unnecessary operations. Also, as
mentioned in the previous section, the knots and basis functions are
obtained separately for the treated and nontreated. Therefore, we can
see the model from Equation \eqref{cslseReg} as two semiparametric LSE
(SLSE) models, one for the treated and one for the nontreated, and
this is the approach that we take in the package. One benefit of this
approach is to allow an extension to multiple treatment models. For
example, a two treatment model, with two treated groups and one
nontreated group, is like a one treatment model with one more SLSE
model.

Since our causal SLSE model is a collection of SLSE models, we start
by presenting how SLSE models are defined in the package. We ignore
for now that our objective is to estimate causal effects and consider
the following SLSE model:

\begin{equation} \label{slseOne}
\begin{split}
Y & = \beta + \sum_{l=1}^q \tp{\psi_l}U_{l}+  \zeta\\
 & \equiv \beta + \tp{\psi}U+  \zeta\,,
\end{split}
\end{equation}

where $U_l=u_l(X_l)=(u_{j,l}(X_l)~:~ 1 \leq j \leq
p_l)\in\reals^{p_{l}}$, $\psi_{l}\in\reals^{p_{l}}$ is an unknown
vector of regression coefficients, $U=u(X)=(u_{l}(X_l)~:~ 1 \leq l
\leq q)\in\reals^{p}$ and $\psi=(\psi_{l}~:~ 1 \leq l \leq
q)\in\reals^{p}$, with $p=\sum_{l=1}^q p_{l}$. The next section
explains how the knots are determined.

### The starting knots \label{sssec:knots}

The starting knots are automatically generated by the function
`slseKnots`. The following is the list of arguments of the function:

- **form**: A formula with the right-hand side being the list of
  covariates. If a left-hand side is provided, the `slseKnots`
  function will ignore it, because its purpose is only to generate
  the knots.
  
- **data**: A `data.frame` containing all variables included in the
  formula.
  
- **X**: Alternatively, we can input directly the matrix of
  covariates. If a matrix `X` is provided, the arguments `form` and
  `data` are ignored.
  
- **nbasis**: A function that determines the number of basis functions
  as explained in the procedure below. The default is
  `nbasis=function(n) n^0.3`.

- **knots**: This argument is used to set the knots manually. We will
  explain how to use this argument in the next section.

The following is the procedure implemented by the function
`slseKnots`. It explains the procedure for any covariate $X$.

1. The starting number of knots, also equal to the number of basis
   functions minus 1, depends on the type of covariate. Unless it has
   a type that restricts the number of knots, which is explained
   below, it is determined by the argument `nbasis`. This is a
   function of one argument, the same size, and it returns the default
   number of basis functions. This number cannot be smaller than 2 (we
   will see other ways of forcing the number of basis functions to be
   equal to 1 below) and must be an integer. To be more specific, the
   number of basis functions is set to the maximum between 2 and the
   ceiling of what the `nbasis` function returns. For example, if the
   sample size is 500, the default starting number of basis functions
   is 7=`ceiling(500^0.3)`, which implies a starting number of knots
   of 6. It is possible to have a number of knots that does not depend
   on the sample size. All we need is to set the argument `nbasis` to
   a function that returns an integer, e.g., `nbasis=function(n) 4`
   for 4 basis functions or 3 knots.
   
2. Let $(p-1)$ be the number of knots determined in the previous
   step. The default knots are obtained by computing $p+1$ sample
   quantiles of $X$ for equally spaced probabilities from 0 to 1, and
   by dropping the first and last quantiles. For example, if the
   number of knots is 3, then the initial knots are given by quantiles
   for the probabilities 0.25, 0.5 and 0.75.
   
3. We drop any duplicated knots and any knots equal to either the max
   or the min of $X$. If the resulting number of knots is equal to 0,
   the vector of knots is set to `NULL`. When the vector of knots is
   equal to `NULL` for a variable $X$, it means that $u(x)=x$.
   
The last step implies that the number of knots for all categorical
variables with two levels is equal to 0. For nominal variables with a
small number of levels, the number of knots, a subset of the levels,
may be smaller than the ones defined by `nbasis`. For example, when
the number of levels for a nominal variable is 3, the number of knots
cannot exceed 1.

To illustrate how to use the package, we are using the dataset from
@lalonde86. The purpose of this dataset is to estimate the causal effect
of a training program on real income, but we ignore it for the
moment. The dataset is included in the `causalSLSE` package and can be
loaded as follows.

```{r}
library(causalSLSE)
data(nsw)
```

The dependent variable is the real income in 1978 (`re78`) and the
dataset contains the following covariates: the continuous variables
age (`age`), education (`ed`) and the 1975 real income (`re75`), and
the binary variables `black`, `hisp`, `married` and `nodeg`. We start
by considering a model that includes the covariates `age`, `re75`,
`ed`, and `married`. Since we do not need to specify the left-hand
side, we can create the initial knots as follows

```{r}
k <- slseKnots(form = ~ age + re75 + ed + married, data = nsw)
```

The function returns an object of class `slseKnots` and its `print`
method produces a nice display separating confounders with and without
knots. For example, the following are the starting knots:

\newpage

```{r}
print(k)
```

The sample size is equal to `r nrow(nsw)` and the default `nbasis` is
$n^{0.3}$, which implies a default number of starting knots equal to
`r ceiling(nrow(nsw)^.3)-1` = `ceiling(``r nrow(nsw)`$^{0.3}$`)`-1.
This is the number of knots we have for `age`. However, the number of
knots for `ed` is 5 and it is 4 for `re75`. To understand why, the
following shows the 7 default quantiles for `re75` and `ed` (the
`type` argument of the `quantile` function is the same as it is
implemented in the package):

```{r}
p <- seq(0,1,len=9)[c(-1,-9)] # these are the probabilities with 7 knots
quantile(nsw[,'re75'], p, type=1)
quantile(nsw[,'ed'], p, type=1)
```

We can see that the first three quantiles of `re75` are equal to its
minimum, so they are removed. For the `ed` variable, 10 and 11 appear
twice, so one 10 and one 11 must be removed.

Note that each object in the package is S3-class, so the elements can
be accessed using the operator $. For example, we can extract the
knots for `age` as follows:

```{r}
k$age
```

Note that the covariates are listed in the "no knots" section when
their values are set to `NULL`. In the above example, it is the case
of `married` because it is a binary variable. As we can see its list
of knots it set to `NULL`:

```{r}
k$married
```

### Creating a SLSE model \label{sssec:slseModel}

The SLSE model is created by the function `slseModel`. The arguments
of the function are the same as for the `slseKnots` function except
for the argument `X`, which is not needed. The difference is that
`form` must include the left-hand side variable. For example, we can
create a SLSE model using `re78` as dependent variable and the same
covariates used in the previous section as follows:

```{r}
mod1 <- slseModel(form = re78 ~ age + re75 + ed + married, data = nsw)
```

The function returns an object of class `slseModel` and its `print`
method provides a summary of its specification:

\newpage

```{r}
print(mod1)
```

Note that the selection method is set to `Default` when the knots are
selected using the procedure described in the previous section. Also,
the `print` function shows the number of knots for each covariate
inside the parentheses next to its name. In this example, we see that
the initial number of knots for `age`, `re75` and `ed` are
respectively 7, 4 and 5. The function selects the knots using the
`slseKnots` and stores them in the object under `knots`. We can print
them using the \$ operator as follows and compare them with the ones
obtained in the previous section:

```{r}
mod1$knots
```

Note that we can also print the knots by running the command
`print(mod1, which="selKnots")`.

In order to present another example with different types of
covariates, the dataset `simDat4` is included in the package. This is
a simulated dataset which contains special types of covariates. It
helps to further illustrate how the knots are determined. The dataset
contains a continuous variable `X1` with a large proportion of zeros,
the categorical variable `X2` with 3 levels, an ordinal variable `X3`
with 3 levels, and a binary variable `X4`. The levels for `X2` are
{"first","second","third"} and for `X3` the levels are {1,2,3}. 

```{r}
data(simDat4)
mod2 <- slseModel(Y ~ X1 + X2 + X3 + X4, data = simDat4)
print(mod2, which="selKnots")
```

Character-type variables are automatically converted into factors. It
is also possible to define a numerical variable like `X3` as a factor
by using the function `as.factor` in the formula. We see that the 2
binary variables `X2second` and `X2third` are created and `X2first` is
omitted to avoid multicollinearity. For the binary variable `X4`, the
number of knots is set to 0, and for the ordinal variable `X3`, the
number of knots is set to 1 because the min and max values 1 and 3
cannot be selected.


### Selecting the knots manually \label{sssec:slseManKnots}

The user has control over the selection of knots through the argument
`knots`. When the argument is missing (the default), all knots are set
automatically as described above. One way to set the number of knots
to 0 for all variables is to set the argument to `NULL`.

```{r}
slseModel(re78 ~ age + re75 + ed + married, data = nsw, knots = NULL)
```

Notice that the selection method is defined as "User Based" whenever
the knots are provided manually by the user. The other option is to
provide a list of knots. For each variable, we have three options:

- `NA`: The knots are set automatically for this variable only.

- `NULL`: The number of knots is set to 0 for this variable only.

-  A numeric vector: The vector cannot contain missing or duplicated
   values and must be strictly inside the sample range of the
   variable.
   
In the following, we describe all possible formats for the list of
knots.

**Case 1: An unnamed list of length equal to the number of covariates.**

In that case, the knots must be defined in the same order as the order
of the variables implied by the formula. For example, if we want to
set an automatic selection for `age`, no knots for `ed` and the knots
$\{1000, 5000, 10000\}$ for `re75`, we proceed as follows. Note that
setting the value to `NA` or `NULL` has the same effect for the binary
variable `married`.

```{r}
selK <- list(NA, c(1000,5000,10000), NULL, NA)
mod <- slseModel(re78 ~ age + re75 + ed + married, data = nsw,
                 knots = selK)
print(mod, which = "selKnots")
```

**Case 2: A named list of length equal to the number of covariates.**

In that case, the order of the list of variables does not matter. The
`slseModel` function will automatically reorder the variables to match
the order implied by the formula. The names must match perfectly the
variable names generated by R. In the following example, we want to
add the interaction between `ed` and `age`. We want the same set of
knots as in the previous example and no knots for the interaction
term. The name of the interaction depends on how we enter it in the
formula. For example, it is "age:ed" if we enter `age*ed` in the
formula and "ed:age" if we enter `ed*age`. For factors, the names
depend on which binary variable is omitted. Using the above example
with the `simDat4` model, if we interact `X2` and `X4` by adding
`X2*X4` to the formula, the names of the interaction terms are
"X2second:X4" and "X2third:X4". When we are uncertain about the names,
we can print the knots of a model with the default sets of knots. In
the following, we change the order of variables to show that the order
does not matter.

```{r}
selK <- list(married = NA, ed = NULL, 'age:ed' = NULL, re75 = c(1000,5000,10000), age = NA)
model <- slseModel(re78 ~ age * ed + re75 + married, data = nsw, knots = selK)
print(model, which="selKnots")
```

**Case 3: A named list of length strictly less than the number of covariates.**

The names of the selected variables must match perfectly the names
generated by R and the order does not matter. This is particularly
useful when the number of covariates is large. If we consider the
previous example, the knots are set manually only for `age`. By default,
all names not included in the list of knots are set to
`NA`. Therefore, we can create the same model from the previous
example as follows:

```{r}
selK <- list(ed = NULL, 'age:ed' = NULL, re75 = c(1000,5000,10000))
model <- slseModel(re78 ~ age * ed + re75 + married, data = nsw, knots = selK)
print(model, which="selKnots")
```

Note that the previous case offers an easy way of setting the number
of knots to 0 for a subset of the covariates. For example, suppose we
want to add more interaction terms and set the knots to 0 for all of
them. We can proceed as follows.

```{r}
selK <- list('age:ed' = NULL, 'ed:re75' = NULL, 'ed:married' = NULL)
model <- slseModel(re78 ~ age * ed + re75 * ed + married * ed,
                   data = nsw, knots = selK)
```

\newpage

```{r}
model
```

Note also that `slseModel` deals with interaction terms as with any other
variable. For example, `ed:black` is like a continuous variable with a
large proportion of zeros. The following shows the default selected
knots for `ed:black`.

```{r}
model <- slseModel(re78 ~ age + ed * black, data = nsw)       
model$knots[["ed:black"]]
```

We can see that the number of knots is smaller than 7. This is because
`ed:black` has many zeros and the quantiles equal to the minimum value
are removed.


### Methods for `slseModel` objects \label{sssec:slseMethods}

Other methods are registered for `slseModel` objects. For example, we
can estimate `slseModel` objects using the `estSLSE` method and
summarize the results using the `summary` method. The following is an
example using a simpler model:


```{r}
mod2 <- slseModel(form = re78 ~ ed + married, data = nsw)
fit2 <- estSLSE(mod2)
summary(fit2)
```

We can also plot the predicted dependent variable as a function of
education using the `plot` method and add a confidence region:

```{r, fig.align='center', out.width='50%', fig.height=5}
plot(fit2, "ed", interval="confidence", level=0.95)
```

In the next section, we present the `cslseModel` class which
represents the causal-SLSE model of Equation \eqref{cslseReg}. We will
see that it is just a list of `slseModel` objects. Therefore, the
methods registered for `cslseModel` objects are derived from
`slseModel` methods. Since this vignette is about causal-SLSE, we
choose to present these methods through the `cslseModel` object.

## The causal-SLSE model (`cslseModel`) \label{ssec:cslseModel}

The function `cslseModel` returns an object of class `cslseModel` (or
causal-SLSE model). It is a list of `slseModel` objects, one for the
treated and one for the nontreated.  The function has the same
arguments as `slseModel`, plus the argument `groupInd` that specifies
which value of $Z$ is associated with the treated and which one is
associated with the nontreated. The default is `groupInd=c(treated =
1, nontreated = 0)`. It is possible to have other values or even
characters as indicator, but the names must be treated and
nontreated. We will allow more names in future version of the package
once the multiple treatment method is implemented.

The argument `form` must include a formula linking the outcome and the
treatment indicator and a formula listing the confounders, separated
by the operator `|`. In the following example, we see the formula
linking the outcome `re78` and the treatment indicator `treat`, and a
list of confounders:

```{r}
model1 <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw)
```

Its `print` method summarizes the characteristics of the model. It is
like the `slseModel` object, but the information is provided by
treatment group

```{r}
model1
```

The object also contains additional information about the model stored
as attributes:

```{r}
attr(model1, "treatedVar")
attr(model1, "groupInd")
```

This information is needed in order to compute the causal effects. The
object `model1` is a list of 2 elements, `treated` and `nontreated`,
which are `slseModel` objects. Following Section
\ref{sssec:slseModel}, we can therefore access the knots of a specific
group as follows:

```{r}
model1$treated$knots
```

Alternatively, we can use the `print` method for `slseModel` objects
and print the knots of a specific group using the command
`print(model1$treated, which="selKnots")`. To print the list of knots
for both groups, we can print the `cslseModel` object as follows:

```{r}
print(model1, which="selKnots")
```

To understand how to create a `cslseModel` when the treatment
indicator is not binary, consider the dataset `simDat4` that we
described in Section \ref{sssec:slseModel}. The dataset also contains
the variable `treat`, which is a character variable equal to "treat"
when `Z=1` and "notreat" when `Z=0`. We can create a `cslseModel`
object using `treat` instead of `Z` by specifying the value associated
with each group in the argument `groupInd`:

```{r}
model2 <- cslseModel(Y ~ treat | ~ X1 + X2 + X3 + X4, data = simDat4,
                     groupInd = c(treated = "treat", nontreated = "notreat"))
model2
```

If some values of the treatment indicator variable differ from the
values in `groupInd`, the function will return an error message.

### Setting the knots manually \label{sssec:cslseManKnots}

As for SLSE models, we can select the knots using the argument
`knots`. The procedure is the same as in Section
\ref{sssec:slseManKnots} (Cases 1 to 3), but we need to specify the
name of the group associated with the knots. If `knots` is set to
`NULL`, the number of knots is set to 0 for all confounders and all
groups. If we only want the number of knots to be 0 for one group, we
need to specify which group. For example, the number of knots is set
to 0 for the treated only in the following:

```{r}
selK <- list(treated=NULL)
cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw,
           knots = selK)
```

If a group is missing from the argument `knots`, the knots of the
missing group are set automatically. For example, if we want to set
the knots as in Case 1, but only for the nontreated and let
`cslseModel` choose them for the treated, we would proceed as follows:

```{r}
selK <- list(nontreated=list(NA, c(1000,5000,10000), NULL, NA))
model <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw,
                    knots = selK)
print(model, which = "selKnots")
```

### Estimating the model \label{sssec:estSLSE}

Given the set of knots from the model object, the estimation is just a
least squares method applied to the regression model given by:

\[ 
Y = \beta_0 (1-Z) + \beta_1 Z + \tp{\psi_0}U_0(1-Z) + \tp{\psi_1}U_1Z + \zeta\,, 
\] 

where $U_0=u_0(X)$ and $U_1=u_1(X)$ are defined above (which depend on
the knots of the model). The method that estimates the model is
`estSLSE` which has two arguments, but one of them is mainly used
internally by other functions. We present them in case they are
needed. The arguments are:

- `model`: A model created by the function `cslseModel`.

- `selKnots`: It is a list of one or two elements, one for each
  group. Each element is a list of integers to select knots for the
  associated group. For example, suppose we have 2 confounders with 5
  knots each. If we want to estimate the model with only the first
  knot for the first confounder and knots 3 and 5 for the second
  confounder for the treated and all knots for the nontreated, we set
  `selKnots` to `list(treated=list(1L,c(3L, 5L)))`. By default it is
  missing and all the knots from the model are used.
  
We illustrate the use of `estSLSE` with a simple model containing 2
confounders and a maximum of one knot.

```{r}
model <- cslseModel(re78 ~ treat | ~ age + married, data = nsw,
                    nbasis = function(n) 2)
```
  
Note that the `cslseModel` object is a list of `slseModel`
objects. Also, we saw that the above model can be written as two
regression models, one for each group. Therefore, the `estSLSE` method
is simply estimating the `slseModel` objects separately. For example,
we can obtain $\{\hat\beta_1,\hat\psi_1\}$ as follows:

```{r}
estSLSE(model$treated)
```

The output is an object of class `slseFit`. When we apply the method
to `model`, the model is estimated for each group:

```{r}
fit <- estSLSE(model)
fit
```

This is an object of class `cslseFit`, which is a list of
`slseFit`. Like `cslseModel` objects, it also contains the information
about the treatment indicator variable and the value associated with
each group stored as attributes:

```{r}
attr(fit, "treatedVar")
attr(fit, "groupInd")
```

When we print, `fit`, the coefficients are separated by group. The
coefficients for the treated correspond to $\{\hat\beta_1,
\hat\psi_1\}$ and the coefficients for the nontreated correspond to
$\{\hat\beta_0, \hat\psi_0\}$. We can access the estimated SLSE model
for each group using the \$ operator. For example, the following is
the estimated model for the treated:

```{r}
fit$treated
```

A more detailed presentation of the results can be obtained using the
`summary` method. The only arguments of `summary` are the `cslseFit`
object and `vcov.`. The latter is a function that returns the
estimated covariance matrix of the LSE. By default, it is equal to the
`vcovHC` function of the `sandwich` package (@sandwich) with its
default `type="HC3"`. The following is an example with the previous
model using the HC0 type:

```{r}
s <- summary(fit, type="HC0")
```

The object `s` is an object of class `summary.cslseFit`, which is a
list of objects of class `summary.slseFit`, one for each group. By
default, if we print `s`, we will see the two LSE summary tables, one
for each group. Alternatively, we can print the result for one group
using the \$ operator. For example, the following is the result for
the nontreated:

```{r}
s$nontreated
```

We can see that the only knot for `age` in the nontreated group is 23:

```{r}
model$nontreated$knots$age
```

Therefore, the coefficient of `U.age_1` is the effect of age for
the nontreated on `re78` when `age`$\leq 23$ and `U.age_2` is the effect
when `age`$> 23$. 

**The `extract` method**

The package comes with an `extract` method for `slseFit` objects,
which allows to present estimated SLSE models in a LaTeX table using
the `texreg` package of @texreg. There is no `extract` method for
`cslseFit` objects, but we can still use `texreg` by converting
`cslseFit` objects into lists using the method `as.list`. Here is an
example (the argument `table=FALSE` is used to avoid having a floating
table):

```{r, results='asis', message=FALSE}
library(texreg)
texreg(as.list(fit), table=FALSE)
```

### The `predict` method \label{sssec:pred}

The `predict` method is very similar to `predict.lm`. We use the same
arguments: `object`, `interval`, `se.fit`, `newdata` and `level`. The
difference is that it returns the predicted outcome for the treated
and nontreated separately, and the argument `vcov.` provides a way of
changing how the least squares covariance matrix is computed. By
default, it is computed using `vcovHC` from the `sandwich`
package. The function returns a list of 2 elements, `treated` and
`nontreated`. By default (`se.fit=FALSE` and `interval="none"`), each
element contains a vector of predictions. Here is an example with the
previously fitted model `fit`:

```{r}
predict(fit, newdata = data.frame(treat = c(1,1,0,0),age = 20:23, married = 1))
```

If `interval` is set to "confidence", but `se.fit` remains equal to
`FALSE`, each element contains a matrix containing the prediction, and
the lower and upper confidence limits, with the confidence level
determined by the argument `level` (set to 0.95 by default). Here is
an example with the same fitted model:

```{r}
predict(fit,  newdata = data.frame(treat = c(1,1,0,0),age = 20:23, married = 1),
        interval = "confidence")
```

If `se.fit` is set to `TRUE`, each element, treated or nontreated, is
a list with the elements `pr`, containing the predictions, and
`se.fit`, containing the standard errors. In the following, we show
the result for the same fitted model:

```{r}
predict(fit, newdata = data.frame(treat = c(1,1,0,0),age = 20:23, married = 1),
        se.fit = TRUE)
```

### The `plot` method \label{sssec:plot}

The `predict` method is called by the `plot` method to visually represent
the predicted outcome for the treated and nontreated with respect to a
given confounder, controlling for the other variables in the
model. Note that this method is very close to the `plot` method for
`slseFit` objects. In fact, the arguments are the same with some
exceptions that we briefly explain below. Since the predicted outcome
is obtained separately for the treated and nontreated the method for
`cslseFit` objects simply applies the method for `slseFit` objects to
each group. The following is the list of arguments:

- **x**: An object of class `cslseFit` (or `slseFit`).

- **y**: An alias for `which` for compatibility with the generic
  `plot` function. 
  
- **which**: confounder to plot against the outcome variable. It can
  be an integer (the position of the confounder) or a character (the
  name of the confounder)

- **interval**: The type of confidence interval to display. The
  default is "none". The alternative is "confidence".
  
- **level**: The confidence level when `interval="confidence"`. The
  default is 0.95.
  
- **fixedCov**: Optional named lists of fixed values for some or all
	other confounders in each group. The values of the confounders not
	specified are determined by the argument `FUN`. To fix some
	confounders for both groups, `fixedCov` is just a named list with
	the names being the variable names. To fix them to different
	values for the treated and nontreated, `fixedCov` is a named list
	of 1 or 2 elements (for the treated, nontreated or both), each
	element being a named list of values for the covariates.  See the
	examples below. When applied to `slseFit` objects, it is just a
	named list with the variables names.
  
- **vcov.**: An optional function to compute the estimated matrix of
  covariance of the least squares estimators. This argument only
  affects the confidence intervals. The default is `vcovHC` from the
  `sandwich` package with `type="HC3"`.

- **add**: Should the curves be added to an existing plot? The
  default is `FALSE`.

- **addToLegend**: An optional character string to add to the legend
  next to "treated" and "nontreated". Note that a legend is not added
  when applied to `slseFit` objects, so this argument has no effect in
  that case.

- **addPoints**: Should we include the scatterplot of the outcome and
  confounder to the graph? The default is `FALSE`.
  
- **FUN**: A function to determine how the other confounders are
  fixed. The default is `mean`. Note that the function is applied to
  each group separately.

- **plot**: By default, the method produces a graph. Alternatively, we
  can set this argument to `FALSE` and it returns one `data.frame` per
  group with the variable selected by `which` and the prediction. This
  could be useful if one wants to design the graphs differently.

- **graphPar**: A list of graphical parameters if not satisfied with
  the default ones.

- **...**: Other arguments are passed to the `vcov.` function. 

The default set of graphical parameters can be obtained by running the
function `causalSLSE:::.initParCSLSE()` (or
`causalSLSE:::.initParSLSE()` for `slseFit` objects). The function
returns a list of four elements: `treated`, `nontreated`, `common`,
`legend`. The first two are lists of two elements: `points` for the
list of parameters of the scatterplot produced when `addPoints=TRUE`
and `lines` for the line parameters. For example, we can see that the
type of points for the treated is initially set to `pch=21` and their
colour to `2`:

```{r}
causalSLSE:::.initParCSLSE()$treated$points
```

The element `common` is for parameters not specific to a group like
the main title or the axis labels, and `legend` are the parameters that
control the legend (for `cslseFit` only). Note, however, that the
colour and line shapes for the legend are automatically determined by
the lines and points parameters of the `treated` and `nontreated`
elements.

The default parameters can be modified by the argument
`graphPar`. This argument must follow the structure of
`causalSLSE:::.initParCSLSE()` (or `causalSLSE:::.initParSLSE()` for
`slseFit` objects). For example, if we want a new title, new x-axis
label, new type of lines for the treated, new type of points for the
nontreated and a different position for the legend, we create the
following `graphPar`:

```{r}
graphPar <- list(treated = list(lines = list(lty=5, col=4)),
                 nontreated = list(points = list(pch=25, col=3)),
                 common = list(xlab = "MyNewLab", main="My New Title"),
                 legend = list(x = "top"))
```

In the following, we illustrate some examples. 

**Example 1:** 

Consider the model:

```{r}
model1 <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw)
fit1 <- estSLSE(model1)
```

Suppose we want to compare the predicted income between the two
treatment groups with respect to age or education, holding the other
variables fixed to their group means (the default). The following are
two examples with some of the default arguments modified. Note that
`vcov.lm` is used in the first plot function and `vcovHC` (the
default) of type HC1 in the second plot.

```{r, fig.show='hold', out.width='50%'}
library(sandwich)
arg1 <- list(treated = list(lines = list(col = "darkred", lty = 4)),
             nontreated = list(lines = list(col = "darkgreen", lty = 2)),
             legend = list(x = "topleft"))
arg2 <- list(legend = list(x = "top", cex=0.7))
plot(fit1, "ed", vcov. = vcov, graphPar=arg1, interval = 'confidence')
plot(fit1, "age", interval = 'confidence', level = 0.9, type = "HC1", graphPar=arg2)
```

**Example 2:**

If we want to fix the other confounders using another function, we can
change the argument `FUN`. The new function must be a function of one
argument. For example, if we want to fix the other confounders to
their group medians, we set `FUN` to `median` (no quotes). We proceed
the same way for any function that requires only one argument. If the
function requires more than one argument, we have to create a new
function. For example, if we want to fix them to their group 20\%
empirical quantiles, we can set the argument to `function(x)
quantile(x, .20)`. The following illustrates the two cases:

```{r, fig.show='hold', out.width='50%'}
plot(fit1, "age", FUN = median)
plot(fit1, "age", FUN = function(x) quantile(x, 0.20))
```

**Example 3:**

It is also possible to set some of the other confounders to a specific
value by changing the argument `fixedCov`. To fix some variables to
the same values for both groups, `fixedCov` must be a named list with
the names corresponding to the variables you want to fix. You can also
add a description to the legend with the argument `addToLegend`. In
the following `re75` is fixed at 10,000 and we compare the predicted
outcome for the married individuals (the left graph) with the
non-married ones (the right graph)

```{r, fig.show='hold', out.width='50%'}
arg2 <- list(legend = list(cex = 0.8), common=list(ylim=c(4000,9000)))
plot(fit1, "age", fixedCov = list(married = 1, re75 = 10000),
     addToLegend = "married", graphPar = arg2)
plot(fit1, "age", fixedCov = list(married = 0, re75 = 10000),
     addToLegend = "non-married", graphPar = arg2)
```

**Example 4:**

To better compare the two groups, it is also possible to have them
plotted on the same graph by setting the argument `add.` to `TRUE`. We
just need to adjust some of the arguments to better distinguish the
different curves. In the following example, we set the colors and line
shapes to different values and change the position of the legend for
the second set of lines. 

```{r, fig.align='center', out.width='70%'}
arg3 <- list(legend = list(cex = 0.7),
             common = list(ylim = c(3000, 10000)))
plot(fit1, "age", fixedCov = list(married = 1, re75 = 10000),
     addToLegend = "married", graphPar = arg3)
arg4 <- list(treated = list(lines = list(col = "darkred", lty = 5)),
             nontreated = list(lines = list(col = "darkgreen", lty = 4)),
             legend = list(x = "topleft", cex = 0.7))
plot(fit1, "age", fixedCov = list(married = 0, re75 = 10000),
     addToLegend = "non-married", add = TRUE, graphPar = arg4)
```

**Example 5:**

Finally, it is also possible to add the actual observations to the
graph.


```{r, fig.show='hold', out.width='50%'}
arg5 <- list(treated = list(lines = list(col = "darkred", lty = 4)),
             nontreated = list(lines = list(col = "darkgreen", lty = 2)),
             legend = list(x = "topleft"))
plot(fit1, "ed", addPoints = TRUE, graphPar = arg5)
plot(fit1, "re75", addPoints = TRUE)
```

### Factors, interactions and functions of confounders \label{sssec:factor}

The package allows some of the confounders to be factors, functions of
other confounders or interactions. For example, the dataset `simDat4`
includes one factor, `X2`, with levels equal to "first", "second" and
"third". We can include this confounder directly to the list of
confounders. For example,

```{r}
data(simDat4)
mod <- cslseModel(Y ~ Z | ~ X1 + X2 + X4, data = simDat4)
mod
```

We see that R has created 2 binary variables, one for `X2="second"`
and one for `X2="third"`. These two variables are automatically
included in the group of confounders not approximated by SLSE because
they are binary variables like `X4`. If we want to plot `Y` against
`X1`, the binary variables `X2second`, `X2third` and `X4` are fixed to
their group averages which, in case of binary variables, represent the
proportions of ones in each group. 

For interaction terms or functions of confounders, `FUN` is applied to
the functions of confounders. This is how we have to proceed to obtain
the average prediction in regression models. For example, if we
interact `X2` and `X4`, we obtain:

```{r}
data(simDat4)
mod <- cslseModel(Y ~ Z | ~ X1 + X2 * X4, data = simDat4)
mod
```

In this case, when `FUN=mean`, `X2second:X4` is replaced by the
proportion of ones in `X2second`$\times$`X4` for each group. It is not
replaced by the proportion of ones in `X2second` times the proportion
of ones in `X4`. The same applies to functions of confounders. For
functions of confounders, which can be defined in the formula using a
built-in function like `log` or using the identity function `I()`
(e.g. we can interact `X1` and `X4` by using `I(X1*X4)`), `FUN` is
applied to the function (e.g. the average `log(X)` or the average
`I(X1*X4)`).

To fix a factor to a specific level, we just set its value in the
`fixedCov`. In the following example, we fix `X2` to "first", so
`X2second` and `X2third` are set to 0.

```{r, eval=FALSE}
fit <- estSLSE(mod)
plot(fit, "X1", fixedCov = list(X2 = "first"))
```

Note that if a function of confounders (or an interaction) involves the
confounder we want to plot the outcome against, we factorize the
confounder out, apply `FUN` to the remaining of the function and add
the confounder back. For example, if we interact `X1` with `X4` and
`FUN=mean`, `X1:X4` is replaced by `X1` times the proportion of ones
in `X4` for each group.

## Optimal selection of the knots \label{ssec:select}

We have implemented two methods for selecting the knots: the backward
semiparametric LSE (BLSE) and the forward semiparametric LSE (FLSE)
methods. For each method, we have 3 criteria: the p-value threshold
(PVT), the Akaike Information criterion (AIC), and the Bayesian
Information criterion (BIC). Note that the consistency of the causal
effect estimators has only been proved for the last two criteria in
@cslse23. The two selection methods can be summarized as follows:

> We first compute one p-value per knot using either the BLSE or
> FLSE method:
>
> **BLSE**: 
>
> 1. We estimate the model with all knots included in the model.
>
> 2. For each knot, we test if the slopes of the basis functions
>    adjacent to the knot are the same, and return the p-value.
>
> 
> **FLSE**: 
>
> 1. We estimate the model by including a subset of the knots, one
>    variable at the time. When we test a knot for one confounder, the
>    number of knots is set to 0 for all other variables. 
>
> 2. For each knot, we test if the adjacent slopes to the knot are the
>    same, and return the p-value. The set of knots used for each test
>    depends on the following:
>
>    - Variables with 1 knot: we return the p-value of the test of
>      equality of the slopes adjacent to the knot.
>
>    - Variables with 2 knots: we include the two knots and return the
>      p-values of the test of equality of the slopes adjacent to each knot.
>
>    - Variables with $p$ knots ($p>2$): We test the equality of the
>      slopes adjacent to knot $i$, for $i=1,...,p$, using the sets of
>      knots $\{1,2\}$, $\{1,2,3\}$, $\{2,3,4\}$, ..., $\{p-2,p-1,p\}$
>      and $\{p-1,p\}$ respectively.
>
> Once we have the p-values, we proceed to step 3:
>
> 3. The knots are selected using one of the following criteria:
>
>    - **PVT**: We remove all knots with a p-value greater than a
>      specified threshold.
>
>    - **AIC** or **BIC**: We order the p-values in ascending
>      order. Then, starting with a model with no knots and going from
>      the smallest to the highest p-value, we add the knot associated
>      with the smallest remaining p-value one by one, estimate the
>      model and return the information criterion. We select the model
>      with the smallest value of the information citerion. 
>
> Note that the SLSE models for the treatment groups contained in
> `cslseModel` objects are estimated separately. However, the AIC and
> BIC are computed as if the they were estimated jointly as in
> Equation \eqref{cslseApp}. As we will see below, a joint selection
> does not necessarily correspond to a selection done group by group.

The knot selection is done using the `selSLSE` method. The arguments are:

- **model**: An object of class `cslseModel`.  The method also exists
  for `slseModel` objects. We will discuss it briefly at the end of
  this section.

- **selType**: This is the selection method. We have the choice between
  "FLSE" and "BLSE" (the default).

- **selCrit**: This is the criterion used by the selection method. We
  have the choice between "AIC" (the default), "BIC" or "PVT".
  
- **pvalT**: This is a function that returns the p-value threshold. It
  is a function of one argument, the average number of basis functions
  per confounder. The default is `function(p) 1/log(p)` and it is
  applied to each group separately. Therefore, the threshold may be
  different for the treated and non-treated. It is also possible to
  set it to a fixed threshold. For example, `function(p) 0.20` sets the
  threshold to 0.2. Note that when the function returns a value
  greater than 1, all knots are kept. This argument affects the result
  only when `selCrit` is set to "PVT".
  
- **vcovType**: The type of LSE covariance matrix used to compute the
  p-values. The options are "HC0" (the default), "HC1", "HC2", "HC3"
  and "Classical" (for the homoskedastic case). Using a
  heteroskedasticity robust covariance matrix is recommended, but
  there is not need to choose the usually recommended HC3 for the
  selection. The reason is that HC3 requires the hat values, which
  slows down the process, especially for FLSE, and simulations from
  @cslse23 show that the choice of HC has very little effect on the
  selection. Note that the model is estimated separately for the
  treated and nontreated. Therefore, we assume a different variance of
  the residuals even when `vcovType` is set to "Classical".
  
- **reSelect**: Should we recompute the optimal knots or use the ones
  already saved in the object. See below for more details.
  
The function returns a model of class `cslseModel` with the optimal
selection of knots. For example, we can compare the starting knots of
the following model, with the ones selected by the default arguments.

```{r}
model1 <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw)
model1
model2 <- selSLSE(model1)
model2
```

For example, the BLSE-AIC method has kept all knots from `re75`
for the treated and kept two knots for the nontreated. The print
method indicates which method was used to select the knots. It is
possible to recover the p-values of all original knots by setting the
argument `which` to `Pvalue`.

```{r}
print(model2, which="Pvalues")
```

In the following example, we use BLSE as selection method and BIC as
criterion. Note that the BIC selects 0 knots for all confounders.

```{r}
model3 <- selSLSE(model1, selType = "BLSE", selCrit = "BIC")
model3
```

Since the `selSLSE` method returns a new model, we can apply the
`estSLSE` to it:

```{r}
estSLSE(selSLSE(model1, selType = "FLSE", selCrit = "BIC"))
```

### Selection for `slseModel` versus `cslseModel` objects

As mentioned in the previous section, the information criteria for
`cslseModel` objects are computed as if the SLSE models of the
treatment groups were estimated using Equation \eqref{cslseApp}. This
approach may lead to a selection different from what we would obtain
by selecting the knots group by group. To see this, consider the
following model and joint selection based on BLSE-AIC:

```{r}
model1 <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw)
model2 <- selSLSE(model1, selType="BLSE", selCrit="AIC")
model2
```

We could also apply the same selection method, but group by group. In
the following, we show how to select the knots group by group by using
the `selSLSE` method for `slseModel` objects. Note that the Selection
method is set to `BLSE-AIC (Sep.)` to indicate that it was done
separately.

```{r}
model3 <- model1
model3$treated <- selSLSE(model3$treated, selType="BLSE", selCrit="AIC")
model3$nontreated <- selSLSE(model3$nontreated, selType="BLSE", selCrit="AIC")
model3
```

We can see that the selected knots are quite different. For example,
the number of knots of `age` for the treated is 2 when selected
jointly and 3 when selected separately. Also, the number of knots of
all confounders for the nontreated are set to 0 when the selection is
done separately. This is very different from the joint
selection. Since the joint selection is the one developed and studied by
@cslse23, it is the one we recommend.

Note that the PVT approach leads to identical selection whether it is
done separately or jointly. The reason is that both approaches produce
identical p-values, and the thresholds depend on the number of knots
in each group.

### Selections saved in `slseModel` objects.

Optimal selection of knots can be time consuming, especially for large
sample sizes. To avoid having to recompute the selection each time we
change the selection method and/or the criterion, every new selection
is saved in the `slseModel` object. Let's consider the following
model:

```{r}
model <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data = nsw)
model
```

Suppose we select the knots using BLSE and AIC. We can replace the
object by the new one.

```{r}
model <- selSLSE(model, selType="BLSE", selCrit="AIC")
```

The main source of computing time for the selection comes from the
number of regressions we need to estimate. Once estimated, all tests
are simple operations, so the proportion of time allocated to testing
the knots is negligible. Once we have the p-values, no additional
regressions are needed for the `PVT` criterion, but we need one
regression for each p-value for the AIC and BIC (plus 2 if we count
the model with no knots). It is therefore worth saving the selection
somewhere. This information is saved separately in each `slseModel`
object under `selections`.

```{r}
names(model$treated$selections)
names(model$treated$selections)
```

We see that the model keeps the original knots, so we loose no
information by replacing the model object with the new one. The other
element of the list is `BLSE` which contains information about the
selection. If we want to compare `BLSE` with `FLSE`, we can call the
`selSLSE` once more and replace `model` with the new one:

```{r}
model <- selSLSE(model, selType="FLSE", selCrit="AIC")
names(model$treated$selections)
names(model$treated$selections)
```

The new selections are added to the model without deleting the
previous ones. The following is what we can find in the element `BLSE`
(or `FLSE`) of a given group:

```{r}
names(model$treated$selections$BLSE)
```

The `pval` element is an object of class `pvalSLSE`. If we print it,
we obtain the same output as when we print the model with the argument
`which="Pvalues"`. The element `Threshold` is the p-value threshold
used for the `PVT` criterion, `JIC` is a matrix of BIC and AIC values,
ordered from the smallest to the highest p-value, and the selections
are saved in the elements `PVT`, `JAIC` and `JBIC`. The `J` in front
of IC, AIC and BIC means that the selection is based on the joint
estimation of the SLSE models. There is no `J` in front of `PVT`
because the selection using this criterion is identical if we proceed
jointly or separately. Note that we do not see the `J` when we print
the object, because it is assumed that AIC and BIC are obtained
jointly for `cslseModel` objects. If we select the knots of one of the
SLSE model separately, we will see the difference when we print the
model. For example, the following replace the SLSE model of the
treated with the SLSE model selected by AIC (separately):

```{r}
model$treated <- selSLSE(model$treated, "BLSE", "AIC")
model
```

The print output indicates that the selection was done jointly for the
nontreated (just `AIC`) and separately for the treated. Now, the SLSE
model for the treated has three more selection elements: `IC`, `AIC` and
`BIC`:

```{r}
names(model$treated$selections$BLSE)
```

With all selections saved in the model, how do we know which knots are
used? We know it by printing the model or by printing the following
knots attribute:

```{r}
attr(model$nontreated$knots, "curSel")
```

It tells us that the current set of knots for that model is BLSE with
a joint AIC. 

The elements `PVT`, `JAIC`, `JBIC`, `AIC` and `BIC` are lists of knots
selection vectors in the same format as the argument `selKnots` of the
`estSLSE` method. For example, the following is the list of selection
vectors using BLSE with AIC for the treated:

```{r}
model$treated$selections$BLSE$JAIC
```

For example, we see that the knots 1 and 3 of `age` are selected by
AIC. We could use this selection list to estimate the model:

```{r}
estSLSE(model$treated, model$treated$selections$BLSE$JAIC)
```

Note that we never explicitly selected the BIC criterion above, but it
was added to the model anyway. The reason is that the cost of
computing the selection based on BIC is negligible once we do it for
AIC. Since it is also negligible to add the selection by `PVT`, if the
`selCrit` argument is set to `"AIC"` or `"BIC"`, all three selections
are computed simultaneously and stored in the model. It is only when
`selCrit` is set to `"PVT"` that the selection is done for this
criterion only.

### Select a saved selection with `update`

Now that we understand that all previous selections are saved in
`slseModel` objects, how do we select them? This is done with the
`update` method. The method is registered for `slseModel` and
`cslseModel` objects and works very similarly. The arguments are
`selType`, `selCrit` and `selKnots`. The latter is used to select the
knots manually as explained in Section \ref{sssec:estSLSE}. The first
two are as in the `selSLSE` method. The purpose of `update` is to
replace the current knots by a selection saved in the model. If the
selection we want is not in the model, `update` will return an error
message. For example, the object `model` from the previous section has
all selections saved, but the current is a mixture of FLSE-AIC and
BLSE-AIC (Sep.).

```{r}
model
```

We can update the model to the FLSE with joint AIC selection as
follows (AIC means joint AIC for `cslseModel`)

```{r}
model <- update(model, selType="FLSE", selCrit="AIC")
model
```

This is done without any computation. The `update` method simply
replaces the knots using the appropriate list of selection vectors. If
we want to recover the model with the initial set of knots, we simply
set the argument `selType` to `"None"`. In the following, we see that
the selection method is back to its initial set of knots:

```{r}
update(model, selType="None")
```

As a last application of the `update` method, the following shows that
the joint AIC and the one applied to each group separately produces
different results. Since we just set the `model` object to FLSE-AIC,
the current treated model in `model` is based on the joint AIC
criterion:

```{r}
model$treated
```

Since we have the AIC computed separately in the treated model, we can
use `update` for `slseModel` to compare the two:

```{r}
update(model$treated, "BLSE", "AIC")
```

The selection for `age` and `ed` are the same, but the AIC selects 3
knots for `re75` and the JAIC selects only 1 knot. 

Note that the `selSLSE` method computes the selection only if the
requested selection is not saved in the model, unless the argument
`reSelect` is set to `TRUE`. Therefore, `selSLSE` is not different
from `update` when the requested selection is saved in the model.

## The `causalSLSE` method for `cslseFit` objects \label{ssec:causalForFit}

The method `causalSLSE` estimates the causal effects from `cslseFit`
objects using the knots included in the estimated model. The arguments
of the method are:

- **object**: An object of class `cslseFit`.
 
- **causal**: What causal effect measure should the function compute? We
  have the choice between "ALL" (the default), "ACE", "ACT" or "ACN".
   
- **vcov.**: An alternative function used to compute the covariance
  matrix of the least squares estimates. This is the
  $\hat\Sigma_{\hat\theta}$ defined in the Introduction section. By
  default, `vcovHC` is used with `type="HC3"`. Simulations from
  @cslse23 show that using `vcovHC` with `type="HC3"` produces the
  most accurate estimate of the variance of ACE, ACT and ACN in small
  and large samples.
  
- **...**: This is used to pass arguments to the `vcov.` function.

In the following example, we estimate the causal effect with the
initial knots (without selection).

```{r}
model1 <- cslseModel(re78 ~ treat | ~ age + re75 + ed + married, data=nsw)
fit1 <- estSLSE(model1)
ce <- causalSLSE(fit1)
ce
```

The method returns an object of class `cslse` and its `print` method
only prints the causal effect estimates. We can extract any causal
estimate and its standard error by using the \$ operator followed by
the type of causal effect estimate. For example, the following is the
ACE:

```{r}
ce$ACE
```

The object also contains the two `slseFit` objects and the model
attributes that specify the name of the treatment indicator variable
and the value of the indicator associated with each group. 

For more details about the estimation, which includes standard errors
and significance tests, we can use the `summary` method:

```{r}
sce <- summary(ce)
sce
```

The `summary` method returns an object of class `summary.cslse` and
the above output is produced by its `print` method. If needed, we can
extract the above table using `$causal`. The summary tables for the
treated and nontreated LSE can be extracted using `$lse`.

The `cslse` object inherits from the class `cslseFit`, so we can apply
the `plot` (or the `predict`) method directly on this object as shown below:

```{r, fig.align='center', out.width='65%'}
plot(ce, "re75")
```

### The `extract` method  \label{sssec:extract}

The package also comes with an `extract` method for objects of class
`cslse`. For example, we can compare different methods in a single
table. In the following example, we compare the SLSE, BLSE-AIC and
FLSE-AIC:

```{r, eval=FALSE}
library(texreg)
c1 <- causalSLSE(fit1)
fit2 <- estSLSE(selSLSE(model1, selType="BLSE"))
fit3 <- estSLSE(selSLSE(model1, selType="FLSE"))
c2 <- causalSLSE(fit2)
c3 <- causalSLSE(fit3)
texreg(list(SLSEC=c1, BLSE_AIC=c2, FLSE_AIC=c3), table=FALSE, digits=4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
library(texreg)
c1 <- causalSLSE(fit1)
fit2 <- estSLSE(selSLSE(model1, selType="BLSE"))
fit3 <- estSLSE(selSLSE(model1, selType="FLSE"))
c2 <- causalSLSE(fit2)
c3 <- causalSLSE(fit3)
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=4)
```
\end{center}

The arguments of the `extract` methods, which control what is printed
and can be modified through the `texreg` function, are:

- **include.nobs**: Include the number of observations. The
  default is `TRUE`.
  
- **include.nknots**: Include the number of knots.  The default is
  `TRUE`.

- **include.rsquared**: Include the $R^2$.  The default is `TRUE`.
  
- **include.adjrs**: Include the adjusted $R^2$.  The default is
  `TRUE`.
  
- **separated.rsquared** Should we return one $R^2$ per `slseModel`?
  By default it is set to `FALSE` and one $R^2$ is computed for the
  joint estimation of Equation \eqref{cslseReg}. This argument applies
  also to the adjusted $R^2$.

- **which**: Which causal effects should be printed? The options are
  "ALL" (the default), "ACE", "ACT", "ACN", "ACE-ACT", "ACE-ACN" or
  "ACT-ACN".

Here is one example on how to change some arguments:

```{r, eval=FALSE}
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=3,
       which="ACE-ACT", include.adjrs=FALSE, separated.rsquared=TRUE)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=3, 
       which="ACE-ACT", include.adjrs=FALSE, separated.rsquared=TRUE)
```
\end{center}

## The `causalSLSE` method for `cslseModel` objects  \label{ssec:causalForMod}

When applied directly to `cslseModel` objects, the `causalSLSE` method
offers the possibility to select the knots and estimate the causal
effects all at once. The method also returns an object of class
`cslse`. The arguments are the same as for the method for `cslseFit`
objects, plus the necessary arguments for the knots selection. The
following are the arguments not already defined for objects of class
`cslseFit`. The details of these arguments are presented in Section
\ref{ssec:select}.

- **object**: An object of class `cslseModel`.

- **selType**: This is the selection method. We have the choice
  between "SLSE" (the default), "FLSE" and "BLSE". The SLSE method
  performs no selection, so all knots from the model are kept. 

- **selCrit**: This is the criterion used by the selection method when
  `selType` is set to "FLSE" or "BLSE". The default is "AIC".
  
- **pvalT**: This is a function that returns the p-value threshold. We
  explained this argument when we presented the `selSLSE` method.
  
For example, we can generate the previous table as follows:

```{r, eval=FALSE}
c1 <- causalSLSE(model1, selType="SLSE")
c2 <- causalSLSE(model1, selType="BLSE")
c3 <- causalSLSE(model1, selType="FLSE")
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(model1, selType="SLSE")
c2 <- causalSLSE(model1, selType="BLSE")
c3 <- causalSLSE(model1, selType="FLSE")
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=4)
```
\end{center}

Note that this `causalSLSE` method calls `selSLSE` for the knot
selection. Therefore, the selection is done without any computation if
the selection is already saved in the model object. It would therefore
be inefficient to compare FLSE with AIC, BIC and PVT using the
following, because it would involve recomputing all FLSE selections
3 times:

```{r, eval=FALSE}
c1 <- causalSLSE(model1, selType="FLSE", selCrit="AIC")
c2 <- causalSLSE(model1, selType="FLSE", selCrit="BIC")
c3 <- causalSLSE(model1, selType="FLSE", selCrit="PVT")
```

It is better to add all `FLSE` selections in the model first and then
call the `causalSLSE` method 3 times as follows:

```{r}
model1 <- selSLSE(model1, selType="FLSE")
c1 <- causalSLSE(model1, selType="FLSE", selCrit="AIC")
c2 <- causalSLSE(model1, selType="FLSE", selCrit="BIC")
c3 <- causalSLSE(model1, selType="FLSE", selCrit="PVT")
```


## The `causalSLSE` method for `formula` objects  \label{ssec:causalForForm}

This last method, offers an alternative way of estimating the causal
effects. It allows the estimation in one step without having to first
create a model. The arguments are the same as for the `cslseModel`
function and the `causalSLSE` method for `cslseModel` objects. It
creates the model, selects the knots and estimates the causal effects
in one step. For example, we can create the previous table as follows:

```{r, eval=FALSE}
c1 <- causalSLSE(re78 ~ treat | ~ age + re75 + ed + married, data=nsw,
                 selType="SLSE")
c2 <- causalSLSE(re78 ~ treat | ~ age + re75 + ed + married, data=nsw,
                 selType="BLSE")
c3 <- causalSLSE(re78 ~ treat | ~ age + re75 + ed + married, data=nsw,
                 selType="FLSE")
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table=FALSE, digits=4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(model1, selType = "SLSE")
c2 <- causalSLSE(model1, selType = "BLSE")
c3 <- causalSLSE(model1, selType = "FLSE")
texreg(list(SLSE=c1, BLSE=c2, FLSE=c3), table = FALSE, digits = 4)
```
\end{center}

Note that this method calls `cslseModel`, `selSLSE`, `estSLSE` and the
method `causalSLSE` for `cslseFit` objects sequentially. It is easier
to simply work with this method, but manually going through all steps
may be beneficial to better understand the procedure. Also, it is more
convenient to work with a model when we want to compare the different
selection methods, or if we want to compare estimations with different
types of standard errors. In particular, this approach does not offer
the more efficient option of computing all selections once and save
them in the model as explained at the end of the previous section.

# Examples \label{sec:example}


## A simulated data set from Model 1  \label{ssec:causalEx1}

In the package, the data set `datSim1` is generated using the
following data generating process with a sample size of 300.

\begin{eqnarray*}
Y(0) &=& 1+X+X^2+\epsilon(0)\\
Y(1) &=& 1-2X+\epsilon(1)\\
Z &=& \mathrm{Bernoulli}[\Lambda(1+X)]\\
Y &=& Y(1)Z + Y(0) (1-Z)\,
\end{eqnarray*}

where $Y(0)$ and $Y(1)$ are the potential outcomes, $X$,
$\epsilon(0)$ and $\epsilon(1)$ are independent standard normal and
$\Lambda(x)$ is the CDF of the standard logistic distribution. The
causal effects ACE, ACT and ACN are approximately equal to -1, -1.6903
and 0.5867 (estimated using a sample size of $10^7$). We can start by
building the starting model:

```{r}
data(simDat1)
mod <- cslseModel(Y ~ Z | ~ X, data = simDat1)
mod <- selSLSE(mod, "BLSE") ## Let's save them all first
mod <- selSLSE(mod, "FLSE")
```

Then we can compare three different methods:

```{r, eval=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "BLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "FLSE", selCrit = "BIC")
texreg(list(SLSE = c1, BLSE = c2, FLSE = c3), table = FALSE, digits = 4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "BLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "FLSE", selCrit = "BIC")
texreg(list(SLSE = c1, BLSE = c2, FLSE = c3), table = FALSE, digits = 4)
```
\end{center}

We see that both selection methods choose to assign 0 knots for the
treated group, which is not surprising since the true $f_1(x)$ is
linear. We can compare the different fits. 

```{r, out.width='50%', fig.show='hold', fig.height=4}
list(common = list(main = "Y vs X using BLSE-BIC"))
plot(c1, "X")
curve(1 - 2 * x, -3, 3, col = "darkgreen", lty = 4, lwd = 3, add = TRUE)
curve(1 + x + x^2, -3, 3, col = "darkorange", lty = 4, lwd = 3, add = TRUE)
legend("bottomleft", c("True-treated", "True-nontreated"),
       col=c("darkgreen", "darkorange"), lty = 4, lwd = 3, bty = 'n')
plot(c2, "X", graphPar = list(common = list(main = "Y vs X using BLSE-BIC")))
curve(1 - 2 * x, -3, 3, col="darkgreen", lty = 4, lwd = 3, add = TRUE)
curve(1 + x + x^2, -3, 3, col = "darkorange", lty = 4, lwd = 3, add = TRUE)
legend("bottomleft", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 4, lwd = 3, bty = 'n')
plot(c3, "X", graphPar = list(common = list(main = "Y vs X using FLSE-BIC")))
curve(1 - 2 * x, -3, 3, col="darkgreen", lty = 4, lwd = 3, add = TRUE)
curve(1 + x + x^2, -3, 3, col = "darkorange", lty = 4, lwd = 3, add = TRUE)
legend("bottomleft", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 4, lwd = 3, bty = 'n')
```

We see that the piecewise polynomials are very close to the true
$f_0(x)$ and $f_1(x)$ for SLSE, BLSE and FLSE. We can see from the
folllowing graph how the lines are fit through the observations by
group.

```{r, out.width='60%', fig.align='center'}
plot(c1, "X", addPoints=TRUE)
```

## A simulated data set from Model 2 \label{ssec:causalEx2}

The dataset `datSim2` is a change point regression model (with unknown
location of change points) defined as follows:

\begin{eqnarray*}
Y(0) &=& (1+X)I(X\leq-1) + (-1-X)I(X>-1) + \epsilon(0)\\
Y(1) &=& (1-2X)I(X\leq 0) + (1+2X)I(X>0) + \epsilon(1)\\
Z &=& \mathrm{Bernoulli}[\Lambda(1+X)]\\
Y &=& Y(1)Z + Y(0) (1-Z)\,
\end{eqnarray*}

where $Y(0)$ and $Y(1)$ are the potential outcomes, $I(A)$ is the
indicator function equal to 1 if $A$ is true, and $X$, $\epsilon(0)$
and $\epsilon(1)$ are independent standard normal. The causal effects
ACE, ACT and ACN are approximately equal to 3.763, 3.858 and 3.545
(estimated with a sample size of $10^7$). We can compare the SLSE,
BLSE-AIC and BLSE-BIC.

```{r}
data(simDat2)
mod <- cslseModel(Y~Z | ~X, data=simDat2)
mod <- selSLSE(mod, "BLSE") ## We just add BLSE because we do not use FLSE
```

```{r, eval=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "BLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "BLSE", selCrit = "AIC")
texreg(list(SLSE = c1, BLSE.BIC = c2, BLSE.AIC = c3), table = FALSE, digits = 4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "BLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "BLSE", selCrit = "AIC")
texreg(list(SLSE = c1, BLSE.BIC = c2, BLSE.AIC = c3), table = FALSE, digits = 4)
```
\end{center}

The following shows the fit of BLSE-AIC with the true $f_1(x)$ and
$f_0(x)$, and the observations.

```{r, out.width='50%', fig.show='hold'}
arg <-  list(common = list(main = "Y vs X using BLSE-AIC"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X", graphPar = arg)
curve((1 -2 * x) * (x <= 0) + (1 + 2 * x) * (x > 0), -3, 3,
      col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve((1 + x) * (x <= -1) + (-1 - x) * (x > -1),
      -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("left", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 3, lwd = 3, bty = 'n', cex = .8)
arg$legend$x <- "topleft"
plot(c2, "X", addPoints = TRUE, graphPar = arg)
```

## A simulated data set from Model 3 \label{ssec:causalEx3}

The data set `datSim3` is generated from model with multiple
confounders defined as follows:

\begin{eqnarray*}
Y(0) &=& [1+X_1+X_1^2] + [(1+X_2)I(X_2\leq-1) + (-1-X_2)I(X_2>-1)] + 
\epsilon(0)\\
Y(1) &=& [1-2X_1] + [(1-2X_2)I(X_2\leq 0) + (1+2X_2)I(X_2>0)] + 
\epsilon(1)\\
Z &=& \mathrm{Bernoulli}[\Lambda(1+X_1+X_2)]\\
Y &=& Y(1)Z + Y(0) (1-Z)\,,
\end{eqnarray*}

where $Y(0)$ and $Y(1)$ are the potential outcomes, $X_1$, $X_2$,
$\epsilon(0)$ and $\epsilon(1)$ are independent standard normal. The
causal effects ACE, ACT and ACN are approximately equal to 2.762,
2.204 and 3.922 (estimated with a sample size of $10^7$). We can
compare the SLSE, FLSE with AIC and FLSE with BIC.

```{r}
data(simDat3)
mod <- cslseModel(Y ~ Z | ~ X1 + X2, data = simDat3)
mod <- selSLSE(mod, "FLSE") ## We just add FLSE because we do not use BLSE
```

```{r, eval=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "FLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "FLSE", selCrit = "AIC")
texreg(list(SLSE = c1, FLSE.BIC = c2, FLSE.AIC = c3), table = FALSE, digits = 4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "FLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "FLSE", selCrit = "AIC")
texreg(list(SLSE = c1, FLSE.BIC = c2, FLSE.AIC = c3), table = FALSE, digits = 4)
```
\end{center}

To illustrate the method, since we have two confounders, we need to
plot the outcome against one confounder holding the other fixed. The
default is to fix it to its sample mean. For the true curve, we fix it
to its population mean, which is 0. We first look at the outcome
against $X_1$. By fixing $X_2$ to 0, the true curve is $X_1+X_1^2$ for
the untreated and $2-2X_1$ for the treated. The following graphs show how
the FLSE-BIC method fits the curves.

```{r, out.width='50%', fig.show='hold'}
arg <-  list(common = list(main = "Y vs X1 using FLSE-AIC"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X1", graphPar = arg)
curve(x + x^2, -3, 3, col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(2 - 2 * x, -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("topleft", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 3, lwd = 3, bty = 'n', cex = .8)
arg$legend$x <- "topleft"
plot(c2, "X1", addPoints = TRUE, graphPar = arg)
```

If we fix $X_1$ to 0, the true curve is $1+[(1+X_2)I(X_2\leq-1) +
(-1-X_2)I(X_2>-1)]$ for the nontreated and $1+[(1-2X_2)I(X_2\leq 0) +
(1+2X_2)I(X_2>0)]$ for the treated. The following graphs illustrates
how these curves are approximated by FLSE-AIC.

```{r, out.width='50%', fig.show='hold'}
arg <-  list(common = list(main = "Y vs X2 using FLSE-AIC"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X2", graphPar = arg)
curve(1 + (1 - 2 * x) * (x <= 0) + (1 + 2 * x) * (x > 0), -3, 3,
      col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(1 + (1 + x) * (x <= -1) + (-1 - x) * (x > -1),
      -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("left", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 3, lwd = 3, bty = 'n', cex = .8)
arg$legend$x <- "topleft"
plot(c2, "X2", addPoints = TRUE, graphPar = arg)
```

## A simulated data set with product terms \label{ssec:causalEx4}


The data set `datSim5` is generated using the following data
generating process with a sample size of 300.

\begin{eqnarray*}
Y(0) &=& [1+X_1+X_1^2] + [(1+X_2)I(X_2\leq-1) + (-1-X_2)I(X_2>-1)] \\
&& + [1+X_1 X_2 + (X_1X_2)^2] + \epsilon(0)\\
Y(1) &=& [1-2X_1] + [(1-2X_2)I(X_2\leq 0) + (1+2X_2)I(X_2>0)] \\
&&+ [1-2X_1X_2] + \epsilon(1)\\
Z &=& \mathrm{Bernoulli}[\Lambda(1+X_1+X_2+X_1X_2)]\\
Y &=& Y(1)Z + Y(0) (1-Z)\,,
\end{eqnarray*}

where $Y(0)$ and $Y(1)$ are the potential outcomes, and $X_1$, $X_2$,
$\epsilon(0)$ and $\epsilon(1)$ are independent standard normal. The
causal effects ACE, ACT and ACN are approximately equal to 1.763,
0.998 and 3.194 (estimated with a sample size of $10^7$). We can
compare the SLSE, FLSE-AIC and FLSE-BIC.

```{r, echo=FALSE}
sim_model <- function(n)
{
	X1 <- rnorm(n)
    X2 <- rnorm(n)
	probX <- plogis(1 + X1 + X2 + X1*X2) 
	Z <- rbinom(n, 1, probX)
	er0 <- rnorm(n)
	er1 <- rnorm(n)
	Y0 <-  (1 + X1 + X1^2) + (1 + X2)  * (X2 <= -1) + 
		(-1 - X2) * (X2 > -1) + (1+X1*X2 + X1^2*X2^2) + er0
	Y1 <- (1 - 2*X1) + (1 - 2*X2)*(X2 <= 0) + 
		(1 + 2*X2) * (X2 > 0) + (1-2*X1*X2) + er1
	Y <- Y0 * (1 - Z) + Y1 * Z
	data.frame(Y=Y, X1=X1, X2=X2, Z=Z, Y1=Y1, Y0=Y0)
}
```

```{r}
data(simDat5)
mod <- cslseModel(Y ~ Z | ~ X1 * X2, data = simDat5)
mod <- selSLSE(mod, "FLSE") ## We just add FLSE because we do not use BLSE
```

```{r, eval=FALSE}
c1 <- causalSLSE(mod, selType = "SLSE")
c2 <- causalSLSE(mod, selType = "FLSE", selCrit = "BIC")
c3 <- causalSLSE(mod, selType = "FLSE", selCrit = "AIC")
texreg(list(SLSE = c1, FLSE.BIC = c2, FLSE.AIC = c3), table = FALSE, digits = 4)
```

\begin{center}
```{r, results='asis', message=FALSE, echo=FALSE}
c1 <- causalSLSE(mod, selType="SLSE")
c2 <- causalSLSE(mod, selType="FLSE", selCrit="BIC")
c3 <- causalSLSE(mod, selType="FLSE", selCrit="AIC")
texreg(list(SLSE=c1, FLSE.BIC=c2, FLSE.AIC=c3), table=FALSE, digits=4)
```
\end{center}

In the case of multiple confounders with interactions, the shape of the
fitted outcome with respect to one confounder depends on the value of
the other confounders. Without interaction, changing the value of the
other confounders only shifts the fitted line without changing its
shape. The following graphs compare the estimated relationship between
$Y$ and $X_1$ for $X_2$ equal to the group means (left graph) and 1
(right graph). Using a sample of $10^7$, we obtain that $\Ex(X_2|Z=1)$
and $\Ex(X_2|Z=0)$ are approximately equal to 0.1982 and -0,3698,
respectively. Therefore, the true curves are
$(1.3698+0.6302x+1.1368x^2)$ for the nontreated and $(3.3964-2.3964x)$ for the
treated. If $X_2=1$, the true curves become $2x+2x^2$ for the treated
and $(5-4x)$ for the nontreated.

```{r, out.width='50%', fig.show='hold'}
x20 <- mean(subset(simDat5, Z == 0)$X2)
x21 <- mean(subset(simDat5, Z == 1)$X2)
arg <-  list(common = list(main = "Y vs X1 (X2 = sample mean for each group)"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X1", fixedCov = list(nontreated = list(X2 = x20), treated = list(X2 = x21)),
     graphPar = arg)
curve(1.3698 + 0.6302 * x + 1.1368 * x^2, -3, 3,
      col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(3.3964 - 2.3964 * x, -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("top", c("True-treated", "True-nontreated"),
       col=c("darkorange", "darkgreen"), lty = 3, lwd = 3, bty = 'n', cex = .8)
arg <-  list(common = list(main = "Y vS X1 (X2 = 1 for each group)"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X1", fixedCov = list(X2 = 1), graphPar = arg)
curve(2 * x + 2 * x^2, -3, 3, col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(5 - 4 * x, -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("top", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 3, lwd = 3, bty = 'n', cex = .8)
```

The following graphs illustrate the relationship between $Y$ and $X_2$
for a given $X_1$. When $X_1$ is equal to its population group means
(they are equal to the population means of $X_2$), the true curves are
$[1.6036-0.3964x)(x\leq 0)+(1+2x)(x>0)]$ for the treated and
$[(1.767-0.3698x+0.1368x^2)+ (1+x)(x\leq-1)+(-1-x)(x>-1)]$ for the
nontreated. If $X_1=1$, the true curves become $[-2x+(1-2x)(x\leq
0)+(1+2x)(x>0)]$ for the treated and
$[(4+x+x^2)+(1+x)(x\leq-1)+(-1-x)(x>-1)]$ for the nontreated.


```{r, out.width='50%', fig.show='hold'}
x10 <- mean(subset(simDat5, Z == 0)$X1)
x11 <- mean(subset(simDat5, Z == 1)$X1)
arg <-  list(common = list(main = "Y vs X2 (X1 = sample mean for each group)"),
             legend = list(x = "right", cex = 0.8))
plot(c2, "X2", fixedCov = list(nontreated = list(X1 = x10), treated = list(X1 = x11)),
     graphPar = arg)
curve(1.603900 - .3964 * x + (1 - 2 * x) * (x <= 0) + (1 + 2 * x) * (x > 0), -3, 3,
      col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(1.767 - 0.3698 * x + 0.1368 * x^2 + (1 + x) * (x <= -1) + (-1 - x) * (x > -1),
      -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("top", c("True-treated", "True-nontreated"),
       col = c("darkorange", "darkgreen"), lty = 3, lwd = 3, bty = 'n', cex = .8)
arg$common$main <- "Y vS X2 (X1 = 1 for each group)"
plot(c2, "X2", fixedCov = list(X1 = 1), graphPar = arg)
curve(-2 * x + (1 - 2 * x) * (x <= 0) + (1 + 2 * x) * (x > 0), -3, 3,
      col = "darkgreen", lty = 3, lwd = 3, add = TRUE)
curve(4 + (1 + x) * (x <= -1) + (-1 - x) * (x > -1) + x + x^2,
      -3, 3, col = "darkorange", lty = 3, lwd = 3, add = TRUE)
legend("top", c("True-treated", "True-nontreated"),
       col = c("darkgreen", "darkorange"), lty = 3, lwd = 3, bty = 'n', cex = .8)
```

\newpage 

# Summary of methods and objects \label{sec:summary}

The following is a list of all objects from the package. For each
object, we explain how it is constructed and give a list of the
registered methods. For more details about the arguments of the
different methods, see the help files. Note, however, that no help
files exist for non-exported methods and the latter must be called
using `causalSLSE:::` before the method names.

- **slseKnots**: The object is created by the function `slseKnots` and
  the only exported registered methods are `print` and `update`. 

- **slseModel** and **cslseModel**: The objects are respectively
  created by the functions `slseModel` and `cslseModel`. The exported
  registered methods are `print`, `estSLSE` (estimate the regression
  model), `selSLSE` (optimal selection of knots), `update` (select
  knots from saved selections) and `causalSLSE` (to estimate the
  causal effects). There are two non-exported methods: `pvalSLSE`
  (used to compute the p-values) and `model.matrix` (to extract the
  matrix of confounders).
  
- **slseFit** and **cslseFit**: The objects are created by the method
  `estSLSE` (the former when applied to `slseModel` objects and the
  latter when applied to `cslseModel` objects) and the exported
  registered methods are `print`, `causalSLSE` (to compute the causal
  effects), `predict` (to predict the outcome), `plot` (to plot the
  outcome as a function of one confounder) and `summary` (to give more
  details about the least squares estimation). There is one
  non-exported method, `as.model`, which extracts the model object from
  the `slseFit` or `cslseFit` object.
  
- **slseFit**: One registered method that only applies to this class
	is `extract`. It is needed to generate LaTeX table with `texreg`.
	
- **cslseFit**: One registered method that only applies to this class
	is `as.list`. It converts the objects into a list of `slseFit`
	objects.
  
- **summary.slseFit** and **summary.cslseFit**: The objects are
  created by the `summary` method for `slseFit` and `cslseFit`
  objects. The only exported registered method is `print`.

- **cslse**: The object is created by any `causalSLSE` method. It
  inherits from `cslseFit` objects. The methods that are common through
  this inheritance are `plot` and `predict`. The exported registered
  methods specific to `cslse` objects are `print`, `summary` (to give
  more details about the causal effect estimation) and `extract` (a
  method needed for `texreg`). There is one non-exported method,
  `as.model`, which extracts the model object.

Note that the method `causalSLSE` is also registered for objects of
class `formula`. 

# Experiments (to be removed before publishing)

```{r}
data(simDat5)
mod <- cslseModel(Y ~ Z | ~ X1 * X2, data = simDat5)
getAlt <- function(which=c("ipw","matching"), ...)
{
    which <- match.arg(which)
    met <- c("ACE","ACT","ACN")
    l <- lapply(met, function(mi)
    {
        arg <- list(...)
        arg$type <- mi
        res <- do.call(which, arg)
    })
    if (length(l[[1]]$estim)==2)
    {
        l <- lapply(1:2, function(i) {
            obj <- lapply(l, function(li) {
                li$estim <- li$estim[i]
                li$se <- li$se[i]
                li})
            names(obj) <- met
            class(obj) <- "allAlt"
            obj})
        names(l) <- c("Matching", "BC_Matching")
    } else {
        names(l) <- met
        class(l) <- "allAlt"
    }
    l
}

setMethod("extract", signature = className("allAlt", package='causalSLSE'),
          definition = function (model, include.nobs = TRUE, ...) 
          {
              se <- sapply(model, function(li) li$se)
              co <- sapply(model, function(li) li$estim)
              pval <- 2*pnorm(-abs(co/se))
              names(co) <- names(se) <- names(pval) <- names(model)
              gof <- numeric()
              gof.names <- character()
              gof.decimal <- logical()
              if (isTRUE(include.nobs)) {
                  n <- nrow(model[[1]]$data)
                  gof <- c(gof, n)
                  gof.names <- c(gof.names, "Num. obs.")
                  gof.decimal <- c(gof.decimal, FALSE)
              }
              tr <- createTexreg(coef.names = names(co), coef = co, se = se, 
                                 pvalues = pval, gof.names = gof.names,
                                 gof = gof, gof.decimal = gof.decimal)
              return(tr)
          })

res <- getAlt("ipw", form=Y~Z, psForm=Z~X1*X2, data=simDat5, normalized=TRUE)
res2 <- getAlt("matching", form=Y~Z, psForm=Z~X1*X2, data=simDat5, bcForm=~X1*X2,
               balm=~X1*X2)
texreg(list(IPW=res, Matching=res2[[1]], BC.Matching=res2[[2]]))
```

# References
