---
title: "Examples in R"
output: 
  rmarkdown::html_vignette:
    toc: true
bibliography: references.bib
link-citations: yes
csl: american-medical-association.csl
vignette: >
  %\VignetteIndexEntry{Examples in R}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  message = TRUE,
  warning = TRUE,
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

Here I show how to produce _P_-value, _S_-value, likelihood, and deviance functions with the `concurve` package using fake data and data from real studies. Simply put, these functions are rich sources of information for scientific inference and the image below, taken from Xie & Singh, 2013[@xie2013isr] displays why. 

<img src="densityfunction.png" align="center" width="750" />

For a more extensive discussion of these concepts, see the following references. [@birnbaum1961ams; @chow2019asb; @fraser2017arsa; @fraser2019as; @Poole1987-nb; @poole1987ajph; @Schweder2002-vh; @schweder2016; @Singh2007-zr; @Sullivan1990-ha; @whitehead1993sm; @xie2013isr; @rothman2008me]

# Simple Models

To get started, we could generate some normal data and combine two vectors in a dataframe

```{r echo=TRUE, fig.height=4.5, fig.width=6}
library(concurve)
set.seed(1031)
GroupA <- rnorm(500)
GroupB <- rnorm(500)
RandomData <- data.frame(GroupA, GroupB)
```

and look at the differences between the two vectors. We'll plug these vectors and the dataframe they're in inside of the `curve_mean()` function. Here, the default method involves calculating CIs using the Wald method.  

``` {r}
intervalsdf <- curve_mean(GroupA, GroupB,
  data = RandomData, method = "default"
)
```

Each of the functions within `concurve` will generally produce a list with three items, and the first will usually contain the function of interest. 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
head(intervalsdf[[1]], 10)
```

We can view the function using the `ggcurve()` function. The two basic arguments that must be provided are the data argument and the "type" argument. To plot a consonance function, we would write "`c`". 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(function1 <- ggcurve(data = intervalsdf[[1]], type = "c", nullvalue = TRUE))
```

We can see that the consonance "curve" is every interval estimate plotted, and provides the _P_-values, CIs, along with the **median unbiased estimate** It can be defined as such,

$$C V_{n}(\theta)=1-2\left|H_{n}(\theta)-0.5\right|=2 \min \left\{H_{n}(\theta), 1-H_{n}(\theta)\right\}$$

Its information counterpart, the surprisal function, can be constructed by taking the $-log_{2}$ of the _P_-value.[@chow2019asb; @greenland2019as; @Shannon1948-uq]

To view the surprisal function, we simply change the type to "`s`". 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(function1 <- ggcurve(data = intervalsdf[[1]], type = "s"))
```


We can also view the consonance distribution by changing the type to "`cdf`", which is a cumulative probability distribution. The point at which the curve reaches 50% is known as the "**median unbiased estimate**". It is the same estimate that is typically at the peak of the _P_-value curve from above. 


```{r echo=TRUE, fig.height=4.5, fig.width=6}
(function1s <- ggcurve(data = intervalsdf[[2]], type = "cdf", nullvalue = TRUE))
```

We can also get relevant statistics that show the range of values by using the `curve_table()` function. There are several formats that can be exported such as .docx, .ppt, and TeX. 

```{r echo=TRUE, fig.height=2, fig.width=4}
(x <- curve_table(data = intervalsdf[[1]], format = "image"))
```

# Comparing Functions

If we wanted to compare two studies to see the amount of "consonance", we could use the `curve_compare()` function to get a numerical output. 

First, we generate some more fake data

```{r echo=TRUE, fig.height=4.5, fig.width=6}
GroupA2 <- rnorm(500)
GroupB2 <- rnorm(500)
RandomData2 <- data.frame(GroupA2, GroupB2)
model <- lm(GroupA2 ~ GroupB2, data = RandomData2)
randomframe <- curve_gen(model, "GroupB2")
```

Once again, we'll plot this data with `ggcurve()`. We can also indicate whether we want certain interval estimates to be plotted in the function with the "`levels`" argument. If we wanted to plot the **50**%, **75**%, and **95**% intervals, we'd provide the argument this way: 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(function2 <- ggcurve(type = "c", randomframe[[1]], levels = c(0.50, 0.75, 0.95), nullvalue = TRUE))
```

