rloopsif-statementiterationoutliers

How to delete rows of outliers rom a nested dataset via an iterative method or pipe operator

I'm trying removing outliers from this nested dataset

df_join
# A tibble: 12 × 2
# Groups:   signals [12]
   sls     data             
   <chr>       <list>           
 1 alphaX1     <tibble [75 × 5]>
 2 betaY2      <tibble [75 × 5]>
 3 gammaZ3     <tibble [75 × 5]>
 4 deltaA4     <tibble [75 × 5]>
 5 epsilonB5   <tibble [75 × 5]>
 6 zetaC6      <tibble [75 × 5]>
 7 etaD7       <tibble [75 × 5]>
 8 thetaE8     <tibble [75 × 5]>
 9 iotaF9      <tibble [75 × 5]>
10 kappaG10    <tibble [75 × 5]>
11 lambdaH11   <tibble [75 × 5]>
12 muI12       <tibble [75 × 5]>

for instance, the first element of it contains this series of variable:

    # A tibble: 75 × 5
   Var1  Var2  Var3  Var4      value
   <fct> <fct> <fct> <fct>     <dbl>
 1 A1    G1    X     A-B       -8.52
 2 A1    G1    X     A-C       -4.77
 3 A1    G1    X     B-C       -3.25
 4 B2    G1    X     A-B        2.13
 5 B2    G1    X     A-C        1.85
 6 B2    G1    X     B-C        3.92
 7 C3    G1    X     A-B       -0.33
 8 C3    G1    X     A-C       -2.10
 9 C3    G1    X     B-C       -1.46
10 D4    G1    X     A-B        0.51

the entire content of it is the following one:

  df_join <- tibble::tibble(
  channel = c("chA", "chB", "chC", "chD", "chE", "chF"),
  recordings = replicate(
    6,
    tibble::tibble(
      subject_id = factor(rep(sprintf("P%02d", 1:25), each = 3)),
      group = factor("T1"),
      sxs = factor("X"),
      clc = factor(rep(c("cond1", "cond2", "cond3"), times = 25)),
      sls_strength = rnorm(75, mean = 0, sd = 5)
    ),
    simplify = FALSE
  )
) |>
  dplyr::group_by(channel)

I've tried to check for the presence of outliers as follows:

 outliers_table <- df_join %>%
  tidyr::unnest(recordings) %>%
  dplyr::select(clc, channel, sls_strength) %>%
  dplyr::group_by(clc) %>%
  rstatix::identify_outliers(sls_strength)

That turns

    # A tibble: 30 × 5
   clc     channel sls_strength is.outlier is.extreme
   <fct>     <chr>             <dbl> <lgl>      <lgl>     
 1 cond1     chA              -10.5  TRUE       FALSE     
 2 cond1     chB               14.2  TRUE       FALSE     
 3 cond1     chC               16.0  TRUE       FALSE     
 4 cond1     chC               15.3  TRUE       FALSE     
 5 cond1     chC               22.8  TRUE       TRUE      
 6 cond1     chC               13.9  TRUE       FALSE     
 7 cond1     chC               12.5  TRUE       FALSE     
 8 cond1     chD              -12.1  TRUE       FALSE     
 9 cond1     chE               18.6  TRUE       FALSE     
10 cond1     chF               24.0  TRUE       TRUE

If I'm interested in delete all of those values that are TRULY EXTREME, how could do I do by using some iterative function orr some if statment?? Please just consider also other alternative in case it is easier (also to keep on the command I've written by adding another %>% command row) that scripring down a for loop or some other function.

Since I'm at the very beginning I've coded the failing code I've created:

outliers_bale <- df_join %>%
  unnest() %>% 
  dplyr::select(clc, sls, value) %>% 
  group_by(clc) %>%  #it is the equivalent to use as grouping variable the time
  identify_outliers(value) %>% 
  filter(is.outlier & is.extreme)

values <- outliers_table$value

df_join[!(df_join$data %in% values), ]

And I am not able to figure out whether it worked or not.

Thanks in advance

Solution

  • All right. Let's do it together step by step. As I understand it, you have serious concerns that in your data (I keep it in the variable df) there are outliers and even extreme values. First, we will extract from your data only one grouped tibble and filter for COND ==" NEG-NOC "

    library(tidyverse)
