Too high statistical significance on random data

Question

I am trying to understand the stats problem that I encountered. Use reproducible R code below. All libraries enlisted are needed.

What I do is: - create two tables $101 \times 1000$, where there are $100$ inputs (a vector of delayed timeseries points), and $1$ output (a one step ahead series point).

The tables are independent, and each row in a table is independent.

• Cross join the tables: a $1,000,000$ row table appears.

• Select $10,000$ random rows from resulting $1,000,000$-row table.

• Calculate $1000^2$ Euclidean distance values for input vectors, and another $1000^2$ Euclidean distances for output. When calculating the input vector distances make feature selection of input space dimensions, where the best set of input dimensions corresponds to a highest positive correlation between: input Euclidean distance and output Euclidean distance.

• Analyze the correlation significance, using adjusted alpha.

Problem: Even when I use rnorm to generate all independent values, I always get significant correlations, much more significant than my significance level.

Question: I wonder what went wrong here. What assumptions were misfit?

rm(list = ls()); gc()

library(data.table)
library(FSelector)
library(magrittr)
library(ggplot2)
library(fNonlinear)

val_numb <- 202000

x <- rnorm(val_numb)

#x <- sin(seq(0.1, val_numb/10, 0.1))

#x <- as.numeric(tentSim(n = val_numb, n.skip = 0, parms = c(a = 2), start = runif(1), doplot = FALSE))

#x <- henonSim(n = val_numb, n.skip = 0, parms = c(a = 1.4, b = 0.3), start = runif(2), doplot = FALSE); x <- as.numeric(as.matrix(x)[, 1]); plot(x[1:500], type = 'l')

#x <- ikedaSim(n = 1000, n.skip = 100, parms = c(a = 0.4, b = 6.0, c = 0.9), start = runif(2), doplot = FALSE)

#x <- as.numeric(logisticSim(n = val_numb, n.skip = 0, parms = c(r = 3.9), start = runif(1), doplot = F))

#x <- lorentzSim(times = seq(0.1, 2200, by = 0.1), parms = c(sigma = 16, r = 45.92, b = 4), start = c(-14, -13, 47), doplot = F); x <- as.numeric(as.matrix(x)[, 3]); plot(x[1:5000], type = 'l')

#x <- roesslerSim(times = seq(0, 100, by = 0.01), parms = c(a = 0.2, b = 0.2, c = 8.0), start = c(-1.894, -9.920, 0.0250), doplot = F)


# fill arrays

inputs <- 100

a <- as.data.table(t(matrix(head(x, length(x) / 2), nrow = inputs + 1, ncol = length(x) / 2 / (inputs + 1))))

colnames(a) <- c(paste0('input_', 1:inputs), 'output')

b <- as.data.table(t(matrix(tail(x, length(x) / 2), nrow = inputs + 1, ncol = length(x) / 2 / (inputs + 1))))

colnames(b) <- c(paste0('input_', 1:inputs), 'output')


CJ.dt = function(X,Y) {
    stopifnot(is.data.table(X),is.data.table(Y))
    k = NULL
    X = X[, c(k=1, .SD)]
    setkey(X, k)
    Y = Y[, c(k=1, .SD)]
    setkey(Y, NULL)
    X[Y, allow.cartesian=TRUE][, k := NULL][]
}

ab <- CJ.dt(a, b)

ab[, output_dist:= abs(output - i.output)]


###

rows <- sample(nrow(ab), 10000, replace = F)
ab <- ab[rows, ]

cor_method <- 'pearson' # 'kendall'

global_alpha <- 0.01

n_tests <- 0

corr_func <- function(subset){

    ab[, input_dist := 
        Map(
            function(x, y) (x - y) ^ 2, 
            .SD[, subset, with = F], 
            .SD[, paste0('i.', subset), with = F]
        ) %>% 
        Reduce(`+`, .) %>% 
        sqrt
       ]

    corr <- cor(x = ab[, input_dist], 
               y = ab[, output_dist],
               method = cor_method)

    n_tests <<- n_tests + 1
    print(subset)
    print(corr)
    return(corr)

}

subset <- forward.search(attributes = paste0('input_', 1:inputs), eval.fun = corr_func)

ab[, input_dist := 
    Map(
        function(x, y) (x - y) ^ 2, 
        .SD[, subset, with = F], 
        .SD[, paste0('i.', subset), with = F]
    ) %>% 
    Reduce(`+`, .) %>% 
    sqrt
   ]

print(
    cor.test(x = ab[, input_dist], 
          y = ab[, output_dist],
          method = cor_method,
          alternative = "greater",
          exact = NULL, 
          continuity = FALSE)$estimate
)

print(
    cor.test(x = ab[, input_dist], 
      y = ab[, output_dist],
      method = cor_method,
      alternative = "greater",
      exact = NULL, 
      continuity = FALSE)$p.value
)

local_alpha <- global_alpha / n_tests

print(local_alpha)

ggplot(ab, aes(x = input_dist, y = output_dist)) +
    geom_point(alpha = 0.05, size = 2) +
    geom_smooth(method = 'lm', level = 0.000001)

Answer 1

I will try to give a possible answer to my question. If you feel I am misled, correct me.

Given that I run a feature selection algorithm, I assume it works perfectly and always converges. That means that I get the highest correlation coefficient based on the selected dimensions my data are projected to.

The total number of possible combinations of dimensions: without repetitions and without misplacements (order does not matter) would be n! / (k! (n - k)!), which yelds 1.267651e+30 for 1 to 100 combinations from 100 available values. This is the actual number of tests that I should use, not the actual number I get after feature selection iterations.

Remember that we assume the most conservative case, where the feature selection reaches the global optimum. Then I make sure the best correlation on random data exceeds a critical value with probability global alpha - hence the local alpha is in my case 0.01 / 1.267651e+30.

After calculating the critical value of the Fisher-transformed value of the Pearson correlation coefficient corresponding to set alpha, and transforming it back to an original scale, I get 0.359. The following simulations on the random data that I run suggest that I can typically expect 0.2-0.3 observed correlation.

Thus the answer is: I got the wrong number of tests to adjust the alpha for multiplicity of the testing.

A code snippet that fills in the gap:

global_alpha <- 0.01

n_tests <- round(sum(factorial(inputs)/(factorial(1:inputs) * factorial(inputs - 1:inputs))), 0)

local_alpha <- global_alpha / n_tests

cor_se <- 1 / sqrt(nrow(ab) - 3)

crit_cor <- tanh(qnorm(p = local_alpha, mean = 0, sd = cor_se, lower.tail = F, log.p = F))

I apologize for formatting and possible typos in advance.

Alexey