###################################################
# This code has been extracted from jaNein.Rnw
# as a supplement to our paper
# Will You Accept Without Knowing What?-A Thuringian Newspaper Experiment of the Yes-No Game"
# by Werner Gueth and Oliver Kirchkamp
# http://www.kirchkamp.de/research/yesNoGame.html
###################################################
#
# please note that in our paper we use data from
# Gehrig, T., Gueth, W., Levati, V., Levinsky, R., Ockenfels, A., Uske,
# T. Weiland, T. (2007). Buying a pig in a poke: An experimental study
# of unconditional veto power. Journal of Economic Psychology, 28(6),
# pp. 692--703.
# Their data is not included in our on-line data set.
# Please contact Gehrig et.al. to obtain their data.


###################################################
### code chunk number 1: init
###################################################
library(Hmisc)
library(xtable)
library(memisc)
library(lme4)
library(geepack)
library(lattice)
library(latticeExtra)
options(scipen=8)

###################################################
### code chunk number 2: init2
###################################################
bootstrapsize<-1000
load("jaNein.Rdata")
levels(jn$type)<-sub("lab large stakes","lab high stakes",levels(jn$type))
pdata=subset(jnL100,!is.na(age) & sex!=""  )
pdata<-within(pdata,{
  age2<-age*age
  a1000<-ifelse(amount==1000,1,0)
  sex<-as.factor(as.character(sex))
})


###################################################
### code chunk number 3: charPart
###################################################
tab<-aggregate(c(length(type),100*mean(!is.na(sex)),100*mean(sex=="female",na.rm=TRUE),100*mean(!is.na(age)),mean(age,na.rm=TRUE),median(age,na.rm=TRUE),100*mean(!is.na(whiteCollar)),100*mean(!whiteCollar,na.rm=TRUE)) ~ type ,data=jn)
tab<-rbind(tab,NA)
dd<-dim(tab)
levels(tab$type)<-c(levels(tab$type),"all")
tab[dd[1],2:dd[2]]<-with(jn,c(length(type),100*mean(!is.na(sex)),100*mean(sex=="female",na.rm=TRUE),100*mean(!is.na(age)),mean(age,na.rm=TRUE),median(age,na.rm=TRUE),100*mean(!is.na(whiteCollar)),100*mean(!whiteCollar,na.rm=TRUE)))
tab[dd[1],1]<-"all"
digits<-c(0,0,0,rep(1,dim(tab)[2]-2))
colnames(tab)<-c("experiment","\\rrr{num. of}{participants}","\\rrr{known}{sex [\\%]}","\\rr{females [\\%]}","\\rrr{known}{age [\\%]}","\\rrr{mean age}{[years]}","\\rrr{median age}{[years]}","\\rrr{known pro-}{fession [\\%]}","\\rrr{blue collar}{[\\%]}")
print(xtable(tab,digits=digits),include.rownames=FALSE,floating=FALSE,hline.after=c(-1,0,nrow(tab)-1,nrow(tab)),sanitize.colnames.function=function(x) x)


###################################################
### code chunk number 4: ageDist
###################################################
xx <- subset(jn,!is.na(age))
xx<-within(xx,{type<-factor(type);type<-reorder(type,age)})
myKey<-list(x=1,y=0,corner=c(1,0),border=TRUE,points=FALSE,lines=TRUE)
mySettings<-list(axis.text=list(cex=1.4),par.xlab.text=list(cex=1.4))
#plot(ecdfplot(~ age ,group= type,data=xx,auto.key=myKey))
plot(bwplot(type~age,data=xx,par.settings=mySettings))


###################################################
### code chunk number 5: offersBoxFig
###################################################
plot(bwplot(offer ~ as.factor(amount) | reorder(place,-offer,length),data=jnL10,layout=c(4,1),scales=list(x=list(rot=45,relation="free")),ylab="relative offer"))


###################################################
### code chunk number 6: rejectRatesFig
###################################################
xx <- with(jnL10,aggregate(c(reject=mean(reject,na.rm=TRUE)) ~ type + amount))
xx <- within(xx,{levels(type)[grep("^lab",levels(type))]<-"lab";type<-factor(type)})
xx<-subset(xx,!is.nan(reject))
myKey<-list(x=0,y=.5,border=TRUE,points=TRUE,lines=TRUE)
myScale<-list(x=list(log=TRUE,at=c(20,100,1000)))
plot(xYplot(reject ~ amount, group=reorder(type,-reject),label.curves=list(tilt=TRUE,cex=1),data=xx,type="b",ylab="freq. of rejection",xlab="amount [Euro]",scales=myScale))


