How would I run multiple paired Wilcoxon signed rank tests over many pairs of vectors in an R data frame?

Question

I have a dataset of 161 immune markers, each a vector in a data frame. Using R, I want to compare 78 pairs of these vectors using the Wilcoxon signed rank (paired) test. The immune markers are distinguished in their names by "_MOM" or "_CB."

Here's a "toy" dataset with example variable names:


# Create toy data frame
toydata = data.frame(CCBB_dyad_number=c(1,2,3,4,5,6,7,8,9,10),
                cCMV_status = c("cCMV+", "cCMV-", "cCMV-", 
                                "cCMV+", "cCMV+", "cCMV-",
                                "cCMV-", "cCMV+", "cCMV+",
                                "cCMV+"),
                maternal_CMV_IgM_status = c("negative", "negative", "positive", 
                                            "negative", "positive", "negative",
                                            "positive", "positive", "positive",
                                            "negative"),
                TB40E_conc_CB = c(1.954727, NA, 1.992956,
                                1.831331, 1.905936, 2.053446,
                                2.055809, 1.739377, 2.052576,
                                1.961838),
                AD169r_conc_CB = c(5.86714, 6.469020, 9.387268,
                                   5.733174, 6.480673, 5.176167,
                                   7.548077, 7.209173, 4.944089,
                                   9.667219),
                TB40E_conc_MOM = c(7.389400, 5.917861, 7.022016,
                                 8.017846, 10.046830, 7.503896,
                                 6.427719, 9.498801, 7.351678,
                                 6.050478),
                AD169r_conc_MOM = c(7.011906, 6.506734, 9.986478,
                                    5.673412, 3.825439, 5.795331,
                                    7.082124, 6.810222, 5.54213,
                                    8.271366)
                )

With some help, I have written code to loop through all 161 vectors and produce a new data frame with p-values and type of test using lapply:



# Pull actual names of variables, not just numbers

excluded_vars <- toydata %>%
  select(., c(CCBB_dyad_number,
              cCMV_status,
              maternal_CMV_IgM_status)) %>%
  names(.)

var_list <- toydata %>%
  select(., -any_of(excluded_vars)) %>%
  names(.)


out = lapply(var_list, function(v){
  #cat(paste0("Wilcox: ", v, "
")) #Loop message for checking
  fmla <- formula(paste(v, " ~ cCMV_status"))
  wilcox.test(fmla, data = toydata, paired = FALSE) %>%
    purrr::flatten() %>% #Unnest/convert to plain list
    as.data.frame(stringsAsFactors=FALSE) %>% #Set as data frame
    mutate(Variable = v) %>% #add new variable column (could also get it from data.name)
    select(Variable, W.statistic=W, P.value=p.value, Method=method) %>%
    mutate(P.value=scientific(P.value, digits=2, format="e"))
}) %>% #%T>% { names(out) <- var_list } %>%  #Didn't actually need this, but could if wanted a named list
  purrr::compact() %>% #Remove any empty data frames/list elements (NULL)
  dplyr::bind_rows() #Bind list of data frames into single data frame

out$FDR_P.value <- p.adjust(out$P.value, method="fdr", n=length(out$P.value)) %>%
  scientific(., digits = 2, format = "e")

col_order <- c("Variable", "W.statistic", "P.value", # Reorder columns for tabling
               "FDR_P.value", "Method")

out <- out[, col_order]
  
kable(out, "html", booktabs = T) %>%
  kable_styling(latex_options = c("striped", "scale_down")) # Print output as a nice table

However, I'm having trouble thinking through how to write code to loop the signed rank test through multiple different pairs of vectors. I'm thinking I would pull the vectors (or just vector names?), like so:



toy_cCMV_pos <- toydata %>%
  filter(cCMV_status == 'cCMV+') %>%
  select(., -any_of(excluded_vars))


variable.set1 <- toy_cCMV_pos %>%
  select(., ends_with("_MOM"))


variable.set2 <- toy_cCMV_pos %>%
  select(., ends_with("_CB"))

