1
$\begingroup$

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)
$\endgroup$
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.