2gcomp2.Rmd

# G-computation using ML

```{r setup01gm, include=FALSE}
require(knitr)
require(glmnet)
require(kableExtra)
require(dplyr)
require(xgboost)
require(caret)
require(SuperLearner)
options(knitr.kable.NA = '')
cachex=TRUE
```

```{block, type='rmdcomment'}
G-computation is highly sensitive to **model misspecification**; and when model is **not correctly specified**, result is subject to bias. 
```

- Therefore, it can be a good idea to use **machine learning** methods, that are more flexible, than parametric methods to estimate the treatment effect.
- Although ML methods are powerful in point estimation, the **coverage probabilities** are usually poor when more flexible methods are used, if inference is one of the goals. Hence we are focusing on **point estimation** here.

## G-comp using Regression tree

```{r reg2sl, cache=cachex, echo = TRUE}
# Read the data saved at the last chapter
ObsData <- readRDS(file = "data/rhcAnalytic.RDS")
baselinevars <- names(dplyr::select(ObsData, !A))
out.formula <- as.formula(paste("Y~ A +",
                               paste(baselinevars,
                                     collapse = "+")))
```

### A tree based algorithm

XGBoost is a fast version of gradient boosting algorithm. Let us use this one to fit the data first. We follow the exact same procedure that we followed in the parametric G-computation setting.

```{r ML1, cache=cachex, echo = TRUE, warning= FALSE}
require(xgboost)
Y <-ObsData$Y 
ObsData.matrix <- model.matrix(out.formula, data = ObsData)
```


```{r ML1fit, cache=cachex, echo = TRUE}
fit3 <- xgboost(data = ObsData.matrix, 
                label = Y,
                max.depth = 10, 
                eta = 1, 
                nthread = 15, 
                nrounds = 100, 
                alpha = 0.5,
                objective = "reg:squarederror", 
                verbose = 0)
```


```{r ML1fit2, cache=cachex, echo = TRUE}
predY <- predict(fit3, newdata = ObsData.matrix)
plot(density(Y), 
     col = "red", 
     main = "Predicted and observed Y",
     xlim = c(1,100))  
legend("topright", 
       c("Y","Predicted Y"), 
       lty = c(1,2), 
       col = c("red","blue"))
lines(density(predY), col = "blue", lty = 2)
caret::RMSE(predY,Y)
```

- What we have done here is we have used the `ObsData.matrix` data to train our model, and we have used `newdata = ObsData.matrix` to obtain prediction. 

```{block, type='rmdcomment'}
When we use same **data** for training and obtaining prediction, often the predictions are highly optimistic (RMSE is unrealistically low for future predictions), and we call this a **over-fitting** problem.
```

- One way to deal with this problem is called Cross-validation.

### Cross-validation

Cross-validation means

- splitting the data into
  - training data
  - testing data


```{block, type='rmdcomment'}
In each iteration: (1) Fitting models in training data (2) obtaining prediction $\hat{Y}$ in test data (3) obtain all RMSEs from each iteration, and (4) average all RMSEs.
```

```{r cvpic, echo = FALSE, out.width = "650px", fig.cap="Cross validation from [wiki](https://en.wikipedia.org/wiki/Cross-validation_(statistics)); training data = used for building model; test data = used for prediction from the model that was built using training data; each iteration = fold"}
knitr::include_graphics("images/CV.png")
```

#### Cross-validation using caret

We use `caret` package to do cross-validation. 

```{block, type='rmdcomment'}
`caret` is a general framework package for machine learning that can also incorporate other ML approaches such as `xgboost`.
```

```{r ML1car, cache=cachex, echo = TRUE}
require(caret)
set.seed(123)
X_ObsData.matrix <- xgb.DMatrix(ObsData.matrix)
Y_ObsData <- ObsData$Y
```

Below we define $K = 3$ for cross-validation. Ideally for a sample size close to $n=5,000$, we would select $K=10$, but for learning / demonstration / computational time-saving purposes, we just use $K = 3$. 

```{r ML1carx, cache=cachex, echo = TRUE}
xgb_trcontrol = trainControl(
  method = "cv",
  number = 3,  
  allowParallel = TRUE,
  verboseIter = FALSE,
  returnData = FALSE
)
```

#### Fine tuning

```{block, type='rmdcomment'}
One of the advantages of `caret` framework is that, it also allows checking the impact of various parameters (can do **fine tuning**).  
```

For example, 

- for interaction depth, we previously use `max.depth = 10`. That means $covariate^{10}$ polynomial.
- We could also check if other interaction depth choices (such as $covariate^{2}$ or $covariate^{4}$) would be better in terms of honest predictions.

