R - T-тест между несколькими колонками, кондиционированными на 2 других колонках

Question

Я пытаюсь вычислить различия в SEX для каждого вида в моем df по 5 переменным (WING; WINGPRI; WEIGHT, BEAK и TARSUS) Я пробовал десятки кодов (я начинающий) и продолжаю получатьошибки или коды не дают мне желаемого результата и т. д.

Это примерно так же близко, как я получил:

library(tidyverse)
library(broom)

df %>% 
      select(SPECIES, SEX, WING, WINGPRI, WEIGHT, BEAK, TARSUS) %>% 
      gather(key = variable, value=value, -SEX, -SPECIES) %>% 
      group_by(SEX, SPECIES, variable) %>%
      summarise(value = list(value)) %>%
      na.rm=TRUE %>%
      tibble::rowid_to_column() %>%
      spread(SPECIES, value) %>%
      group_by(variable) %>% 
      mutate(p_value = t.test(unlist("F"), unlist("M"))$p.value,
             t_value = t.test(unlist("F"), unlist("M"))$statistic)

Но после устранения многих ошибок я не могу получитьмимо этого Error: is.data.frame(df) is not TRUE

Любые идеи о том, как решить этот код или совершенно другой код, который будет делать то, что я пытаюсь сделать, высоко ценится

> dput(sample)
structure(list(RING = c("A264874", "A432586", "O92477", "B9124", 
"C95571", "A395011", "C88213", "C58443", "A95422", "C58409"), 
    SPECIES = c("CARDUELIS CARDUELIS", "SYLVIA ATRICAPILLA", 
    "ESTRILDA ASTRILD", "ALCEDO ATTHIS", "CHLORIS CHLORIS", "FRINGILLA COELEBS", 
    "SYLVIA ATRICAPILLA", "CHLORIS CHLORIS", "SYLVIA ATRICAPILLA", 
    "PARUS MAJOR"), SEX = c("U", "M", "F", "F", "F", "F", "F", 
    "M", "F", "M"), AGE = c(2L, 3L, 3L, 3L, 4L, 2L, 4L, 4L, 6L, 
    3L), FAT = c(0L, 0L, 0L, NA, 0L, 0L, 0L, 0L, 1L, 0L), WEIGHT = c(148, 
    185, 85, 32, 225, 20, 245, 22, 197, 19), WING = c(775, 69, 
    45, 76, 82, 84, 77, 83, 69, 72), WINGPRI = c(58L, NA, 32L, 
    NA, NA, NA, 57L, 64L, 52L, 54L), BEAK = c(156, 132, 86, NA, 
    NA, 138, 125, 13, 145, 125), TARSUS = c(148, 199, 146, NA, 
    NA, 178, 18, 177, 207, 205), BROODPATCH = c(0L, NA, 0L, 0L, 
    0L, NA, 0L, 0L, 0L, 0L), MUSCLE = c(2L, 3L, 2L, 2L, 2L, 2L, 
    2L, 1L, 2L, 2L), PROGRAM = c("MAI", "MAI", "MAI", "MIGRA<c7><c3>O", 
    "PEEC", "MAI", "MAI", "PEEC", "MAI", "MIGRA<c7><c3>O")), .Names = c("RING", 
"SPECIES", "SEX", "AGE", "FAT", "WEIGHT", "WING", "WINGPRI", 
"BEAK", "TARSUS", "BROODPATCH", "MUSCLE", "PROGRAM"), row.names = c(NA, 
-10L), class = c("tbl_df", "tbl", "data.frame"))

Answer 1

Рассмотрим пакет base, который поставляется с каждой установкой R и загружается в каждом сеансе R, не требуя base::library().В частности, используйте его by (объектно-ориентированная оболочка tapply) для подстановки фрейма данных на SPECIES и построения матрицы t-stat и p-value с помощью sapply.

Поскольку в данных OP недостаточно наблюдений для t.test, ниже генерируется набор данных

set.seed(10102018)

