### R commands to replicate results for ### Fabian Bornhorst, Andrea Ichino, Oliver Kirchkamp, Karl H. Schlag, Eyal Winter: ### "How do People Play a Repeated Trust Game? Experimental Evidence" ### Encoding: UTF-8 ### Note: The paper was published in 2010. ### For our calculations we also use the mcmcsamp command from the lme4 library which was part of lme4 at that time. ### Since lme4 1.0-4 (September 2013) the mcmcsamp command is no longer part of lme4. ### ### You have two options: ### ### Either you want to exactly replicate our results. ### Then you should use the mcmcsamp command from lme4 before 1.0-4. ### Please download the library from https://cran.r-project.org/src/contrib/Archive/lme4/lme4_0.999999-2.tar.gz ### Install the library with a command like the following: ### R CMD INSTALL -l . lme4_0.999999-2.tar.gz ### Load this version of the library from within R like this one: ### library(lme4,lib.loc=".") ### ### Alternatively you can use the bootMer command which is part of the modern lme4 packages. ### The command requires a different syntax and you may then get slightly different results. bsize=1000 bseed=123 set.seed(bseed) library(foreign) library(relaimpo) library(xtable) library(lme4,lib.loc=".") ################################################### ### some functions ################################################### mer2ctab <- function (mer,relimp=TRUE,type="pmvd",boot=NULL) { ctab <- summary(mer)@coefs if (!is.null(boot)) { ctab[,"Std. Error"] <- apply(boot@fixef,1,sd) ctab[,"t value"] <- ctab[,"Estimate"] / ctab[,"Std. Error"] } if (length(colnames(ctab))==4) { ctab <- cbind(ctab, ctab[,"Estimate"] + cbind(ctab[,"Std. Error"]*abs(qnorm(.025))) %*% c(-1,1)) } else { df <- mer@dims["n"]-mer@dims["q"]-mer@dims["p"]+1 ctab <- cbind(ctab,2*pt(-abs(ctab[,"t value"]),df=df)) ctab <- cbind(ctab, ctab[,"Estimate"] + cbind(ctab[,"Std. Error"]*abs(qt(.025,df=df))) %*% c(-1,1)) } colnames(ctab)[4:6]=c("p value","confsub","confsup") ctab<-as.data.frame(ctab) if (relimp) { rownames(ctab)<-sub("TRUE","",rownames(ctab)) lmformula <- paste(sapply(c(2,3),function(i) attr(mer@frame,"terms")[[i]]),collapse="~") mer.lm <- lm(lmformula,data=mer@frame) relimp <- eval(parse(text=sprintf("calc.relimp(mer.lm,type='%s',rela=TRUE)@%s",type,type))) ctab[[type]]<-NA for(name in names(relimp)) ctab[name,type]<-relimp[name] } ctab } mer2la <- function (mer=NULL,ctab=NULL,trans=NULL,relimp=TRUE,type="pmvd",strip=NULL,boot=NULL) { if (is.null(ctab)) ctab <- mer2ctab(mer,relimp=relimp,type=type,boot=boot) relimp=dim(ctab)[2]==7 if(relimp) type=colnames(ctab)[7] cc <- ctab if (!is.null(strip)) { delVec <- grep(strip,rownames(ctab)) delVec2 <- delVec[ctab[delVec,type]<.05] if (length(delVec2)>0) { ctab <- ctab[-delVec2,] if(relimp) ctab["other",type]<-1-sum(ctab[,type],na.rm=TRUE) } } colnames(ctab)[1:2]=c("$\\beta$","$\\sigma$") colnames(ctab)[4:6]=c("$p$ value","95\\% conf","interval") colnames(ctab)[colnames(ctab)=="t value"]<-"$t$" colnames(ctab)[colnames(ctab)=="z value"]<-"$z$" if (length(dim(trans))==2) { rn <- rownames(ctab) rn <- sub("TRUE$","",rn) apply(trans,1,function(x) {rn<<-sub(paste("^",x[1],"$",sep=""),x[2],rn);}) rownames(ctab)<-rn } if (relimp) print(xtable(ctab,digits=c(1,3,3,3,4,3,3,3),display=c("fg","fg","fg","fg","f","fg","fg","f"),align="l||c|c|c|c|cc|c"),floating=FALSE,sanitize.colnames.function = function(x) x,sanitize.rownames.function = function(x) x) else print(xtable(ctab,digits=c(1,3,3,3,4,3,3),display=c("fg","fg","fg","fg","f","fg","fg"),align="l||c|c|c|c|cc|"),floating=FALSE,sanitize.colnames.function = function(x) x,sanitize.rownames.function = function(x) x) cc } creaInt <- function(var) { tg4[[paste(var,"1",sep="")]]<<-tg4[[var]] * tg4[["period1"]] tg4[[paste(var,"6",sep="")]]<<-tg4[[var]] * tg4[["period6"]] tg4[[paste(var,"25",sep="")]]<<-tg4[[var]] * tg4[["period25"]] tg4[[paste(var,"T",sep="")]]<<-tg4[[var]] * tg4[["treatnum"]] tg4[[paste(var,"N",sep="")]]<<-tg4[[var]] * tg4[["north"]] tg4[[paste(var,"P",sep="")]]<<-tg4[[var]] * tg4[["period"]] } trans <- rbind(c("\\(Intercept\\)","$1$"), c("north","$\\\\north$"), c("north1","$\\\\northi$"), c("north25","$\\\\northiiv$"), c("north6","$\\\\northvi$"), c("northT","$\\\\northT$"), c("male","$\\\\male$"), c("male6","$\\\\malevi$"), c("maleN","$\\\\maleN$"), c("ownAge","$\\\\ownAge$"), c("ownAge6","$\\\\ownAgevi$"), c("ownAgeN","$\\\\ownAgeN$"), c("shownAge","$\\\\shownAge$"), c("shownAge6","$\\\\shownAgevi$"), c("shownAgeN","$\\\\shownAgeN$"), c("randomAge","$\\\\randomAge$"), c("randomAge6","$\\\\randomAgevi$"), c("randomAgeN","$\\\\randomAgeN$"), c("oneSibling","$\\\\oneSibling$"), c("oneSibling6","$\\\\oneSiblingvi$"), c("oneSiblingN","$\\\\oneSiblingN$"), c("manySiblings","$\\\\manySiblings$"), c("manySiblings6","$\\\\manySiblingsvi$"), c("manySiblingsN","$\\\\manySiblingsN$"), c("treatnum","$\\\\round$"), c("treatnumN","$\\\\roundN$"), c("treatnum6","$\\\\round$"), c("lat","$\\\\lat$"), c("lon","$\\\\lon$"), c("transfer","$t^\\\\sent$"), c("transfer6","$\\\\transferiv$"), c("transfer25","$\\\\transferiiv$"), c("transfer1","$\\\\transferi$"), c("transferT","$\\\\transferT$"), c("transferN","$\\\\transferN$"), c("transferT","$\\\\transfer_{\round}$"), c("transferN1","$\\\\transfer_{1N}$"), c("transferN25","$\\\\transfer_{2-5N}$"), c("transferN6","$\\\\transfer_{6N}$"), c("transferNT","$\\\\transfer_{TN}$"), c("retRatio","$\\\\retRatio$"), c("return","$r^\\\\received$"), c("return6","$r_6^\\\\received$"), c("returnN","$r^\\\\received_N$"), c("totalget","$t^\\\\received$"), c("totalget6","$t_6^\\\\received$"), c("totalget25","$t_{2-5}^\\\\received$"), c("totalget1","$t_1^\\\\received$"), c("totalgetT","$t_T^\\\\received$"), c("totalgetN","$t^\\\\received_N$"), c("totalgetN1","$t^\\\\received_{1N}$"), c("totalgetN25","$t^\\\\received_{2-5N}$"), c("totalgetN6","$t^\\\\received_{6N}$"), c("totalgetNT","$t^\\\\received_{TN}$"), c("period1:transfer","$t\\\\cdot d_1$"), c("northP","$P\\\\cdot d_N$"), c("period","$P$"), c("period1","$d_1$"), c("period25","$d_{2-5}$"), c("period6","$d_6$"), c("transfer:period6","$t\\\\cdot d_6$"), c("transfer:treatnum","$t\\\\cdot t$")) #-------------------------------------------------------------------------------- btest <- function (bs,expression,pos=NULL,res=c("mean"),lmer=NULL,digits=3) { results<-"" if ("expr" %in% res) results <- paste(results,expression) if ("calc" %in% res) { calc <- expression vname <- names(cc.lmer@fixef) for(n in vname[order(-sapply(vname,nchar))]) { number <- signif(cc.lmer@fixef[n],3) if (number<0) number<- paste("(",number,")") calc <- gsub(n,number,calc) } calc <- gsub("\\*","\\\\\\\\cdot",calc) calc <- gsub("\\+ *-","-",calc) calc <- gsub("\\- *-","+",calc) calc <- gsub("^ *\\+","",calc) results <- paste(results,calc,"=") } x <- with(as.data.frame(t(bs@fixef)),eval(parse(text=expression))) mmean <- mean(x) msd <- sd(x) if(!is.null(lmer)) { mmean <- with(as.data.frame(t(lmer@fixef)),eval(parse(text=expression))) msd <- sqrt(mean((x-mmean)^2)) } cinf <- mmean + qnorm(c(.025,.975)) * msd names(cinf)<-c("confint.lower","confint.upper") if ("mean" %in% res) results <- paste(results,signif(mmean,digits),sep="") if ("confint" %in% res) results <- paste(results,", \\\\mathrm{CI}_{95}=[",paste(signif(cinf,digits),collapse=","),"]",sep="") if ("pvalue" %in% res) { if(is.null(pos)) pval <- 2*pnorm(abs(mmean/msd),lower.tail=FALSE) else if(pos) pval <- pnorm(mmean/msd,lower.tail=FALSE) else pval <- pnorm(-mmean/msd,lower.tail=FALSE) pval<-if (pval<.00001) "p<0.0001" else sprintf("p=%.4f",pval) results <- paste(results,pval) } results } bootmerge <- function (b1,b2,prefix1="",prefix2="ii") { rownames(b1@fixef) <- sub("^",prefix1,rownames(b1@fixef)) rownames(b2@fixef) <- sub("^",prefix2,rownames(b2@fixef)) b1@fixef <- rbind(b1@fixef,b2@fixef) b1 } lmermerge <- function (l1,l2,prefix1="",prefix2="ii") { names(l1@fixef) <- sub("^",prefix1,names(l1@fixef)) names(l2@fixef) <- sub("^",prefix2,names(l2@fixef)) l1@fixef <- c(l1@fixef,l2@fixef) l1 } cReturn <- function (period=1,which,ratio=FALSE,res=c("mean","confint")) { Ret <- "" if ("north" %in% which) { Ret <- "(north%d + period%d + 2.5 * (treatnum + northT) + (te.period%d + 2.5 * te.treatnum)*3*((totalget%d+totalgetN%d)+2.5*(totalgetNT+totalgetT)))" if(ratio) Ret <- paste(Ret,"/3/(te.period%d + 2.5 * te.treatnum)") } if ("south" %in% which) { Ret <- "(period%d + 2.5 * treatnum + (te.period%d + 2.5 * te.treatnum)*3*(totalget%d+2.5*(totalgetT)))" if(ratio) Ret <- paste(Ret,"/3/(te.period%d + 2.5 * te.treatnum)") } if ("trans" %in% which) Ret <- "(tr.north%d + 2.5 * tr.northT) * (3 * ret.totalget%d - 1)" cc.lmer@fixef if ("return" %in% which) Ret <- paste(Ret,"- totalgetN%d * te.period%d * 3 - north%d") Out <- "" for (p in period) { if (length(period)==3 & p=="25") Out <- paste(Out,"+ 4*(",gsub("%d",p,Ret),")") else Out <- paste(Out,"+",gsub("%d",p,Ret)) } if (length(period)==3) Out <- paste("(",Out,")/6") btest(cc.boot,Out,res=res,lmer=cc.lmer) } ################################################### mix<-"(1|session)+(1|subject)" mixfix<-"as.factor(session)+as.factor(subject)" mixfix2<-"as.factor(subject)+(1|session)" ### tgns.ori <- read.dta("../NS/trustgame_ns.dta") ### we only need the first four treatments: ### write.csv(subset(tgns.ori,as.numeric(treatment)<5),file="trust.csv",row.names=FALSE) tgns<-read.csv("trust.csv") tgns <- within(tgns, lat <- round(lat_p0,2)) lon<-c(13.20,4,-2,10,5.45,26,2,9,22,-8,12.5,20,-8,-4,15) lat<-c(47.2,50.5,54.0,56.0,52.3,64.0,46.0,51.0,39.0,53.0,42.5,52.0,39.3,40.0,62.0) latlon<-as.data.frame(cbind(lon,lat)) tgns<-merge(tgns,latlon,by="lat",all.x=TRUE) sessionPar <- aggregate(tgns$subject,list(tgns$session),function(x) length(unique(x))) tg4 <- subset(tgns,as.integer(tgns$treatment)<5) tg4 <- within(tg4,{ treatnum <- as.numeric(treatment) ownRetRatio <- getback/transfer/3 treat<-as.numeric(treatment) period1=ifelse(period==1,1,0) period6=ifelse(period==6,1,0) period25=ifelse(period>1 & period<6,1,0) period1m6=period1-period6 period1p6=period1+period6 transfer6<-transfer*period6 transfer1<-transfer*period1 transfer25<-transfer*period25 transferT<-transfer*as.integer(treatment) young <- age < mean(age) old <- !young north<- dn_p0 male <- ifelse(sex=="male",1,0) ownAge <- rage shownAge <- age randomAge <- age-rage oneSibling <- ifelse(sib==1,1,0) manySiblings <- ifelse(sib>1,1,0) payoff <- with(tg4,100 - transfer + getback + totalget - totalpayback) sess<-as.numeric(session) nat=as.integer(nat_p0) profitI <- totalget + getback - transfer - totalpayback subjectid <- as.factor(subject) north <- ifelse(lat>mean(lat),1,0) totalgetN<-with(tg4,totalget*north) transferN<-with(tg4,transfer*north) }) creaInt("north") creaInt("totalget") creaInt("transfer") creaInt("totalgetN") creaInt("transferN") ################################################### ### Figure 1 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent middle <- subset(tg4,as.integer(treatment)<5,drop=TRUE) middle[["tt"]]<-as.integer(middle[["treatment"]]) xx <- with(middle,by(middle,list(period,tt), function(x) {with(x,rbind(period[1], tt[1], mean(transfer), sum(getback)/sum(transfer)/3, sum(totalpayback)/sum(totalget), mean(ownRetRatio,na.rm=TRUE), mean(transfer==0), mean(getback==0)))})) xx <- as.data.frame(t(sapply(xx,c))) colnames(xx)<-c("period","treatment","transfer","return ratio","return ratio2","avg. of indiv. return ratios","freq. of zero transfer","freq. of zero return") overSession <- function (var,legpos=NULL) { ylim <- c(min(xx[[var]]),max(xx[[var]])) plot(NULL,xlim=c(1,6),ylim=ylim,xlab="Period",ylab=var,main="") qq <- by(xx,xx[["treatment"]],function(z) lines(z[[var]] ~ z[["period"]],lty=median(z[["treatment"]]))) if(!is.null(legpos)) legend(legpos,sprintf("repetition %d",1:4),lty=1:4) } par(mfrow = c(1, 2),mar=c(4,5,0,0),mex=.63); overSession("transfer","bottom") overSession("return ratio","bottomleft") ################################################### ### code chunk number 4: transPeriodC ################################################### creaInt("north") tg4[["totalgetN"]]<-with(tg4,totalget*north) tg4[["transferN"]]<-with(tg4,transfer*north) creaInt("totalget") creaInt("transfer") creaInt("totalgetN") creaInt("transferN") te.lmer <- lmer (sprintf("transfer ~ period1 + period25 -1 + period6 + treatnum + %s",mix),data=tg4) set.seed(bseed) te.boot <- mcmcsamp(te.lmer,n=bsize) ################################################### ### code chunk number 5: returnPeriodC ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent ret16.lmer <- lmer (sprintf("totalpayback ~ period1 + period25 + period6 + treatnum + totalget1 + totalget25 + totalget6 + totalgetT + %s -1 ",mix),data=tg4) set.seed(bseed) ret16.boot <- mcmcsamp(ret16.lmer,n=bsize) ################################################### ### code chunk number 6: meanTransReturnC ################################################### middle <- subset(tg4,period>1 & period<6) #myCols <- c("ownRetRatio","transfer","ret_ratio","session","subject") avTrans <- as.data.frame(t(sapply(by(middle,middle[["subject"]],function(x) {with(x,c(sum(totalpayback)/sum(totalget),mean(ownRetRatio,na.rm=TRUE),mean(transfer),sum(getback)/sum(transfer),session[1],subject[1],dn_p0[1]))}),c))) colnames(avTrans)=c("ownRetRatio","ownRetRatio2","transfer","ret_ratio","session","subject","north") avT.lmer <- lmer(ownRetRatio ~ transfer + (1|session),data=avTrans) ownRetTransP <- sprintf("%.4f",mer2ctab(avT.lmer,relimp=FALSE)["transfer","p value"]) ################################################### ### Figure 2: ################################################### par(mar=c(4,4,0,0),mex=0.63) plot(ownRetRatio ~ transfer,data=avTrans,pch=1,ylab="average return ratio",xlab="average transfer") abline(summary(avT.lmer)@coefs[,"Estimate"],lty=2) ################################################### ### Figure 3 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent par(mar=c(4,5,0,0),mex=0.63) psd <- with(tg4,aggregate(cbind(profit,profitI,payoff,north),list(subject=subject),mean)) myECDF <- function (data) { nprof <- with(data,profitI - mean(profitI)) ecdf(nprof) } plot(myECDF(psd),verticals = TRUE, do.points = FALSE,main="",xlab="average participant profit - average total profit") abline(v=0,lty=3) ################################################### ### Figure 4 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent collapse=c("payoff","north","lat","lon","male","ownAge","randomAge","shownAge","oneSibling","manySiblings","sib","sess","subject") subjectSummary <- with(tg4,aggregate(tg4[,collapse],list(subject),mean)) subjectSummary <- within(subjectSummary, old <- ownAge>median(ownAge)) myECDF <- function (var,label) { lst<-c("solid","dotted","dotdash") plot(NULL,xlim=range(subjectSummary$payoff),ylim=c(0,1),xlab="average payoff",ylab="F(x)",main="",cex.lab=1.5) i<-1 qq <- by(subjectSummary,list(subjectSummary[[var]]), function(x) {lines(ecdf(x[["payoff"]]), lty=lst[i],do.points=FALSE,verticals=TRUE);i<<-i+1}) legend("bottomright",label,lty=lst[c((i-1):1)],bg="white",cex=1.5) } par(mfrow = c(1, 3),mar=c(4.5,5,0,0),mex=.63) myECDF("north",c("North","South")) myECDF("male",c("male","female")) subjectSummary[subjectSummary$sib>2,"sib"]<-2 myECDF("sib",c("0 siblings","1 sibling",expression(symbol(">")*1*plain(sibling)))) ################################################### ### Table 1 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent pay.lmer <- lmer (payoff ~ north + male + ownAge + oneSibling + manySiblings + (1|sess),data=subjectSummary) set.seed(bseed) pay.boot <- mcmcsamp(pay.lmer,n=bsize) xx<-sort(pay.lmer@resid) iqAge<-quantile(tg4[["ownAge"]],c(.25,.75)) %*% c(-1,1) pay.la<-mer2la(pay.lmer,trans=trans,boot=pay.boot) ################################################### ### Figure 5 ################################################### par(mex=0.63) avPayPeriod <- with(tg4,aggregate(payoff,list(north=north,period=period),mean,na.rm=TRUE)) range <- with(avPayPeriod,c(min(x),max(x))) plot(x ~ period,ylab="average payoff",data=avPayPeriod,subset=north==1,ylim=range,t="l") lines(x ~ period,data=avPayPeriod,subset=north!=1,lty=2) legend("bottomleft",c("North","South"),lty=1:2) par(mex=0.63) avPayPeriod2 <- within(reshape(avPayPeriod,timevar="north",idvar="period",direction="wide"),{period;diff <- x.1-x.0}) with(avPayPeriod2,plot(diff ~ period,ylim=c(0,max(diff)),t="l",ylab="average payoff north - south")) ################################################### ### Table 2 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent #=== transfer.lmer <- lmer ( transfer ~ period1 + period25 + period6 -1 + treatnum + north1 + north25 + north6 + northT + (1|session)+(1|subject),data=tg4) set.seed(bseed) transfer.boot <- mcmcsamp(transfer.lmer,n=bsize) transfer.lmr <- mer2la(transfer.lmer,relimp=FALSE,type="lmg",trans=trans,boot=transfer.boot) ################################################### ### Table 3 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent #=== totalpayback.lmer <- lmer ( sprintf("totalpayback ~ period1 + period25 + period6 -1 + treatnum + north1 + north25 + north6 + northT + totalget1 + totalget25 + totalget6 + totalgetT + totalgetN1 + totalgetN25 + totalgetN6 + totalgetNT + %s",mix),data=tg4) set.seed(bseed) totalpayback.boot <- mcmcsamp(totalpayback.lmer,n=bsize) cc.boot <- bootmerge(totalpayback.boot,te.boot,prefix2="te.") cc.boot <- bootmerge(cc.boot,transfer.boot,prefix2="tr.") cc.boot <- bootmerge(cc.boot,ret16.boot,prefix2="ret.") cc.lmer <- lmermerge(totalpayback.lmer,te.lmer,prefix2="te.") cc.lmer <- lmermerge(cc.lmer,transfer.lmer,prefix2="tr.") cc.lmer <- lmermerge(cc.lmer,ret16.lmer,prefix2="ret.") Avtrans1 <- btest(te.boot,"period1 + 2 * treatnum",res=c("mean"),lmer=te.lmer) Avtrans1.ci <- btest(te.boot,"period1 + 2 * treatnum",res=c("mean","confint"),lmer=te.lmer) lmr.totalpayback <- mer2la(totalpayback.lmer,relimp=FALSE,type="lmg",trans=trans,boot=totalpayback.boot) ################################################### ### Table 4 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent #=== totalget.lmer <- lmer ( sprintf("totalget ~ period1 + period25 + period6 -1 + treatnum + north1 + north25 + north6 + northT + %s",mix) ,data=tg4) set.seed(bseed) totalget.boot <- mcmcsamp(totalget.lmer,n=bsize) totalget.lmr <- mer2la(totalget.lmer,relimp=FALSE,type="lmg",trans=trans,boot=totalget.boot) ################################################### ### Table 5 ################################################### # transfer = t^\sent # getback = r^\received # totalget = t^\received # totalpayback = r^\sent getback.lmer <- lmer (sprintf("getback ~ period1 + period25 + period6 -1 + treatnum + north1 + north25 + north6 + northT + transfer1 + transfer25 + transfer6 + transferT + transferN1 + transferN25 + transferN6 + transferNT + %s",mix) ,data=tg4) set.seed(bseed) getback.boot <- mcmcsamp(getback.lmer,n=bsize) getback.lmr <- mer2la(getback.lmer,relimp=FALSE,type="lmg",trans=trans,boot=getback.boot) ################################################### ### Tables 6-9: ### The dynamic choice model has been done with Stata ### For the code please contact Fabian Bornhorst ################################################### ################################################### ### Table 10 from Online Appendix ################################################### te.cc<-mer2la(te.lmer,trans=trans,relimp=FALSE,boot=te.boot) ## d1/d6 p-value: ################################################### ### Table 11 from Online Appendix ################################################### te2.lmer <- lmer (sprintf("transfer ~ period1p6 + period1m6 + treatnum + %s",mix),data=tg4) ret.cc<-mer2la(ret16.lmer,trans=trans,relimp=FALSE,boot=ret16.boot) ################################################### ### Table 12 from Online Appendix ################################################### pay2.lmer <- lmer (payoff ~ lat + lon + male + ownAge + oneSibling + manySiblings + (1|sess),data=subjectSummary) set.seed(bseed) pay2.boot <- mcmcsamp(pay2.lmer,n=bsize) pay2.la<-mer2la(pay2.lmer,trans=trans,boot=pay2.boot) ################################################### ### Figure 10 from Online Appendix ################################################### par(mar=c(4.5,5.5,4,0),mex=0.63) ss <- subjectSummary n <- dim(ss)[1] simul <- function (x) { ss[["north"]]<-subjectSummary[order(runif(n)),"north"] xx<-lmer (payoff ~ north + male + ownAge + randomAge + oneSibling + manySiblings + (1|sess),data=ss) xx@fixef["north"] } set.seed(123) sim <- sapply(1:1000,simul) qqnorm(sim,ylab=expression(hat(beta)[d[n]])) qqline(sim) abline(h=0,v=0) psim <- sprintf("%.4f",pnorm(-(pay.lmer@fixef["north"]-mean(sim))/sd(sim))) ################################################### ### Table 13 from Online Appendix ################################################### pay.lmer2 <- lmer(sprintf("payoff ~ north + period + northP + period6 + %s",mix),data=tg4) set.seed(bseed) pay.boot2<-mcmcsamp(pay.lmer2,n=bsize) qq<-mer2la(pay.lmer2,trans=trans,boot=pay.boot2)