# Too high statistical significance on random data

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)


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