```{r ML1carxy, cache=cachex, echo = TRUE}
xgbGrid <- expand.grid(
  nrounds = 100, 
  max_depth = seq(2,10,2),
  eta = 1,
  gamma = 0,
  colsample_bytree = 0.1,
  min_child_weight = 2,
  subsample = 0.5 
)
```


#### Fit model with CV

once we set

- resampling or cross-validation settings
- parameter grid

we can fit the model:

```{r ML1car2, cache=cachex, echo = TRUE}
fit.xgb <- train(
  X_ObsData.matrix, Y_ObsData,  
  trControl = xgb_trcontrol,
  method = "xgbTree",
  tuneGrid = xgbGrid,
  verbose = FALSE
)
fit.xgb
```

Based on the loss function (say, RMSE) it automatically chose the best tuning parameter set:

```{r ML1car3, cache=cachex, echo = TRUE}
fit.xgb$bestTune$max_depth
``` 

```{r}
predY <- predict(fit.xgb, newdata = ObsData.matrix)
plot(density(Y), 
     col = "red", 
     main = "Predicted and observed Y",
     xlim = c(1,100))  
legend("topright", 
       c("Y","Predicted Y"), 
       lty = c(1,2), 
       col = c("red","blue"))
lines(density(predY), col = "blue", lty = 2)
caret::RMSE(predY,Y)
```


### G-comp step 2: Extract outcome prediction as if everyone is treated

```{r ML12, cache=cachex, echo = TRUE}
ObsData.matrix.A1 <- ObsData.matrix 
ObsData.matrix.A1[,"A"] <- 1
ObsData$Pred.Y1 <- predict(fit.xgb, newdata = ObsData.matrix.A1)
summary(ObsData$Pred.Y1)
```

### G-comp step 3: Extract outcome prediction as if everyone is untreated


```{r ML13, cache=cachex, echo = TRUE}
ObsData.matrix.A0 <- ObsData.matrix
ObsData.matrix.A0[,"A"] <- 0
ObsData$Pred.Y0 <- predict(fit.xgb, newdata = ObsData.matrix.A0)
summary(ObsData$Pred.Y0)
```

### G-comp step 4: Treatment effect estimate


```{r ML13b, cache=cachex, echo = TRUE}
ObsData$Pred.TE <- ObsData$Pred.Y1 - ObsData$Pred.Y0  
```


Mean value of predicted treatment effect 

```{r reg2acnx1b, cache=cachex, echo = TRUE}
TE1 <- mean(ObsData$Pred.TE)
TE1
summary(ObsData$Pred.TE)
```

Notice that the mean is slightly different than the parametric G-computation method.

## G-comp using regularized methods

### A regularized model

```{block, type='rmdcomment'}
LASSO is a regularized method. One of the uses of these methods is "variable selection" or addressing concerns of multicollinearity. 
```

Let us use this method to fit our data. 

- We are again using cross-validation here, and we chose $K=3$.

```{r ML1r, cache=cachex, echo = TRUE, warning=FALSE}
require(glmnet)
Y <-ObsData$Y
ObsData.matrix <- model.matrix(out.formula, data = ObsData)
fit4 <-  cv.glmnet(x = ObsData.matrix, 
                y = Y,
                alpha = 1,
                nfolds = 3,
                relax=TRUE)
```


### G-comp step 2: Extract outcome prediction as if everyone is treated

```{r ML12r, cache=cachex, echo = TRUE}
ObsData.matrix.A1 <- ObsData.matrix 
ObsData.matrix.A1[,"A"] <- 1
ObsData$Pred.Y1 <- predict(fit4, newx = ObsData.matrix.A1,
                           s = "lambda.min")
summary(ObsData$Pred.Y1)
```

### G-comp step 3: Extract outcome prediction as if everyone is untreated


```{r ML13r, cache=cachex, echo = TRUE}
ObsData.matrix.A0 <- ObsData.matrix
ObsData.matrix.A0[,"A"] <- 0
ObsData$Pred.Y0 <- predict(fit4, newx = ObsData.matrix.A0,
                           s = "lambda.min")
summary(ObsData$Pred.Y0)
```

### G-comp step 4: Treatment effect estimate


```{r ML13br, cache=cachex, echo = TRUE}
ObsData$Pred.TE <- ObsData$Pred.Y1 - ObsData$Pred.Y0  

```


Mean value of predicted treatment effect 

```{r reg2acnx1br, cache=cachex, echo = TRUE}
TE2 <- mean(ObsData$Pred.TE)
TE2
summary(ObsData$Pred.TE)
```

Notice that the mean is very similar to the parametric G-computation method.

## G-comp using SuperLearner

```{block, type='rmdcomment'}
SuperLearner is an ensemble ML technique, that uses **cross-validation** to find a weighted combination of estimates provided by different **candidate learners** (that help predict).
```