###################################################
### code chunk number 7: wtest
###################################################
require(exactRankTests,quietly=TRUE)
require(coin)
moffer100 <-  (mean(jn$offer100,na.rm=TRUE))
moffer1000<- (mean(jn$offer1000,na.rm=TRUE))
(tt <- with(jn,t.test(offer1000,offer100,paired=TRUE)))
(wt <- with(jn,wilcox.exact(offer100, offer1000,paired=TRUE,exact=TRUE)))
rejRates<-round(with(jnL100,aggregate(reject,list(lab),mean,na.rm=TRUE))*100,1)$x


###################################################
### code chunk number 8: demAgeMF
###################################################
myKey<-list(x=.5,y=.9,background="white",corner=c(.5,1),border=TRUE,points=FALSE,lines=TRUE)
myScale<-list(x=list(log=TRUE,at=c(10,100,1000)))
xx<-subset(jnL100,!is.na(sex))
plot(xyplot(offer + expOffer~ age ,xx, ylab="relative offer + expOffer",ylim=c(.22,.51),groups=factor(sex),panel=function(...) {panel.plsmo(...,label.curves=list(tilt=TRUE))},type="l",auto.key=FALSE) + layer(panel.smoother(...,lty=0,col.se=gray(.9),alpha.se=1),under=TRUE))
#


###################################################
### code chunk number 9: acceptAgeMF
###################################################
rPanel <- function(group,family='gaussian',ylim=c(0,.3)) {
  xyplot(reject ~ age ,data=xx,groups= group,ylim=ylim,ylab='reject + expReject',panel=function(...) {
    panel.plsmo(...,label.curves=list(tilt=TRUE))},type="l",auto.key=FALSE) + 
  layer(panel.smoother(...,lty=0,col.se=gray(.9),alpha.se=1),under=TRUE)
}
ePanel <- function(group,family='gaussian',ylim=c(0,.3)) {
  xyplot(expReject ~ age ,data=xx,groups= group,ylim=ylim,panel=function(...) {
    panel.plsmo(...,family='symmetric',label.curves=list(tilt=TRUE))},type="l",auto.key=FALSE) + 
  layer(panel.smoother(...,family='symmetric',lty=0,col.se=gray(.9),alpha.se=1),under=TRUE)
}
with(xx,plot(c(reject=rPanel(sex),expReject=ePanel(sex))))


###################################################
### code chunk number 10: demAge
###################################################
noage<-with(data=subset(jnL10,executive | gehrig),aggregate(c(expOffer=mean(expOffer,na.rm=TRUE),offer=mean(offer,na.rm=TRUE)) ~ type))
noage<-merge(noage,as.data.frame(list(age=20:40)))
noage<-within(noage, age[type=="executives"]<-age[type=="executives"]+40)
xx<-subset(jnL100,!is.na(age))
xx<-within(xx,{type<-factor(type);levels(type)[grep("lab",levels(type))]<-"labLargeSt";a2<-as.factor(amount);levels(a2)<-sprintf("%s Euro",levels(a2))})
qqLab<-quantile(subset(xx,lab)[["age"]],c(0,1))
x1<-merge(subset(xx,!lab | (age>=qqLab[1] &  age<=qqLab[2])),noage,all=TRUE)
x1<-within(x1,{levels(type)[grep("lab small",levels(type))]<-"labSmallSt";type<-reorder(factor(type),-offer)})
myKey<-list(x=.5,y=.9,background="white",corner=c(.5,1),border=TRUE,points=FALSE,lines=TRUE)
myScale<-list(x=list(log=TRUE,at=c(10,100,1000)))
plot(xyplot(offer + expOffer ~ age,groups="type",ylab="relative offer + expOffer",ylim=c(.22,.51),data=xx,panel=function(...) panel.smoother(...,lty=0,col.se=gray(.9),alpha.se=1),type="l") + xyplot(offer + expOffer~ age , groups=type,panel=function(...) {panel.plsmo(...,label.curves=list(tilt=TRUE))},x1,type="l",auto.key=FALSE))
#


