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, 20131 displays why.

For a more extensive discussion of these concepts, see the following references.113

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.

tibble::tibble(intervalsdf[[1]])
#> # A tibble: 10,000 x 1
#>    `intervalsdf[[1… $upper.limit $intrvl.width $intrvl.level  $cdf $pvalue
#>               <dbl>        <dbl>         <dbl>         <dbl> <dbl>   <dbl>
#>  1           -0.113       -0.113     0              0        0.5     1    
#>  2           -0.113       -0.113     0.0000154      0.0001   0.500   1.00 
#>  3           -0.113       -0.113     0.0000309      0.0002   0.500   1.00 
#>  4           -0.113       -0.113     0.0000463      0.000300 0.500   1.00 
#>  5           -0.113       -0.113     0.0000617      0.0004   0.500   1.00 
#>  6           -0.113       -0.113     0.0000772      0.0005   0.500   1.00 
#>  7           -0.113       -0.113     0.0000926      0.000600 0.500   0.999
#>  8           -0.113       -0.113     0.000108       0.0007   0.500   0.999
#>  9           -0.113       -0.112     0.000123       0.0008   0.500   0.999
#> 10           -0.113       -0.112     0.000139       0.0009   0.500   0.999
#> # … with 9,990 more rows, and 1 more variable: $svalue <dbl>

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.3,14,15

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

Lower Limit

Upper Limit

Interval Width

Interval Level (%)

CDF

P-value

S-value (bits)

-0.132

-0.093

0.039

25.0

0.625

0.750

0.415

-0.154

-0.071

0.083

50.0

0.750

0.500

1.000

-0.183

-0.042

0.142

75.0

0.875

0.250

2.000

-0.192

-0.034

0.158

80.0

0.900

0.200

2.322

-0.201

-0.024

0.177

85.0

0.925

0.150

2.737

-0.214

-0.011

0.203

90.0

0.950

0.100

3.322

-0.233

0.008

0.242

95.0

0.975

0.050

4.322

-0.251

0.026

0.276

97.5

0.988

0.025

5.322

-0.271

0.046

0.318

99.0

0.995

0.010

6.644

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
))
#> [1] "AUC = Area Under the Curve"
#> [[1]]
#> 
#> 
#>  AUC 1   AUC 2   Shared AUC   AUC Overlap (%)   Overlap:Non-Overlap AUC Ratio
#> ------  ------  -----------  ----------------  ------------------------------
#>  0.098   0.073        0.024            16.309                           0.195
#> 
#> [[2]]

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

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
))
#> [1] "AUC = Area Under the Curve"
#> [[1]]
#> 
#> 
#>  AUC 1   AUC 2   Shared AUC   AUC Overlap (%)   Overlap:Non-Overlap AUC Ratio
#> ------  ------  -----------  ----------------  ------------------------------
#>  3.947   1.531        1.531            38.801                           0.634
#> 
#> [[2]]

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

Survival Modeling

Here, we’ll look at how to create consonance functions from the coefficients of predictors of interest in a Cox regression model.

We’ll use the carData package for this. Fox & Weisberg, 2018 describe the dataset elegantly in their paper,

The Rossi data set in the carData package contains data from an experimental study of recidivism of 432 male prisoners, who were observed for a year after being released from prison (Rossi et al., 1980). The following variables are included in the data; the variable names are those used by Allison (1995), from whom this example and variable descriptions are adapted:

week: week of first arrest after release, or censoring time.

arrest: the event indicator, equal to 1 for those arrested during the period of the study and 0 for those who were not arrested.

fin: a factor, with levels “yes” if the individual received financial aid after release from prison, and “no” if he did not; financial aid was a randomly assigned factor manipulated by the researchers.

age: in years at the time of release.

race: a factor with levels “black” and “other”.

wexp: a factor with levels “yes” if the individual had full-time work experience prior to incarceration and “no” if he did not.

mar: a factor with levels “married” if the individual was married at the time of release and “not married” if he was not.

paro: a factor coded “yes” if the individual was released on parole and “no” if he was not.

prio: number of prior convictions.

