---
title: "Plot rms model summaries and predictions"
author: "Richard Meitern"
date: "06/11/2019"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Plot rms model summaries and predictions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warnings = FALSE,
  fig.width = 6,
  fig.height = 4
)
library(rms)
```



### Load the required libraries
```{r setup}
#the regression modelling by F. Harrell
library(rms)

#the extension for plotting the rms summary objects
library(ormPlot)

#to modify plots
library(ggplot2)
```

### Looking at the bundled data
OrmPlot has `educ_data` bundled. `educ_data` is used in this vignette
to show the functionality of the package. To familiarize with the data:
```{r eval=FALSE}
#show first 6 rows
head(educ_data)
```


```{r echo=FALSE, results="asis"}
pander::pandoc.table(head(educ_data), split.tables=Inf)

```
```{r, eval = FALSE}
#show variable explanation
help(educ_data)
```

### Setting up rms datadist
The rms package requires that a `datadist` object is set up properly. According
to `datadist` documentation:
`q.effect` is a set of two quantiles for computing the range of continuous
variables to use in estimating regression effects.
Defaults are `c(.25,.75)`, which yields inter-quartile-range odds ratios
```{r}
#q.effect determines for what range the odds ratios are given on plots
dd <- datadist(educ_data, q.effect = c(0.5, 0.75))
#set it also to options
options(datadist="dd")
```

### Creating a model

```{r}
#see help(orm) for further info
orm_model<-orm(educ_3 ~ Rural + sex + n_siblings + cran_rzs + height_rzs + 
    FW_rzs + YOBc + (YOBc * sex) + (YOBc * Rural), data = educ_data)
```

### Plotting the model predictions with CI

The main advantage of using ormPlot is that you get plots with confidence intervals
shown on the plot.

The simplest way is to predict for only one value. Plotting returns a 
customizable ggplot object.
```{r}
plot(orm_model, cran_rzs)
```

For more complex models specify facet column and rows.
```{r}
plot(orm_model, cran_rzs,  Rural, sex)

```


You can easily set custom labels.
```{r}

p<-plot(orm_model, cran_rzs,  Rural, sex,
        xlab = "Cranial volume (residuals to age an birth date)",
        facet_labels = list(Rural = c("Urban", "Rural"),sex=c("Male","Female")))



colors <- c("#4a9878", "#0a191e", "#d8b65c")
educ_names <- c("Primary", "Secondary", "Tertiary")

# further modifing like any other ggplot
final_plot<-p + labs(color = "Education", fill = "Education") + 
  scale_color_manual(values = colors, labels = educ_names) +
  scale_fill_manual(values = colors, labels = educ_names) 
 
final_plot 
```


Save like any ggplot graph.
```{r eval = FALSE}
ggsave("educ_cran.svg",final_plot, height = 8 ,width = 8)
```


### Plotting the model summary
The easiest way is to just plot the summary object
```{r fig.width=7, fig.height=4 }
forestplot(summary(orm_model))
```

If this does not look nice enough you can also get `ggplot2` objects to customize to your needs.
The best way to get customizable plots is to specify `return_ggplots=TRUE`
```{r fig.width=6}
# you can use also use plot instead of forestplot
plots<-forestplot(summary(orm_model), return_ggplots=T)
plots[[1]]
plots[[2]]
```


These can be joined using the join_ggplots command. You can edit the plots as 
any ggplot plot
```{r, fig.width=7}

p1 <- plots[[1]] + 
  scale_x_discrete(labels=c("Mean", "Lower CI", "Upper CI"), 
                   position = "top",
                   name = NULL)

# the x axis is actually y axis because the cordinates are flipped with coord_flip()
p2 <- plots[[2]] + scale_y_continuous(breaks = c(0.5, 0.7, 0.9, 1.1),
                                position = "right")
  

forestplot<-join_ggplots(p1,p2)

```

To save as svg fur further editing just use `ggsave` from `ggplot2`
```{r, eval=FALSE}
ggsave("forestplot.svg",forestplot)
```