Функция для запуска линейной модели и anova для нескольких переменных и сбора значений p в кадре данных

Question

Из переменных точка A, точка B, точка C запускаются в линейной модели, и я хотел бы собрать p-значения во фрейме данных. В моем реальном наборе данных есть около 30 переменных, поэтому я ищу функцию, выполняющую их и собирающую p-значения во фрейме данных. Более эффективный способ, чем выполнять каждый из них и вручную объединять их во фрейм данных.

set.seed(1)
id <- rep(1:3,each=4)
trt <- rep(c("A","OA", "B", "OB"),3)
pointA <- sample(1:10,12, replace=TRUE)
pointB<- sample(1:10,12, replace=TRUE)
pointC<- sample(1:10,12, replace=TRUE)
df <- data.frame(id,trt,pointA, pointB,pointC)
df

id trt pointA pointB pointC
1   1   A      8      8     10
2   1  OA      2      7      3
3   1   B      8      5      5
4   1  OB      5      9      4
5   2   A      9      5      7
6   2  OA      7      3      3
7   2   B      8      1      5
8   2  OB      6      1      8
9   3   A      6      4      1
10  3  OA      8      6      9
11  3   B      1      7      4
12  3  OB      5      5      9


data <- function(i){
  lmdf <- lm(df[,i]~trt, data=df)
  anv <- anova(lmdf)
  pvalue <- anv$`Pr(>F)`
  return(pvalue)
  }
data(5)

Мне бы хотелось, чтобы это выглядело примерно так:

 variable pvalue
1   pointA  0.714
2   pointB  0.949
3   pointC  0.080

Answer 1

Вы можете рассматривать столбцы данных как список и использовать sapply для перебора столбцов, которые вы хотите использовать в качестве результатов.

Получение p-значения из каждого ANOVA немного возможно, кто-то знает лучший способ, но это работает:

pvalues <- sapply( df[,3:5] , FUN = function(x) 
  summary(aov(x~df$trt))[[1]]$`Pr(>F)`[1]  
  )
data.frame(pvalues)

         pvalues
pointA 0.9737895
pointB 0.2482931
pointC 0.7660808

Я действительно думал об использовании многомерной регрессии, но было непросто получить p-значения для каждой выходной переменной таким образом.

Answer 2

Вы можете попробовать использовать пакет метлы, потому что он имеет гибкость более чем на 1 семестр в anova, а также, если вам нужны другие характеристики.

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

library(dplyr)
library(tidyr)
library(broom)

df %>% pivot_longer(-c(id,trt)) 
# A tibble: 36 x 4
      id trt   name   value
   <int> <fct> <chr>  <int>
 1     1 A     pointA     9
 2     1 A     pointB     6
 3     1 A     pointC     9
 4     1 OA    pointA     4
 5     1 OA    pointB    10
 6     1 OA    pointC     1
 7     1 B     pointA     7
 8     1 B     pointB     7
 9     1 B     pointC     4
10     1 OB    pointA     1

С этого момента она группируется по имени столбца, выполняя анову с тем же самым формула и сбор статистики (используя tidy из broom):

res = df %>% pivot_longer(-c(id,trt)) %>% 
group_by(name) %>% do(tidy(anova(lm(value ~ trt,data=.)))) 
# A tibble: 6 x 7
# Groups:   name [3]
  name   term         df  sumsq meansq statistic p.value
  <chr>  <chr>     <int>  <dbl>  <dbl>     <dbl>   <dbl>
1 pointA trt           3   2.67  0.889    0.0711   0.974
2 pointA Residuals     8 100.   12.5     NA       NA    
3 pointB trt           3  27.7   9.22     1.68     0.248
4 pointB Residuals     8  44.0   5.50    NA       NA    
5 pointC trt           3  14.0   4.67     0.386    0.766
6 pointC Residuals     8  96.7  12.1     NA       NA   

Вам просто нужно отфильтровать результаты из trt:

res %>% filter(term=="trt")
# A tibble: 3 x 7
# Groups:   name [3]
  name   term     df sumsq meansq statistic p.value
  <chr>  <chr> <int> <dbl>  <dbl>     <dbl>   <dbl>
1 pointA trt       3  2.67  0.889    0.0711   0.974
2 pointB trt       3 27.7   9.22     1.68     0.248
3 pointC trt       3 14.0   4.67     0.386    0.766

Вы можете сделать это в baseR:

library(broom)
res = lapply(colnames(df)[-c(1:2)],function(i){
this_formula = as.formula(paste(i,"~ trt"))
anova_fit = anova(lm(this_formula, data=df))
data.frame(name=i,statistic=anova_fit[1,4],p.value=anova_fit[1,5])
})
res = do.call(rbind,res)

    name     statistic   p.value
1 pointA 0.07111111 0.9737895
2 pointB 1.67676768 0.2482931
3 pointC 0.38620690 0.7660808