educ: education, a categorical variable coded numerically, with codes 2 (grade 6 or less), 3 (grades 6 through 9), 4 (grades 10 and 11), 5 (grade 12), or 6 (some post-secondary).

emp1–emp52: factors coded “yes” if the individual was employed in the corresponding week of the study and “no” otherwise.

We read the data file into a data frame, and print the first few cases (omitting the variables emp1 – emp52, which are in columns 11–62 of the data frame):

library(carData)
Rossi[1:5, 1:10]
#>   week arrest fin age  race wexp         mar paro prio educ
#> 1   20      1  no  27 black   no not married  yes    3    3
#> 2   17      1  no  18 black   no not married  yes    8    4
#> 3   25      1  no  19 other  yes not married  yes   13    3
#> 4   52      0 yes  23 black  yes     married  yes    1    5
#> 5   52      0  no  19 other  yes not married  yes    3    3

Thus, for example, the first individual was arrested in week 20 of the study, while the fourth individual was never rearrested, and hence has a censoring time of 52. Following Allison, a Cox regression of time to rearrest on the time-constant covariates is specified as follows:

library(survival)
mod.allison <- coxph(Surv(week, arrest) ~
fin + age + race + wexp + mar + paro + prio,
data = Rossi
)
mod.allison
#> Call:
#> coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp + 
#>     mar + paro + prio, data = Rossi)
#> 
#>                    coef exp(coef) se(coef)      z       p
#> finyes         -0.37942   0.68426  0.19138 -1.983 0.04742
#> age            -0.05744   0.94418  0.02200 -2.611 0.00903
#> raceother      -0.31390   0.73059  0.30799 -1.019 0.30812
#> wexpyes        -0.14980   0.86088  0.21222 -0.706 0.48029
#> marnot married  0.43370   1.54296  0.38187  1.136 0.25606
#> paroyes        -0.08487   0.91863  0.19576 -0.434 0.66461
#> prio            0.09150   1.09581  0.02865  3.194 0.00140
#> 
#> Likelihood ratio test=33.27  on 7 df, p=2.362e-05
#> n= 432, number of events= 114

Now that we have our Cox model object, we can use the curve_surv() function to create the function.

If we wanted to create a function for the coefficient of prior convictions, then we’d do so like this:

z <- curve_surv(mod.allison, "prio")

Then we could plot our consonance curve and density and also produce a table of relevant statistics. Because we’re working with ratios, we’ll set the measure argument in ggcurve() to “ratio”.

ggcurve(z[[1]], measure = "ratio", nullvalue = TRUE)

ggcurve(z[[2]], type = "cd", measure = "ratio", nullvalue = TRUE)

curve_table(z[[1]], format = "image")

Lower Limit

Upper Limit

Interval Width

Interval Level (%)

CDF

P-value

S-value (bits)

1.086

1.106

0.020

25.0

0.625

0.750

0.415

1.075

1.117

0.042

50.0

0.750

0.500

1.000

1.060

1.133

0.072

75.0

0.875

0.250

2.000

1.056

1.137

0.080

80.0

0.900

0.200

2.322

1.052

1.142

0.090

85.0

0.925

0.150

2.737

1.045

1.149

0.103

90.0

0.950

0.100

3.322

1.036

1.159

0.123

95.0

0.975

0.050

4.322

1.028

1.168

0.141

97.5

0.988

0.025

5.322

1.018

1.180

0.162

99.0

0.995

0.010

6.644

We could also construct a function for another predictor such as age

x <- curve_surv(mod.allison, "age")
ggcurve(x[[1]], measure = "ratio")

ggcurve(x[[2]], type = "cd", measure = "ratio")

curve_table(x[[1]], format = "image")

Lower Limit

Upper Limit

Interval Width

Interval Level (%)

CDF

P-value

S-value (bits)

0.938

0.951

0.013

25.0

0.625

0.750

0.415

0.930

0.958

0.028

50.0

0.750

0.500

1.000

0.921

0.968

0.048

75.0

0.875

0.250

2.000

0.918

0.971

0.053

80.0

0.900

0.200

2.322

0.915

0.975

0.060

85.0

0.925

0.150

2.737

0.911

