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

References