library(rstatix)
library(outliers)

data = df$data[[1]] %>% filter(COND=="NEG-NOC")

    Now let's consider what method of outlier identification we will use. We can use the boxplot function for this.

    boxplot.stats(data$value)$out
#[1] 8.164181

    This is fine, but it only gives us outliers in vector form. The second way is to use identify_outliers. This gives us a tibble but still only with those lines that have these outlier values.

    data %>% identify_outliers(variable = "value")
# # A tibble: 1 x 7
# ID    GR    SES   COND    value is.outlier is.extreme
# <fct> <fct> <fct> <fct>   <dbl> <lgl>      <lgl>     
#   1 11    RP    V     NEG-NOC  8.16 TRUE       FALSE

    Well, let's use the outlier function from the outliers package. This can give us a logic vector.

    outlier(data$value, opposite = T)
#[1] 8.164181
outlier(data$value, opposite = T, logical = T)
# [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#[22] FALSE FALSE FALSE FALSE

    However, neither of these methods will assist you in deciding what to do with these outliers. Please read this carefully . As you can see, you have three options to choose from: Imputation, Capping, Prediction. Which one will you choose? I chose Capping. So I wrote a tiny function that identifies outliers, extreme values and additionally returns your values after Capping.

    fOutCapp = function(data){
  x = data$value
  qnt = quantile(x, probs=c(.25, .75), na.rm = T)
  caps = quantile(x, probs=c(.05, .95), na.rm = T)
  H = 1.5 * IQR(x, na.rm = T)
  He = 3 * IQR(x, na.rm = T)
  is.outlier = (x < (qnt[1] - H)) | (x > (qnt[2] + H))
  x[x < (qnt[1] - H)] <- caps[1]
  x[x > (qnt[2] + H)] <- caps[2]
  data %>% group_by(COND) %>% 
    mutate(
      is.outlier = is.outlier,
      is.extreme = (x < (qnt[1] - He)) | (x > (qnt[2] + He)),
      cap.value = x
    )
}

    Let's see if it works

    data %>% fOutCapp() %>% filter(is.outlier)
# A tibble: 1 x 8
# ID    GR    SES   COND    value is.outlier is.extreme cap.value
# <fct> <fct> <fct> <fct>   <dbl> <lgl>      <lgl>          <dbl>
#   1 11    RP    V     NEG-NOC  8.16 TRUE       FALSE           4.95
data %>% fOutCapp()
# A tibble: 25 x 8
# ID    GR    SES   COND      value is.outlier is.extreme cap.value
# <fct> <fct> <fct> <fct>     <dbl> <lgl>      <lgl>          <dbl>
#   1 01    RP    V     NEG-NOC -11.1   FALSE      FALSE        -11.1  
# 2 04    RP    V     NEG-NOC   0.239 FALSE      FALSE          0.239
# 3 06    RP    V     NEG-NOC  -2.96  FALSE      FALSE         -2.96 
# 4 07    RP    V     NEG-NOC   1.09  FALSE      FALSE          1.09 
# 5 08    RP    V     NEG-NOC   2.99  FALSE      FALSE          2.99 
# 6 09    RP    V     NEG-NOC   5.42  FALSE      FALSE          5.42 
# 7 10    RP    V     NEG-NOC  -2.83  FALSE      FALSE         -2.83 
# 8 11    RP    V     NEG-NOC   8.16  TRUE       FALSE          4.95 
# 9 12    RP    V     NEG-NOC  -9.83  FALSE      FALSE         -9.83 
# 10 13    RP    V     NEG-NOC   2.12  FALSE      FALSE          2.12 
# ... with 15 more rows

    Note, however, that your data inside the variable data is grouped after the variable COND. So let's write one more tiny function that will do our fOutCapp on each of the groups.

    fOutCappGroup = function(data) data %>% group_by(COND) %>% 
  group_modify(~fOutCapp(.x))