0.979

0.068

90.0

0.950

0.100

3.322

0.904

0.986

0.081

95.0

0.975

0.050

4.322

0.899

0.992

0.093

97.5

0.988

0.025

5.322

0.892

0.999

0.107

99.0

0.995

0.010

6.644

That’s a very quick look at creating functions from Cox regression models.

Meta-Analysis

Here, we’ll use an example dataset taken from the metafor website, which also comes preloaded with the metafor package.

library(metafor)
#> Loading required package: Matrix
#> Loading 'metafor' package (version 2.1-0). For an overview 
#> and introduction to the package please type: help(metafor).
dat.hine1989
#>   study         source n1i n2i ai ci
#> 1     1  Chopra et al.  39  43  2  1
#> 2     2       Mogensen  44  44  4  4
#> 3     3    Pitt et al. 107 110  6  4
#> 4     4   Darby et al. 103 100  7  5
#> 5     5 Bennett et al. 110 106  7  3
#> 6     6 O'Brien et al. 154 146 11  4

I will quote Wolfgang here, since he explains it best,

"As described under help(dat.hine1989), variables n1i and n2i are the number of patients in the lidocaine and control group, respectively, and ai and ci are the corresponding number of deaths in the two groups. Since these are 2×2 table data, a variety of different outcome measures could be used for the meta-analysis, including the risk difference, the risk ratio (relative risk), and the odds ratio (see Table III). Normand (1999) uses risk differences for the meta-analysis, so we will proceed accordingly. We can calculate the risk differences and corresponding sampling variances with:

dat <- escalc(measure = "RD", n1i = n1i, n2i = n2i, ai = ai, ci = ci, data = dat.hine1989)
dat
#>   study         source n1i n2i ai ci     yi     vi 
#> 1     1  Chopra et al.  39  43  2  1 0.0280 0.0018 
#> 2     2       Mogensen  44  44  4  4 0.0000 0.0038 
#> 3     3    Pitt et al. 107 110  6  4 0.0197 0.0008 
#> 4     4   Darby et al. 103 100  7  5 0.0180 0.0011 
#> 5     5 Bennett et al. 110 106  7  3 0.0353 0.0008 
#> 6     6 O'Brien et al. 154 146 11  4 0.0440 0.0006

"Note that the yi values are the risk differences in terms of proportions. Since Normand (1999) provides the results in terms of percentages, we can make the results directly comparable by multiplying the risk differences by 100 (and the sampling variances by \(100^{2}\)):

dat$yi <- dat$yi * 100
dat$vi <- dat$vi * 100^2

We can fit a fixed-effects model with the following

fe <- rma(yi, vi, data = dat, method = "FE")

Now that we have our metafor object, we can compute the consonance function using the curve_meta() function.

fecurve <- curve_meta(fe)

Now we can graph our function.

ggcurve(fecurve[[1]], nullvalue = TRUE)

We used a fixed-effects model here, but if we wanted to use a random-effects model, we could do so with the following, which will use a restricted maximum likelihood estimator for the random-effects model

re <- rma(yi, vi, data = dat, method = "REML")

And then we could use curve_meta() to get the relevant list

recurve <- curve_meta(re)

Now we can plot our object.

ggcurve(recurve[[1]], nullvalue = TRUE)

We could also compare our two models to see how much consonance/overlap there is

curve_compare(fecurve[[1]], recurve[[1]], plot = TRUE)
#> [1] "AUC = Area Under the Curve"
#> [[1]]
#> 
#> 
#>  AUC 1   AUC 2   Shared AUC   AUC Overlap (%)   Overlap:Non-Overlap AUC Ratio
#> ------  ------  -----------  ----------------  ------------------------------
#>  2.085   2.085        2.085               100                             Inf
#> 
#> [[2]]

The results are practically the same and we cannot actually see any difference, and the AUC % overlap also indicates this.

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 studies16,17 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")
))

The Bootstrap and Consonance Functions

Some authors have shown that the bootstrap distribution is equal to the confidence distribution because it meets the definition of a consonance distribution.1,18,19 The bootstrap distribution and the asymptotic consonance distribution would be defined as:

