#-------------------------------------------------------------------------------
# R script
# Study 3 - Data analysis
# Goswami and Urminsky "When should the ask be a nudge?" 
# Journal of Marketing Research 2016
#
#-------------------------------------------------------------------------------

#R library needed for lasso variable selection
library(glmnet)
#remove other variables from memory
rm(list=ls())  

##############################################################################
#data reading and top-line numbers
#check VariableMeaning.xls for variable definition

#read the data file as a .csv
dataOW <- read.csv("Study3Data-Pub.csv", header =T)
#number of observations
nrow(dataOW)
#overall donation rate
mean(dataOW$donated)*100 



##############################################################################
#data preparation


#convert dollar amounts to log scale
dataOW$LnDonAmt <- log(dataOW$donAmt + 1)
dataOW$LnSchoolLifetimeGiving <- log(dataOW$SchoolLifetimeGivingZ + 1)
dataOW$LnAnnualFundLifetimeGiving <- log(dataOW$AnnualFundLifetimeGivingZ + 1)
dataOW$LnMenu2Num <- log(dataOW$Menu2Num)

#convert categorical variables into factors in R
dataOW$defLvl <- factor(dataOW$defLvl, c("NO", "LO", "MD", "HI"))
dataOW$Segment <- factor(dataOW$Segment, c( "Other", "Retention", "1-2 SYBUNT"))
dataOW$numOpt <- factor(dataOW$numOpt, c("2", "5"))

#recode default-levels for comparing control vs high plus medium
dataOW$defLvl2[dataOW$defLvl == "NO"] <- "NO"
dataOW$defLvl2[dataOW$defLvl == "LO"] <- "LO"
dataOW$defLvl2[dataOW$defLvl %in% c("MD", "HI")] <- "HI"
dataOW$defLvl2 <- factor(dataOW$defLvl2, c("NO", "LO", "HI"))

#recode default-levels for comparing control vs high plus medium (ignore control)
dataOW$defLvl6[dataOW$defLvl == "NO"] <- NA
dataOW$defLvl6[dataOW$defLvl == "LO"] <- "LO"
dataOW$defLvl6[dataOW$defLvl %in% c("MD", "HI")] <- "HI"
dataOW$defLvl6 <- factor(dataOW$defLvl6, c("LO", "HI"))


#dummy code for the Regention segment 
dataOW$LYBUNT <- ifelse ( dataOW$Segment == "Retention", 1,
                          (ifelse (dataOW$Segment %in% 
                                     c( "1-2 SYBUNT", "Other"), 0, NA)))
dataOW$LYBUNT <- factor(dataOW$LYBUNT, c(0,1))


################################################################################
#Table 2


#Table 2, column 2
table(dataOW$defLvl)

#Table 2, column 3
table(dataOW$donated, dataOW$defLvl)

#Table 2, column 4 manually computed from column 2 and column 3

#Table 2, column 5 
with(dataOW[dataOW$donated==1,], tapply(donAmt, list(defLvl), mean))

#Table 2, column 6 
with(dataOW, tapply(donAmt, list(defLvl), mean))



################################################################################
#Table 3


### Regression Estimates: Log of Revenue per Person ###


#1. preping data for lasso
d1 <- dataOW[c("LnDonAmt", 
               "AgeZ_adj", 
               "YrsOfAssosZ_adj", 
               "AFConsYoGZ", 
               "LnAnnualFundLifetimeGiving",
               "LnSchoolLifetimeGiving",  
               "LnMenu2Num", 
               "defLvl" ,
               "numOpt",  
               "prevRemind", 
               "Segment", 
               "Gender")]
d1 <- na.omit(d1)


#2. preping columns, categorical and then continuous
xfactors <- model.matrix(LnDonAmt ~ defLvl + 
                           prevRemind + 
                           numOpt + 
                           Segment + 
                           Gender, 
                         data=d1)[,-1]
x <- with(d1, as.matrix(data.frame(AgeZ_adj, YrsOfAssosZ_adj, AFConsYoGZ, 
                                   LnAnnualFundLifetimeGiving,  
                                   LnSchoolLifetimeGiving, 
                                   LnMenu2Num, xfactors)))
y <- d1$LnDonAmt


#3. run lasso

nvars <- 14 # nunber of RHS variables excluding intercept
nobs <- nrow(d1)
p.fac <- rep(1, nvars)
p.fac[c(7, 8, 9, 10, 11)] <- 0 #retain experimental factor

lasso <- glmnet(x,y=y, penalty.factor = p.fac)


#4. cross validation

