---
title: "Quickstart"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{quickstart}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center"
)
library(loewesadditivity)
```

# Installation

To install `loewesadditivity`, run the following code.

```{r eval = FALSE}
if(!require(loewesadditivity)){
  library(devtools)
  devtools::install_github("shannong19/loewesadditivity")
}
library(loewesadditivity)
```

```{r echo = F, message = F, warning = F}
# hidden load currently
#library(devtools)
devtools::load_all()
```

# Quickstart

```{r setup, message = FALSE, warning = FALSE}
library(ggplot2)
library(dplyr)
library(knitr)
library(kableExtra)

results <- rh5_ama1ron2 %>% 
  rename(dose_A = "RH5", dose_B = "AMA1RON2")  %>% 
  fortify_gia_data() %>%
  estimate_params() 

results$params_est


```


## Look at the raw data


This is a raw set of GIA data that has been collected and uploaded to this package.  We measure the GIA% for different experiments, plates (a or b) and repetitions.  We are using the compounds AMA1RON2 and RH5.  

The variables iRBC and uRBC are optional rows which measure the minimal and maximum GIA respectively.

```{r data, message = FALSE, warning = FALSE}
data(rh5_ama1ron2)

rh5_ama1ron2 %>% head() %>% knitr::kable() %>%
  kableExtra::kable_styling(full_width = T, font_size = 7)

```


## Fortify data

The standard output we use is that of the long, data.frame with the columns

1. `dose_A`
2. `dose_B`
3. `GIA`

where `dose_A` and `dose_B` are the unscaled doses (e.g. mg/mL) and GIA is the percent GIA.

We provide the function `fortify_gia_data()` which works when the data is in the above format with column headings of  `exp{x}{y}rep{z}` along with the two doses (must have names `dose_A` and `dose_B`.



```{r data_f, message = FALSE, warning = FALSE}

rh5_ama1ron2_f <-
  rh5_ama1ron2 %>% rename(dose_A = 'RH5', dose_B = 'AMA1RON2') %>%
  fortify_gia_data()
  
  rh5_ama1ron2_f %>% head() %>%  kable(align = "c") %>%
  kableExtra::kable_styling(full_width = F)
```

The fortification has averaged over the repetitions but has kept the plates and the experiments.


## Estimate parameters

You can use the function `estimate_params()` to estimate the parameters.  The parameters and their meaning can be found in the documentation of `?base_GIA`.


```{r est, message = FALSE, warning = FALSE}
## Set up initial guesses and parameters
model_params <- c(
  "beta_A" = .5,
  "beta_B" = .5,
  "gamma_A" = .5,
  "gamma_B" = .5,
  "tau_1" = 0,
  "tau_2" = 0
  )
  n_boot <- 10
  fn_list <- NULL
  alpha <- .05
  verbose <- TRUE
  
  
out <- estimate_params(
    data = rh5_ama1ron2_f,
    init_params = model_params,
    n_boot = n_boot,
    alpha = alpha,
    verbose = verbose
  )
```


Hewlett's S can be found below.

```{r S, message = FALSE, warning = FALSE}
out$S_est  %>%  kable(align = "c", digits = 3) %>%
  kableExtra::kable_styling(full_width = F)
```


The parameter estimates from the model are shown here:

```{r par_est, message = FALSE, warning = FALSE}
out$params_est %>%  kable(align = "c", digits = 3) %>%
  kableExtra::kable_styling(full_width = F)
```

Finally, the (first six) individual observation estimates and 95% CIs can be found here:

```{r GIA_est, message = FALSE, warning = FALSE}
out$GIA_est %>% head()  %>%  kable(align = "c", digits = 3) %>%
  kableExtra::kable_styling(full_width = F)
```


## Plot results

We can plot the estimated surface, individual curves, and the isobologram.

```{r, message = FALSE, warning = FALSE, fig.width=8, fig.height = 6}
# Note xlab is dose_A and ylab is dose_B
g1 <- plot_surface(out, xlab = "RH5", ylab = "AMA1RON2")
```

```{r, message = FALSE, warning = FALSE, fig.width=8, fig.height = 12}
g2 <- plot_curves(out, dose_A = "RH5", dose_B = "AMA1RON2" )
```

```{r, message = FALSE, warning = FALSE, fig.width=8, fig.height = 6}
g3 <- plot_isobologram(out, dose_A = "RH5", dose_B = "AMA1RON2" )
```


## Simulation of coverage

We provide the function `simulate_coverage()` as a way to get an estimate of the power of the experiment to detect synergy or antagonism.  For given model parameters, an experimental grid of dose combinations, and an assumed noise structure for GIA, we can determine 1) percent of times we expect 0 to in the 95% CI of the interaction parameter $\tau_1$ and 2) percent of times we expect the true given value of the model parameters to be in the 95% CI.

**WARNING** It is advised to use at least 100 bootstraps for each of 10 simulations but is not shown here.

```{r coverage, message = FALSE, warning = FALSE}
  model_params <- c("beta_A" = .250, "beta_B" = .146,
                    "gamma_A" = .527, "gamma_B" = .921,
                    "tau_1" = -.058, "tau_2" = -.289)
  experimental_grid <- make_grid(par = model_params, 
                                 n = 5)
  n_boot <- 3
  n_sims <- 5
  GIA_fn <- base_GIA
  S_fn <- calc_S_base
  fn_list <- NULL
  alpha <- .05
  verbose <- TRUE
  out <- simulate_coverage(n_sims = n_sims,
                         n_boot = n_boot,
                         verbose = FALSE,
                         experimental_grid = experimental_grid,
                         model_par = model_params,
                         alpha = .05,
                         noise_par = c("a0" = 3, "a1" = .01),
                         GIA_fn = base_GIA,
                         fn_list = NULL)
  
  out

```


## Design an experiment

A unique feature of `loewesadditivity` is the ability to generate code to use to create your own coverage simulation in conjunction with a (hopefully) intuitive shiny App.  Run the below code, copy and paste the results into R and see what happens!

Ex.

```{r}
 out <- design_experiment(n_rep = 2)
```