Someone suggested looping through the vectors like this. However, I keep getting an "undefined columns selected" error, and because I don't quite understand what the code below is doing, I can't troubleshoot.

for (a in variable.set1) {
  groups = unique(toy_cCMV_pos[,a])
  for (b in variable.set2) {
    wilcox.test(x=toy_cCMV_pos[which(toy_cCMV_pos[a]==groups[1]),b], 
                y=toy_cCMV_pos[which(toy_cCMV_pos[a]==groups[2]),b], 
                paired=TRUE)
  }
}

# Keep getting error "undefined columns selected"

I would like to be able to pull results, including p-values, into a new data frame as with the rank sum tests.

Could anyone help me think through how to do run these paired tests?

Claudio · Accepted Answer

EDIT: the original solution deleted missing values row-wise, so some valid data were delete too leading with results inconsistent with other methods.

Here is a more correct approach:

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

toydata = data.frame(CCBB_dyad_number=c(1,2,3,4,5,6,7,8,9,10),
                cCMV_status = c("cCMV+", "cCMV-", "cCMV-", 
                                "cCMV+", "cCMV+", "cCMV-",
                                "cCMV-", "cCMV+", "cCMV+",
                                "cCMV+"),
                maternal_CMV_IgM_status = c("negative", "negative", "positive", 
                                            "negative", "positive", "negative",
                                            "positive", "positive", "positive",
                                            "negative"),
                TB40E_conc_CB = c(1.954727, NA, 1.992956,
                                1.831331, 1.905936, 2.053446,
                                2.055809, 1.739377, 2.052576,
                                1.961838),
                AD169r_conc_CB = c(5.86714, 6.469020, 9.387268,
                                   5.733174, 6.480673, 5.176167,
                                   7.548077, 7.209173, 4.944089,
                                   9.667219),
                TB40E_conc_MOM = c(7.389400, 5.917861, 7.022016,
                                 8.017846, 10.046830, 7.503896,
                                 6.427719, 9.498801, 7.351678,
                                 6.050478),
                AD169r_conc_MOM = c(7.011906, 6.506734, 9.986478,
                                    5.673412, 3.825439, 5.795331,
                                    7.082124, 6.810222, 5.54213,
                                    8.271366))

toydata |> 
  select(ends_with("MOM"), ends_with("CB")) |> 
  pivot_longer(everything(),
               names_to=c(".value", "group"),
               names_sep="_(?!.*_)") |> 
  pivot_longer(-group,
               names_to="variable",
               values_to="value") |>  
  group_by(variable) |> 
  do(broom::tidy(wilcox.test(.$value ~ .$group, paired=TRUE, na.action=na.pass)))
#> # A tibble: 2 × 5
#> # Groups:   variable [2]
#>   variable    statistic p.value method                          alternative
#>                                                   
#> 1 AD169r_conc        28 1       Wilcoxon signed rank exact test two.sided  
#> 2 TB40E_conc          0 0.00391 Wilcoxon signed rank exact test two.sided

^{Created on 2021-09-09 by the reprex package (v2.0.1)}

The results match those of the individual calculations:

> wilcox.test(toydata$TB40E_conc_CB, toydata$TB40E_conc_MOM, paired=TRUE)

    Wilcoxon signed rank exact test

data:  toydata$TB40E_conc_CB and toydata$TB40E_conc_MOM
V = 0, p-value = 0.003906
alternative hypothesis: true location shift is not equal to 0

And

> wilcox.test(toydata$AD169r_conc_CB, toydata$AD169r_conc_MOM, paired=TRUE)

    Wilcoxon signed rank exact test

data:  toydata$AD169r_conc_CB and toydata$AD169r_conc_MOM
V = 28, p-value = 1
alternative hypothesis: true location shift is not equal to 0

The result of the proposed solution is a tibble/dataframe, so you can modify it selecting only the needed columns.

How would I run multiple paired Wilcoxon signed rank tests over many pairs of vectors in an R data frame?

Answers (2)

Related Questions