cv.glmmod <- cv.glmnet(x,y=y,penalty.factor = p.fac)
cv.glmmod$lambda.1se; log(cv.glmmod$lambda.1se)
coef(cv.glmmod, s=(cv.glmmod$lambda.1se))  


#5. covariate selection based on cross-validation results

#lasso model
m1 <-lm(formula = LnDonAmt ~ defLvl + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)



### Regression Estimates: Donation date  ###
### similar steps as Log of Revenue per Person 


d1 <- dataOW[c("donated", "AgeZ_adj", "YrsOfAssosZ_adj", 
               "AFConsYoGZ", "LnAnnualFundLifetimeGiving",
               "LnSchoolLifetimeGiving",  "LnMenu2Num", 
               "defLvl" , "numOpt",  "prevRemind", 
               "Segment", "Gender")]
d1 <- na.omit(d1)
d1$donated <- factor(d1$donated, c(0,1))


xfactors <- model.matrix(donated ~ defLvl + prevRemind + numOpt +
                           Segment + Gender, data=d1)[,-1]
x <- with(d1, as.matrix(data.frame(AgeZ_adj, YrsOfAssosZ_adj,
                            AFConsYoGZ, LnAnnualFundLifetimeGiving,  
                            LnSchoolLifetimeGiving, LnMenu2Num, xfactors)))
y <- d1$donated

p.fac <- rep(1, 14)
p.fac[c(7, 8, 9, 10, 11)] <- 0
lasso <- glmnet(x, y=y, penalty.factor = p.fac, family="binomial");


cv.glmmod <- cv.glmnet(x, y=y, penalty.factor = p.fac, family="binomial")
coef(cv.glmmod, s=(cv.glmmod$lambda.1se)) 


m1 <- glm(donated ~ defLvl + numOpt + prevRemind +
            AgeZ_adj + AFConsYoGZ + LYBUNT, data=dataOW, family = "binomial" )
summary(m1)



### Regression Estimates: Log of AverAgeZ_adj Donation  ###
### similar steps as Log of Revenue per Person/Donation rate



d1 <- dataOW[dataOW$donated ==1, ][c("LnDonAmt", "AgeZ_adj", "YrsOfAssosZ_adj", 
                                     "AFConsYoGZ", "LnAnnualFundLifetimeGiving",
                                     "LnSchoolLifetimeGiving",  "LnMenu2Num", 
                                     "defLvl" , "numOpt",  "prevRemind", 
                                     "Segment", "Gender")]
d1 <- na.omit(d1)

xfactors <- model.matrix(LnDonAmt ~ defLvl + prevRemind + numOpt +
                           Segment + Gender, data=d1)[,-1]
x <- with(d1, as.matrix(data.frame(AgeZ_adj, YrsOfAssosZ_adj, AFConsYoGZ, 
                                   LnAnnualFundLifetimeGiving,  
                                   LnSchoolLifetimeGiving, 
                                   LnMenu2Num, xfactors)))
y=d1$LnDonAmt


p.fac <- rep(1, 14)
p.fac[c(7, 8, 9, 10, 11)] <- 0
lasso<-glmnet(x,y=y, penalty.factor = p.fac);

cv.glmmod <- cv.glmnet(x,y=y,penalty.factor = p.fac)
coef(cv.glmmod, s=(cv.glmmod$lambda.1se))  

m1 <-lm(formula = LnDonAmt ~ defLvl + numOpt + prevRemind
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)




################################################################################
# stats in text, in sequential order of their appearance
# some stats are taken from Table 2 and Table 3


##############################
## Choice of default option ##


#identifying those selecting low, medium, high default amounts 
#on their own in control condition

