<jupyter_start><jupyter_text>Cross Validation[<NAME>](https://natematias.com), April 21, 2020This notebook is part of [COMM 4940: The Design & Governance of Field Experiments](https://natematias.com/courses/comm4940/)In today's class, we're discussing the idea of cross-validation, where we confirm the reproducibility of statistical tests by holding back data. The reading for today is here:* <NAME>., <NAME>., & <NAME>. (2018). [Cross-validation approaches for replicability in psychology](https://www.frontiersin.org/articles/10.3389/fpsyg.2018.01117/full). Frontiers in Psychology, 9, 1117.Note that for this class, we're going to be thinking about cross-validation differently from prediction tasks (which is how it's often discussed in machine learning). Here, we're focused on meta-analysis, leading into analysis of the Upworthy Research Archive:* With prediction, you test the error in the overall model* With meta-analysis, it's fine if the model predicts things badly, so long as we learn the average treatment effect. In fact, adding multiple predictors might lead us to mis-estimate the average treatment effect.<jupyter_code>options("scipen"=9, "digits"=4)
library(MASS)
library(ggplot2)
library(rlang)
library(gmodels)  #includes CrossTable
library(corrplot)
library(stringr)
library(plm)      #Fixed effects models

## DOCUMENTATION AT: https://cran.r-project.org/web/packages/DeclareDesign/DeclareDesign.pdf
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
options(repr.plot.width=7, repr.plot.height=3.5)

## Set Random Seed:
set.seed(35985)<jupyter_output><empty_output><jupyter_text>Minimizing False Positives with Cross-ValidationCross-validation can help us minimize false positives in research.In the following example, we create a dataset with 100 normal distributions that aren't correlated in any way. We then see how many of those correlations survive cross validation<jupyter_code># total sample size
n.size = 1000

## create dataframe
nocorr.df <- data.frame(
    id = seq(1, n.size),
    a = rnorm(n.size),
    b = rnorm(n.size),
    c = rnorm(n.size),
    d = rnorm(n.size),
    e = rnorm(n.size),
    f = rnorm(n.size),
    g = rnorm(n.size),
    h = rnorm(n.size),
    i = rnorm(n.size),
    j = rnorm(n.size)
)<jupyter_output><empty_output><jupyter_text>Creating the Exploratory and Confirmatory (holdout) DatasetsHere, I create two separate, independent datasets from the larger data:* `nocorr.exploratory.df` is a 30% dataset for exploratory analysis* `nocorr.confirmatory.df` is a 70% holdout dataset for confirmatory analysisThere's no hard and fast rule for how large the datasets should be, though there are some tradeoffs, as I illustrate at the end of the notebook.<jupyter_code>## the proportion of observations
## to include in the exploratory dataset
exp.proportion = 0.3

nocorr.exploratory.df <- nocorr.df[
    sample(1:nrow(nocorr.df), n.size*exp.proportion, replace=FALSE),]

nocorr.confirmatory.df <- subset(nocorr.df, (id %in% nocorr.exploratory.df$id)!=TRUE)

print(paste(nrow(nocorr.exploratory.df), "in the exploratory dataset"))
print(paste(nrow(nocorr.confirmatory.df), "in the confirmatory dataset"))<jupyter_output>[1] "300 in the exploratory dataset"
[1] "700 in the confirmatory dataset"<jupyter_text>Developing a Model with the Exploratory Dataset<jupyter_code># Let's imagine I was looking for just any correlation in the dataset, with the goal of
# getting some kind of statistically-significant result. Here, we see that both b and d 
# are associated with g

ndt <- cor.mtest(nocorr.exploratory.df, conf.level = .95)
ndc <- cor(nocorr.exploratory.df)
corrplot(ndc, p.mat=ndt$p, sig.level=0.05)
# in this linear regression model, we see that both b and d are associated with g
summary(lm(g ~ b + d, data=nocorr.exploratory.df))<jupyter_output><empty_output><jupyter_text>Testing the spurious correlations with the confirmatory datasetWhen we test the correlations on a (larger) sample from the population, we fail to reject the null hypothesis that these results are too large to have arisen from chance.<jupyter_code>summary(lm(g ~ b + d, data=nocorr.confirmatory.df))<jupyter_output><empty_output><jupyter_text>Handling Overfitting with Cross-ValidationSometimes the best model that fits the data is so specific that it doesn't apply to anything beyond the data. That's what researchers call overfitting (and it's what happened in the above example too).By giving researchers a chance to test their model on more than one dataset, cross-validation enables researchers to confirm whether their model is overfit.For this example, we are using the [titanic dataset](http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic.html), a record of who survived a profound tragedy in 1912 when 1,500 out of 2,208 people [died when the H.M.S. Titanic sunk in the Atlantic](https://en.wikipedia.org/wiki/Passengers_of_the_RMS_Titanic). This dataset includes information about 1309 passengers ([description here](http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3info.txt)).In this case, I want to know if people who paid more money for their ticket (`fare`) had a higher chance of survival, holding all else equal.<jupyter_code>## Reset Random Seed:
set.seed(69982)

titanic.df <- read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.csv")

titanic.df$is.minor <- as.numeric(titanic.df$age < 18)

t.exploratory.df <- titanic.df[
    sample(1:nrow(titanic.df), nrow(titanic.df)*exp.proportion, replace=FALSE),]

t.confirmatory.df <- subset(titanic.df, (name %in% t.exploratory.df$name)!=TRUE)

print(paste(nrow(t.exploratory.df), "in the exploratory dataset"))
print(paste(nrow(t.confirmatory.df), "in the confirmatory dataset"))
colnames(t.exploratory.df)<jupyter_output>[1] "392 in the exploratory dataset"
[1] "917 in the confirmatory dataset"<jupyter_text>First,  Develop an Exploratory ModelUsing the exploratory dataset, I can create a model that tests for the correlation I'm curious about. Here, I use the classic exploratory technique of fitting a sequence of models until I find the one that has the highest R2.The key hypothesis test will be whether `fare` is a statistically-significant predictor of survival, accounting for all other variables.<jupyter_code>summary(lm(survived ~ fare, data=t.exploratory.df))
summary(lm(survived ~ fare + pclass, data=t.exploratory.df))
summary(lm(survived ~ fare + pclass + sex, data=t.exploratory.df))
summary(lm(survived ~ fare + pclass + sex + age, data=t.exploratory.df))<jupyter_output><empty_output><jupyter_text>Here's the final model, which adjusts for passenger class, sex, age, and the number of siblings or parents on the boat.<jupyter_code>summary(lm(survived ~ fare + pclass + sex + age + sibsp, data=t.exploratory.df))<jupyter_output><empty_output><jupyter_text>Next test the hypothesis in the confirmatory datasetAt this point, there are two approaches that you could take when doing hypothesis tests (in contrast with prediction tasks, where we are testing the overall error of the model). * Test the hypotheses on the confirmatory dataset alone* Test the hypotheses on the combined confirmatory + exploratory dataset (the full dataset) (this is the approach we will take here)Here, when we look for a correlation within the full dataset, we fail to reject the null hypothesis. Accounting for other factors, in this model, there was no relationship between chance of survival and the fare that someone paid.<jupyter_code>summary(lm(survived ~ fare + pclass + sex + age + sibsp, data=titanic.df))<jupyter_output><empty_output><jupyter_text>So which is the case? Is there a relationship or not? In one sample, I fount support for my hypothesis. In the other, I did not. Which is right? My failure to reject the null hypothesis in the confirmatory dataset might be because there's no true relationship in the population, or it might be the luck of the draw. Cross-validation can also yield a false negative, where you miss an important relationship that is actually there. Imagine that the exploratory dataset didn't yield a statistically-significant relationship but there was genuinely a relationship in the population of Titanic passengers. And imagine that I hadn't gone on to test that hypothesis because it wasn't statistically-significant in the exploratory dataset. I might have missed an important relationship.In the case of the Titanic and the Upworthy Research Archive, we can't go back and increase the sample. If prediction were our goal, we could calculate root mean square error using the following techniques:* K-fold cross validation, where we create a series of exploratory sets and then test them against the confirmatory sets* Leave-one-out cross-validation, where we iterate through all possible samples and test the prediction against the left-out observationIn all of these cases, we would use measures like root mean square error to evaluate the goodness of fit of the model.  Cross-Validation of a Hypothesis from the Upworthy Research ArchiveIn this example, I use the upper bound exploratory and confirmatory datasets of experiments naming notable people to:* develop a hypothesis with the exploratory dataset* test it with the full dataset<jupyter_code>data_dir="../assignments/upworthy-archive-project/"

max.exp.df <- read.csv(paste(data_dir,"upworthy_archive_exploratory_max_effect_size_dataset.csv", sep=""))
max.conf.df <- read.csv(paste(data_dir,"upworthy_archive_confirmatory_max_effect_size_dataset.csv", sep=""))

headlines.df <- read.csv(paste(data_dir,"headlines.csv", sep=""))<jupyter_output><empty_output><jupyter_text>Fit Exploratory ModelHere, I fit a fixed effects model on the hypothesis that including a notable person's name in a headline increases the chance that someone will click on the headline. In the upper bound dataset, across 40 tests, including a notable person's name in the headline increases the chance of someone clicking by 3 tenths of a percentage point.<jupyter_code>summary(plm(clicked ~ has_treatment,
            index="clickability_test_id",
            model="within",
            data=max.exp.df))<jupyter_output><empty_output><jupyter_text>What does exploration look like for the Upworthy Research Archive?In this case, how might we adjust our model based on this exploratory data analysis, since we only have one predictor?* Maybe headline selection isn't precise enough:    * If we don't have enough names in our list, maybe we're inadvertently comparing two headlines that both have notable names    * Maybe the comparison headlines (without notable names) are very different in some other way that skews the results, for example have a different length or more positive/negative phrasing* Maybe there's a `mediating` factor, such as how famous the person is. Maybe we need to add another variable for how popular they were at the timeThese are all areas where further adjustment and exploration could improve the model or could lead to overfitting.  Testing the Confirmatory HypothesisIn the next step, it's important to be very precise about what we are confirming. Options include:* a hypothesis about the existence of an effect (the hypothesis that has_treatment is statistically-significant)* a hypothesis about the sign of the effect (the hypothesis that has_treatment is greater than zero and statistically-significant)* a prediction about magnitude of the estimate (the hypothesis that has_treatment is at least X and statistically-significant)In this case, I'm choosing the following hypothesis test (which I did actually decide on before testing against the confirmatory dataset):**Hypothesis:** including a notable person's name in a headline as defined in `selecting_upworthy_archive_packages.py`, has an upper bound effect that is positive and statistically significant.Here, since we're doing meta-analysis and want to include all possible information, I'm *combining* the exploratory and confirmatory datasets into one final model.**Result:** the hypothesis is confirmed. On average across the upper-bound dataset of 297 A/B tests (and 2.1 million viewers), including a notable person's name in a headline increased the chance that someone clicked by 4 tenths of a percentage point, an effect that is positive and statistically-significant.<jupyter_code>max.df <- rbind(max.exp.df, max.conf.df)
summary(plm(clicked ~ has_treatment,
            index="clickability_test_id",
            model="within",
            data=max.df))<jupyter_output><empty_output><jupyter_text>Other ways to cross-validateI used this technique because it is relatively straightforward and involves just expanding the available data. This is what I will be doing for the final.For meta-analysis, I'll be honest that I'm not (yet) sure how to conduct the other techniques, since the exploratory process is focused on defining variables and selecting the tests that go into the dataset.I'm not currently aware of anyone else who has used cross-validation to test hypotheses from a large dataset of experiments (though I bet this has been done in medicine).For correlational analyses like the Titanic example, k-fold and leave-one-subject-out are simpler because you can compare the Root Mean Square Error of a model with the key predictor (`fare`) to one without the predictor. If your preferre model is a better fit, then you uphold the finding from that model.  Software EnvironmentThis example was generated using the following R configurations and libraries.<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07782272894446488, 'fraction_numerical': 0.01147480691430673, 'mean_word_length': 5.359214218896165, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 55, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '345520', 'license': 'MIT', 'n_tokens_mistral': 3913, 'n_tokens_neox': 3630, 'n_words': 1622, 'path': '2020/lecture-code/lecture-21-cross-validation.ipynb', 'repo_name': 'natematias/Cornell-Experimenting-Society', 'size': '13614'}