- There exists many candidate learners. Here we are using a combination of
  - linear regression
  - Regularized regression (lasso)
  - gradient boosting (tree based)
  
### Steps  

|  |  |
|-|-|
|Step 1| Identify candidate learners|
|Step 2| Choose Cross-validation K|
|Step 3| Select loss function for meta learner|
|Step 4| Find SL prediction: (1) Discrete SL (2) Ensamble SL |


#### Identify candidate learners

- Choose variety of candidate learners 
   - parametric (linear or logistic regression)
    - regularized (LASSO, ridge, elasticnet)
    - stepwise
   - non-parametric 
    - transformation (SVM, NN)
    - tree based (bagging, boosting)
   - smoothing or spline (gam)
- tune the candidate learners for better performance
  - tree depth
  - tune regularization parameters
  - variable selection

```{r ML1s00, cache=cachex, echo = TRUE}
SL.library.chosen=c("SL.glm", "SL.glmnet", "SL.xgboost")
```

  
**SuperLearner** is an ensemble learning method. Let us use this one to fit the data first.

#### Choose Cross-validation K

To combat against optimism, we use cross-validation. SuperLearner first **splits** the data according to chosen $K$ fold for the cross-validation. 

```{r ML1s000, cache=cachex, echo = TRUE}
cvControl.chosen = list(V = 3)
```

#### Select loss function for meta learner and estimate risk

```{block, type='rmdcomment'}
The goal is to minimize the estimated risk (i.e., minimize the difference of $Y$ and $\hat{Y}$) that comes out of a model. 
```

<!---
For each fold, estimate a **measure of performance** (could be RMSE) in test sets based on models that was built using training sets

$RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^n (Y - \hat{Y})^2}$ for continuous $Y$

- we obtain risk estimate in each fold (from test data)
- we average all the estimates risks
---->

We can chose a (non-negative) least squares loss function for the meta learner (explained below):

```{r ML1s0, cache=cachex, echo = TRUE}
loss.chosen = "method.NNLS"
```

#### Find SL prediction

We first fit the super learner:

```{r ML1s, cache=cachex, echo = TRUE, results='hide', warning = FALSE}
require(SuperLearner)
ObsData.noY <- dplyr::select(ObsData, !Y)
fit.sl <- SuperLearner(Y=ObsData$Y, 
                       X=ObsData.noY, 
                       cvControl = cvControl.chosen,
                       SL.library=SL.library.chosen,
                       method=loss.chosen,
                       family="gaussian")
```

We can also obtain the predictions from each candidate learners.

```{r ML12stest0, cache=cachex, echo = TRUE}
all.pred <- predict(fit.sl, type = "response")
Yhat <- all.pred$library.predict
head(Yhat)
```  

We can obtain the $K$-fold cross-validated risk estimates for each candidate learners.

```{r ML12stestrisk, cache=cachex, echo = TRUE}
fit.sl$cvRisk
```

Once we have the performance measures and predictions from candidate learners, we could go one of **two routes** here

##### Discrete SL

```{block, type='rmdcomment'}
Get measure of performance from all folds are averaged, and choose the **best** one. The prediction from the chosen learners are then used.
```

`glmnet` has the lowest cross-validated risk

```{r ML12stest, cache=cachex, echo = TRUE}
lowest.risk.learner <- names(which(
  fit.sl$cvRisk == min(fit.sl$cvRisk)))
lowest.risk.learner
as.matrix(head(Yhat[,lowest.risk.learner]), 
          ncol=1)
```  

##### Ensamble SL 

Here are the first 6 rows from the candidate learner predictions:

```{r ML12stestx, cache=cachex, echo = TRUE}
head(Yhat)
``` 

```{block, type='rmdcomment'}
fit a **meta learner** (optimal weighted combination; below is a simplified description)
```

- using
  - linear regression (without intercept, but could produce -ve coefs) or 
  - preferably non-negative least squares for 
  
  $Y_{obs}$ $\sim$ $\hat{Y}_{SL.glm}$ + $\hat{Y}_{SL.glmnet}$ + $\hat{Y}_{SL.xgboost}$. 

- Obtain the regression coefs $\mathbf{\beta}$ = ($\beta_{SL.glm}$, $\beta_{SL.glmnet}$, $\beta_{SL.xgboost}$) for each $\hat{Y}$, 
- scale them to 1 
  - $\mathbf{\beta_{scaled}}$ = $\mathbf{\beta}$ / $\sum_{i=1}^3{\mathbf{\beta}}$; 
  - so that the  sum of scaled coefs =  1 
- Scaled coefficients $\mathbf{\beta_{scaled}}$ represents the **value / importance of the corresponding candidate learner**. 

