Analyzing Experimental Data¶
%%R
suppressMessages({
# install.packages(c('dplyr', 'tidyr', 'ggplot2', 'broom'))
library(dplyr)
library(tidyr)
library(ggplot2)
library(broom)
})
Exposure logs¶
First we logged everytime a user loaded up the user interface. Often we want to only take the first exposure by a user.
%%R
src2ideology <- c(
# <=liberal=====moderate====conservative=>
'msnbc' = -1, 'cnn' = 0, 'foxnews' = 1)
Exp <- read.csv('data/exposure_mturk2_log.csv', stringsAsFactors = FALSE) %>%
group_by(userid, story_key) %>%
arrange(time) %>%
summarise(source = first(source), position = first(position), time = first(time)) %>%
mutate(source.ideology = src2ideology[source]) %>%
ungroup
cat('We have', nrow(Exp), 'exposures.', '\n')
cat('We have', with(Exp, length(unique(userid))), 'potential subjects.', '\n\n')
head(Exp)
%%R
Exp %>% group_by(source.ideology) %>% summarise(n = n())
Events logs¶
The events log is produced by user actions that take place sometime after the exposure. So in this experiment we logged two kinds of these: summaries and survey responses.
%%R
Sum <- read.csv('data/summary_mturk2_log.csv', stringsAsFactors = FALSE) %>%
group_by(userid, story_key) %>%
arrange(time) %>%
summarise(summary = first(summary)) %>%
ungroup
cat('We have', nrow(Sum), 'summaries.', '\n')
cat('We have', with(Sum, length(unique(userid))), 'subjects who wrote summaries.', '\n')
cat(with(Sum, round(mean(summary), 2)), 'stories were summarized.', '\n\n')
head(Sum)
Survey responses¶
After the experiment, we survey the subjects about their political ideology.
%%R
Survey <- read.csv('data/survey_mturk2_log.csv', stringsAsFactors = FALSE) %>%
filter(!is.na(response)) %>%
group_by(userid) %>%
arrange(time) %>%
summarise(response = first(response), reversed_scale = first(reversed_scale)) %>%
# Original scale is [0, 4] (from conservative to liberal)
# let's convert it to [-2, 2] (liberal to conservative)
mutate(subject.ideology = 2 - response) %>%
ungroup
cat('We have ', nrow(Survey), 'survey responses.', '\n')
head(Survey)
Joining exposures to events¶
Typically we need to perform joins from the exposures to the outcomes in order to analyze the data.
%%R
Results <- Exp %>%
inner_join(Sum, c('userid', 'story_key')) %>%
inner_join(Survey, 'userid') %>%
mutate(opposed = as.numeric((subject.ideology * source.ideology) < 0)) %>%
group_by(userid) %>%
mutate(num_stories = n(), num_summaries = sum(summary)) %>%
ungroup %>%
# We only study the people who followed instructions and have a political affiliation.
filter(num_stories == 3, subject.ideology != 0)
with(Results, summary(factor(opposed)))
EDA¶
It's good to do some EDA on your data to make sure all the data we collected looks good.
%%R
# When did people show up to do the experiment?
ggplot(Exp, aes(x = time)) +
geom_histogram() +
theme_bw()
%%R
# Check that we have balanced exposures.
Results %>%
group_by(source.ideology) %>%
summarise(n = n(), num.users = length(unique(userid))) %>%
ggplot(aes(x = source.ideology, y = n)) +
geom_bar(stat = 'identity') +
xlab('(Manipulated) Source Ideology') +
ylab('Number of observations') +
theme_bw()
%%R
# Summary rate by story to see if any stories are weird.
Sum %>%
group_by(story_key) %>%
summarise(k = sum(summary), n = n(),
p = mean(summary),
lcl = binom.test(k, n)$conf.int[1],
ucl = binom.test(k, n)$conf.int[2]) %>%
ggplot(aes(x = factor(story_key), y = p, ymin = lcl, ymax = ucl)) +
geom_pointrange() +
xlab('Story ID') +
ylab('Probability of Summary') +
theme_bw()
%%R
# Make sure our survey results look sane.
ideology.labels <- c('very liberal', 'liberal', 'moderate', 'conservative', 'very conservative')
Survey %>%
ggplot(aes(x = factor(subject.ideology, labels = ideology.labels), fill = subject.ideology)) +
geom_histogram() +
scale_fill_gradient(low = "blue", high = "red") + # haha
theme_bw() +
theme(legend.position = 'none') +
xlab('Political Affiliation')
Measuring the ATE¶
Now we can tabulate the selections by source ideology to measure the ATE.
%%R
Results %>%
group_by(opposed) %>%
summarise(p.summary = mean(summary),
measurement = 'p.summary') %>%
mutate(opposed = factor(opposed, labels = c('not_opposed', 'opposed'))) %>%
spread(opposed, p.summary) %>%
mutate(ate = opposed - not_opposed)
Modeling the Results¶
We usually prefer modeling the data using a regression. We'll show how to use linear regression to perform correct inference on data with dependence, control variables, and larger than you can fit into memory.
Outline¶
- i.i.d. case (with and without covariates)
- bootstrapping to deal with dependence
- using data reduction with weighted regressions
I.i.d. case¶
If your data are i.i.d (which is pretty rare), then you can use an OLS regression and trust the standard errors on the treatment variable that R produces.
This analysis is inappropriate for our experiment because of user- and story-level dependence across observations. Later we'll show how we account for this using a bootstrap.
%%R
lm(summary ~ opposed, data = Results) %>%
tidy %>%
ggplot(aes(x = term, y = estimate,
ymin = estimate - 1.96*std.error,
ymax = estimate + 1.96*std.error)) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
I.i.d. case with control variables¶
Including control variables in your regression shouldn't change the ATE estimate, but it might help with precision (giving you smaller confidence intervals).
%%R
# In this case we think the position should help explain variation in the outcome (which is plausible).
lm(summary ~ opposed + factor(position), data = Results) %>%
tidy %>%
ggplot(aes(x = term, y = estimate,
ymin = estimate - 1.96*std.error,
ymax = estimate + 1.96*std.error)) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
Bootstrapping to get confidence intervals¶
If your data exhibits dependence, then you want to use a bootstrap to get correct inference. In our experiment we resample both users and stories, so the true model has two way dependence. We'll start with just user-level dependence:
%%R
fitter <- function(.data) {
lm(summary ~ opposed, data = .data, weights = .weights) %>%
tidy
}
iid.replicates <- iid.bootstrap(Results, fitter, .combine = bind_rows)
oneway.replicates <- clustered.bootstrap(Results, c('userid'), fitter, .combine = bind_rows)
twoway.replicates <- clustered.bootstrap(Results, c('userid', 'story_key'), fitter, .combine = bind_rows)
head(iid.replicates)
%%R
replicates <- bind_rows(
iid.replicates %>% mutate(bootstrap = 'iid'),
oneway.replicates %>% mutate(bootstrap = 'oneway'),
twoway.replicates %>% mutate(bootstrap = 'twoway')
)
replicates %>%
group_by(term, bootstrap) %>%
summarise(y = mean(estimate), lcl = quantile(estimate, 0.025), ucl = quantile(estimate, 0.975)) %>%
ggplot(aes(x = bootstrap, y = y, ymin = lcl, ymax = ucl)) +
facet_grid(term ~ .) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
Using Data Reduction¶
If you have a lot of data, sometimes it's a lot easier to group it and fit weighted models. We show how to do that here.
%%R
Results.small <- Results %>%
group_by(summary, opposed) %>%
summarise(cases = n())
Results.small
Fitting Reduced Data¶
%%R
lm(summary ~ opposed, data = Results.small, weights = cases) %>%
coeftest(vcov = sandwich.lm) %>%
tidy %>%
ggplot(aes(x = term, y = estimate,
ymin = estimate - 1.96*std.error,
ymax = estimate + 1.96*std.error)) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
Reducing data when you have covariates¶
If your covariates are discrete and low-cardinality, you can still use data reduction.
%%R
Results.medium <- Results %>%
group_by(summary, opposed, position) %>%
summarise(cases = n())
lm(summary ~ opposed + factor(position), data = Results.medium, weights = cases) %>%
coeftest(vcov = sandwich.lm) %>%
tidy %>%
ggplot(aes(x = term, y = estimate,
ymin = estimate - 1.96*std.error,
ymax = estimate + 1.96*std.error)) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
Bootstrapping Reduced Data¶
If your data are dependent and large, you can combine data reduction with the bootstrap:
- Generate a bootstrap replicate
- Reduce that replicate
- Fit a model on the reduced replicate
%%R
reduced.fitter <- function(.data) {
reduced.data <- .data %>%
group_by(summary, opposed, position) %>%
summarise(cases = sum(.weights))
lm(summary ~ opposed + factor(position), data = reduced.data, weights = cases) %>%
tidy
}
twoway.replicates <- clustered.bootstrap(Results, c('userid', 'story_key'), reduced.fitter, .combine = bind_rows)
head(twoway.replicates)
%%R
twoway.replicates %>%
group_by(term) %>%
summarise(y = mean(estimate), lcl = quantile(estimate, 0.025), ucl = quantile(estimate, 0.975)) %>%
ggplot(aes(x = term, y = y, ymin = lcl, ymax = ucl)) +
geom_pointrange() +
geom_hline(yintercept = 0.0, linetype = 'dashed') +
coord_flip() +
theme_bw()