df$data[[1]] %>% fOutCappGroup()
# # A tibble: 75 x 8
# # Groups:   COND [3]
# COND    ID    GR    SES     value is.outlier is.extreme cap.value
# <fct>   <fct> <fct> <fct>   <dbl> <lgl>      <lgl>          <dbl>
#   1 NEG-CTR 01    RP    V     -11.6   FALSE      FALSE        -11.6  
# 2 NEG-CTR 04    RP    V      -0.314 FALSE      FALSE         -0.314
# 3 NEG-CTR 06    RP    V      -0.214 FALSE      FALSE         -0.214
# 4 NEG-CTR 07    RP    V      -2.83  FALSE      FALSE         -2.83 
# 5 NEG-CTR 08    RP    V       4.24  FALSE      FALSE          4.24 
# 6 NEG-CTR 09    RP    V       9.57  FALSE      FALSE          9.57 
# 7 NEG-CTR 10    RP    V      -6.13  FALSE      FALSE         -6.13 
# 8 NEG-CTR 11    RP    V       0.529 FALSE      FALSE          0.529
# 9 NEG-CTR 12    RP    V      -7.74  FALSE      FALSE         -7.74 
# 10 NEG-CTR 13    RP    V       1.27  FALSE      FALSE          1.27 
# # ... with 65 more rows

    Bingo. Everything works great. Now we only needs to do one simple mutation.

    df %>% group_by(signals) %>% 
  mutate(data = map(data, ~fOutCappGroup(.x))) %>% 
  unnest(data)

    output

    # A tibble: 450 x 9
# Groups:   signals [6]
   signals COND    ID    GR    SES     value is.outlier is.extreme cap.value
   <chr>   <fct>   <fct> <fct> <fct>   <dbl> <lgl>      <lgl>          <dbl>
 1 P3FCz   NEG-CTR 01    RP    V     -11.6   FALSE      FALSE        -11.6  
 2 P3FCz   NEG-CTR 04    RP    V      -0.314 FALSE      FALSE         -0.314
 3 P3FCz   NEG-CTR 06    RP    V      -0.214 FALSE      FALSE         -0.214
 4 P3FCz   NEG-CTR 07    RP    V      -2.83  FALSE      FALSE         -2.83 
 5 P3FCz   NEG-CTR 08    RP    V       4.24  FALSE      FALSE          4.24 
 6 P3FCz   NEG-CTR 09    RP    V       9.57  FALSE      FALSE          9.57 
 7 P3FCz   NEG-CTR 10    RP    V      -6.13  FALSE      FALSE         -6.13 
 8 P3FCz   NEG-CTR 11    RP    V       0.529 FALSE      FALSE          0.529
 9 P3FCz   NEG-CTR 12    RP    V      -7.74  FALSE      FALSE         -7.74 
10 P3FCz   NEG-CTR 13    RP    V       1.27  FALSE      FALSE          1.27 
# ... with 440 more rows

    This is how your sentence has been completed. Not only did we identify outliers, but we also applied capping to them. Now decide whether to use the value variable or the cap.value variable for further analysis. The decision is yours.

    A small update for a @little_statistician

    First, we will load all your data.

    #Loading libraries
library(tidyverse)
library(rstatix)
library(ggpubr)
library(readxl)

#Upload data
df_join <- read_excel("df_join.xlsx")

df = df_join  %>%
  mutate_at(vars(ID:COND), factor) %>%
  pivot_longer(P3FCz:LPP2Pz, names_to = "signals") %>%
  group_by(signals) %>%
  nest()

    Now let's define the fOutCapp and fOutCappGroup functions once again. Note, in the original version of fOutCapp there is no need for the group_by function.

    fOutCapp = function(data){
  x = data$value
  qnt = quantile(x, probs=c(.25, .75), na.rm = T)
  caps = quantile(x, probs=c(.05, .95), na.rm = T)
  H = 1.5 * IQR(x, na.rm = T)
  He = 3 * IQR(x, na.rm = T)
  is.outlier = (x < (qnt[1] - H)) | (x > (qnt[2] + H))
  x[x < (qnt[1] - H)] <- caps[1]
  x[x > (qnt[2] + H)] <- caps[2]
  data %>%  
    mutate(
      is.outlier = is.outlier,
      is.extreme = (x < (qnt[1] - He)) | (x > (qnt[2] + He)),
      cap.value = x
    )
}