species_df <- data.frame(
  SPECIES = sample(unique(df$SPECIES), 500, replace=TRUE),
  SEX = sample(c("F", "M"), 500, replace=TRUE),
  WING = runif(500) * 100,
  WINGPRI = runif(500) * 100,
  WEIGHT = runif(500) * 100,
  BEAK = runif(500) * 100,
  TARSUS =  runif(500) * 100,
  stringsAsFactors = FALSE
)

# NAMED LIST OF MATRICES
mat_list <- by(species_df, species_df$SPECIES, function(sub){     
   sapply(c("WING", "WINGPRI", "WEIGHT", "BEAK", "TARSUS"), function(col) {
      tryCatch({
        output <- t.test(sub[sub$SEX=="F", col], sub[sub$SEX=="M", col])          
        c(p_value = output$p.value, t_value = output$statistic)
      }, warning = function(w) c(p_value=NA, t_value=NA)
       , error = function(e) c(p_value=NA, t_value=NA))
    })      
})

Выход

mat_list

# species_df$SPECIES: ALCEDO ATTHIS
#                WING    WINGPRI     WEIGHT      BEAK     TARSUS
# p_value   0.7273165  0.8382798  0.3180979 0.6450270  0.3856571
# t_value.t 0.3501749 -0.2048995 -1.0055505 0.4629014 -0.8733496
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: CARDUELIS CARDUELIS
#                WING    WINGPRI    WEIGHT      BEAK    TARSUS
# p_value   0.5200729  0.8520463 0.3370721 0.8189008 0.1212502
# t_value.t 0.6470729 -0.1873091 0.9678003 0.2299977 1.5716422
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: CHLORIS CHLORIS
#                 WING    WINGPRI     WEIGHT       BEAK    TARSUS
# p_value    0.1115453  0.5689228 0.94825726  0.5989776 0.9108546
# t_value.t -1.6129915 -0.5725928 0.06514506 -0.5284384 0.1124033
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: ESTRILDA ASTRILD
#                 WING   WINGPRI    WEIGHT      BEAK    TARSUS
# p_value   0.09291222 0.7700545 0.6859697 0.1958938 0.6452502
# t_value.t 1.70719717 0.2935269 0.4062293 1.3054498 0.4624954
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: FRINGILLA COELEBS
#                 WING   WINGPRI    WEIGHT      BEAK   TARSUS
# p_value   0.06157204 0.8636649 0.2183259 0.4757378 0.274626
# t_value.t 1.89924201 0.1723255 1.2416417 0.7170863 1.101813
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: PARUS MAJOR
#                 WING   WINGPRI    WEIGHT       BEAK    TARSUS
# p_value   0.96688923 0.5857059 0.1140328  0.5055508 0.5747242
# t_value.t 0.04168846 0.5481212 1.6046303 -0.6694396 0.5643418
# ---------------------------------------------------------------------------------------- 
# species_df$SPECIES: SYLVIA ATRICAPILLA
#                 WING   WINGPRI    WEIGHT      BEAK     TARSUS
# p_value    0.4350621 0.5446387 0.7073097 0.3911381  0.7631614
# t_value.t -0.7851506 0.6091449 0.3770283 0.8628441 -0.3024993

Answer 2

Здравствуйте, земляк,

Вы пытаетесь запустить код в списке, предполагая, что это фрейм данных.Я не пытаюсь сначала преобразовать df в фрейм данных, и они запускают код.Я бы посоветовал вам попробовать этот пример и посмотреть, как он работает:

# Creates single vectors for test
RING = c("A264874", "A432586", "O92477", "B9124", 
         "C95571", "A395011", "C88213", "C58443", "A95422", "C58409")
SPECIES = c("CARDUELIS CARDUELIS", "SYLVIA ATRICAPILLA", 
            "ESTRILDA ASTRILD", "ALCEDO ATTHIS", "CHLORIS CHLORIS", "FRINGILLA COELEBS", 
            "SYLVIA ATRICAPILLA", "CHLORIS CHLORIS", "SYLVIA ATRICAPILLA", 
            "PARUS MAJOR")
SEX = c("U", "M", "F", "F", "F", "F", "F", 
        "M", "F", "M")
