Apply function to each level of a factor (participant) to remove outliers based on distance from mean in standard deviation in R

Viewed 234

I have a data.frame in R with a column indicating participants subject in an experiment, another column indicating condition, trial_type, and a final column with my numeric dependent variable rt.

Here's a reproducible example of my data generated with the dput() function for the first 64 rows of my data.frame:

 structure(list(subject = structure(c(21L, 21L, 21L, 21L, 21L, 
21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 24L, 24L, 
24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 
24L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 
27L, 27L, 27L, 27L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 
47L, 47L, 47L, 47L, 47L, 47L, 47L), .Label = c("1p12", "1p13", 
"1p15", "1p30", "1p36", "1p39", "1p43", "1p46", "1p49", "1p59", 
"1p60", "1p67", "1p69", "1p79", "1p80", "1p81", "1p84", "1p85", 
"1p88", "1p9", "2p1", "2p11", "2p18", "2p2", "2p22", "2p25", 
"2p3", "2p31", "2p33", "2p42", "2p44", "2p5", "2p50", "2p58", 
"2p63", "2p72", "2p76", "2p78", "2p8", "2p83", "3p10", "3p16", 
"3p20", "3p28", "3p32", "3p34", "3p4", "3p41", "3p47", "3p54", 
"3p55", "3p56", "3p61", "3p64", "3p66", "3p7", "3p77", "3p82", 
"3p86", "3p87", "4p14", "4p17", "4p19", "4p21", "4p24", "4p26", 
"4p27", "4p29", "4p35", "4p37", "4p38", "4p48", "4p51", "4p57", 
"4p6", "4p62", "4p68", "4p70", "4p74", "4p75"), class = "factor"), 
    rt = c(4303L, 5616L, 1317L, 1663L, 1353L, 645L, 648L, 457L, 
    2359L, 2497L, 832L, 523L, 1427L, 511L, 483L, 1300L, 873L, 
    1185L, 1752L, 2037L, 4849L, 2975L, 1621L, 1235L, 3008L, 1560L, 
    1075L, 4596L, 1129L, 1093L, 1302L, 1414L, 5542L, 2369L, 4944L, 
    2338L, 1274L, 1837L, 3384L, 1338L, 2002L, 1756L, 2516L, 1868L, 
    2017L, 1337L, 1106L, 1388L, 6812L, 5579L, 1695L, 1976L, 1897L, 
    4484L, 3095L, 1865L, 2283L, 1659L, 1328L, 1882L, 1483L, 1993L, 
    1776L, 2256L), condition = structure(c(1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("reliable", 
    "unreliable"), class = "factor"), trial_type = structure(c(2L, 
    2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 
    2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 
    2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 
    1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 
    1L, 2L, 2L), .Label = c("same", "switch"), class = "factor"), 
    accuracy = c(0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 
    1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 
    0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 
    1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 
    1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L)), row.names = c(NA, -64L), class = c("tbl_df", 
"tbl", "data.frame"))

I want to remove outliers based on the distance of each score in rt from the mean (in standard deviation units) per participant. Each participant has 16 rows, and there are 80 participants. R correctly interprets each of the 16 rows for a given participant as a single level of a factor.

I want to use this to remove outliers from rt separately for each level in subject and for each level in trial_type. The formula I used to remove outliers across all participants was:

# Calculate mean and `sd` for each level of `trial_type`

# For 'same':

MeanSame <- mean(RTs$rt[RTs$trial_type == "same"])
SDSame <- sd(RTs$rt[RTs$trial_type == "same"])

# For 'switch':

MeanSwitch <- mean(RTs$rt[RTs$trial_type == "switch"])
SDSwitch <- sd(RTs$rt[RTs$trial_type == "switch"])

# Create upper and lower cut for level 'same' of 'trial_type':
UpperSame <- MeanSame + 2.5*SDSame
LowerSame <- MeanSame - 2.5*SDSame

# Create upper and lower cut for level 'switch' of 'trial_type':
UpperSwitch <-2.5*SDSwitch + MeanSwitch
LowerSwitch <- MeanSwitch - 2.5*SDSwitch


#Identify Outliers in Same

OutliersSameUpper <- which(RTs$rt > UpperSame & RTs$trial_type == "same") 

OutliersSameLower <- which(RTs$rt < LowerSame & RTs$trial_type == "same")


# Identify Outliers in Switch

OutliersSwitchUpper <- which(RTs$rt > UpperSwitch  & RTs$trial_type == "switch") 

OutliersSwitchLower <- which(RTs$rt < LowerSwitch & RTs$trial_type == "switch")


# Create new data.frame without the identified outliers:

RTsClean <- RTs[-c(OutliersSameUpper,OutliersSameLower,OutliersSwitchUpper,OutliersSwitchLower),]

The way I did this, for all participants, was to calculate the mean and SD for each condition and then find the rows in rt that were above or beyond the cutting points. However, I cannot figure out how I could do this for each level of the subject vector. the apply functions are not working for me, because it is not enough to apply a single function. I need to keep track of other variables as well (the mean and sd, and the upper and lower cut points).

It seems to me the way to go would be to create a function, or a for loop, or both. But that is beyond my skill level in R.

I would appreciate anyone who could help me find the best way to apply the method for outlier removal I specified above, but in a way that can be applied separately to each level of the subject vector. This means that the mean and sd variables need to be specified for each level of subject.

If possible, it would be even better to create a new column in the data.frame that specified, for each row in rt, whether that row was considered to be an outlier or not. But I don't have a clue how I could achieve that.

Thanks in advance for the help.

3 Answers

based on this requirement of OP

I would appreciate anyone who could help me find the best way to apply the method for outlier removal I specified above, but in a way that can be applied separately to each level of the subject vector. This means that the mean and sd variables need to be specified for each level of subject.

It removes three rows as outliers


library(dplyr)

RTs %>% group_by(subject) %>%
  filter(rt <= mean(rt) + (2.5 * sd(rt)),  rt >= mean(rt) - (2.5 * sd(rt)))

#> # A tibble: 61 x 5
#> # Groups:   subject [4]
#>    subject    rt condition trial_type accuracy
#>    <fct>   <int> <fct>     <fct>         <int>
#>  1 2p1      4303 reliable  switch            0
#>  2 2p1      1317 reliable  switch            0
#>  3 2p1      1663 reliable  same              1
#>  4 2p1      1353 reliable  switch            1
#>  5 2p1       645 reliable  same              1
#>  6 2p1       648 reliable  same              1
#>  7 2p1       457 reliable  same              1
#>  8 2p1      2359 reliable  switch            0
#>  9 2p1      2497 reliable  switch            0
#> 10 2p1       832 reliable  same              1
#> # ... with 51 more rows

Based on this further requirement of OP

If possible, it would be even better to create a new column in the data.frame that specified, for each row in rt, whether that row was considered to be an outlier or not. But I don't have a clue how I could achieve that.

RTs %>%
  group_by(subject) %>%
  mutate(OUTLIER = rt >= mean(rt) + (2.5 * sd(rt)) | rt <= mean(rt) - (2.5 * sd(rt)))

# A tibble: 64 x 6
# Groups:   subject [4]
   subject    rt condition trial_type accuracy OUTLIER
   <fct>   <int> <fct>     <fct>         <int> <lgl>  
 1 2p1      4303 reliable  switch            0 FALSE  
 2 2p1      5616 reliable  switch            0 TRUE   
 3 2p1      1317 reliable  switch            0 FALSE  
 4 2p1      1663 reliable  same              1 FALSE  
 5 2p1      1353 reliable  switch            1 FALSE  
 6 2p1       645 reliable  same              1 FALSE  
 7 2p1       648 reliable  same              1 FALSE  
 8 2p1       457 reliable  same              1 FALSE  
 9 2p1      2359 reliable  switch            0 FALSE  
10 2p1      2497 reliable  switch            0 FALSE  
# ... with 54 more rows

Clearly TRUE means it is outlier whereas FALSE indicate otherwise


BaseR way

RTs$outlier <- as.logical(ave(RTs$rt, RTs$subject, 
               FUN = function(.x) (.x >= mean(.x) + 2.5 * sd(.x)) | 
                 (.x <= mean(.x) - 2.5 * sd(.x))))
RTs

# A tibble: 64 x 6
   subject    rt condition trial_type accuracy outlier
   <fct>   <int> <fct>     <fct>         <int> <lgl>  
 1 2p1      4303 reliable  switch            0 FALSE  
 2 2p1      5616 reliable  switch            0 TRUE   
 3 2p1      1317 reliable  switch            0 FALSE  
 4 2p1      1663 reliable  same              1 FALSE  
 5 2p1      1353 reliable  switch            1 FALSE  
 6 2p1       645 reliable  same              1 FALSE  
 7 2p1       648 reliable  same              1 FALSE  
 8 2p1       457 reliable  same              1 FALSE  
 9 2p1      2359 reliable  switch            0 FALSE  
10 2p1      2497 reliable  switch            0 FALSE  
# ... with 54 more rows

Here you can get your results with group_by and summarise. You can adapt -> whatever you need:

library(dplyr)
RTs %>% 
  group_by(subject, trial_type) %>% 
  summarise(mean= mean(rt), sd = sd(rt), Upper = 2.5*sd + mean, Lower = mean - 2.5*sd,
            OutliersUpper <- rt > Upper, OutliersLower <- rt < Lower)

Output:

   subject trial_type  mean    sd Upper  Lower `OutliersUpper <- rt > Upper` `OutliersLower <- rt < Lower`
   <fct>   <fct>      <dbl> <dbl> <dbl>  <dbl> <lgl>                         <lgl>                        
 1 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 2 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 3 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 4 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 5 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 6 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 7 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 8 2p1     same        720.  400. 1720.  -280. FALSE                         FALSE                        
 9 2p1     switch     2522. 1616. 6562. -1519. FALSE                         FALSE    

You can adapt: for example

RTs %>% 
  group_by(subject, trial_type, condition, accuracy) %>% 
  summarise(mean= mean(rt), sd = sd(rt), Upper = 2.5*sd + mean, Lower = mean - 2.5*sd,
            OutliersUpper <- rt > Upper, OutliersLower <- rt < Lower)

A generalised version of the other answers.

First, convert your existing filtering to a function:

dropOutliers <- function(.x, .y, multiplier=2.5) {
  limits <- .x %>% 
              group_by(trial_type) %>% 
              summarise(
                SD=sd(rt), 
                Mean=mean(rt), 
                Lower=Mean - multiplier * SD, 
                Upper=Mean + multiplier * SD
              )
  .x %>% 
    left_join(
      limits, 
      by="trial_type"
    ) %>% 
    filter(rt >= Lower && rt <= Upper) %>% 
    select(-SD, -Mean, -Upper, -Lower)
}

I've used .x and .y as parameter names as they match the documentation of group_map(), which I will use in a moment. multiplier is for testing.

Now apply the function to the data frame grouped in any way you wish. group_map returns a list of tibbles, so row_bind() the return value to get a single, combined tibble as required.

d %>% group_by(subject) %>% group_map(dropOutliers) %>% bind_rows()

This doesn't seem to drop any rows, so apply more agressive filtering to check:

d %>% group_by(subject) %>% group_map(dropOutliers, multiplier=1) %>% bind_rows()
# A tibble: 16 x 4
      rt condition trial_type accuracy
   <int> <fct>     <fct>         <int>
 1   873 reliable  switch            1
 2  1185 reliable  same              1
 3  1752 reliable  same              1
 4  2037 reliable  same              1
 5  4849 reliable  switch            1
 6  2975 reliable  switch            0
 7  1621 reliable  switch            0
 8  1235 reliable  same              1
 9  3008 reliable  switch            0
10  1560 reliable  same              1
11  1075 reliable  switch            0
12  4596 reliable  same              1
13  1129 reliable  same              1
14  1093 reliable  switch            0
15  1302 reliable  same              1
16  1414 reliable  switch            0
Related