Now that we have two datasets and two functions, we can compare them using the `curve_compare()` function.

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(curve_compare(
  data1 = intervalsdf[[1]], data2 = randomframe[[1]], type = "c",
  plot = TRUE, measure = "default", nullvalue = TRUE
))
```

This function will provide us with the area that is shared between the curve, along with a ratio of overlap to non-overlap. 

Another way to compare the functions is to use the `cowplot`  `plot_grid()` function. 

```{r echo=TRUE}
cowplot::plot_grid(function1, function2)
```

We can also do this for the surprisal function simply by changing `type` to "`s`".

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(curve_compare(
  data1 = intervalsdf[[1]], data2 = randomframe[[1]], type = "s",
  plot = TRUE, measure = "default", nullvalue = FALSE
))
```

It's clear that the outputs have changed and indicate far more overlap than before.


# Constructing Functions From Single Intervals

We can also take a set of confidence limits and use them to construct a consonance, surprisal, likelihood or deviance function using the `curve_rev()` function. This method is computed from the approximate normal distribution. 

Here, we'll use two epidemiological studies[@brown2017j; @brown2017jcp] that studied the impact of SSRI exposure in pregnant mothers, and the rate of autism in children. 

Both of these studies suggested a null effect of SSRI exposure on autism rates in children. 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
curve1 <- curve_rev(point = 1.7, LL = 1.1, UL = 2.6, type = "c", measure = "ratio", steps = 10000)
(ggcurve(data = curve1[[1]], type = "c", measure = "ratio", nullvalue = TRUE))
curve2 <- curve_rev(point = 1.61, LL = 0.997, UL = 2.59, type = "c", measure = "ratio", steps = 10000)
(ggcurve(data = curve2[[1]], type = "c", measure = "ratio", nullvalue = TRUE))
```

The null value is shown via the red line and it's clear that a large mass of the function is away from it. 

We can also see this by plotting the likelihood functions via the `curve_rev()` function. 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
lik1 <- curve_rev(point = 1.7, LL = 1.1, UL = 2.6, type = "l", measure = "ratio", steps = 10000)
(ggcurve(data = lik1[[1]], type = "l1", measure = "ratio", nullvalue = TRUE))
lik2 <- curve_rev(point = 1.61, LL = 0.997, UL = 2.59, type = "l", measure = "ratio", steps = 10000)
(ggcurve(data = lik2[[1]], type = "l1", measure = "ratio", nullvalue = TRUE))
```

We can also view the amount of agreement between the likelihood functions of these two studies. 

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(plot_compare(
  data1 = lik1[[1]], data2 = lik2[[1]], type = "l1", measure = "ratio", nullvalue = TRUE, title = "Brown et al. 2017. J Clin Psychiatry. vs. \nBrown et al. 2017. JAMA.",
  subtitle = "J Clin Psychiatry: OR = 1.7, 1/6.83 LI: LL = 1.1, UL = 2.6 \nJAMA: HR = 1.61, 1/6.83 LI: LL = 0.997, UL = 2.59", xaxis = expression(Theta ~ "= Hazard Ratio / Odds Ratio")
))
```

and the consonance functions

```{r echo=TRUE, fig.height=4.5, fig.width=6}
(plot_compare(
  data1 = curve1[[1]], data2 = curve2[[1]], type = "c", measure = "ratio", nullvalue = TRUE, title = "Brown et al. 2017. J Clin Psychiatry. vs. \nBrown et al. 2017. JAMA.",
  subtitle = "J Clin Psychiatry: OR = 1.7, 1/6.83 LI: LL = 1.1, UL = 2.6 \nJAMA: HR = 1.61, 1/6.83 LI: LL = 0.997, UL = 2.59", xaxis = expression(Theta ~ "= Hazard Ratio / Odds Ratio")
))
```


* * * 

# References

* * * 
