normalize subject data

normalize.R

# install.packages('dplyr')
# install.packages('tidyr')
# install.packages('ggplot2')
library(dplyr)
library(tidyr)
library(ggplot2)

# make some trials
data = expand.grid(
  subjectID = seq(1, 1.5, length.out=4), 
  condition = 1:3,
  trial = 1:20
)

# response is random but dependent on the condition and the subject
data$response = rnorm(
  n = nrow(data), 
  mean = data$condition*2 * data$subjectID, 
  sd = 1
)

########## Normalize ###########

# get the average result for each condition per subject
dataBaseline = data %>%
  # uncomment this line to use one condition as the basedline
  #filter(condition == 1) %>%
  group_by(subjectID, condition) %>%
  summarize(responseBaseline = mean(response))

# get the average across conditions per subject
dataBaseline = dataBaseline %>%
  group_by(subjectID) %>%
  summarize(responseBaseline = mean(responseBaseline))

# distribute the aggregated value to each row and normalize
data = data %>% 
  left_join(dataBaseline) %>%
  mutate(responseNormalized = response / responseBaseline)


############ scalar x-axis (smooth takes care of stats) ################

ggplot(data, aes(x=condition, color=factor(subjectID), fill=factor(subjectID))) +
  geom_smooth(method='lm', aes(y=response)) +
  labs(title = 'not normalized')

ggplot(data, aes(x=condition, color=factor(subjectID), fill=factor(subjectID))) +
  geom_smooth(method='lm', aes(y=responseNormalized)) +
  labs(title = 'normalized')

############ categorical x-axis (aggregate yourself) ################

# categorize conditions
data$condition = LETTERS[data$condition]

# r doesn't have a standard error function
stderr = function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))
ci95 = function(x) stderr(x) * 1.96

# get the average per condition per subject
dataAggregated = data %>%
  group_by(condition, subjectID) %>%
  summarize(response = mean(response),
            responseNormalized = mean(responseNormalized))

# get the average and CI per condition (across subjects)
dataPerCondition = dataAggregated %>%
  group_by(condition) %>%
  summarize(response_ci95 = ci95(response),
            response_mean = mean(response),
            responseNormalized_ci95 = ci95(responseNormalized),
            responseNormalized_mean = mean(responseNormalized))

ggplot(dataPerCondition, aes(x=condition)) +
  geom_pointrange(size=1, aes(
    y=response_mean, 
    ymin=response_mean-response_ci95, 
    ymax=response_mean+response_ci95)) +
  labs(title = 'not normalized')

ggplot(dataPerCondition, aes(x=condition)) +
  geom_pointrange(size=1, aes(
    y=responseNormalized_mean, 
    ymin=responseNormalized_mean - responseNormalized_ci95, 
    ymax=responseNormalized_mean + responseNormalized_ci95)) +
  labs(title = 'normalized')

####### Different approach: Cousineau within-subject error bars #########

# get average across all conditions per subject
dataAggregated = dataAggregated %>%
  group_by(subjectID) %>%
  mutate(averagePerID = mean(response))

# get the average of those subbject averages
bigAverage = mean(dataAggregated$averagePerID)

# scale based on those values
dataAggregated = dataAggregated %>%
  mutate(response_Cousineau = response - averagePerID + bigAverage)

# get the average and CI per condition (across subjects)
dataPerCondition = dataAggregated %>%
  group_by(condition) %>%
  summarize(response_ci95 = ci95(response),
            response_mean = mean(response))

ggplot(dataPerCondition, aes(x=condition)) +
  geom_pointrange(size=1, aes(
    y=response_mean, 
    ymin=response_mean-response_ci95, 
    ymax=response_mean+response_ci95)) +
  labs(title = 'Within Subject Error Bars')