####replication code
library(haven)
library(dplyr)
library(tidyr)
library(ggplot2)
library(nnet)
library(knitr)
library(broom)


#1. Import wave 4 data 
fp <-"C://TempTY/Manuscripts/NSF_Panel_Proposal/Analysis/FinalDataSets_for_ICPSR/"

wave4<-read_sas(paste0(fp, "wave4_final.sas7bdat"))

#2. Create new variables for dates

wave4a<-wave4 %>% 
  mutate(day_start = as.Date(StartDate)) %>% 
  mutate(day_complete = as.Date(CompleteDate))

#3. Import wave 1 clean data and subset demographic variables

wave1<-read_sas(paste0(fp, "wave1_final.sas7bdat"))

demographics<-wave1 %>% 
  select(SampleID,DEMO_1,DEMO_2,DEMO_6,DEMO_7,DEMO_8_1,DEMO_8_2,DEMO_8_3,DEMO_8_4,DEMO_8_5,DEMO_8_6,DEMO_14, 
         wave4_contact_grp, wave4_invited) %>% 
  mutate(MALE = if_else(DEMO_2 == 1,1,0)) %>% 
  mutate(AGE = if_else(DEMO_1<25,1,
                       if_else(DEMO_1<35,2,
                               if_else(DEMO_1<50,3,4)))) %>% 
  mutate(AGE_rec = if_else(DEMO_1<50,1,2)) %>% 
  mutate(HISPANIC = if_else(DEMO_7 == 1,0,1)) %>% 
  mutate(MULTIRACE = 0) %>%
  mutate(MULTIRACE = if_else(DEMO_8_1 == "1",MULTIRACE+1,MULTIRACE)) %>% 
  mutate(MULTIRACE = if_else(DEMO_8_2 == "1",MULTIRACE+1,MULTIRACE)) %>%
  mutate(MULTIRACE = if_else(DEMO_8_3 == "1",MULTIRACE+1,MULTIRACE)) %>%
  mutate(MULTIRACE = if_else(DEMO_8_4 == "1",MULTIRACE+1,MULTIRACE)) %>%
  mutate(MULTIRACE = if_else(DEMO_8_5 == "1",MULTIRACE+1,MULTIRACE)) %>%
  mutate(MULTIRACE = if_else(DEMO_8_6 == "1",MULTIRACE+1,MULTIRACE)) %>%
  mutate(RACE = if_else(HISPANIC == 1,1,0)) %>% 
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_1 == 1,2,RACE)) %>% #NH-WHITE
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_2 == 1,3,RACE)) %>% #NH-BLACK
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_3 == 1,4,RACE)) %>% #NH-OTHER
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_4 == 1,4,RACE)) %>% #NH-OTHER
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_5 == 1,4,RACE)) %>% #NH-OTHER
  mutate(RACE = if_else(MULTIRACE == 1 & HISPANIC == 0 & DEMO_8_6 == 1,4,RACE)) %>% #NH-OTHER
  mutate(RACE = if_else(MULTIRACE != 1 & HISPANIC == 0,4,RACE)) %>% #NH-OTHER
  mutate(EDUCATION = if_else(DEMO_6 == 6,2,1)) %>%
  mutate(INCOME = if_else(DEMO_14 == 8,2,1)) %>%
  mutate(sampleid=SampleID)

#4. Import Wave 2, Wave 3 data
wave2<-read_sas(paste0(fp,"wave2_final.sas7bdat")) %>% select(sampleid, wave2_complete)

wave3<-read_sas(paste0(fp,"wave3_final.sas7bdat")) %>% select(sampleid, wave3_complete)


#5 Merge wave 4 full with demographics

wave4_sample <- wave4a %>% 
  right_join(demographics, by = c("sampleid")) %>% filter(wave4_invited == 1)


wave4_full<-wave4a %>% 
  left_join(demographics, by = c("sampleid")) %>%
  left_join(wave2, by = c("sampleid")) %>%
  left_join(wave3, by = c("sampleid")) %>%
  mutate(COMPLETE3 = if_else(wave2_complete + wave3_complete == 2,1,0)) 

wave4_full$COMPLETE3[is.na(wave4_full$COMPLETE3)] <- 0

#5. Analysis-Table 1

wave4_sample %>% 
  count(wave4_contact_grp.y)

wave4_full %>% 
  count(wave4_contact_grp.y)

# Logistic regression

log_model<-wave4_sample %>%
  filter(!is.na(wave4_contact_grp.y)) %>% 
  mutate(second = if_else(wave4_contact_grp.y == 3 | wave4_contact_grp.y == 4,1,0)) %>% 
  mutate(amazon = if_else(wave4_contact_grp.y == 1 | wave4_contact_grp.y == 3,1,0)) 

log_model$wave4_complete[is.na(log_model$wave4_complete)] <- 0

model<-glm(wave4_complete ~ second + amazon + second*amazon, data = log_model, family = "binomial")

