stargazer.Rmd

---
title: "Document to Produce Stargazer Table and analysis on data combined from all sources"
author: "Group"
output:
  pdf_document: default
  html_document: default
---

```{r}
# Load all data
dff <- read.csv("./cleaneddata/ff.csv")
dl <- read.csv("./cleaneddata/lucid.csv")
dmt <- read.csv("./cleaneddata/mturk.csv")

```

```{r}
# Adding columns for control, treatment and clusters (M, F)
t = c("FTreatment", "MTreatment")
s = c("FTreatment", "FControl")
dmt$treat = ifelse(dmt$FL_3_DO %in% t, 1, 0 )
dmt$female = ifelse(dmt$FL_3_DO %in% s, 1, 0)
dmt$trust = ifelse(dmt$Q7 == "I trust what the speaker was saying.", 1,
                   ifelse (dmt$Q6 != "Pritikin" | dmt$Q7 == "", NA, 0))

dl$treat = ifelse(dl$FL_3_DO %in% t, 1, 0 )
dl$female = ifelse(dl$FL_3_DO %in% s, 1, 0)
dl$trust = ifelse(dl$Q7 == "I trust what the speaker was saying.", 1,
                   ifelse (dl$Q6 != "Pritikin" | dl$Q7 == "", NA, 0))

dff$treat = ifelse(dff$FL_3_DO %in% t, 1, 0 )
dff$female = ifelse(dff$FL_3_DO %in% s, 1, 0)
dff$trust = ifelse(dff$Q7 == "I trust what the speaker was saying.", 1,
                   ifelse (dff$Q6 != "Pritikin" | dff$Q7 == "", NA, 0))
```

```{r}
dall <- rbind(dff, dl, dmt)
```

```{r}
# Models with female covariate

lmff = lm(trust ~ treat + female * treat, data = dff)
lml = lm(trust ~ treat + female * treat, data = dl)
lmmt = lm(trust ~ treat + female * treat, data = dmt)
lmall = lm(trust ~ treat + female * treat, data = dall)
```
```{r}
# Results from the regression with all subjects
summary(lmall)
ate = lmall$coefficients['treat']
paste0('The ATE of the dissonace video is: ', ate)
```
** The ATE is -0.22262 (0.07585) with a p value of 0.00364 ****

```{r star, results = 'asis', warning=FALSE, message=FALSE}
# Stargazer

library(stargazer)

stargazer(lmff, lml, lmmt, lmall, title="Experimental Results", column.labels=c("FriendsFamily","Lucid","MTurk","AllSources"))
```
```{r}
# Randomization inference among all the subjects
# Not worrying about M & F distributions in this one
dall$treat_rnd = sample(dall$treat)
head(dall)
summary(dall)
```
```{r}
# Generic function to randomly(pseudo) pick the treatment assignments
coeff_treat_rnd <- function(d) {
  d$treat_rnd = sample(d$treat)
  lmr = lm(trust ~ treat_rnd + female * treat, data=d, na.action = na.omit)
  return (lmr$coefficients['treat_rnd'])
}

# Trying another randomization
coeff_treat_rnd(dall)
```
```{r}
# 100K Randomizations
h.distribution.under.sharp.null <- replicate(100000, coeff_treat_rnd(dall))
h.ate.mean <- mean(h.distribution.under.sharp.null)
paste0("Mean fromm 1k randomizations: ", h.ate.mean)

```
```{r}
# Graph for the estimates
plot(density(h.distribution.under.sharp.null), 
     main = "Density under Sharp Null")

# Adding our original ATE to the plot
abline(v = ate, col = "blue")

# Adding the mean to the plot
abline(v = mean(h.distribution.under.sharp.null), col = 'orange')
```
```{r}
# Values close to what we get
pn = sum(abs(h.distribution.under.sharp.null) >= abs(ate))
paste0('Values bigger than the ATE: ', pn)

# 2 tailed p value
pv = mean(abs(h.distribution.under.sharp.null) >= abs(ate))
paste0('p value: ', pv)
```
```{r}
# RI by randomizing within male and female groups
dallf = dall[which(dall$female == 1), ]
dallf$treat_rnd = sample(dallf$treat)
```
```{r}
# RI by randomizing within male and female groups
dallm = dall[which(dall$female == 0), ]
dallm$treat_rnd = sample(dallm$treat)
```
```{r}
dall2 = rbind(dallf, dallm)

```

```{r}
# Generic function to randomly(pseudo) pick the treatment assignments
# This randomizes within the M and F groupos
coeff_treat_rnd_s <- function(d) {
  
  df = d[which(d$female == 1), ]
  df$treat_rnd = sample(df$treat)
  
  dm = d[which(d$female == 0), ]
  dm$treat_rnd = sample(dm$treat)
  
  da = rbind(df, dm)
  
  lms = lm(trust ~ treat_rnd + female * treat, data=da, na.action = na.omit)
  return (lms$coefficients['treat_rnd'])
}

# Trying another randomization
coeff_treat_rnd_s(dall2)
```

```{r}
# 50K Randomizations
h.distribution.under.sharp.null.2 <- replicate(50000, coeff_treat_rnd_s(dall2))
h.ate.mean.2 <- mean(h.distribution.under.sharp.null.2)
paste0("Mean fromm 1k randomizations: ", h.ate.mean.2)

```
```{r}
# Graph for the estimates
plot(density(h.distribution.under.sharp.null.2), 
     main = "Density under Sharp Null")

# Adding our original ATE to the plot
abline(v = ate, col = "blue")

# Adding the mean to the plot
abline(v = mean(h.distribution.under.sharp.null.2), col = 'orange')
```
```{r}
# Values close to what we get
pn2 = sum(abs(h.distribution.under.sharp.null.2) >= abs(ate))
paste0('Number of observations bigger than the experiment: ', pn2)

# 2 tailed p value
pv2 = mean(abs(h.distribution.under.sharp.null.2) >= abs(ate))
paste0('p value: ', pv2)
```