fOutCappGroup = function(data) data %>% group_by(COND) %>% 
  group_modify(~fOutCapp(.x))

    Now is the time to mutate.

    df = df %>% group_by(signals) %>% 
  mutate(data = map(data, ~fOutCappGroup(.x))) %>% 
  unnest(data) %>% # step 1
  mutate(old.value = value,
         value = cap.value) %>% #Step 2
  nest(data=COND:old.value)  #Step 3

    It is very important that you understand what is really going on here. So in step 1 we group your tibble by the signals variable. It is simple and you certainly understand it. In step 2 we mutate the data variable, which is a list consisting of data for individual signals.

    output after step 2

    # A tibble: 12 x 2
# Groups:   signals [12]
   signals     data                 
   <chr>       <list>               
 1 P3FCz       <grouped_df [75 x 8]>
 2 P3Cz        <grouped_df [75 x 8]>
 3 P3Pz        <grouped_df [75 x 8]>
 4 LPPearlyFCz <grouped_df [75 x 8]>
 5 LPPearlyCz  <grouped_df [75 x 8]>
 6 LPPearlyPz  <grouped_df [75 x 8]>
 7 LPP1FCz     <grouped_df [75 x 8]>
 8 LPP1Cz      <grouped_df [75 x 8]>
 9 LPP1Pz      <grouped_df [75 x 8]>
10 LPP2FCz     <grouped_df [75 x 8]>
11 LPP2Cz      <grouped_df [75 x 8]>
12 LPP2Pz      <grouped_df [75 x 8]>

    This way your inner tibbles have gained new variables. You will see it after the unnest in step 3.

    output after step 3

    # A tibble: 900 x 9
# Groups:   signals [12]
   signals COND    ID    GR    SES     value is.outlier is.extreme cap.value
   <chr>   <fct>   <fct> <fct> <fct>   <dbl> <lgl>      <lgl>          <dbl>
 1 P3FCz   NEG-CTR 01    RP    V     -11.6   FALSE      FALSE        -11.6  
 2 P3FCz   NEG-CTR 04    RP    V      -0.314 FALSE      FALSE         -0.314
 3 P3FCz   NEG-CTR 06    RP    V      -0.214 FALSE      FALSE         -0.214
 4 P3FCz   NEG-CTR 07    RP    V      -2.83  FALSE      FALSE         -2.83 
 5 P3FCz   NEG-CTR 08    RP    V       4.24  FALSE      FALSE          4.24 
 6 P3FCz   NEG-CTR 09    RP    V       9.57  FALSE      FALSE          9.57 
 7 P3FCz   NEG-CTR 10    RP    V      -6.13  FALSE      FALSE         -6.13 
 8 P3FCz   NEG-CTR 11    RP    V       0.529 FALSE      FALSE          0.529
 9 P3FCz   NEG-CTR 12    RP    V      -7.74  FALSE      FALSE         -7.74 
10 P3FCz   NEG-CTR 13    RP    V       1.27  FALSE      FALSE          1.27 
# ... with 890 more rows

    And since you already have a very nice function that generates beautiful boxplot-violin plots with different stats, let's do one small mutation (step 4) replacing value with cap.value.

    output after step 4

    # A tibble: 900 x 10
# Groups:   signals [12]
   signals COND    ID    GR    SES     value is.outlier is.extreme cap.value old.value
   <chr>   <fct>   <fct> <fct> <fct>   <dbl> <lgl>      <lgl>          <dbl>     <dbl>
 1 P3FCz   NEG-CTR 01    RP    V     -11.6   FALSE      FALSE        -11.6     -11.6  
 2 P3FCz   NEG-CTR 04    RP    V      -0.314 FALSE      FALSE         -0.314    -0.314
 3 P3FCz   NEG-CTR 06    RP    V      -0.214 FALSE      FALSE         -0.214    -0.214
 4 P3FCz   NEG-CTR 07    RP    V      -2.83  FALSE      FALSE         -2.83     -2.83 
 5 P3FCz   NEG-CTR 08    RP    V       4.24  FALSE      FALSE          4.24      4.24 
 6 P3FCz   NEG-CTR 09    RP    V       9.57  FALSE      FALSE          9.57      9.57 
 7 P3FCz   NEG-CTR 10    RP    V      -6.13  FALSE      FALSE         -6.13     -6.13 
 8 P3FCz   NEG-CTR 11    RP    V       0.529 FALSE      FALSE          0.529     0.529
 9 P3FCz   NEG-CTR 12    RP    V      -7.74  FALSE      FALSE         -7.74     -7.74 
