A partir de las variables punto A, punto B, punto C se ejecutan en un modelo lineal y anova, me gustaría recopilar los valores p en un marco de datos. En mi conjunto de datos reales hay alrededor de 30 variables, por lo que estoy buscando una función que las ejecute y recopile los valores p en un marco de datos. Una forma más eficiente que ejecutar cada uno y unirlos manualmente en un marco de datos.

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)

Me gustaría que se parezca a algo como:

 variable pvalue
1   pointA  0.714
2   pointB  0.949
3   pointC  0.080
r
1
user11916948 27 feb. 2020 a las 22:49

2 respuestas

La mejor respuesta

Puede tratar las columnas de un marco de datos como una lista y usar sapply para iterar sobre las columnas que desea usar como resultados.

Obtener el valor p de cada ANOVA es un poco complicado, tal vez alguien conozca una mejor manera, pero esto funciona:

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

Pensé en usar una regresión multivariada, pero no fue más fácil obtener los valores p para cada variable de resultado de esa manera.

2
George Savva 27 feb. 2020 a las 23:01

Puede intentar usar el paquete de escoba, porque tiene la flexibilidad para más de 1 término en anova, y también si necesita otras estadísticas.

Primero pivotamos más para que su punto de prueba ahora sea una variable de columna:

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

A partir de este momento, se agrupa por el nombre de la columna, realiza un anova con la misma fórmula y recopila las estadísticas (usando tidy de 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   

Solo necesita filtrar sus resultados de 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

Puedes hacerlo en 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
1
StupidWolf 28 feb. 2020 a las 09:21