library(foreign) Mdata <- read.dta("marbles.dta") attach(Mdata) Half <- rep(c(-1,1),c(150,150)) Ss <- unique(subject)[-1] # first subject is blank Session <- unique(cbind(subject,session))[-1,2] Ss <- Ss[-127] # subject 127 was an outlier, <.4 Session <- Session[-127] Etable <- matrix(NA,161,6) # table for Epsison effect Btable <- matrix(NA,161,4) # table for Base Epsilon <- Base <- Good <- Win <- matrix(NA,161,300) for (i in 1:161) {w1 <- which(subject==Ss[i]) Good[i,round[w1]] <- goodchoi[w1] Base[i,round[w1]] <- round(10*baseodds[w1]-2.5) Win[i,round[w1]] <- win[w1] Epsilon[i,round[w1]] <- epsilon[w1] Etable[i,] <- tapply(Good[i,],Epsilon[i,],mean,na.rm=T) Btable[i,] <- tapply(Good[i,],Base[i,],mean,na.rm=T) } Opt1 <- rowMeans(Good[,1:150],na.rm=T) Opt2 <- rowMeans(Good[,151:300],na.rm=T) Opt <- (Opt1+Opt2)/2 Eps <- sign(Epsilon) # size doesn't matter #Coef1 <- Coef2 <- matrix(NA,162,3) #Coef <- matrix(NA,162,4) #colnames(Coef) <- c("Int","Base","Eps","Half") #colnames(Coef1) <- colnames(Coef2) <- c("Int","Base","Eps") #for (i in 1:162) # {Coef[i,] <- summary(glm(Good[i,]~Base[i,]+Eps[i,]+Half,family=binomial(link="probit")))$coef[,1] # Coef1[i,] <- summary(glm(Good[i,Half<0]~Base[i,Half<0]+Eps[i,Half<0],family=binomial(link="probit")))$coef[,1] # Coef2[i,] <- summary(glm(Good[i,Half>0]~Base[i,Half>0]+Eps[i,Half>0],family=binomial(link="probit")))$coef[,1]} # Inc10h1 <- (Session==3) + (Session==5) - (Session==1) - (Session==4) - (Session==6) # Inc10h1[Inc10h1==0] <- NA # Inc10h2 <- (Session==3) + (Session==4) - (Session==1) - (Session==5) - (Session==6) # Inc10h2[Inc10h2==0] <- NA # between-s tests of incentive #t.test(Coef1[,1]~Inc10h1,var.eq=T) # not close #t.test(Coef2[,1]~Inc10h2,var.eq=T) # wrong way but n.s. # within-s test of incentive #t.test(Coef[Session==4,1],Coef[Session==5,1],var.eq=T) Bias1 <- rowMeans(Good[,Half<0] * Eps[,Half<0],na.rm=T) Bias2 <- rowMeans(Good[,Half>0] * Eps[,Half>0],na.rm=T) Bias <- (Bias1+Bias2)/2 #t.test(Bias[Session==1],Bias[Session==3]) # n.s. t.test((Opt1-Opt2)[Session==4],(Opt1-Opt2)[Session==5],var.eq=T) # t.test((Bias1-Bias2)[Session==4],(Bias1-Bias2)[Session==5],var.eq=T) # not quite signif one tailed Opt1s <- Opt1[Session==4 | Session==5] Opt2s <- Opt2[Session==4 | Session==5] Bias1s <- Bias1[Session==4 | Session==5] Bias2s <- Bias2[Session==4 | Session==5] S45 <- Session[Session==4 | Session==5] Opt12s <- Opt2s-Opt1s # note reverse Bias12s <- Bias1s-Bias2s #summary(lm(Bias12s~S45+Opt12s)) se <- function(x) {sqrt(var(x)/length(x))} B1 <- tapply(Bias1,Session,mean) B2 <- tapply(Bias2,Session,mean) M1 <- tapply(Opt1,Session,mean) M2 <- tapply(Opt2,Session,mean) B <- rbind("1st"=B1,"2nd"=B2) M <- rbind("1st"=M1,"2nd"=M2) B1se <- tapply(Bias1,Session,se) B2se <- tapply(Bias2,Session,se) M1se <- tapply(Opt1,Session,se) M2se <- tapply(Opt2,Session,se) Bse <- rbind("1st"=B1se,"2nd"=B2se) Mse <- rbind("1st"=M1se,"2nd"=M2se) stop() library(gplots) postscript(file="fig3.eps",family="NimbusSan",width=8,height=8, horiz=F,onefile=F,pointsize=20,paper="special") barplot2(B,beside=T,col=gray(c(0,0,.5,.5,1,1,0,1,1,0,0,0)), plot.ci=T,ci.u=Bse+B,ci.l=B,xlab="Session",ylab="Mean bias") dev.off() postscript(file="fig2.eps",family="NimbusSan",width=8,height=8, horiz=F,onefile=F,pointsize=20,paper="special") barplot2(M,beside=T,col=gray(c(0,0,.5,.5,1,1,0,1,1,0,0,0)), plot.ci=T,ci.u=Mse+M,ci.l=M,xlab="Session",ylab="Proportion optimal") dev.off() #### sequential effects table(Good[,2:300],Good[,1:299],Win[,1:299])