###################################################
### code chunk number 11: acceptAge
###################################################
with(xx,plot(c(reject=rPanel(a2),expReject=ePanel(a2))))


###################################################
### code chunk number 12: offExp
###################################################
expect <- with(subset(jnL100,!executive),aggregate(offer,list(offer=offer,expectation=expOffer),length))
plot(with(expect,xyplot(expectation ~ offer,cex=.22*sqrt(x),xlab="own relative offer",ylab="expected relative offer"))+layer(panel.abline(lm(expectation ~ offer,weight=x,data=expect),lty=1))+layer(panel.abline(a=0,b=1,lty=2))+layer(panel.key(text=(c("OLS",expression(45^o))),corner=c(1,1),background="white",border=TRUE,points=FALSE,lines=TRUE)))


###################################################
### code chunk number 13: rejExp
###################################################
par(mar=c(4.5,4.6,1.3,0),cex=1,mex=0.7)
expect2 <- with(subset(jnL100,!executive),table(reject,expReject))
rownames(expect2)=c("accept","reject")
expect2<-expect2[,dim(expect2)[2]:1]
colnames(expect2)[c(2:5,8)]<-""
mosaicplot(expect2,off=0,sort=2:1,dir="h",main="",xlab="own choice",ylab="expected rejection rate",las=1,col=gray(c(1,.5)),cex.axis=1)


###################################################
### code chunk number 14: HistD2
###################################################
type=""
jnL10<-within(jnL10,{a2<-as.factor(amount);levels(a2)<-sprintf("%s Euro",levels(a2))})
plot(densityplot(~offer ,groups=a2,subset=amount>20,data=jnL10,cut=0,auto.key=TRUE,type=type))


###################################################
### code chunk number 15: HistD1
###################################################
myKey<-list(x=1,y=1,background="white",corner=c(1,1),border=TRUE,points=FALSE,lines=TRUE)
plot(densityplot(~offer ,groups=type,subset=!gehrig,data=jnL10,cut=0,auto.key=TRUE,type=type))
# plot(ecdfplot(~offer ,groups=type,data=jnL10,cut=0,auto.key=TRUE))


###################################################
### code chunk number 16: HistD3
###################################################
plot(densityplot(~offer ,groups=sex,data=jnL10,cut=0,auto.key=TRUE,type=type))


###################################################
### code chunk number 17: offer1
###################################################
offerws.lme <- lmer (offer ~ a1000 + age*sex + internet + lab  + sex + whiteCollar + (1|id),data=pdata)
offerw.lme <- lmer (offer ~ a1000 + age + age2 + internet + lab  + sex + whiteCollar + (1|id),data=pdata)
offer1a.lme <- lmer (offer ~ a1000 + age + age2 + internet + lab  + sex + (1|id),data=pdata)
offer1.lme <- lmer (offer ~ a1000 + age + internet + lab + sex  + (1|id),data=pdata)
offer0.lme <- lmer (offer ~ a1000 + age + internet + sex  + (1|id),data=pdata)
set.seed(123)
mtable("1"=offer0.lme,"2"=offer1.lme,"3"=offer1a.lme,"4"=offerw.lme,"5"=offerws.lme,coef.style="pci",summary.stats=c("indep.obs.","N"))


###################################################
### code chunk number 19: reject1
###################################################
reject1.lme <- lmer (reject ~ lab*a1000 + expOffer + age + internet + sex + (1|id),data=pdata,family=binomial(link = "logit"))


###################################################
### code chunk number 21: offerE
###################################################
pdata<-within(pdata,rexpoffer<-expOffer/amount)
offerE.lme <- lmer (expOffer ~  a1000 + age + internet + lab + sex + (1|id),data=pdata)
offerEo.lme <- lmer (expOffer ~ offer + a1000 + age + internet + lab + sex + (1|id),data=pdata)
set.seed(123)
mtable("1"=offerE.lme,"2"=offerEo.lme,coef.style="pci",summary.stats=c("indep.obs.","N"))



###################################################
### code chunk number 23: rejectE
###################################################
rejectE.lme <- glmer(expReject ~  age + internet + lab*a1000 + sex + (1|id),data=pdata,family=binomial(link = "logit"))