summary(model)
exp(coef(model))

# Logistic regression # 1 day before sending 2nd incentive

log_model2<-log_model %>% 
  select(sampleid,wave4_contact_grp.y,wave4_complete, day_start, day_complete, second, amazon) %>% 
  mutate(daystocomplete = as.numeric(day_complete)) %>%
  mutate(start = as.Date("2023-02-13", origin = "2023-02-13")) %>% 
  mutate(start = as.numeric(start)) %>% 
  mutate(daystocomplete = abs(start - daystocomplete)) %>% 
  mutate(complete_2 = if_else(daystocomplete<14,1,0)) %>% 
  mutate(complete_2 = if_else(is.na(complete_2),0,complete_2))

model2<-glm(complete_2 ~ second + amazon + second*amazon, data = log_model2, family = "binomial")

summary(model2)
exp(coef(model2))

###Table 2
#8.2. Table 3.c. Subgroups and chi square test

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(MALE)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(AGE)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(AGE_rec)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(RACE)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(COMPLETE3)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(EDUCATION)

wave4_full %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(INCOME)


chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$MALE, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$AGE, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$AGE_rec, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$RACE, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$COMPLETE3, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$EDUCATION, correct = FALSE)
chisq.test(wave4_full$wave4_contact_grp.y, wave4_full$INCOME, correct = FALSE)


##Figure 2

cumulativeRR<-wave4_full %>% 
  mutate(daystocomplete = as.numeric(day_complete)) %>%
  mutate(start = as.Date("2023-02-13", origin = "2023-02-13")) %>% 
  mutate(start = as.numeric(start)) %>% 
  mutate(daystocomplete = abs(start - daystocomplete)) %>% 
  group_by(wave4_contact_grp.y) %>% 
  count(daystocomplete)

cumulativeRR2<-cumulativeRR %>% 
  mutate(RR = if_else(wave4_contact_grp.y==1,n/332,
                      if_else(wave4_contact_grp.y==2,n/332,
                              if_else(wave4_contact_grp.y==3,n/332,n/331))))

cumulativeRR3<-cumulativeRR2 %>%
  mutate(cumulativeRR = if_else(wave4_contact_grp.y==1, cumsum(RR),
                                if_else(wave4_contact_grp.y==2,cumsum(RR),
                                        if_else(wave4_contact_grp.y==3,cumsum(RR),cumsum(RR)))))


cumulativeRR3 %>% 
  mutate(daystocomplete = as.numeric(daystocomplete)) %>% 
  mutate(cumulativeRR = as.numeric(cumulativeRR)) %>% 
  mutate(Contact_Group = as.factor(wave4_contact_grp.y)) %>% 
  ggplot(aes(x = daystocomplete, y = cumulativeRR, colour = Contact_Group, linetype = Contact_Group, shape = Contact_Group))+
  geom_line(size = 1.5)+
  geom_point(size = 5)+
  scale_linetype_manual(values =c("solid","dashed","dotted","dotted", "dotdash"),breaks=c("1", "2", "3", "4"), labels=c("Amazon ($5/$0)","Cash ($5/$0)","Amazon ($2/$3)","Cash ($2/$3)"))+
  scale_colour_manual(values=c("#FF3300","#66cc00","#000FFF","#CC9900"),breaks=c("1", "2", "3", "4"), labels=c("Amazon ($5/$0)","Cash ($5/$0)","Amazon ($2/$3)","Cash ($2/$3)"))+
  scale_shape_manual(values=c(19,15,19,15),breaks=c("1", "2", "3", "4"), labels=c("Amazon ($5/$0)","Cash ($5/$0)","Amazon ($2/$3)","Cash ($2/$3)"))+
  scale_x_continuous(name="Days")+
  scale_y_continuous(name="Cumulative RR")+
  theme(legend.position="bottom")+
  theme(legend.title = element_text(size=18), #change legend title font size,
        legend.text = element_text(size=23))+ #change legend text font size
  theme(axis.text=element_text(size=20),
        axis.title=element_text(size=24,face="bold"))+
  theme(axis.text.x = element_text(face="bold", 
                                   size=24),
        axis.text.y = element_text(face="bold", 
                                   size=24))+
  #scale_color_discrete("Contact Group", labels=c("Amazon ($5/$0)","Cash ($5/$0)","Amazon ($2/$3)","Cash ($2/$3)"))+
  geom_vline(xintercept = 7, linetype="dashed", 
             color = "black", size=0.6)+
  geom_vline(xintercept = 10, linetype="dashed", 
             color = "black", size=0.6)+
  geom_vline(xintercept = 14, linetype="dashed", 
             color = "black", size=0.6)+
  geom_vline(xintercept = 21, linetype="dashed", 
             color = "black", size=0.6)+
  geom_vline(xintercept = 28, linetype="dashed", 
             color = "black", size=0.6)+
  annotate("text", x=6.3, y=0.16, label="Email/Text reminder", angle=90, size =6.5)+
  annotate("text", x=9.3, y=0.14, label="Email reminder", angle=90, size =6.5)+
  annotate("text", x=13.3, y=0.17, label="Letter (2nd incentive)", angle=90, size =6.5)+
  annotate("text", x=20.3, y=0.17, label="Email/Text reminder*", angle=90, size =6.5)+
  annotate("text", x=27.3, y=0.15, label="Email reminder**", angle=90, size =6.5)