dataOW$ControlLO <- 0
dataOW$ControlLO[dataOW$cond == 7 & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$ControlLO[dataOW$cond == 8 & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$ControlLO[dataOW$cond == 15 & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$ControlLO[dataOW$cond == 16 & dataOW$donAmt == dataOW$Menu3Num] <- 1

dataOW$ControlMD <- 0
dataOW$ControlMD[dataOW$cond == 7 & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$ControlMD[dataOW$cond == 8 & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$ControlMD[dataOW$cond == 15 & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$ControlMD[dataOW$cond == 16 & dataOW$donAmt == dataOW$Menu2Num] <- 1


dataOW$ControlHI <- 0
dataOW$ControlHI[dataOW$cond == 7 & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$ControlHI[dataOW$cond == 8 & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$ControlHI[dataOW$cond == 15 & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$ControlHI[dataOW$cond == 16 & dataOW$donAmt == dataOW$Menu1Num] <- 1


dataOW$ControlDefSel <- ifelse(dataOW$ControlLO ==1, "LO",
                               (ifelse (dataOW$ControlMD ==1, "MD",
                               (ifelse(dataOW$ControlHI ==1, "HI", 
                               (ifelse(dataOW$donAmt != 0 & 
                                         dataOW$cond %in% c( 7, 8, 15, 16 ) &
                                         dataOW$ControlLO ==0 &
                                         dataOW$ControlMD ==0 & 
                                         dataOW$ControlHI ==0, "Other", NA)))))))

table(dataOW$ControlDefSel)
prop.table(table(dataOW$ControlDefSel))

#identifying those selecting low, medium, high default amounts 
#in the LOW defaults condition

dataOW$LowLO <- 0
dataOW$LowLO[dataOW$defLvl == "LO" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$LowLO[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$LowLO[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$LowLO[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu3Num] <- 1

dataOW$LowMD <- 0
dataOW$LowMD[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$LowMD[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$LowMD[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$LowMD[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu2Num] <- 1


dataOW$LowHI <- 0
dataOW$LowHI[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$LowHI[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$LowHI[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$LowHI[dataOW$defLvl ==  "LO" & dataOW$donAmt == dataOW$Menu1Num] <- 1


dataOW$LowDefSel <- ifelse(dataOW$LowLO ==1, "LO",
                               (ifelse (dataOW$LowMD ==1, "MD",
                                        (ifelse(dataOW$LowHI ==1, "HI", 
                                        (ifelse(dataOW$donAmt != 0 & 
                                          dataOW$defLvl ==  "LO" &
                                          dataOW$LowLO ==0 &
                                          dataOW$LowMD ==0 & 
                                          dataOW$LowHI ==0, "Other", NA)))))))

table(dataOW$LowDefSel)
prop.table(table(dataOW$LowDefSel))


#identifying those selecting low, medium, high default amounts 
#in the MEDIUM defaults condition

dataOW$MedLO <- 0
dataOW$MedLO[dataOW$defLvl == "MD" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$MedLO[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$MedLO[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$MedLO[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu3Num] <- 1

dataOW$MedMD <- 0
dataOW$MedMD[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$MedMD[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$MedMD[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$MedMD[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu2Num] <- 1


dataOW$MedHI <- 0
dataOW$MedHI[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$MedHI[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$MedHI[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$MedHI[dataOW$defLvl ==  "MD" & dataOW$donAmt == dataOW$Menu1Num] <- 1


dataOW$MedDefSel <- ifelse(dataOW$MedLO ==1, "LO",
                           (ifelse (dataOW$MedMD ==1, "MD",
                                    (ifelse(dataOW$MedHI ==1, "HI", 
                                    (ifelse(dataOW$donAmt != 0 & 
                                    dataOW$defLvl ==  "MD" &
                                    dataOW$MedLO ==0 &
                                    dataOW$MedMD ==0 & 
                                    dataOW$MedHI ==0, "Other", NA)))))))

table(dataOW$MedDefSel)
prop.table(table(dataOW$MedDefSel))

#identifying those selecting low, medium, high default amounts 
#in the HIGH defaults condition

dataOW$HighLO <- 0
dataOW$HighLO[dataOW$defLvl == "HI" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$HighLO[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$HighLO[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu3Num] <- 1
dataOW$HighLO[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu3Num] <- 1

dataOW$HighMD <- 0
dataOW$HighMD[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$HighMD[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$HighMD[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu2Num] <- 1
dataOW$HighMD[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu2Num] <- 1


dataOW$HighHI <- 0
dataOW$HighHI[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$HighHI[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$HighHI[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu1Num] <- 1
dataOW$HighHI[dataOW$defLvl ==  "HI" & dataOW$donAmt == dataOW$Menu1Num] <- 1


dataOW$HighDefSel <- ifelse(dataOW$HighLO ==1, "LO",
                           (ifelse (dataOW$HighMD ==1, "MD",
                                    (ifelse(dataOW$HighHI ==1, "HI", 
                                    (ifelse(dataOW$donAmt != 0 & 
                                    dataOW$defLvl ==  "HI" &
                                    dataOW$HighLO ==0 &
                                    dataOW$HighMD ==0 & 
                                    dataOW$HighHI ==0, "Other", NA)))))))

table(dataOW$HighDefSel)
prop.table(table(dataOW$HighDefSel))



################################################################
## Effect of defaults on revenue per person (i.e. net effect) ##


with(dataOW, tapply(donAmt, list(defLvl), mean, na.rm=T))


#covariates based on losso for predicting LnDonAmt

#Table A81
#def present vs not
m1 <-lm(formula = LnDonAmt ~ defPres + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)


#Table 3
#low/medium or high default vs no-default (3 regression coefs)
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "NO")
m1 <-lm(formula = LnDonAmt ~ defLvl + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)


#Table A82
#high/medium defaut vs low default
dataOW$defLvl2 <- relevel(dataOW$defLvl2, ref = "LO")
m1 <-lm(formula = LnDonAmt ~ defLvl2 + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)


#interactions

#Table A83
#AgeZ_adj
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "NO")
m1 <-lm(formula = LnDonAmt ~ defLvl*AgeZ_adj + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)


#Table A84
#years of consecutive giving
m1 <-lm(formula = LnDonAmt ~ defPres*AFConsYoGZ + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)


#Table A85
#donated in the prior campaign
m1 <-lm(formula = LnDonAmt ~ defPres*LYBUNT + numOpt + prevRemind + AFConsYoGZ, 
        data = dataOW)
summary(m1)



###############################################################
## Effect of defaults on participation (i.e. donation rate) ##


#low vs no default
with(dataOW, tapply(donated, list(defLvl), mean, na.rm=T))*100
dtmp <- dataOW[dataOW$defLvl %in% c("NO", "LO"),]
dtmp$defLvl <- factor(dtmp$defLvl, c("NO", "LO"))
with(dtmp, prop.test(table(defLvl, donated), correct=FALSE)) 


#low vs high/med combined
with(dataOW, tapply(donated, list(defLvl2), mean, na.rm=T))*100
dtmp <- dataOW[dataOW$defLvl2 %in% c("HI", "LO"),]
dtmp$defLvl2 <- factor(dtmp$defLvl2, c("HI", "LO"))
with(dtmp, prop.test(table(defLvl2, donated), correct=FALSE))


#donation rates, default present vs absent
with(dataOW, tapply(donated, list(defPres), mean, na.rm=T))*100
chisq.test(table(dataOW$defPres, dataOW$donated), correct=F)


#covariates based on losso for donation rate (i.e. donated)

#Table A86
#def present vs not
m1 <- glm(donated ~ defPres + numOpt + prevRemind +
            AgeZ_adj + AFConsYoGZ + LYBUNT, data=dataOW, family = "binomial" )
summary(m1)


#low vs control, 3 regression coefs
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "NO")
m1 <- glm(donated ~ defLvl + numOpt + prevRemind +
            AgeZ_adj + AFConsYoGZ + LYBUNT, data=dataOW, family = "binomial" )
summary(m1)


#low vs combining HI and MED
dataOW$defLvl2 <- relevel(dataOW$defLvl2, ref = "LO")
m1 <- glm(donated ~ defLvl2 + numOpt + prevRemind +
            AgeZ_adj + AFConsYoGZ + LYBUNT, data=dataOW, family = "binomial" )
summary(m1)



#########################################################################
## Effect of defaults on conditional contribution (i.e. avg. donation) ##


with(dataOW[dataOW$donated==1,], tapply(donAmt, list(defPres), mean))

# some descriptives based on table 2

#covariates based on losso for predicting LnDonAmt (for participant only)

#low vs no defaults, from Table 3
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "NO")
m1 <-lm(formula = LnDonAmt ~ defLvl + numOpt + prevRemind
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)


#Table A87
#same comparison with low defaults as reference
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "LO")
m1 <-lm(formula = LnDonAmt ~ defLvl + numOpt + prevRemind
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)


#Table A88
#presence of defaults
m1 <-lm(formula = LnDonAmt ~ defPres + numOpt + prevRemind
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)


#effect of experimental factors

#Table A89
#effect of reminders * default levels (last coef of HI used)
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "LO")
m1 <-lm(formula = LnDonAmt ~ defLvl*prevRemind + numOpt
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)


#Table A90
#effect of reminders * default levels, different coding (not used in main paper)
dataOW$defLvl2 <- relevel(dataOW$defLvl2, ref = "LO")
m1 <-lm(formula = LnDonAmt ~ defLvl2*prevRemind + numOpt
        + LnAnnualFundLifetimeGiving + LnMenu2Num,  data = dataOW[dataOW$donated==1,])
summary(m1)



#Table A91
#middle option interaction on conditional contribution
dataOW$defLvl <- relevel(dataOW$defLvl, ref = "NO")
m1 <-lm(formula = LnDonAmt ~ defLvl*LnMenu2Num + numOpt + prevRemind
        + LnAnnualFundLifetimeGiving ,  data = dataOW[dataOW$donated==1,])
summary(m1)


#Table A92 (not used in main paper)
m1 <-lm(formula = LnDonAmt ~ defPres*LnSchoolLifetimeGiving +
          LnAnnualFundLifetimeGiving + LnMenu2Num + 
          numOpt + prevRemind ,  data = dataOW[dataOW$donated==1,])
summary(m1)




