---
title: "Workflow to Develop and Validate a Prediction model in Multiply Imputed data"
author: "Martijn W Heymans"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Workflow to Develop and Validate a Prediction model in Multiply Imputed data}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

# Introduction {-}

This vignette shows how to develop, internally and externally validate
a (logistic) regression prediction model with `mice` and `psfmi`.

```{r, include = FALSE, message = FALSE}
library(knitr)
opts_chunk$set(tidy = FALSE, cache = FALSE)
library(mice)
library(psfmi)
```

# Steps {-}

   + [Step 1 - Install the psfmi and mice package]
   + [Step 2 - Impute the missing data with mice]
   + [Step 3 - Develop the model]
   + [Step 4 - Determine model performance]
   + [Step 5 - Internally validate the model]
   + [Step 6 - Shrink pooled coefficients and determine new intercept]
   + [Step 7 - Externally validate the model]
   

## Step 1 - Install the psfmi and mice package {-}

You can install the released version of psfmi with:

``` r
install.packages("psfmi")
```
And the development version from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
devtools::install_github("mwheymans/psfmi")
```

You can install the released version of mice with:

``` r
install.packages("mice")

``` 

Back to [Steps]

## Step 2 - Impute the missing data with mice {-}

I generated 5 imputed datasets to handle the missing values in the lbp_orig dataset that is
included in the `psfmi` package. Depending on the amount of missing data, the number of imputed 
datasets may be increased (see on-line book: 
[Applied Missing Data analysis](https://bookdown.org/mwheymans/bookmi/multiple-imputation.html#number-of-imputed-datasets-and-iterations)). 

```{r}

  imp <- mice(lbp_orig, m=5, maxit=10, seed = 750) 
  
```

Use the `complete` function of the `mice` package to extract the completed datasets. By 
setting `action = "long"` and `include = FALSE`, the imputed datasets are stacked under
each other to form one long dataset and the original dataset (that included missing 
values) is not included in that long dataset.

```{r}
  data_comp <- complete(imp, action = "long", include = FALSE)
```

Back to [Steps]

## Step 3 - Develop the model {-}

Use the `psfmi_lr` function in the `psfmi` package to perform backward selection over
the 5 multiply imputed datasets with a p-value of 0.05 and as selection method D1.  
**Note** that with this function it is possible to include cubic spline coefficients in case
of non-linear relationships during backward selection and that predictors may be forced 
in the model during backward selection by including them at the setting `keep.predictors`.
When the setting `direction` is changed to FW, forward selection is done.

```{r}

  pool_lr <- psfmi_lr(data=data_comp, nimp=5, impvar=".imp", 
              formula=Chronic ~ Gender + Age + JobControl + Tampascale + 
              Pain + Radiation + JobDemands + SocialSupport + Smoking + 
              factor(Satisfaction) + factor(Carrying) + rcs(Function, 3), 
              keep.predictors = "Radiation", 
              p.crit = 0.157, method="D1", direction = "BW")
  
  pool_lr$RR_model_final
  pool_lr$multiparm_final

```

Back to [Steps]

## Step 4 - Determine model performance {-}

The performance of the model selected at step 3 can be determined by using the
`pool_performance` function.

```{r plot, fig.height = 3, fig.width = 5}

  perf <- pool_performance(data=data_comp, nimp=5, impvar=".imp", 
    formula = pool_lr$formula_final[[1]], 
    cal.plot=TRUE, plot.method="mean", groups_cal = 10)
  perf

```

## Step 5 - Internally validate the model {-}

To internally validate the model we use the `psfmi_` function. With this function five different
methods can be used to internally validate models in MI data [see these Vignettes](https://mwheymans.github.io/psfmi/articles/),
Three methods use cross-validation and two bootstrapping in combination with MI.

We will use cross-validation and more specific the method [cv_MI_RR](https://mwheymans.github.io/psfmi/articles/cv_MI_RR.html). With this method it is possible
to integrate backward selection into the cross-validation procedure. The last model that is selected 
by the psfmi_lr function is internally validated. So, if we want to apply backward selection during 
cross-validation from the full model we first have to apply the psfmi_lr function without variable 
selection. That is what we apply here, because we know that variable selection 
is one of the main reasons that prediction models are over-fitted.

```{r}

  pool_val <- psfmi_lr(data=data_comp, formula = Chronic ~ Gender + Age + JobControl + Tampascale + 
                      Pain + Radiation + JobDemands + SocialSupport + Smoking + factor(Satisfaction) + 
                      factor(Carrying) + rcs(Function, 3), p.crit = 1, direction="BW",
                      nimp=5, impvar=".imp", method="D1")

  set.seed(200)
  res_cv <- psfmi_validate(pool_val, val_method = "cv_MI_RR", data_orig = lbp_orig, folds = 5,
                     p.crit=0.05, BW=TRUE, nimp_mice = 10, miceImp = miceImp, printFlag = FALSE)

  res_cv

```

Back to [Steps]

## Step 6 - Shrink pooled coefficients and determine new intercept {-}

We can use the slope value of 0.8449148 that was estimated at the previous step as a 
shrinkage factor to prevent our model from being overfitted in new data. 
We do this by multiplying the pooled coefficients with the shrinkage factor and 
also to determine a new intercept value that is aligned with the shrunken coefficients.

We use the `pool_intadj` function for this that will provide the adjusted coefficients and 
new intercept value.

```{r}

  pool_select <- psfmi_lr(data=data_comp, nimp=5, impvar=".imp", 
  formula = pool_lr$formula_final[[1]],  
  p.crit = 1, method="D1")

  res <- pool_intadj(pool_select, shrinkage_factor = 0.8449148)

  res$int_adj
  res$coef_shrink_pooled

```

The last step is to externally validate this adjusted model.

Back to [Steps]

## Step 7 - Externally validate the model {-}

We validate the model in an external dataset that is imputed five times
due to missing data. 

```{r plot2, fig.height = 3, fig.width = 5}

  res_extval <- mivalext_lr(data.val = lbpmi_extval, nimp = 5, impvar = "Impnr",
    formula = pool_lr$formula_final[[1]], lp.orig = c(res$int_adj, 
    res$coef_shrink_pooled), cal.plot = TRUE, plot.method = "mean")
 
  res_extval$ROC
 
  res_extval$R2
 
  res_extval$HLtest
 
  res_extval$LP_pooled_ext
 
```

In the external dataset the AUC is 0.8482 and the R-squared 0.45587. The
Hosmer and Lemeshow test has a p-value of 0.67820, which means that
the fit is satisfactory. This is also confirmed by the calibration plot
where the calibration curves averaged across multiply imputed datasets is 
near the optimal (dashed) line. The slope (regression coefficients) has a 
value of 0.88741 and slightly deviates from the value of 1, which means 
that the coefficient values differ between the development and external dataset.

Back to [Steps]