##Nonresponse Bias

NR_BIAS <- log_model %>% select(SampleID,wave4_complete) %>%
  left_join(wave1, by = c("SampleID")) 

NR_BIAS_v2<-NR_BIAS %>% 
  mutate(EXP2_1 = if_else(EXP2_1 == 1,1,0)) %>% 
  mutate(EXP2_2 = if_else(EXP2_2 == 1,1,0)) %>% 
  mutate(EXP2_3 = if_else(EXP2_3 == 1,1,0)) %>% 
  mutate(EXP2_4 = if_else(EXP2_4 == 1,1,0)) %>% 
  mutate(EXP2_5 = if_else(EXP2_5 == 1,1,0)) %>% 
  mutate(EXP2b_1 = if_else(EXP2b_1 == 1,1,0)) %>% 
  mutate(EXP2b_2 = if_else(EXP2b_2 == 1,1,0)) %>% 
  mutate(EXP2b_3 = if_else(EXP2b_3 == 1,1,0)) %>% 
  mutate(EXP2b_4 = if_else(EXP2b_4 == 1,1,0)) %>% 
  mutate(EXP1_1a = if_else(EXP1_1a == 1,1,0)) %>% 
  mutate(EXP1_1b = if_else(EXP1_1b == 1,1,0)) %>% 
  mutate(EXP1_1c = if_else(EXP1_1c == 1,1,0)) %>% 
  mutate(EXP1_1d = if_else(EXP1_1d == 1,1,0)) %>% 
  mutate(EXP1_1e = if_else(EXP1_1e == 1,1,0)) %>% 
  mutate(EXP1_1f = if_else(EXP1_1f == 1,1,0)) %>% 
  mutate(EXP1_1g = if_else(EXP1_1g == 1,1,0)) %>% 
  mutate(EXP1_2a = if_else(EXP1_2a == 1,1,0)) %>% 
  mutate(EXP1_2b = if_else(EXP1_2b == 1,1,0)) %>% 
  mutate(EXP1_3a = if_else(EXP1_3a == 1,1,0)) %>% 
  mutate(EXP1_3b = if_else(EXP1_3b == 1,1,0)) %>% 
  mutate(EXP1_3c = if_else(EXP1_3c == 1,1,0)) %>% 
  mutate(EXP1_3d = if_else(EXP1_3d == 1,1,0)) %>% 
  mutate(EXP1_3e = if_else(EXP1_3e == 1,1,0)) %>% 
  mutate(EXP1_3f = if_else(EXP1_3f == 1,1,0)) 
 

NR_BIAS_v2 %>%
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP2_1),mean(EXP2_2),mean(EXP2_3),mean(EXP2_4),mean(EXP2_5))

NR_BIAS_v2 %>%
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP2b_1),mean(EXP2b_2),mean(EXP2b_3),mean(EXP2b_4))

NR_BIAS_v2 %>%
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_1a),mean(EXP1_1b),mean(EXP1_1c),mean(EXP1_1d),mean(EXP1_1e),mean(EXP1_1f),mean(EXP1_1g))

NR_BIAS_v2 %>%
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_2a),mean(EXP1_2b))

NR_BIAS_v2 %>%
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_3a),mean(EXP1_3b),mean(EXP1_3c),mean(EXP1_3d),mean(EXP1_3e),mean(EXP1_3f))

NR_BIAS_v2 %>%
  filter(wave4_complete==1) %>% 
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP2_1),mean(EXP2_2),mean(EXP2_3),mean(EXP2_4),mean(EXP2_5))

NR_BIAS_v2 %>%
  filter(wave4_complete==1) %>% 
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP2b_1),mean(EXP2b_2),mean(EXP2b_3),mean(EXP2b_4))

NR_BIAS_v2 %>%
  filter(wave4_complete==1) %>% 
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_1a),mean(EXP1_1b),mean(EXP1_1c),mean(EXP1_1d),mean(EXP1_1e),mean(EXP1_1f),mean(EXP1_1g))

NR_BIAS_v2 %>%
  filter(wave4_complete==1) %>% 
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_2a),mean(EXP1_2b))

NR_BIAS_v2 %>%
  filter(wave4_complete==1) %>% 
  group_by(wave4_contact_grp) %>% 
  summarise(mean(EXP1_3a),mean(EXP1_3b),mean(EXP1_3c),mean(EXP1_3d),mean(EXP1_3e),mean(EXP1_3f))