\[H_{n}(\theta)=1-P\left(\hat{\theta}-\hat{\theta}^{*} \leq \hat{\theta}-\theta | \mathbf{x}\right)=P\left(\hat{\theta}^{*} \leq \theta | \mathbf{x}\right)\]

Certain bootstrap methods such as the BCa method and t-bootstrap method also yield second order accuracy of consonance distributions.

\[H_{n}(\theta)=1-P\left(\frac{\hat{\theta}^{*}-\hat{\theta}}{\widehat{S E}^{*}\left(\hat{\theta}^{*}\right)} \leq \frac{\hat{\theta}-\theta}{\widehat{S E}(\hat{\theta})} | \mathbf{x}\right)\]

Here, I demonstrate how to use these particular bootstrap methods to arrive at consonance curves and densities.

We’ll use the Iris dataset and construct a function that’ll yield a parameter of interest.

The Nonparametric Bootstrap

iris <- datasets::iris
foo <- function(data, indices) {
  dt <- data[indices, ]
  c(
    cor(dt[, 1], dt[, 2], method = "p")
  )
}

We can now use the curve_boot() method to construct a function. The default method used for this function is the “Bca” method provided by the bcaboot package.19

I will suppress the output of the function because it is unnecessarily long. But we’ve placed all the estimates into a list object called y.

The first item in the list will be the consonance distribution constructed by typical means, while the third item will be the bootstrap approximation to the consonance distribution.

ggcurve(data = y[[1]], nullvalue = TRUE)

ggcurve(data = y[[3]], nullvalue = TRUE)

We can also print out a table for TeX documents

(gg <- curve_table(data = y[[1]], format = "image"))

Lower Limit

Upper Limit

Interval Width

Interval Level (%)

CDF

P-value

S-value (bits)

-0.142

-0.093

0.048

25

0.625

0.75

0.415

-0.169

-0.067

0.102

50

0.750

0.50

1.000

-0.205

-0.031

0.174

75

0.875

0.25

2.000

-0.214

-0.021

0.194

80

0.900

0.20

2.322

-0.266

0.031

0.296

95

0.975

0.05

4.322

-0.312

0.077

0.389

99

0.995

0.01

6.644

More bootstrap replications will lead to a smoother function. But for now, we can compare these two functions to see how similar they are.

plot_compare(y[[1]], y[[3]])

If we wanted to look at the bootstrap standard errors, we could do so by loading the fifth item in the list

knitr::kable(y[[5]])
theta sdboot z0 a sdjack
est -0.1175698 0.0755961 0.0576844 0.0304863 0.075694
jsd 0.0000000 0.0010234 0.0274023 0.0000000 0.000000

where in the top row, theta is the point estimate, and sdboot is the bootstrap estimate of the standard error, sdjack is the jacknife estimate of the standard error. z0 is the bias correction value and a is the acceleration constant.

The values in the second row are essentially the internal standard errors of the estimates in the top row.


The densities can also be calculated accurately using the t-bootstrap method. Here we use a different dataset to show this

library(Lock5Data)
dataz <- data(CommuteAtlanta)
func <- function(data, index) {
  x <- as.numeric(unlist(data[1]))
  y <- as.numeric(unlist(data[2]))
  return(mean(x[index]) - mean(y[index]))
}

Our function is a simple mean difference. This time, we’ll set the method to “t” for the t-bootstrap method

z <- curve_boot(data = CommuteAtlanta, func = func, method = "t", replicates = 2000, steps = 1000)
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
ggcurve(data = z[[1]], nullvalue = FALSE)

ggcurve(data = z[[2]], type = "cd", nullvalue = FALSE)

The consonance curve and density are nearly identical. With more bootstrap replications, they are very likely to converge.

(zz <- curve_table(data = z[[1]], format = "image"))

Lower Limit

Upper Limit

Interval Width

Interval Level (%)

CDF

P-value

S-value (bits)

-39.400

-39.075

0.325

25.0

0.625

0.750

0.415

-39.611

-38.876

0.735

50.0

0.750

0.500

1.000

-39.873

-38.608

1.265

75.0

0.875

0.250