10 P3FCz   NEG-CTR 13    RP    V       1.27  FALSE      FALSE          1.27      1.27 
# ... with 890 more rows

    Finally, let's roll it all back to its original form with the variable data in step 5.

    output after step 5

    # A tibble: 12 x 2
# Groups:   signals [12]
   signals     data             
   <chr>       <list>           
 1 P3FCz       <tibble [75 x 9]>
 2 P3Cz        <tibble [75 x 9]>
 3 P3Pz        <tibble [75 x 9]>
 4 LPPearlyFCz <tibble [75 x 9]>
 5 LPPearlyCz  <tibble [75 x 9]>
 6 LPPearlyPz  <tibble [75 x 9]>
 7 LPP1FCz     <tibble [75 x 9]>
 8 LPP1Cz      <tibble [75 x 9]>
 9 LPP1Pz      <tibble [75 x 9]>
10 LPP2FCz     <tibble [75 x 9]>
11 LPP2Cz      <tibble [75 x 9]>
12 LPP2Pz      <tibble [75 x 9]>

    Well now let's make a graph!

    #Function to special boxplot3
SpecBoxplot3 = function(data, signal, parametric = FALSE, autor = "G. Anonim"){
  if(parametric) {
    pwc = data %>%
      pairwise_t_test(value~COND, paired = TRUE,
                      p.adjust.method = "bonferroni") %>%
      add_xy_position(x = "COND") %>%
      mutate(COND="NEG-CTR",
             lab = paste(p, " - ", p.adj.signif))
    res.test = data %>% anova_test(value~COND)
  } else {
    pwc = data %>% pairwise_wilcox_test(value~COND) %>%
      add_xy_position(x = "COND") %>%
      mutate(COND="NEG-CTR",
             lab = paste(p, " - ", p.adj.signif))
    res.test = data %>% kruskal_test(value~COND)
  }
  
  data %>% ggplot(aes(COND, value, fill=COND))+
    geom_violin(alpha=0.2)+
    geom_boxplot(outlier.shape = 23,
                 outlier.size = 3,
                 alpha=0.6)+
    geom_jitter(shape=21, width =0.1)+
    stat_pvalue_manual(pwc, step.increase=0.05, label = "lab")+
    ylab(signal)+
    labs(title = get_test_label(res.test, detailed = TRUE),
         subtitle = get_pwc_label(pwc),
         caption = autor)
}


#special boxplot for the P3FCz signal
df$data[[1]] %>% SpecBoxplot3("P3FCz", TRUE)
df$data[[1]] %>% SpecBoxplot3("P3FCz", FALSE)

    enter image description here

    As you can see on the chart, there are no outliers anymore!

    Now we're ready to plot each signal!

    #A function that creates a special boxplot3 and adds it to a data frame
AddSignalBoxplot3 = function(df, signal, printPlot=TRUE) {
  plot1 = SpecBoxplot3(df$data[[1]], signal, TRUE)
  plot2 = SpecBoxplot3(df$data[[1]], signal, FALSE)
  if(printPlot) print(plot1)
  if(printPlot) print(plot2)
  df %>% mutate(boxplot1 = list(plot1),
                boxplot2 = list(plot2),
  )
}

#Added special boxplot3
df %>% group_by(signals) %>%
  group_modify(~AddSignalBoxplot3(.x, .y))

    enter image description here

    Good luck on your further analysis !!

    Last update

    create.plot2 = function(df, group){
  data = df$data[[1]]
  minv = min(data$value)
  maxv = max(data$value)
  df.stat = data %>% group_by(COND) %>% 
    summarise(
      n = n(),
      mean = mean(value),
      sd = sd(value),
      min = minv,
      max = maxv,
      x = seq(min, max, length.out = n*100),
      value = dnorm(x, mean, sd) 
    ) 
  data %>% ggplot(aes(value, fill=COND))+
    geom_histogram(aes(y=..density..), colour="black", fill="white", bins = 30)+
    geom_density(alpha=.2, fill="red", col="red")+
    geom_line(aes(x, value), data=df.stat, col="blue")+
    xlab(group)+
    facet_grid(cols = vars(COND))
}

df %>% group_by(signals) %>% 
  group_map(create.plot2)

    enter image description here