# Examples in R

```
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.

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

```
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.

```
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.

`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`

".

`(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`

".

`(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.

`(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.

`(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

```
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:

`(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.

```
(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.

`cowplot::plot_grid(function1, function2)`

We can also do this for the surprisal function simply by changing `type`

to "`s`

".

```
(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.

```
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.

```
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.

```
(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

```
(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")
))
```