2.000

-39.932

-38.530

1.402

80.0

0.900

0.200

2.322

-40.026

-38.456

1.570

85.0

0.925

0.150

2.737

-40.118

-38.354

1.763

90.0

0.950

0.100

3.322

-40.294

-38.174

2.120

95.0

0.975

0.050

4.322

-40.442

-38.026

2.416

97.5

0.988

0.025

5.322

-40.636

-37.806

2.830

99.0

0.995

0.010

6.644

The Parametric Bootstrap

For the examples above, we mainly used nonparametric bootstrap methods. Here I show an example using the parametric Bca bootstrap and the results it yields.

First, we’ll load our data again and set our function.

data(diabetes, package = "bcaboot")
X <- diabetes$x
y <- scale(diabetes$y, center = TRUE, scale = FALSE)
lm.model <- lm(y ~ X - 1)
mu.hat <- lm.model$fitted.values
sigma.hat <- stats::sd(lm.model$residuals)
t0 <- summary(lm.model)$adj.r.squared
y.star <- sapply(mu.hat, rnorm, n = 1000, sd = sigma.hat)
tt <- apply(y.star, 1, function(y) summary(lm(y ~ X - 1))$adj.r.squared)
b.star <- y.star %*% X

Now, we’ll use the same function, but set the method to “bcapar” for the parametric method.

df <- curve_boot(method = "bcapar", t0 = t0, tt = tt, bb = b.star)

Now we can look at our outputs.

ggcurve(df[[1]], nullvalue = FALSE)

ggcurve(df[[3]], nullvalue = FALSE)

We can compare the functions to see how well the bootstrap approximations match up

plot_compare(df[[1]], df[[3]])

We can also look at the density function

ggcurve(df[[5]], type = "cd", nullvalue = FALSE)

That concludes our demonstration of the bootstrap method to approximate consonance functions.

Using Profile Likelihoods

For this last example, we’ll explore the curve_lik() function, which can help generate profile likelihood functions, and deviance statistics with the help of the ProfileLikelihood package.

library(ProfileLikelihood)
#> Loading required package: nlme
#> Loading required package: MASS

We’ll use a simple example taken directly from the ProfileLikelihood documentation where we’ll calculate the likelihoods from a glm model

data(dataglm)
xx <- profilelike.glm(y ~ x1 + x2,
  data = dataglm, profile.theta = "group",
  family = binomial(link = "logit"), length = 500, round = 2
)
#> Warning message: provide lo.theta and hi.theta

Then, we’ll use curve_lik() on the object that the ProfileLikelihood package created.

lik <- curve_lik(xx, dataglm)
tibble::tibble(lik[[1]])
#> # A tibble: 500 x 1
#>    `lik[[1]]`$values $likelihood $loglikelihood  $support $deviancestat
#>                <dbl>       <dbl>          <dbl>     <dbl>         <dbl>
#>  1             -1.41    9.26e-21          -9.79 0.0000560          9.79
#>  2             -1.40    1.00e-20          -9.71 0.0000606          9.71
#>  3             -1.39    1.08e-20          -9.63 0.0000655          9.63
#>  4             -1.38    1.17e-20          -9.56 0.0000708          9.56
#>  5             -1.37    1.26e-20          -9.48 0.0000765          9.48
#>  6             -1.35    1.37e-20          -9.40 0.0000826          9.40
#>  7             -1.34    1.47e-20          -9.32 0.0000892          9.32
#>  8             -1.33    1.59e-20          -9.25 0.0000963          9.25
#>  9             -1.32    1.72e-20          -9.17 0.000104           9.17
#> 10             -1.31    1.85e-20          -9.10 0.000112           9.10
#> # … with 490 more rows

Next, we’ll plot three functions, the relative likelihood, the log-likelihood, the likelihood, and the deviance function.

ggcurve(lik[[1]], type = "l1", nullvalue = TRUE)

ggcurve(lik[[1]], type = "l2")

ggcurve(lik[[1]], type = "l3")

ggcurve(lik[[1]], type = "d")

The obvious advantage of using reduced likelihoods is that they are free of nuisance parameters

