0
$\begingroup$

Originally, I started out with two .csv files, each 7 GB in size. Once read in into R they each hold about 20 GB of RAM. That is why I calculated the files df1 and df2, which only contain a mean and the standard deviation. Both files have about 8 million rows. The sample size for every mean/standard deviation pair is 157.

df1:

ID       Mean         Standard deviation
16091071 0.1551044586 0.0120334914
16091086 0.1528095541 0.0125274201
16091147 0.3395656051 0.0186907447
16091154 0.2788547771 0.017261902
16091227 0.250456051 0.0176726877
16091236 0.2776751592 0.0175430484

df2:

ID       Mean         Standard deviation
16091071 0.0999431847 0.0106193089
16091086 0.3864509554 0.0181384583
16091147 0.3864509554 0.0181384583
16091154 0.0999431847 0.0106193089
16091227 0.0999431847 0.0106193089
16091236 0.0999431847 0.0106193089

The question now is how to do inferential statistics in R only based on mean and standard deviation data row by row while simultaneously still considering the multiple testing problem.

Improve this question
$\endgroup$
14
  • $\begingroup$ What is it that you want to test, and do you have the sample sizes that yielded each mean and standard deviation? $\endgroup$
    – Dave
    Commented Jul 6, 2020 at 10:28
  • $\begingroup$ Oh, sure, of course. I should have given the sample size. For both, df1 and df2, sample size is 157 for every row. At the end, I would like to have 8 million p-values. So whether or not the data sets are significantly different row by row. $\endgroup$
    – PolII
    Commented Jul 6, 2020 at 10:34
  • $\begingroup$ Significantly different how? Also, your problem of dealing with eight-million p-values warrants its own question. Always remember the XKCD on jelly beans: xkcd.com/882. $\endgroup$
    – Dave
    Commented Jul 6, 2020 at 10:44
  • $\begingroup$ I am actually well aware of the jelly bean analogy. In this case though the question behind the statistics should be treated as it is done with the p-values of a genome-wide association study (GWAS) with post-hoc Bonferroni correction. $\endgroup$
    – PolII
    Commented Jul 6, 2020 at 11:20
  • $\begingroup$ Bonferroni on eight million p-values? That means you’re testing at $\alpha= 0.00000000625$ for the usual $5\%$ type I error rate. You’re sure you want to do that? $\endgroup$
    – Dave
    Commented Jul 6, 2020 at 11:25

1 Answer 1

Reset to default
0
$\begingroup$

Here is my solution with R using your attached dataframes df1 and df2.

library(dplyr)
#> 
#> Attache Paket: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# Dataframes: 
df1 <- data.frame(ID = c(16091071, 16091086, 16091147, 16091154, 16091227, 16091236), 
                  Mean = c(0.1551044586, 0.1528095541, 0.3395656051, 0.2788547771, 
                           0.250456051, 0.2776751592), 
                  Standard_deviation = c(0.0120334914, 0.0125274201, 0.0186907447, 
                                         0.017261902, 0.0176726877, 0.0175430484)
)

df2 <- data.frame(ID = c(16091071, 16091086, 16091147, 16091154, 16091227, 16091236), 
                  Mean = c(0.0999431847, 0.3864509554, 0.3864509554, 0.0999431847,
                           0.0999431847, 0.0999431847), 
                  Standard_deviation = c(0.0106193089, 0.0181384583, 0.0181384583,
                                         0.0106193089, 0.0106193089, 0.0106193089)
)

# Merging: 
df_all <- dplyr::left_join(df1, df2, by = "ID")

# Function to perform t.test for given means and standard deviations: 
# mean_x: mean of sample x
# mean_y: mean of sample y
# sd_x: standard deviation of sample x
# sd_y: standard deviation of sample y 
# n_x: sample size of sample x
# n_y: sample size of sample y
# mu: difference in means 
# alpha: significance level; default is 0.05.. important for F Test of variance comparison. 
t.test_new <- function(mean_x, mean_y, sd_x, sd_y, n_x, n_y, mu = 0, alpha = 0.05)
{
  
  # test for equal variances: 
  F_test <- sd_x^2 / sd_y^2
  F_crit <- qf(p = c(alpha / 2, 1 - alpha / 2), df1 = n_x - 1, df2 = n_y - 1)
  
  if(dplyr::between(F_test, F_crit[1], F_crit[2])) var_equal <- TRUE else var_equal <- FALSE 
  
  if (var_equal) {
   sd_test <- sqrt( (((n_x - 1) * sd_x^2 + (n_y - 1) * sd_y^2) / (n_x + n_y - 2)) * ((n_x + n_y) / (n_x * n_y)) ) 
   dof <- n_x + n_y - 2
  } else {
    sd_test <- sqrt((sd_x^2 / n_x) + (sd_y^2 / n_y))
    c_welch <- (sd_x^2 / n_x) / ((sd_x^2 / n_x) + (sd_y^2 / n_y))
    dof <- 1 / ( (c_welch^2 / (n_x - 1)) + ((1 - c_welch)^2 / (n_y - 1)) )
  }

  t_test <- ((mean_x - mean_y) - mu) / sd_test
  output <- c(
    statistic = t_test, 
    parameter = dof, 
    p.value = 2 * pt(q = abs(t_test), df = dof, lower.tail = FALSE), 
    stderr = sd_test
  )
  output
}


test_mat <- t(sapply(1:nrow(df_all), FUN = function(i) {
  t.test_new(mean_x = df_all[i, ]$Mean.x, mean_y = df_all[i, ]$Mean.y, 
                 sd_x = df_all[i, ]$Standard_deviation.x, 
             sd_y = df_all[i, ]$Standard_deviation.y, 
                 n_x = 157, n_y = 157
  )}
))

final_df <- dplyr::bind_cols(df_all, as.data.frame(test_mat))
final_df
#>         ID    Mean.x Standard_deviation.x     Mean.y Standard_deviation.y
#> 1 16091071 0.1551045           0.01203349 0.09994318           0.01061931
#> 2 16091086 0.1528096           0.01252742 0.38645096           0.01813846
#> 3 16091147 0.3395656           0.01869074 0.38645096           0.01813846
#> 4 16091154 0.2788548           0.01726190 0.09994318           0.01061931
#> 5 16091227 0.2504561           0.01767269 0.09994318           0.01061931
#> 6 16091236 0.2776752           0.01754305 0.09994318           0.01061931
#>    statistic parameter       p.value      stderr
#> 1   43.06580  312.0000 2.469975e-133 0.001280860
#> 2 -132.80314  277.2394 5.436359e-253 0.001759306
#> 3  -22.55594  312.0000  1.686821e-67 0.002078626
#> 4  110.61225  259.2849 3.304864e-220 0.001617466
#> 5   91.47052  255.6602 2.418289e-197 0.001645479
#> 6  108.59699  256.7911 1.257999e-216 0.001636620

Created on 2020-07-07 by the reprex package (v0.3.0)

Improve this answer
$\endgroup$
1
  • $\begingroup$ Wow, thank you! That is a great solution! $\endgroup$
    – PolII
    Commented Jul 15, 2020 at 12:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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