library(lattice) library(latticeExtra) library(lme4,lib.loc="~/R/lme4") # we use lme4_0.999999-2 here ! options(scipen=5) bootSize<-10000 library(languageR) library(xtable) library(plyr) library(car) aggregate.formula<-getFromNamespace("aggregate.formula","stats") print.xtable<-function(...) xtable:::print.xtable(floating=FALSE,...) testMER <- function(mcmc,r,tag="",mean) { x <- as.matrix(mcmc) %*% r stats<-cbind(mean,HPDinterval(x),2*min(mean(x<0),mean(x>0))) colnames(stats)<-c("Estimate","HPD95lower","HPD95upper","pMCMC") rownames(stats)<-tag stats } mSegplot<-function(...) { segplot(...,colorkey=FALSE,draw.bands=FALSE,segments.fun = panel.arrows,angle=90,length=1,unit='mm',ends="both") } share1a<-100*(1 - 0.8)/(1-.8*.9) share1Asym<-100*(1 - 0.95)/(1-.9*.95) share2a<-100-100*(1 - 0.9)/(1-.8*.9) share2Asym<-100-100*(1 - 0.9)/(1-.9*.95) # the data from Boehm and Vollstaedt is not included here load("asymBargain.Rdata") subGlob<-within(subGlob[grep("^Ind",subGlob$Treatmenttype,invert=TRUE),], Treatmenttype<-factor(Treatmenttype)) subGlob[["GroupInd"]]<-FALSE subGlob[with(subGlob,grep("^Group",Treatmenttype)),"GroupInd"]<-TRUE # subOffer<-within(subOffer[grep("^Ind",subOffer$Treatmenttype,invert=TRUE),], Treatmenttype<-factor(Treatmenttype)) subOffer[["GroupInd"]]<-FALSE subOffer[with(subOffer,grep("^Group",Treatmenttype)),"GroupInd"]<-TRUE subGlob <- subset(subGlob, Date!="091201_0959" & Date!="091201_1202") subOffer <- subset(subOffer, Date!="091201_0959" & Date!="091201_1202") subGlobAll <- dim(subGlob)[1] with(subGlob,table(Problem,ProblemPreview,useNA="always")) fracProb <- signif(100*(mean(with(subGlob,Problem | ProblemPreview),na.rm=TRUE)),3) subOffer <- subset(subOffer, vorschaurunde.o==-1) xx <- subOffer[,c("Date","Treatment","Treatmenttype","GroupInd","Period","Group","Subject","sid","Type","runde","runde.o","Situation","share_red", "share_red.o","accept", "accept.o", "valid_preview", "last_entry_offers", "vorschaurunde.o")] xx1 <- with(xx, aggregate(xx$valid_preview, list(Date=Date, Treatment=Treatment, Treatmenttype=Treatmenttype, Period=Period, Group=Group, Subject=Subject, runde.o=runde.o), max)) roundOffer <- merge(xx, xx1, by.x=c("Date", "Treatment", "Treatmenttype", "Period", "Group", "Subject", "runde.o", "valid_preview"), by.y=c("Date", "Treatment", "Treatmenttype", "Period", "Group", "Subject", "runde.o", "x")) roundOffer$Broff <- ifelse(roundOffer$Situation %in% c("break_off"), 1, 0) aa <- subset(roundOffer[,c("Period","sid","share_red.o","Type","runde.o")], (Type=="redA"|Type=="redB") & runde.o==1)[,c("Period","sid","share_red.o")] names(aa)[3] <- "share_red.o1" groupOffers <- subset(subOffer, (runde==runde.o | Situation %in% c("break_off", "wait_screen")))[,c("Date","Treatment","Treatmenttype","GroupInd","Period","Group","Subject","sid","Type","runde","runde.o","Situation","share_red")] groupOffers$Broff <- ifelse(groupOffers$Situation %in% c("break_off", "wait_screen"), 1, 0) groupOffers <- with(groupOffers, aggregate(groupOffers[,c("runde","runde.o","Broff")], list(Date=Date,Treatment=Treatment,Treatmenttype=Treatmenttype,Period=Period,Subject=Subject, sid=sid, Type=Type,Group=Group, share_red=share_red), max)) groupOffers$runde.oB <- ifelse(groupOffers$Broff==0, groupOffers$runde.o, groupOffers$runde.o+1) groupOffers$runde.oBln <- log(groupOffers$runde.oB) groupOffers <- merge(groupOffers, aa, all=TRUE) maxR <- max(groupOffers[["runde"]]) xx1 <- with(groupOffers, aggregate(runde, list(runde=runde, broff=Broff), length)) xx.o <- with(groupOffers, aggregate(runde, list(runde=runde.o, broff=Broff),length)) breakFrak <- signif(100*mean(groupOffers$Broff), 3) avRound <- signif(mean(subset(groupOffers,Broff==0)$runde), 3) avRoundS <- signif(mean(subset(groupOffers,Broff==1)$runde.o), 3) avRoundS2 <- signif(mean(subset(groupOffers,Broff==1)$runde), 3) xx<-breakFrakL<-aggregate(Broff~Treatmenttype,FUN=mean,data=subset(groupOffers,substring(Treatmenttype,1,5)!="Group")) breakFrakL<-xx$Broff names(breakFrakL)<-xx$Treatmenttype breakFrakL<-signif(breakFrakL*100,3) xx<-with(subset(roundOffer,Type %in% c("blueA","blueB","blue")),aggregate(runde.o,list(sid=sid,Treatment=Treatment,Treatmenttype=Treatmenttype,Period=Period),min)) acceptance1st<-subset(join(roundOffer,xx),runde.o==x) share1a<-100*(1 - 0.8)/(1-.8*.9) share1Asym<-100*(1 - 0.95)/(1-.9*.95) share2a<-100-100*(1 - 0.9)/(1-.8*.9) share2Asym<-100-100*(1 - 0.9)/(1-.9*.95) xx <- subset(roundOffer, Period<6 & runde.o==1 & ((Treatmenttype=="2playerSingle" & (Type=="redA" | Type=="redB")) | ((Treatmenttype=="2player" | Treatmenttype=="3playerSym" | Treatmenttype=="3playerAsym" | Treatmenttype=="3playerPatient") & Type=="redA"))) xxa<-aggregate(share_red.o ~ Period + Treatmenttype,data=xx,FUN=mean,na.rm=TRUE) xyplot(share_red.o ~ Period | Treatmenttype,data=xxa,layout=c(NA,1),t="l",par.strip.text = list(cex = 0.8),ylab="proposer's share",ylim=extendrange(xxa$share_red.o),between=list(x=c(1,0,0,0))) xx<-subset(roundOffer, runde.o==1 & Period>0 & ((Treatmenttype=="2playerSingle" & (Type=="redA" | Type=="redB")) | ((Treatmenttype %in% c("2player","3playerSym","3playerAsym","3playerPatient")) & Type=="redA") | Treatmenttype %in% c("GGroup0809","GGroup0908") )) # xx<-within(xx,{Treatmenttype<-factor(Treatmenttype);sro<-rank(share_red.o)/length(share_red.o)}) bootSize<-10000 est<-by(xx,list(xx$Treatmenttype),function(x) { mer<-lmer(share_red.o ~ 1 + (1|Date) + (1|sid),data=x) set.seed<-123 # <- this is, of course, not very clever, but it should not hurt, either # it should be set.seed(123) instead :-( mcmc<-mcmcsamp(mer,bootSize) list(mer=mer,mcmc=mcmc,name=x$Treatmenttype[1]) }) linTest <- function(i,j,est) { diff<-est[[i]]$mcmc@fixef["(Intercept)",]-est[[j]]$mcmc@fixef["(Intercept)",] res<-cbind(mean=fixef(est[[i]]$mer)-fixef(est[[j]]$mer), rbind(quantile(diff,c(0.025,.975))),p=mean(diff<0)) rownames(res)<-sprintf("%s-%s",est[[i]]$name,est[[j]]$name) res } as.data.frame(rbind(linTest(2,3,est),linTest(3,5,est),linTest(4,5,est),linTest(3,4,est))) propShare<-as.data.frame(t(sapply(est,function(x) { hh<-HPDinterval(x$mcmc)$fixef["(Intercept)",] c(fixef(x$mer),hh[1],hh[2],Std.Error=sd(x$mcmc@fixef),p=2*(mean(x$mcmc@fixef<0)),sqrt(unlist(VarCorr(x$mer)))) }))) colnames(propShare)[1]<-"mean" propShare<-within(propShare,{name<-reorder(factor(rownames(propShare)),mean)}) mSegplot(name ~ lower + upper,center=mean,data=propShare,xlab=list(label="share of proposer",cex=.7),xlim=range(share1Asym-2,share1a+2))+layer(panel.abline(v=c(share1Asym,share1a),lty=2:1,lwd=3)) diffEst<-as.data.frame(rbind(linTest(2,3,est),linTest(3,5,est),linTest(4,5,est),linTest(3,4,est))) colnames(diffEst)[2:4]<-c("95\\%","\\!\\!\\!\\!\\!-conf.int.","p-value") print(xtable(diffEst,digits=c(0,2,2,2,4)),sanitize.colnames.function=function(x)x) fitBeta <- function(grid,Treatment) { x<-(diag(grid[["shares"]])-propShare[Treatment,"mean"]/100)^2 grid[["a"]][which(x==min(x))] } offers <- subset(roundOffer, runde.o==1 & ((Treatmenttype=="2playerSingle" & (Type=="redA" | Type=="redB")) | ((Treatmenttype %in% c("2player","3playerSym","3playerAsym","3playerPatient")) & Type=="redA") | (Treatmenttype %in% c("Group0809","Group0908") & Type=="red" )))[,c("Date","Treatment","Treatmenttype","Period","Group","share_red.o")] offers<-with(offers,aggregate(share_red.o,list(Treatment=Treatment,Treatmenttype=Treatmenttype,Date=Date,Period=Period,Group=Group),mean)) acceptance<-acceptance1st[,c("Date","Treatment","Treatmenttype","Period","Group","Type","sid","accept.o")] acceptance<-with(acceptance,aggregate(accept.o,list(Treatment=Treatment,Treatmenttype=Treatmenttype,Date=Date,Period=Period,Group=Group,Type=Type,sid=sid),mean)) names(acceptance)[8]<-"accept" names(offers)[6]<-"offer" offerAccept<-subset(join(acceptance,offers),match(Treatmenttype,c("Group0809","Group0908"),nomatch=0)==0) offerAccept<-subset(offerAccept, ! (Treatmenttype %in% c("2player") & Type=="blueB")) offerAccept<-within(offerAccept,{ levels(Treatmenttype)<-c(levels(Treatmenttype),"3pAsym.Strong"); Treatmenttype[Treatmenttype=="3playerAsym" & Type=="blueB"]<-"3pAsym.Strong" levels(Treatmenttype)[grep("3playerAsym",levels(Treatmenttype))]<-"3pAsym.Weak" Treatmenttype<-factor(as.character(Treatmenttype)) Treatmenttype<-relevel(Treatmenttype,"2playerSingle") }) accept.mer <- glmer(accept ~ offer + Treatmenttype -1 + (1|Date) + (1|sid),family=binomial,data=offerAccept) accept.coefs<-summary(accept.mer)@coefs logSlope<- accept.coefs["offer","Estimate"] accept.coefs<-(accept.coefs/-logSlope)[-1,] xx<-as.data.frame(cbind(accept.coefs[,c("Estimate","Std. Error" )]) %*% rbind(1,qnorm(c(.025,.5,.975)))) xx[["name"]]<-reorder(as.factor(sub("^Treatmenttype","",rownames(xx))),-c(1,2,5,4,6,3)) mSegplot(name ~ V1 + V3,center=V2,data=xx,xlab=list(label="share of proposer",cex=.7),xlim=c(share1Asym-3,share1a+3))+layer(panel.abline(v=c(share1Asym,share1a),lty=2:1,lwd=3)) zVal<-sqrt(c( linearHypothesis(accept.mer,"Treatmenttype3playerSym = Treatmenttype3playerPatient")$Chisq[2], linearHypothesis(accept.mer,"Treatmenttype3pAsym.Weak = Treatmenttype3pAsym.Strong")$Chisq[2], linearHypothesis(accept.mer,"Treatmenttype3pAsym.Weak - Treatmenttype3pAsym.Strong = Treatmenttype3playerSym - Treatmenttype3playerPatient")$Chisq[2] )) cbind(accept.mer@fixef) Estimate<-rbind(accept.mer@fixef) %*% cbind(c(0,0,0,0,0,-1,1),c(0,0,0,-1,1,0,0),c(0,0,0,1,-1,-1,1)) xx<-as.data.frame(t(rbind(Estimate,Estimate/zVal)) %*% rbind(1,qnorm(c(.025,.5,.975))) / -logSlope) xx$name<-c("(1) 3plSym-3plPatient","(2) Weak-Strong","(1) - (2)") xx<-within(xx,name<-reorder(as.factor(xx$name),c(3,2,1))) mSegplot(name ~ V1 + V3,center=V2,data=xx,xlab=list(label="difference in demands",cex=.7),xlim=c(-10,-share1Asym+share1a+3)) + layer(panel.abline(v=c(-share1Asym+share1a,0))) numberRounds1a <- glmer(runde.oB ~ Treatmenttype - 1+ (1|Date) + (1|sid), family=poisson, data=subset(groupOffers,Type=="redA")) xx<-as.data.frame(cbind(summary(numberRounds1a)@coefs[,c("Estimate","Std. Error" )]) %*% rbind(1,qnorm(c(.025,.5,.975)))) xx[["name"]]<-reorder(as.factor(sub("^Treatmenttype","",rownames(xx))),c(3,4,2,1,0)) mSegplot(name ~ V1 + V3,center=V2,data=xx,xlab=list(label="$\\beta_t$",cex=.7))+layer(panel.abline(v=c(share1Asym,share1a),lty=2:1,lwd=3)) xx<-with(subset(subOffer,!GroupInd),aggregate(sid,list(Treatmenttype=Treatmenttype,Date=Date,sid=sid),length)) xx<-with(xx,aggregate(xx[,c("Date","sid")],list(Treatmenttype),function(x) length(unique(x)))) rownames(xx)<-xx[,1] colnames(xx)[2:3]<-c("sessions","participants") sessionOverview<-t(xx[,c("sessions","participants")]) xtable(sessionOverview,1,sum) # simpler <- function(xx,name) { zz<-as.data.frame(xx[,2]) rownames(zz)<-xx[,1] colnames(zz)<-name; zz } xx <- simpler(with(subset(roundOffer, runde.o==1 & ((Treatmenttype=="2playerSingle" & (Type=="redA" | Type=="redB")) | (Treatmenttype != "2playerSingle" & Type=="redA"))), aggregate(share_red.o, list(Treatmenttype=Treatmenttype), mean, na.rm=TRUE)),"Mean first demand (%)") # zz<-with(groupOffers, aggregate(runde.oB, list(Treatmenttype=Treatmenttype), mean)) zz<-simpler(zz,"Mean number of rounds") # yy<-simpler(with(subset(roundOffer, runde.o==2 & !GroupInd & ((Treatmenttype=="2player" & Type=="blueA") | ((Treatmenttype!="2player") & (Type=="blueA"|Type=="blueB")))), aggregate(share_red.o, list(Treatmenttype=Treatmenttype), mean, na.rm=TRUE)),"Mean second offer (%)") # bb <- cbind("Stopped participants (%)"=breakFrakL) allX<-t(round(cbind(xx,yy,zz,bb),2)) colnames(allX)<-sub("player","pl.",colnames(allX)) xtable(allX,caption="Summary statistics by treatment") colnames(propShare)<-c("$\\hat{\\beta}_0$","95\\%\\,-","\\!\\!\\!\\!\\!\\!conf.int.","$\\hat{\\sigma}_\\beta$","Pr$(>|\\hat{\\beta}_0|)$","$\\hat{\\sigma}_i$","$\\hat{\\sigma}_k$") print(xtable(propShare[,c(1,4,2,3,5,6,7)],digits=c(0,2,2,2,2,4,2,2)),sanitize.colnames.function=function(x)x) print.mer <- function(xx) { cc<-summary(xx)@coefs rr<-summary(xx)@REmat rownames(cc)<-sub("^Treatmenttype","",rownames(cc)) cat("\\begin{center}\\begin{tabular}{l}") print(xtable(cc,digits=c(2,2,3,2,4))) cat("\\\\\n") rr[,1]<-sub("^sid","i", rr[,1]) rr[,1]<-sub("^Date","k", rr[,1]) rr[,3:4]<-sprintf("%g",as.numeric(rr[,3:4])) print(xtable(rr,digits=c(2,2,2,8,8)),include.rownames=FALSE) cat("\\end{tabular}\\end{center}") } print.mer(accept.mer) print.mer(numberRounds1a)