AGE = c(2L, 3L, 3L, 3L, 4L, 2L, 4L, 4L, 6L,3L)
FAT = c(0L, 0L, 0L, NA, 0L, 0L, 0L, 0L, 1L, 0L)
WEIGHT = c(148,185, 85, 32, 225, 20, 245, 22, 197, 19)
WING = c(775, 69, 45, 76, 82, 84, 77, 83, 69, 72)

#List them all
df <- list(WEIGHT,FAT,AGE,SEX,SPECIES,RING,WING)

# Create dataframe
data.frame <- do.call(cbind.data.frame, df)
colnames(data.frame) <- c("WEIGHT","FAT","AGE","SEX","SPECIES","RING","WING")

Проверьте, работает ли ваш код сейчас

Answer 3

Использование другой библиотеки:

library(matrixTests)

vars <- c("WING", "WINGPRI", "WEIGHT", "BEAK", "TARSUS")
col_t_welch(sample[sample$SEX=="F", vars], sample[sample$SEX=="M", vars])

        obs.x obs.y obs.tot    mean.x    mean.y mean.diff      var.x      var.y    stderr       df  statistic    pvalue   conf.low conf.high alternative mean.null conf.level
WING        6     3       9  72.16667  74.66667  -2.50000   204.5667   54.33333  7.225341 6.873827 -0.3460044 0.7396874  -19.64900  14.64900   two.sided         0       0.95
WINGPRI     3     2       5  47.00000  59.00000 -12.00000   175.0000   50.00000  9.128709 2.985075 -1.3145341 0.2805428  -41.13398  17.13398   two.sided         0       0.95
WEIGHT      6     3       9 134.00000  75.33333  58.66667 10074.4000 9022.33333 68.458097 4.318219  0.8569719 0.4364231 -126.00493 243.33827   two.sided         0       0.95
BEAK        4     3       7 123.50000  90.00000  33.50000   693.6667 4459.00000 40.740029 2.471493  0.8222871 0.4826893 -113.34042 180.34042   two.sided         0       0.95
TARSUS      4     3       7 137.25000 193.66667 -56.41667  6940.9167  217.33333 42.516745 3.247234 -1.3269282 0.2701275 -186.07859  73.24525   two.sided         0       0.95

Как упомянуто @AntoniosK - нельзя разделить по видам с использованием предоставленного вами набора данных.

Answer 4

Вот пример, основанный на данных, которые вы разместили

library(tidyverse)
library(broom)

df %>% 
  select(SPECIES, SEX, WING, WINGPRI, WEIGHT, BEAK, TARSUS) %>% 
  gather(key = variable, value=value, -SEX, -SPECIES) %>% 
  filter(SEX != "U") %>%
  group_by(variable) %>%
  summarise(tt = list(tidy(t.test(value ~ SEX)))) %>%
  unnest()

# # A tibble: 5 x 11
#   variable estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method                  alternative
#   <chr>       <dbl>     <dbl>     <dbl>     <dbl>   <dbl>     <dbl>    <dbl>     <dbl> <chr>                   <chr>      
# 1 BEAK         33.5     124.       90       0.822   0.483      2.47   -113.      180.  Welch Two Sample t-test two.sided  
# 2 TARSUS      -56.4     137.      194.     -1.33    0.270      3.25   -186.       73.2 Welch Two Sample t-test two.sided  
# 3 WEIGHT       58.7     134        75.3     0.857   0.436      4.32   -126.      243.  Welch Two Sample t-test two.sided  
# 4 WING         -2.5      72.2      74.7    -0.346   0.740      6.87    -19.6      14.6 Welch Two Sample t-test two.sided  
# 5 WINGPRI     -12        47        59      -1.31    0.281      2.99    -41.1      17.1 Welch Two Sample t-test two.sided  

Я также не могу группировать по SPECIES, потому что у меня недостаточно наблюдений для t.test для работы.

Обратите внимание , что t.test автоматически игнорирует NA значения, но вы можете добавить ... %>% na.omit() %>% ... в свой код перед выполнением теста.

Примечание что tt = list(tidy(t.test(value ~ SEX))) автоматически использует соответствующий поднабор данных в зависимости от вашей группировки.Однако, если вы используете tt = list(tidy(t.test(value ~ SEX, data= .))), он будет использовать весь набор данных каждый раз.