Pearson similarity and distance from control

suppressPackageStartupMessages(library(magrittr))
suppressPackageStartupMessages(library(tidyverse))

1 Generate data

n <- 5 # observations
d <- 100 # dimensions
eps <- .5 # standard deviation

n_draws <- 10 # number of draws i.e. samples from from d-dim gaussian
n_biases <- 10 # number of shifts i.e. vectors to shift each sample
n_scalings <- 20 # number of scalings i.e. scaling factors for each shift
scaling_max <- 3
scalings <- scaling_max * seq(n_scalings) / n_scalings # scaling factors for each shift

Function to compute vector of cross-correlation values

corvec <- function(m) {cm <- cor(m); cm[upper.tri(cm)]}

1.1 Generate control distribution

a <- matrix(rnorm(n*d*10), n*10, d) * eps # sample from d-dim Gaussian

amean <- colMeans(a)

asd <- apply(a, 2, sd)

1.2 Generate treatment distributions and compute metrics

• Generate n_draws x n_biases x n_scalings sets of samples from a d-dim Gaussian, with n samples per set.
• Compute median replicate correlation for each set
• Compute distance from control for each set - this is just the scaling value

df <-
  map_df(1:n_draws,
         function(draw) {
           b <-
             scale(matrix(rnorm(n * d), n, d)) # sample from d-dim gaussion

           map_df(1:n_biases,
                  function(bias) {
                    bshift <- rnorm(d) # create shift vector

                    bshift <- d * bshift / norm(as.matrix(bshift))
                    
                    map_df(scalings,
                           function(bshift_scale) {
                             # scale the shift vector and then add to sample
                             bi <- scale(b, center = -bshift * bshift_scale, scale = FALSE) 
                             
                             # normalize the samples wrt control
                             bn <- scale(bi, center = amean, scale = asd)
                             
                             data.frame(draw,
                                        bias,
                                        bshift_scale, # distance from control
                                        medrepcor = median(corvec(t(bn))) # median replicate correlation
                                        )
                           })
                  })
         })

2 Plot data

2.1 Plot accross all shifts and draws

ggplot(df, aes(bshift_scale, medrepcor, group = interaction(draw, bias))) + 
  geom_line(alpha = 0.1) + 
  xlab("Distance from control") + 
  ylab("Median replicate correlation")

2.2 Plot accross all shifts grouped by draws

df %>%
  filter(draw %in% sample(1:n_draws, min(12, n_draws))) %>%
  ggplot(aes(bshift_scale, medrepcor, group = bias)) + 
  geom_line(alpha = 0.1) + 
  xlab("Distance from control") + 
  ylab("Median replicate correlation") + 
  facet_wrap(~draw, labeller = "label_both")

2.3 Plot accross all draws grouped by shifts

df %>%
  filter(bias %in% sample(1:n_biases, min(12, n_biases))) %>%
  ggplot(aes(bshift_scale, medrepcor, group = draw)) + 
  geom_line(alpha = 0.1) + 
  xlab("Distance from control") + 
  ylab("Median replicate correlation") + 
  facet_wrap(~bias, labeller = "label_both")

3 Generate data for 2 replicates and scale the shift

d <- 10 # dimensions

b0 <- matrix(rnorm(2 * d), 2, d) 

bshift <- rnorm(d)

bshift <- bshift / norm(as.matrix(bshift))

norm(as.matrix(bshift))

[1] 1

bd <- map_df(seq(20),
             function(bshift_scale) {
               bi <- sweep(b0, 2, bshift * bshift_scale, "+")
               
               bind_rows(
                 as.data.frame(t(b0)) %>% mutate(type = "orig"),
                 as.data.frame(t(bi)) %>% mutate(type = "shifted")
               ) %>%
                 mutate(bshift_scale = bshift_scale)
             })

4 Plot data

ggplot(bd, 
       aes(V1, V2, color = type)) + 
  geom_point() + facet_wrap(~bshift_scale, labeller = "label_both")

bd %>% 
  filter(type == "shifted") %>%
  group_by(bshift_scale) %>% 
  summarise(corval = cor(V1, V2), .groups = "keep") %>%
  ggplot(aes(bshift_scale, corval)) + geom_line() +
  xlab("Distance from control") + 
  ylab("Median replicate correlation")

unit_vec <- function(vec){
   return(vec/(sqrt(sum(vec**2,na.rm=TRUE))))
}

d <- 1000
x <- rnorm(d)
y <- rnorm(d)

df <- 
map_df(1:1000, 
       function(i) {
         b <- unit_vec(rnorm(d)) * sqrt(d)
         data.frame(corval_xy = cor(x+b, y+b), corval_bmeanxy = cor(b, (x+y)/2))
       }
       )

ggplot(df, aes(corval_xy, corval_bmeanxy)) + geom_point()