\[L_{t_{n}}(\theta)=f_{n}\left(F_{n}^{-1}\left(H_{p i v}(\theta)\right)\right)\left|\frac{\partial}{\partial t} \psi\left(t_{n}, \theta\right)\right|=h_{p i v}(\theta)\left|\frac{\partial}{\partial t} \psi(t, \theta)\right| /\left.\left|\frac{\partial}{\partial \theta} \psi(t, \theta)\right|\right|_{t=t_{n}}\] thus, giving summaries of the data that can be incorporated into combined analyses.


References


1. Xie M-g, Singh K. Confidence Distribution, the Frequentist Distribution Estimator of a Parameter: A Review. International Statistical Review. 2013;81(1):3-39. doi:10.1111/insr.12000

2. Birnbaum A. A unified theory of estimation, I. The Annals of Mathematical Statistics. 1961;32(1):112-135. doi:10.1214/aoms/1177705145

3. Chow ZR, Greenland S. Semantic and Cognitive Tools to Aid Statistical Inference: Replace Confidence and Significance by Compatibility and Surprise. arXiv:190908579 [statME]. September 2019. http://arxiv.org/abs/1909.08579.

4. Fraser DAS. P-Values: The Insight to Modern Statistical Inference. Annual Review of Statistics and Its Application. 2017;4(1):1-14. doi:10.1146/annurev-statistics-060116-054139

5. Fraser DAS. The P-value function and statistical inference. The American Statistician. 2019;73(sup1):135-147. doi:10.1080/00031305.2018.1556735

6. Poole C. Beyond the confidence interval. American Journal of Public Health. 1987;77(2):195-199. doi:10.2105/AJPH.77.2.195

7. Poole C. Confidence intervals exclude nothing. American Journal of Public Health. 1987;77(4):492-493. doi:10.2105/ajph.77.4.492

8. Schweder T, Hjort NL. Confidence and Likelihood*. Scand J Stat. 2002;29(2):309-332. doi:10.1111/1467-9469.00285

9. Schweder T, Hjort NL. Confidence, Likelihood, Probability: Statistical Inference with Confidence Distributions. Cambridge University Press; 2016.

10. Singh K, Xie M, Strawderman WE. Confidence distribution (CD) – distribution estimator of a parameter. August 2007. http://arxiv.org/abs/0708.0976.

11. Sullivan KM, Foster DA. Use of the confidence interval function. Epidemiology. 1990;1(1):39-42. doi:10.1097/00001648-199001000-00009

12. Whitehead J. The case for frequentism in clinical trials. Statistics in Medicine. 1993;12(15-16):1405-1413. doi:10.1002/sim.4780121506

13. Rothman KJ, Greenland S, Lash TL. Precision and statistics in epidemiologic studies. In: Rothman KJ, Greenland S, Lash TL, eds. Modern Epidemiology. 3rd ed. Lippincott Williams & Wilkins; 2008:148-167.

14. Greenland S. Valid P-values behave exactly as they should: Some misleading criticisms of P-values and their resolution with S-values. The American Statistician. 2019;73(sup1):106-114. doi:10.1080/00031305.2018.1529625

15. Shannon CE. A mathematical theory of communication. The Bell System Technical Journal. 1948;27(3):379-423. doi:10.1002/j.1538-7305.1948.tb01338.x

16. Brown HK, Ray JG, Wilton AS, Lunsky Y, Gomes T, Vigod SN. Association between serotonergic antidepressant use during pregnancy and autism spectrum disorder in children. JAMA. 2017;317(15):1544-1552. doi:10.1001/jama.2017.3415

17. Brown HK, Hussain-Shamsy N, Lunsky Y, Dennis C-LE, Vigod SN. The association between antenatal exposure to selective serotonin reuptake inhibitors and autism: A systematic review and meta-analysis. The Journal of Clinical Psychiatry. 2017;78(1):e48-e58. doi:10.4088/JCP.15r10194

18. Efron B, Tibshirani RJ. An Introduction to the Bootstrap. CRC Press; 1994.

19. Efron B, Narasimhan B. The automatic construction of bootstrap confidence intervals. October 2018:17.