# Computing plausibility functions In flipr: Flexible Inference via Permutations in R

knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(dplyr) library(ggplot2) library(purrr) library(tibble) library(tidyr) library(flipr) load("../R/sysdata.rda")  ngrid_in <- 10 ngrid_out <- 100 nperms <- 100000 n1 <- 30 n2 <- 30 set.seed(1301) x1 <- rnorm(n1, mean = 0, sd = 1) x2 <- rnorm(n2, mean = 3, sd = 1) y1 <- rnorm(n1, mean = 0, sd = 1) y2 <- rnorm(n2, mean = 0, sd = 2) z1 <- rnorm(n1, mean = 0, sd = 1) z2 <- rnorm(n2, mean = 3, sd = 2)  The concept of plausibility functions pertains to assessing the$p$-value of a set of null hypotheses and to plot this$p$-value surface on the domain defined by the set of null hypotheses. The idea behind is that, if such a plausibility function is available, you can deduce from it point estimates or confidence interval estimates for parameters used to define the nulls or extract a single$p$-value for a specific null of interest [@martin2017; @fraser2019; @infanger2019]. In particular, there is another R package dedicated to plausibility functions called pvaluefunctions. ## Plausibility function for the mean null_spec <- function(y, parameters) { map(y, ~ .x - parameters) } stat_functions <- list(stat_t) stat_assignments <- list(delta = 1) pf <- PlausibilityFunction$new(
null_spec = null_spec,
stat_functions = stat_functions,
stat_assignments = stat_assignments,
x1, x2,
seed = 1234
)
pf$set_nperms(nperms) pf$set_point_estimate(mean(x2) - mean(x1))
pf$set_parameter_bounds( point_estimate = pf$point_estimate,
conf_level = pf$max_conf_level ) pf$set_grid(
parameters = pf$parameters, npoints = ngrid_in ) pf$set_alternative("two_tail")
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- rename(pf$grid, two_tail = pvalue) pf$set_alternative("left_tail")
pf$grid$pvalue <- NULL
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- bind_rows(
df,
rename(pf$grid, left_tail = pvalue) ) pf$set_alternative("right_tail")
pf$grid$pvalue <- NULL
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- bind_rows(
df,
rename(pf$grid, right_tail = pvalue) ) pf$set_grid(
parameters = pf$parameters, npoints = ngrid_out ) df_mean <- tibble( delta = pf$grid$delta, two_tail = approx(df$delta, df$two_tail, delta)$y,
left_tail = approx(df$delta, df$left_tail, delta)$y, right_tail = approx(df$delta, df$right_tail, delta)$y,
) %>%
pivot_longer(-delta)

df_mean %>%
ggplot(aes(delta, value, color = name)) +
geom_line() +
labs(
title = "P-value function for the mean",
subtitle = "t-statistic",
x = expression(delta),
y = "p-value",
color = "Type"
) +
geom_hline(
yintercept = 0.05,
color = "black",
linetype = "dashed"
) +
geom_vline(
xintercept = mean(x2) - mean(x1),
color = "black"
) +
geom_vline(
xintercept = stats::t.test(x2, x1, var.equal = TRUE)$conf.int, color = "black", linetype = "dashed" ) + scale_y_continuous(breaks = seq(0, 1, by = 0.05), limits = c(0, 1))  ## Plausibility function for the variance null_spec <- function(y, parameters) { map(y, ~ .x / parameters) } stat_functions <- list(stat_f) stat_assignments <- list(rho = 1) pf <- PlausibilityFunction$new(
null_spec = null_spec,
stat_functions = stat_functions,
stat_assignments = stat_assignments,
y1, y2,
seed = 1234
)
pf$set_nperms(nperms) pf$set_point_estimate(sd(y2) / sd(y1))
pf$set_parameter_bounds( point_estimate = pf$point_estimate,
conf_level = pf$max_conf_level ) pf$set_grid(
parameters = pf$parameters, npoints = ngrid_in ) pf$set_alternative("two_tail")
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- rename(pf$grid, two_tail = pvalue) pf$set_alternative("left_tail")
pf$grid$pvalue <- NULL
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- bind_rows(
df,
rename(pf$grid, left_tail = pvalue) ) pf$set_alternative("right_tail")
pf$grid$pvalue <- NULL
pf$evaluate_grid( grid = pf$grid,
ncores = 1
)
df <- bind_rows(
df,
rename(pf$grid, right_tail = pvalue) ) pf$set_grid(
parameters = pf$parameters, npoints = ngrid_out ) df_sd <- tibble( rho = pf$grid$rho, two_tail = approx(df$rho, df$two_tail, rho)$y,
left_tail = approx(df$rho, df$left_tail, rho)$y, right_tail = approx(df$rho, df$right_tail, rho)$y,
) %>%
pivot_longer(-rho)

