setwd("~/Documents/Dropbox/Research/Eric/study 2a replication")
all<-read.csv("study1+rep_ALLDVs.csv", header=T, sep=",")
all$betac <- all$beta - mean(all$beta)
all$deltac<- all$delta - mean(all$delta)

all$cond[all$cond1SS2LL3NO ==1]<-"SS"
all$cond[all$cond1SS2LL3NO ==2]<-"LL"
all$cond[all$cond1SS2LL3NO ==3]<-"no"

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

summary(glm(SS~betac, data=all, family=binomial))
## 
## Call:
## glm(formula = SS ~ betac, family = binomial, data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6246  -1.0269  -0.8507   1.1776   2.0078  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.1059     0.0634  -1.670   0.0949 .  
## betac        -4.6126     0.5330  -8.654   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1522.0  on 1100  degrees of freedom
## Residual deviance: 1422.8  on 1099  degrees of freedom
## AIC: 1426.8
## 
## Number of Fisher Scoring iterations: 4
s1<-(glm(SS~betac, data=all, family=binomial))
library(nlme)
## Warning: package 'nlme' was built under R version 3.2.5
logLik(s1)
## 'log Lik.' -711.3866 (df=2)
library(lme4)
## Warning: package 'lme4' was built under R version 3.2.5
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
nest.reg <- glmer(SS ~  (1|cond)+betac, family = binomial, data = all)
coef(nest.reg) # this will give the estimated coeficients by condition
## $cond
##    (Intercept)     betac
## LL -0.38096863 -4.702969
## no -0.07243852 -4.702969
## SS  0.12606050 -4.702969
## 
## attr(,"class")
## [1] "coef.mer"
?glmer
fixef(nest.reg) # this will give the model averaging over all conditions
## (Intercept)       betac 
##  -0.1094669  -4.7029688
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.27150170
## no  0.03702841
## SS  0.23552743
library(arm)
## Loading required package: MASS
## 
## arm (Version 1.8-6, built: 2015-7-7)
## Working directory is /Users/NeuroethicsMacBookPro11/Documents/Dropbox/research/Eric
display(nest.reg)
## glmer(formula = SS ~ (1 | cond) + betac, data = all, family = binomial)
##             coef.est coef.se
## (Intercept) -0.11     0.15  
## betac       -4.70     0.54  
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  cond     (Intercept) 0.23    
##  Residual             1.00    
## ---
## number of obs: 1101, groups: cond, 3
## AIC = 1419.6, DIC = 1398.6
## deviance = 1406.1
logLik(nest.reg)
## 'log Lik.' -706.7798 (df=3)

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                   1100     1522.0              
## betac  1    99.21      1099     1422.8 < 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 ~ (betac|cond), family = binomial, data = all)
coef(nest.reg2) 
## $cond
##        betac (Intercept)
## LL -4.519581  -0.3686238
## no -1.919259  -0.1740988
## SS -8.186344   0.2088735
## 
## attr(,"class")
## [1] "coef.mer"
fixef(nest.reg2) 
## (Intercept) 
##  -0.4199777
ranef(nest.reg2)
## $cond
##    (Intercept)     betac
## LL  0.05135394 -4.519581
## no  0.24587890 -1.919259
## SS  0.62885125 -8.186344
display(nest.reg2)
## glmer(formula = SS ~ (betac | cond), data = all, family = binomial)
## coef.est  coef.se 
##    -0.42     0.28 
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  cond     (Intercept) 0.41           
##           betac       5.67     -0.88 
##  Residual             1.00           
## ---
## number of obs: 1101, groups: cond, 3
## AIC = 1413.2, DIC = 1364.6
## deviance = 1384.9
logLik(nest.reg2)
## 'log Lik.' -702.594 (df=4)

And test these 2 models against each other

anova(nest.reg, nest.reg2, test='Chisq')
## Data: all
## Models:
## nest.reg: SS ~ (1 | cond) + betac
## nest.reg2: SS ~ (betac | cond)
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## nest.reg   3 1419.6 1434.6 -706.78   1413.6                            
## nest.reg2  4 1413.2 1433.2 -702.59   1405.2 8.3716      1   0.003811 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ryan0 <- glm(SS~betac,family=binomial,data=all)
ryan1 <- glm(SS~cond+betac,family=binomial,data=all)
ryan2 <- glm(SS~cond*betac,family=binomial,data=all)
anova(ryan0,ryan1,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: SS ~ betac
## Model 2: SS ~ cond + betac
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1099     1422.8                          
## 2      1097     1405.6  2   17.211 0.0001831 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(ryan1,ryan2,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: SS ~ cond + betac
## Model 2: SS ~ cond * betac
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1097     1405.6                          
## 2      1095     1384.0  2   21.607 2.033e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1