Сохранение всех значений образца rnorm в больших количествах (вектор) для теста Тьюки позже - PullRequest
0 голосов
/ 21 июня 2020

Так что я почти у цели! Что я пытаюсь сделать, так это создать сценарий, в котором вы заполняете «ожидаемые средние (mu)» и «стандартные отклонения (sd)» двух популяций, и в качестве результата вы получаете таблицу, в которой вы показываете долю значительных различий, обнаруженных между этими двумя популяции (определяемые ANOVA и тестом Тьюки) на основе размера выборки. Конечно, ожидая, что чем больше размер выборки, тем выше будет доля значимых результатов. Теперь происходит что-то странное: когда я сравниваю два средних (200 и 260) с sd = 30, кажется, что доля значимых p-значений случайным образом распределяется по размерам выборки (ниже я использовал размер выборки 1,2,3. ..9). Я надеялся обнаружить, что по мере увеличения размера выборки доля значимых различий, обнаруженных с помощью ANOVA и Тьюки, также будет увеличиваться из-за того, как работает ANOVA. Чем больше размер выборки, тем больше solid моделей и, следовательно, больше solid p-значений. Я ожидал что-то вроде этого:

n   fraction significant
1        0.1
2        0.3
3        0.4
4        0.6
5        0.8
6        0.9

Этого не происходит, я получаю что-то вроде этого:


n   fraction significant
1        0.1
2        0.0
3        0.2
4        0.4
5        0.0
6        0.1

Я думаю, это потому, что изначально, когда я запускаю свои циклы и rnorm, он сохраняет только одно значение в качестве «выборки» на начальный размер выборки: например, для начального размера выборки n = 1 и «exp» = 1 (это мои итерации) я получаю 1 случайно сгенерированное значение выборки, что логично. Но для n = 9 и «exp = 1» я тоже получаю только одно значение вместо 9 значений. В этом случае я должен ввести 9 значений в ANOVA и Тьюки, что, вероятно, приведет к большему количеству итераций «exp», что приведет к значительной разнице. Я думаю, что если мне удастся достичь чего-то подобного, это должно сработать, хотя я не знаю, как это сделать, я пробовал 3 дня и застрял:

n    sample pop 1         sample pop 2
1    200.4                270.1           compare by ANOVA and Tukey -> p value
1    190.3                234.7           compare by ANOVA and Tukey -> p value
2.   [200.2 184.6]        [260.2 265.3]   compare by ANOVA and Tukey -> p value
2.   [230.1 222.7]        [280.1 204.6]   compare by ANOVA and Tukey -> p value
3.   [270.3 230.1 201.3]  [232.2 222.1]   compare by ANOVA and Tukey -> p value
#etc... 
#(Right now 10 times per n-value, but ultimately I want to do 10.000 iterations)

Вот сценарий, который я сделал, прошу прощения у кого есть опыт;)

{
  library(tidyverse)
  library(dplyr)
  library(reshape2)
  library(multcomp)
  library(stringr)
}
rm(list=ls())
{
  mu <- 200
  sd <- 30
  n_array_1 <- c()
  exp_array_1 <- c()
  pop_sample_array_1 <- c()
  ID_array_1 <- "population1"
  for (n_1 in c(1:9)) {
    for (exp_1 in c(1:10)) {
      n_array_1 <- c(n_array_1, n_1)
      exp_array_1 <- c(exp_array_1, exp_1)
      population_sample_1 <- rnorm(n_1, mu, sd)
      pop_sample_array_1 <- c(pop_sample_array_1, population_sample_1)
      
    }
}

    mu_2 <- 260
    sd_2 <- 30
    n_array_2 <- c()
    exp_array_2 <- c()
    pop_sample_array_2 <- c()
    ID_array_2 <- "population2"
   
    for (n_2 in c(1:9)) {
      for (exp_2 in c(1:10)) {
        n_array_2 <- c(n_array_2, n_2)
        exp_array_2 <- c(exp_array_2, exp_2)
        population_sample_2 <- rnorm(n_2, mu_2, sd_2)
        pop_sample_array_2 <- c(pop_sample_array_2, population_sample_2)
      
        
      }
}
  
  df1 <- data.frame(n=n_array_1, exp=exp_array_1, sample=pop_sample_array_1, ID=ID_array_1)
  df2 <- data.frame(n=n_array_2, exp=exp_array_2, sample=pop_sample_array_2, ID=ID_array_2)
}

combineddf <- rbind(df1, df2)
view(combineddf)
#
## Running an ANOVA comparing the 2 groups specified
{
  combineddf$IDcombined <- as.factor(paste(combineddf$IDcombined, combineddf$ID, combineddf$n, combineddf$exp,  sep = "_"))
  combineddf$IDcombined <- as.factor(combineddf$IDcombined)
  combineddf
  resaov <- aov(sample ~ IDcombined, data=combineddf)
  test <- TukeyHSD(resaov)
  TK <- test
  TK_data<-as.data.frame(TK[1:1])
  TK_data$IDcombined <- rownames(TK_data)
}
{
LC1 <- separate(TK_data, col=IDcombined, into = c("A", "B"), sep = "-")
LC2 <- separate(LC1, col = A, into=c("ignore", "comparison_factor_1", "sample_size", "iteration"), sep="_")
LC3 <- separate(LC2, col=B, into=c("ignore", "comparison_factor_2", "sample_size_2", "iteration_2"), sep="_")
  }
{
  results <- dplyr::select(LC3, comparison_factor_1, comparison_factor_2, iteration, iteration_2, sample_size, sample_size_2, IDcombined.p.adj)
results <- dplyr::filter(results, comparison_factor_1=="population1" & comparison_factor_2=="population2" |  comparison_factor_1=="population2" & comparison_factor_2=="population1")
results <- dplyr::filter(results, sample_size == sample_size_2 & iteration == iteration_2)
}

