setwd("~/Documents/Dropbox/Research/Eric/study 2a replication")
combined2<-read.csv("combined2.csv", header=T, sep=",")
combined2$bd<-combined2$beta*((1/(1+combined2$delta))^(2/52))*1.5

First use only intercept and time preference to predict choice of SS

summary(glm(SS~bd, data=combined2, family=binomial))
## 
## Call:
## glm(formula = SS ~ bd, family = binomial, data = combined2)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.695  -1.069  -0.653   1.169   1.823  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.8782     0.5150   7.530 5.06e-14 ***
## bd           -2.6548     0.3431  -7.739 1.01e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1509.0  on 1088  degrees of freedom
## Residual deviance: 1432.7  on 1087  degrees of freedom
## AIC: 1436.7
## 
## Number of Fisher Scoring iterations: 3
s1<-(glm(SS~bd, data=combined2, family=binomial))
library(nlme)
logLik(s1)
## 'log Lik.' -716.3642 (df=2)

Then look at the effect of condition as a 3 level factor

summary(glm(SS~cond, data=combined2, family=binomial))
## 
## Call:
## glm(formula = SS ~ cond, family = binomial, data = combined2)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.249  -1.111  -1.100   1.246   1.257  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.18514    0.10171  -1.820   0.0687 .
## condno       0.02587    0.15386   0.168   0.8665  
## condSS       0.35177    0.14314   2.457   0.0140 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1509.0  on 1088  degrees of freedom
## Residual deviance: 1501.7  on 1086  degrees of freedom
## AIC: 1507.7
## 
## Number of Fisher Scoring iterations: 3

Then look at the first nested model

library(lme4)
nest.reg <- glmer(SS ~  (1|cond), family = binomial, data = combined2)
coef(nest.reg) # this will give the estimated coeficients by #condition.
## $cond
##    (Intercept)
## LL -0.13365788
## no -0.11158158
## SS  0.07981324
## 
## attr(,"class")
## [1] "coef.mer"
fixef(nest.reg) # this will give the model averaging over all #conditions
## (Intercept) 
## -0.05519815
ranef(nest.reg) # the errors (specificity) at the condition #level. If you sum fixef with ranef you will get coef results
## $cond
##    (Intercept)
## LL -0.07845973
## no -0.05638342
## SS  0.13501139
library(arm)
display(nest.reg)
## glmer(formula = SS ~ (1 | cond), data = combined2, family = binomial)
## coef.est  coef.se 
##    -0.06     0.09 
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  cond     (Intercept) 0.13    
##  Residual             1.00    
## ---
## number of obs: 1089, groups: cond, 3
## AIC = 1511.3, DIC = 1498.5
## deviance = 1502.9
logLik(nest.reg)
## 'log Lik.' -753.6315 (df=2)

We can then test the simple model against the nested model

anova(s1, nest.reg, test="Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: SS
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                  1088     1509.0              
## bd    1   76.277      1087     1432.7 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now nest time preference in condition

nest.reg2 <- glmer(SS ~  (bd|cond), family = binomial, data = combined2)
coef(nest.reg2) 
## $cond
##           bd (Intercept)
## LL -1.028199    1.307330
## no -3.048204    4.479321
## SS -5.285455    7.992450
## 
## attr(,"class")
## [1] "coef.mer"
fixef(nest.reg2) 
## (Intercept) 
##  -0.3072392
ranef(nest.reg2)
## $cond
##    (Intercept)        bd
## LL    1.614569 -1.028199
## no    4.786560 -3.048204
## SS    8.299689 -5.285455
display(nest.reg2)
## glmer(formula = SS ~ (bd | cond), data = combined2, family = binomial)
## coef.est  coef.se 
##    -0.31     0.13 
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  cond     (Intercept) 5.77           
##           bd          3.67     -1.00 
##  Residual             1.00           
## ---
## number of obs: 1089, groups: cond, 3
## AIC = 1421.3, DIC = 1385.4
## deviance = 1399.3
logLik(nest.reg2)
## 'log Lik.' -706.6602 (df=4)

And test these 2 models against each other

anova(nest.reg, nest.reg2, test="Chisq")
## Data: combined2
## Models:
## nest.reg: SS ~ (1 | cond)
## nest.reg2: SS ~ (bd | cond)
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## nest.reg   2 1511.3 1521.2 -753.63   1507.3                             
## nest.reg2  4 1421.3 1441.3 -706.66   1413.3 93.943      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here is the regular analysis, using the beta-delta combined measure of present value of LL, first for SS+control conditions:

SScontrol<-subset(combined2, cond!="LL")
SScontrol<-droplevels(SScontrol)

LLcontrol<-subset(combined2, cond!="SS")
LLcontrol<-droplevels(LLcontrol)

summary(glm(SS ~ cond * bd, data=SScontrol, family="binomial"))
## 
## Call:
## glm(formula = SS ~ cond * bd, family = "binomial", data = SScontrol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9784  -1.0341   0.1524   1.1550   1.9788  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.7602     1.1065   4.302 1.69e-05 ***
## condSS        3.4452     1.6631   2.072   0.0383 *  
## bd           -3.2706     0.7253  -4.509 6.51e-06 ***
## condSS:bd    -2.1291     1.0951  -1.944   0.0519 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 968.90  on 698  degrees of freedom
## Residual deviance: 867.08  on 695  degrees of freedom
## AIC: 875.08
## 
## Number of Fisher Scoring iterations: 5
f1<-(glm(SS~cond*bd, data=SScontrol, family="binomial"))
library(visreg)
col<-c("black", "blue")
col2<-c("gray80", "cornflowerblue")
col3<-adjustcolor(col2, alpha=.7)
visreg(f1, "bd", by="cond", line=list(col=col), fill=list(col=col3), overlay=TRUE, partial=FALSE,  scale="response", xlab="Present Value of LL", ylab="Percent Choosing SS")

And for LL+control conditions - note here that the pattern is inconsistent with the idea that the default is helping people align their preferences with their choice. The LL default appears to make patient people more likely to choose SS and impatient people less likely to choose SS.

summary(glm(SS~cond*bd, data=LLcontrol,family="binomial"))
## 
## Call:
## glm(formula = SS ~ cond * bd, family = "binomial", data = LLcontrol)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9784  -1.0518  -0.8927   1.2471   1.9788  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.3139     0.6582   1.996  0.04591 * 
## condno        3.4462     1.2875   2.677  0.00743 **
## bd           -1.0276     0.4454  -2.307  0.02105 * 
## condno:bd    -2.2430     0.8512  -2.635  0.00841 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 954.11  on 691  degrees of freedom
## Residual deviance: 919.26  on 688  degrees of freedom
## AIC: 927.26
## 
## Number of Fisher Scoring iterations: 4
f2<-(glm(SS~cond*bd, data=LLcontrol,family="binomial"))
col<-c("red", "black")
col2<-c("indianred1" , "gray80")
col3<-adjustcolor(col2, alpha=.45)
visreg(f2, "bd", by="cond", line=list(col=col), fill=list(col=col3), overlay=TRUE, partial=FALSE,  xlab="Present value of LL", scale="response", ylab="P (choosing SS)")