```{r ML12stestxreg, cache=cachex, echo = TRUE, include= FALSE}
# fit.meta <- lm(Y ~ 0 + Yhat)
# lm produces -ve
# non-negative least squares (NNLS)
fit.meta <- nnls::nnls(A=Yhat,b=Y)
# alternate would be to use glmnet
# fit.meta <- glmnet(X, Yhat, 
#                    lambda = 0, 
#                    lower.limits = 0, 
#                    intercept = FALSE)
fit.meta
coefs <- coef(fit.meta)
coefs
scaled.coefs <- abs(coefs)/sum(abs(coefs))
scaled.coefs
``` 

Scaled coefs

```{r ML12stestcoef, cache=cachex, echo = TRUE}
fit.sl$coef
```

```{r ML12stestcoef2b, cache=cachex, echo = TRUE}
sum(fit.sl$coef)
```

Hence, in creating superlearner prediction column,

a. Linear regression has no contribution
b. lasso has majority contribution
c. gradient boosting of tree has some minimal contribution

- A new prediction column is produced based on the fitted values from this meta regression.

You can simply multiply these coefs to the predictions from candidate learners, and them sum them to get ensable SL. Here are the first 6 values:

```{r ML12stestb0, cache=cachex, echo = TRUE}
SL.ens <- t(t(Yhat)*fit.sl$coef)
head(SL.ens)
as.matrix(head(rowSums(SL.ens)), ncol = 1)
```  

Alternatively, you can get them directly from the package: here are the first 6 values

```{r ML12stestb, cache=cachex, echo = TRUE}
head(all.pred$pred)
```  

The last column is coming from Ensamble SL.

### G-comp step 2: Extract outcome prediction as if everyone is treated

We are going to use **Ensamble SL** predictions in the following calculations. If you wanted to use discrete SL predictions instead, that would be fine too.

```{r ML12s, cache=cachex, echo = TRUE}
ObsData.noY$A <- 1
ObsData$Pred.Y1 <- predict(fit.sl, newdata = ObsData.noY,
                           type = "response")$pred
summary(ObsData$Pred.Y1)
```

### G-comp step 3: Extract outcome prediction as if everyone is untreated


```{r ML13s, cache=cachex, echo = TRUE}
ObsData.noY$A <- 0
ObsData$Pred.Y0 <- predict(fit.sl, newdata = ObsData.noY,
                           type = "response")$pred
summary(ObsData$Pred.Y0)
```

### G-comp step 4: Treatment effect estimate


```{r ML13bs, cache=cachex, echo = TRUE}
ObsData$Pred.TE <- ObsData$Pred.Y1 - ObsData$Pred.Y0  

```


Mean value of predicted treatment effect 

```{r reg2acnx1bs, cache=cachex, echo = TRUE}
TE3 <- mean(ObsData$Pred.TE)
TE3
summary(ObsData$Pred.TE)
```

### Additional details for SL

#### Choice of K

- simplest cross-validation splits the data into $K=2$ parts, but can go higher.
  - select $K$ judiciously
    - large sample size means small $K$ may be adequate
      - for $n \lt 10,000$ consider $K=3$
      - for $n \lt 500$ consider $K=20$
    - smaller sample size means larger $K$ may be necessary
      - for $n \lt 30$ consider leave 1 out 

#### Alternative to CV

- other similar algorithms such as **cross-fitting** had been shown to have better performances

#### Rare outcome

- for rare outcomes, consider using **stratification** to attempt to maintain training and test sample ratios the same

#### Dependant sample

- if data is clustered and not independent and identically distributed, use ID for the **cluster**

#### Choice of meta learner method

It is easy to show that, depending on the choice of meta-learners, the coefficients of the meta learners can be slightly different.

```{r ML12lossf, cache=cachex, echo = TRUE}
fit.sl2 <- recombineSL(fit.sl, Y = Y, 
                       method = "method.NNLS2")
fit.sl2$coef
fit.sl2 <- recombineSL(fit.sl, Y = Y, 
                       method = "method.CC_LS")
fit.sl2$coef
fit.sl4 <- recombineSL(fit.sl, Y = Y, 
                       method = "method.CC_nloglik")
fit.sl4$coef
```

- `method.CC_LS` is [suggested](https://si.biostat.washington.edu/sites/default/files/modules/lab1_0.pdf) as a good method for continuous outcome
- `method.CC_nloglik` is [suggested](https://si.biostat.washington.edu/sites/default/files/modules/lab1_0.pdf) as a good method for binary outcome

```{r, cache=TRUE, echo = TRUE}
saveRDS(TE1, file = "data/gcompxg.RDS")
saveRDS(TE2, file = "data/gcompls.RDS")
saveRDS(TE3, file = "data/gcompsl.RDS")
```