installing necessary libraries

library(dplyr) 
library(texreg)
library(lme4)
library(rmarkdown)
library(ggplot2)
library(GGally)
library(effects)
library(visreg)

reading in and transforming data

all <- read.csv("~/Documents/Dropbox/Research/Eric/Study 2a/2a v2/all conds with beta delta.csv") %>% #reading in
        mutate(cond_label = factor(ifelse(cond == 1, "no_default", #making a condition description column
                           ifelse(cond == 3, "SS_default",
                           ifelse(cond == 2, "LL_default", NA)))),
               DV_label = factor(ifelse(DV == 0, "SS",
                          ifelse(DV == 1, "LL", NA))), 
                  cond = factor(cond))%>% # making condition a factor 
        select(-coding) #dropping coding column
        all<-subset(all, cond!=1) #dropping control condition

centering continuous measures

all$betac <- all$beta - mean(all$beta)
all$deltac<- all$delta - mean(all$delta)
all$betacR<-(all$betac) * (-1)
all$combined <- all$betacR + all$deltac

all$adr<- exp(all$deltac * 365)
all$adrc<-all$adr - mean(all$adr)
all$adrc<- -1 * (all$adrc)

the holy grail continuous analysis (combined = beta + delta)

fit3<-glm(DV ~ combined*cond_label + I(combined^2)*cond_label, data=all, family=binomial)
visreg(fit3, "combined", by="cond_label", overlay=TRUE, partial=FALSE,  xlab="patience", scale="response", ylab="P (choosing SS)")

this works using just beta too

fit4<-glm(DV ~ betac*cond_label + I(betac^2)*cond_label, data=all, family=binomial)
visreg(fit4, "betac", by="cond_label", overlay=TRUE, partial=FALSE,  xlab="beta", scale="response", ylab="P (choosing SS)")

and just delta!

fit5<-glm(DV ~ deltac*cond_label + I(deltac^2)*cond_label, data=all, family=binomial)
visreg(fit5, "deltac", by="cond_label", overlay=TRUE, partial=FALSE,  xlab="delta", scale="response", ylab="P (choosing SS)")

this is the ‘discrete-ized’ version of the above continuous analysis - shows that there is only a significalt default effect for “middle” patient people, using the cut-offs defined by the floodlight analysis (roughly +/-1 sd from mean)

all$cond<-droplevels(all$cond)
patient<-subset(all, combined<=-0.1)
middle<-subset(all, combined > -0.1)
middle<-subset(middle, combined <0.2)
impatient<-subset(all, combined >= 0.2)
prop.test(table(patient$cond, patient$DV)) #no default
## Warning in prop.test(table(patient$cond, patient$DV)): Chi-squared
## approximation may be incorrect
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(patient$cond, patient$DV)
## X-squared = 0.2024, df = 1, p-value = 0.6528
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1943869  0.4296810
## sample estimates:
##    prop 1    prop 2 
## 0.2352941 0.1176471
prop.test(table(impatient$cond, impatient$DV)) #no default
## Warning in prop.test(table(impatient$cond, impatient$DV)): Chi-squared
## approximation may be incorrect
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(impatient$cond, impatient$DV)
## X-squared = 0.2763, df = 1, p-value = 0.5992
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.181475  0.429841
## sample estimates:
##    prop 1    prop 2 
## 0.8888889 0.7647059
prop.test(table(middle$cond, middle$DV)) # yes default
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(middle$cond, middle$DV)
## X-squared = 9.6061, df = 1, p-value = 0.001939
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.06415276 0.28753708
## sample estimates:
##    prop 1    prop 2 
## 0.5562130 0.3803681

and here is the original plot that inspired crystal’s insight (I used SPSS to do the floodlight analysis [identifying the precise points where the effect of patience becomes significant] because I couldn’t figure it out in R)

note this uses p(choosing default) on the y-axis

nest2 <- glm(chosedefault1yes ~ combined* cond_label, data = all, family=binomial)
effect_test <- effect("combined*cond_label", nest2)       
plot(effect_test, multiline=T, ylab="Probability of Choosing DEFAULT OPTION", xlab="Patience")

and using p(choosing SS) on the y axis:

nest3 <- glm(DV ~ combined* cond_label, data = all, family=binomial)
effect_test <- effect("combined*cond_label", nest3)       
plot(effect_test, multiline=T, ylab="Probability of Choosing SS OPTION", xlab="Patience")

EVERYTHING BELOW HERE IS OLDER STUFF

the logistic regression below gives the effects of conditions 2 and 3 relative to condition 1, which I don’t think we want… so I use ANOVA instead to test for an overall effect of condition Note that DV is choice of SS/LL payment option, cond is the default condition (SS / LL / No default)

summary(glm(DV ~ cond, data=all, family="binomial"))  
## 
## Call:
## glm(formula = DV ~ cond, family = "binomial", data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3707  -1.0788   0.9957   0.9957   1.2793  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2364     0.1410  -1.676 0.093654 .  
## cond3         0.6801     0.2030   3.350 0.000807 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 555.00  on 400  degrees of freedom
## Residual deviance: 543.61  on 399  degrees of freedom
## AIC: 547.61
## 
## Number of Fisher Scoring iterations: 4
summary(aov(DV ~ cond, data=all))
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## cond          1   2.83  2.8273   11.61 0.000724 ***
## Residuals   399  97.20  0.2436                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

given the significant effect of default condition revealed by ANOVA, test for differences in the proportions of Ps choosing SS between each condition

first just look at the numbers:

all$cond_label<-droplevels(all$cond_label)
table(all$cond_label, all$DV_label) # columns: 0=SS, 1=LL. rows: 1=no default, 2=SS default, 3=LL default
##             
##               LL  SS
##   LL_default  90 114
##   SS_default 120  77
prop.test(table(all$cond_label, all$DV_label)) #proportions shown = % choosing SS
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(all$cond_label, all$DV_label)
## X-squared = 10.6712, df = 1, p-value = 0.001088
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.26930938 -0.06661179
## sample estimates:
##    prop 1    prop 2 
## 0.4411765 0.6091371

beta and delta (centered) are correlated:

##splom
ggpairs(all %>% select(betac, adr), 
        lower=list(continuous="smooth", params=c(colour="blue")),
        diag=list(continuous="density", params=c(colour="blue")),
        upper=list(params=list(corSize=5)), axisLabels='show',
        title = "beta / delta correlation")

evaluate the nested logit as per eric’s suggestion - this tell us the effect of ADR in each condition:

nest <- glm(formula = DV ~ combined %in% cond_label, family = "binomial", data = all)
summary(nest)
## 
## Call:
## glm(formula = DV ~ combined %in% cond_label, family = "binomial", 
##     data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6995  -1.1319   0.5901   1.0006   2.4668  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    0.03848    0.10931   0.352    0.725    
## combined:cond_labelLL_default -7.43951    1.60164  -4.645  3.4e-06 ***
## combined:cond_labelSS_default -6.47750    1.50142  -4.314  1.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 555.00  on 400  degrees of freedom
## Residual deviance: 495.75  on 398  degrees of freedom
## AIC: 501.75
## 
## Number of Fisher Scoring iterations: 4

we can then plot these three lines

effect_test <- effect("combined*cond_label", nest)       
 plot(effect_test, multiline=T, ylab="Probability of Choosing SS", xlab="Patience", cex=100, main="Effect of Patience on Probability of Choosing SS")

try spotlight analysis - at 1st and 3rd quartiles of continuous delta measure

summary(all$deltac)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.004542 -0.001643 -0.001166  0.000000  0.001606  0.007140
high_D <- subset(all, deltac >=.001571) #impatient
low_D <- subset(all, deltac <= -.001713) #patient
summary(glm(DV ~ cond, data=high_D, family="binomial")) 
## 
## Call:
## glm(formula = DV ~ cond, family = "binomial", data = high_D)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7495  -0.7495  -0.6536  -0.6536   1.8158  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4351     0.3519  -4.079 4.53e-05 ***
## cond3         0.3091     0.4839   0.639    0.523    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 105.88  on 100  degrees of freedom
## Residual deviance: 105.47  on  99  degrees of freedom
## AIC: 109.47
## 
## Number of Fisher Scoring iterations: 4
summary(glm(DV ~ cond, data=low_D, family="binomial"))
## 
## Call:
## glm(formula = DV ~ cond, family = "binomial", data = low_D)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1574   0.4531   0.4531   0.7215   0.7215  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.2130     0.3434   3.532 0.000412 ***
## cond3         1.0116     0.6284   1.610 0.107449    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 80.734  on 88  degrees of freedom
## Residual deviance: 77.888  on 87  degrees of freedom
## AIC: 81.888
## 
## Number of Fisher Scoring iterations: 4