#
## calculations of fractions significant p-values of a tukey test comparing the 2 populations based on the different sample sizes

results$significant <- results$IDcombined.p.adj < 0.05
for (sample_size in unique(results$sample_size)) {
  subset_results <- results[results$sample_size == sample_size, ]
  n_significant <- sum(subset_results$significant)
  n_total <- nrow(subset_results)
  fraction_significant <- n_significant / n_total
  cat(sample_size, " ", fraction_significant, "\n")
}

1 Ответ

0 голосов
/ 24 июня 2020

Добавление еще одного «for» l oop (с «репликами») устранило мои проблемы. Теперь мне просто нужно сохранить вывод функции «cat» в таблице или фрейме данных, и я могу расширить сценарий (добавив больше групп для сравнения). Скрипт ниже:

{
  library(tidyverse)
  library(dplyr)
  library(reshape2)
  library(multcomp)
  library(stringr)
}
rm(list=ls())
#
## Syncom C - (t=10)
{
  {
    mu_1 <- 628
    sd_1 <- 130
    n_array_1 <- c()
    exp_array_1 <- c()
    pop_sample_array_1 <- c()
    ID_array_1 <- "population1"
    for (n_1 in c(1:20)) {
      for (exp_1 in c(1:30)) {
        population_sample_1 <- rnorm(n_1, mu_1, sd_1)
        for (replicate_1 in population_sample_1) {
          n_array_1 <- c(n_array_1, n_1)
          exp_array_1 <- c(exp_array_1, exp_1)
          pop_sample_array_1 <- c(pop_sample_array_1, replicate_1)
        }
      }
    }
  }
}
#
## Control - (t=10)
{
  {
    mu_2 <- 423
    sd_2 <- 96
    n_array_2 <- c()
    exp_array_2 <- c()
    pop_sample_array_2 <- c()
    ID_array_2 <- "population2"
    for (n_2 in c(1:20)) {
      for (exp_2 in c(1:30)) {
        population_sample_2 <- rnorm(n_2, mu_2, sd_2)
        for (replicate_2 in population_sample_2) {
          n_array_2 <- c(n_array_2, n_2)
          exp_array_2 <- c(exp_array_2, exp_2)
          pop_sample_array_2 <- c(pop_sample_array_2, replicate_2)
        }
      }
    }
  }
}
{  
  df1 <- data.frame(n=n_array_1, exp=exp_array_1, sample=pop_sample_array_1, ID=ID_array_1)
  df2 <- data.frame(n=n_array_2, exp=exp_array_2, sample=pop_sample_array_2, ID=ID_array_2)
}

combineddf <- rbind(df1, df2)
#view(combineddf)
#
## Running an ANOVA comparing the 2 groups specified
{
  combineddf$IDcombined <- as.factor(paste(combineddf$IDcombined, combineddf$ID, combineddf$n, combineddf$exp,  sep = "_"))
  combineddf$IDcombined <- as.factor(combineddf$IDcombined)
  combineddf
  resaov <- aov(sample ~ IDcombined, data=combineddf)
  test <- TukeyHSD(resaov)
  TK <- test
  TK_data<-as.data.frame(TK[1:1])
  TK_data$IDcombined <- rownames(TK_data)
}
{
LC1 <- separate(TK_data, col=IDcombined, into = c("A", "B"), sep = "-")
LC2 <- separate(LC1, col = A, into=c("ignore", "comparison_factor_1", "sample_size", "iteration"), sep="_")
LC3 <- separate(LC2, col=B, into=c("ignore", "comparison_factor_2", "sample_size_2", "iteration_2"), sep="_")
  }
{
  results <- dplyr::select(LC3, comparison_factor_1, comparison_factor_2, iteration, iteration_2, sample_size, sample_size_2, IDcombined.p.adj)
results <- dplyr::filter(results, comparison_factor_1=="population1" & comparison_factor_2=="population2" |  comparison_factor_1=="population2" & comparison_factor_2=="population1")
results <- dplyr::filter(results, sample_size == sample_size_2 & iteration == iteration_2)
}

#
## calculations of fractions significant p-values of a tukey test comparing the 2 populations based on the different sample sizes

results$significant <- results$IDcombined.p.adj < 0.05
for (sample_size in unique(results$sample_size)) {
  subset_results <- results[results$sample_size == sample_size, ]
  n_significant <- sum(subset_results$significant)
  n_total <- nrow(subset_results)
  fraction_significant <- n_significant / n_total
  cat(sample_size, " ", fraction_significant, "\n")
}
...