df_sd %>%
ggplot(aes(rho, value, color = name)) +
geom_line() +
labs(
title = "P-value function for the standard deviation",
subtitle = "F-statistic",
x = expression(rho),
y = "p-value",
color = "Type"
) +
geom_hline(
yintercept = 0.05,
color = "black",
linetype = "dashed"
) +
geom_vline(
xintercept = sqrt(stats::var.test(y2, y1)$statistic), color = "black" ) + geom_vline( xintercept = sqrt(stats::var.test(y2, y1)$conf.int),
color = "black",
linetype = "dashed"
) +
scale_y_continuous(breaks = seq(0, 1, by = 0.05), limits = c(0, 1))


## Plausibility function for both mean and variance

Assume that we have two r.v. $X$ and $Y$ that differ in distribution only in their first two moments. Let $\mu_X$ and $\mu_Y$ be the means of $X$ and $Y$ respectively and $\sigma_X$ and $\sigma_Y$ be the standard deviations. We can therefore write

$$Y = \delta + \rho X.$$

In this case, we have

$$\begin{cases} \mu_Y = \delta + \rho \mu_X \ \sigma_Y^2 = \rho^2 \sigma_X^2 \end{cases} \Longleftrightarrow \begin{cases} \delta = \mu_Y - \frac{\sigma_Y}{\sigma_X} \mu_X \ \rho = \frac{\sigma_Y}{\sigma_X} \end{cases}$$

In the following example, we have $\delta = 3$ and $\rho = 2$.

null_spec <- function(y, parameters) {
map(y, ~ (.x - parameters[1]) / parameters[2])
}
stat_functions <- list(stat_t, stat_f)
stat_assignments <- list(delta = 1, rho = 2)

pf <- PlausibilityFunction$new( null_spec = null_spec, stat_functions = stat_functions, stat_assignments = stat_assignments, z1, z2, seed = 1234 ) pf$set_nperms(nperms)
pf$set_point_estimate(c( mean(z2) - sd(z2) / sd(z1) * mean(z1), sd(z2) / sd(z1) )) pf$set_parameter_bounds(
point_estimate = pf$point_estimate, conf_level = pf$max_conf_level
)

# Fisher combining function
pf$set_aggregator("fisher") pf$set_grid(
parameters = pf$parameters, npoints = ngrid_in ) pf$evaluate_grid(grid = pf$grid, ncores = 1) grid_in <- pf$grid
pf$set_grid( parameters = pf$parameters,
npoints = ngrid_out
)
if (requireNamespace("interp", quietly = TRUE)) {
Zout <- interp::interp(
x = grid_in$delta, y = grid_in$log_rho,
z = grid_in$pvalue, xo = sort(unique(pf$grid$delta)), yo = sort(unique(pf$grid$log_rho)) ) pf$grid$pvalue <- as.numeric(Zout$z)
} else
pf$grid$pvalue <- rep(NA, nrow(pf$grid)) df_fisher <- pf$grid

# Tippett combining function
pf$set_aggregator("tippett") pf$set_grid(
parameters = pf$parameters, npoints = ngrid_in ) pf$evaluate_grid(grid = pf$grid, ncores = 1) grid_in <- pf$grid
pf$set_grid( parameters = pf$parameters,
npoints = ngrid_out
)
if (requireNamespace("interp", quietly = TRUE)) {
Zout <- interp::interp(
x = grid_in$delta, y = grid_in$log_rho,
z = grid_in$pvalue, xo = sort(unique(pf$grid$delta)), yo = sort(unique(pf$grid$log_rho)) ) pf$grid$pvalue <- as.numeric(Zout$z)
} else
pf$grid$pvalue <- rep(NA, nrow(pf$grid)) df_tippett <- pf$grid

df_fisher %>%
ggplot(aes(delta, log_rho, z = pvalue)) +
geom_contour_filled(binwidth = 0.05) +
labs(
title = "Contour plot of the p-value surface",
subtitle = "Using Fisher's non-parametric combination",
x = expression(delta),
y = expression(log(rho)),
fill = "p-value"
) +
theme_minimal()

df_tippett %>%
ggplot(aes(delta, log_rho, z = pvalue)) +
geom_contour_filled(binwidth = 0.05) +
labs(
title = "Contour plot of the p-value surface",
subtitle = "Using Tippett's non-parametric combination",
x = expression(delta),
y = expression(log(rho)),
fill = "p-value"
) +
theme_minimal()


## Try the flipr package in your browser

Any scripts or data that you put into this service are public.

flipr documentation built on Aug. 23, 2023, 9:06 a.m.