## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# number of pairs
n <- 1000
# matrix with different pairs in rows 
X <- matrix( rnorm( n * 2 ), nrow = n, ncol = 2 )

## -----------------------------------------------------------------------------
library(corsym)
Xe <- partial_order( X )

## -----------------------------------------------------------------------------
Xr <- randomize_order( Xe )

## -----------------------------------------------------------------------------
data <- data.frame(
  type = c('Random 1', 'Random 2', 'Extreme'),
  x = c(
    sum( X[,1] < X[,2] ),
    sum( Xr[,1] < Xr[,2] ),
    sum( Xe[,1] < Xe[,2] )
  ),
  n = n
)

## -----------------------------------------------------------------------------
order_bias_test( data )

## -----------------------------------------------------------------------------
# Each function returns 3 values:
# estimate, and lower and upper confidence intervals (CIs)
pearson( X )
pearson( Xr )
pearson( Xe )
corsym( X )
corsym( Xr )
corsym( Xe )

## -----------------------------------------------------------------------------
# sequence of forced order proportions between 0 and 1
# (q=0 is original data, q=1 is extreme order)
qs <- ( 0:10 ) / 10

# data frame to grow
data <- NULL
for ( q in qs ) {
    # force a proportion `q` of rows to have extreme order 
    Xp <- partial_order( X, q )
    # pearson estimates
    out <- pearson( Xp )

    # gather into data frame for plotting
    data_q <- data.frame(
        n = n,
	q = q,
        # counts for biased order test
        x = sum( Xp[,1] < Xp[,2] ),
	# Pearson estimates
	r = out[1],
	CIL = out[2],
	CIU = out[3]
    )
    data <- rbind( data, data_q )
}

# perform order bias tests on these data
data <- order_bias_test( data )

# inspect table
data

## -----------------------------------------------------------------------------
library(ggplot2)
ggplot( data, aes( x = q, y = f ) ) +
	geom_point() +
	theme_classic() +
	geom_abline( intercept = 0.5, slope = 0.5, linetype = 'dashed', color = 'gray' )

## -----------------------------------------------------------------------------
ggplot( data, aes( x = q, y = r ) ) +
	geom_point() +
	theme_classic() +
   	geom_errorbar( aes( ymin = CIL, ymax = CIU ) ) +
	geom_hline( yintercept = 0, linetype = 'dashed', color = 'gray' )

## -----------------------------------------------------------------------------
# Bonferroni threshold
pcut <- 0.05 / nrow( data )
pcut

ggplot( data, aes( x = q, y = -log10( pval ) ) ) +
	geom_point() +
	theme_classic() +
	geom_hline( yintercept = -log10( pcut ), linetype = 'dashed', color = 'gray' )