the proportions these tests give are the % of people choosing SS

prop.test(table(high_D$cond, high_D$DV)) #impatient people choose SS often, but default condition didn't make a difference
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(high_D$cond, high_D$DV)
## X-squared = 0.159, df = 1, p-value = 0.6901
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1283874  0.2335679
## sample estimates:
##    prop 1    prop 2 
## 0.8076923 0.7551020
prop.test(table(low_D$cond, low_D$DV)) #patient people choose SS seldom, and default cond makes a difference
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(low_D$cond, low_D$DV)
## X-squared = 1.8745, df = 1, p-value = 0.171
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04062698  0.30383836
## sample estimates:
##     prop 1     prop 2 
## 0.22916667 0.09756098

note: hid (patient people) make a ‘preference-consistent’ choice when they choose the LL option. Since DV is coded as 1=LL 0=SS, the DV variable is also a ‘chose pref-consistent option’variable. lod(impatient people) make a ’pref-consistent’ choice when DV=0, so for these people the proportion choosing the pref-consistent option is 100 - the proportion given by the prop tests above

significance tests for 2 proportions at a time:

hiD12 <- subset(high_D, cond!=3)
hiD12$cond <- droplevels(hiD12$cond)
prop.test(table(hiD12$cond, hiD12$DV))
## 
##  1-sample proportions test with continuity correction
## 
## data:  table(hiD12$cond, hiD12$DV), null probability 0.5
## X-squared = 18.4808, df = 1, p-value = 1.716e-05
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6703039 0.8991531
## sample estimates:
##         p 
## 0.8076923
hiD13 <- subset(high_D, cond!=2)
hiD13$cond <- droplevels(hiD13$cond)
prop.test(table(hiD13$cond, hiD13$DV))
## 
##  1-sample proportions test with continuity correction
## 
## data:  table(hiD13$cond, hiD13$DV), null probability 0.5
## X-squared = 11.7551, df = 1, p-value = 0.0006068
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6082190 0.8619044
## sample estimates:
##        p 
## 0.755102
hiD23 <- subset(high_D, cond!=1)
hiD23$cond <- droplevels(hiD23$cond)
prop.test(table(hiD23$cond, hiD23$DV)) #none of these comparisons are significant - the default condition doesn't seem to shift choice for impatient people
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(hiD23$cond, hiD23$DV)
## X-squared = 0.159, df = 1, p-value = 0.6901
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1283874  0.2335679
## sample estimates:
##    prop 1    prop 2 
## 0.8076923 0.7551020
loD12 <- subset(low_D, cond!=3)
loD12$cond <- droplevels(loD12$cond)
prop.test(table(loD12$cond, loD12$DV))
## 
##  1-sample proportions test with continuity correction
## 
## data:  table(loD12$cond, loD12$DV), null probability 0.5
## X-squared = 13.0208, df = 1, p-value = 0.000308
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.1251445 0.3766604
## sample estimates:
##         p 
## 0.2291667
loD13 <- subset(low_D, cond!=2)
loD13$cond <- droplevels(loD13$cond)
prop.test(table(loD13$cond, loD13$DV))
## 
##  1-sample proportions test with continuity correction
## 
## data:  table(loD13$cond, loD13$DV), null probability 0.5
## X-squared = 24.9756, df = 1, p-value = 5.806e-07
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.03172088 0.24059474
## sample estimates:
##          p 
## 0.09756098
loD23 <- subset(low_D, cond!=1)
loD23$cond <- droplevels(loD23$cond)
prop.test(table(loD23$cond, loD23$DV))  #two of these are at least marginal - the default condition seems capable of shifting choice for patient people
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(loD23$cond, loD23$DV)
## X-squared = 1.8745, df = 1, p-value = 0.171
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04062698  0.30383836
## sample estimates:
##     prop 1     prop 2 
## 0.22916667 0.09756098