#Statistical Analyses for: Adult bonobos show no prosociality in both prosocial choice task and group service paradigm #--------------------------------------------------------------------------------------------------------------------- #Prosocial choice task #--------------------------------------------------- # Import data binom=read.table(file="binomdatapasscriterion.txt", header=T, sep="\t") attach(binom) binom$z.trial=as.vector(scale(binom$trial)) binom$session=as.factor(binom$session) binom$position=as.factor(binom$position) binom$condition=as.factor(binom$condition) contr=glmerControl(optimizer = c("bobyqa"), restart_edge = FALSE, sparseX = FALSE, check.nobs.vs.rankZ = "warningSmall", check.nobs.vs.nlev = "stop", check.nlev.gtreq.5 = "ignore", check.nlev.gtr.1 = "stop", tolPwrss = 1e-07, compDev = TRUE, optCtrl = list(maxfun=1000000000)) #Using binomial GLMM binmodel.PCT <- glmer(pros ~ condition+position+(1|actor),data = binom, family = binomial(link="logit"),control=contr) vif(binmodel.PCT) drop1(binmodel.PCT, test='Chisq') Anova(binmodel.PCT) ##### diagnostics for binmodel.PCT resmodel_PCT <- simulateResiduals(binmodel.PCT, n=250, plot=TRUE) testOutliers(resmodel_PCT, type = 'bootstrap') testUniformity(resmodel_PCT) testResiduals(resmodel_PCT) testQuantiles(resmodel_PCT, predictor = NULL, quantiles = c(0.25, 0.5, 0.75), plot = T) testZeroInflation(resmodel_PCT) #Post-hoc emm.binmodel.cond<-emmeans(binmodel.PCT, "condition",type="response",nesting = NULL) emm.binmodel.posi<-emmeans(binmodel.PCT, "position",type="response",nesting = NULL) pairs(emm.binmodel.cond) pairs(emm.binmodel.posi) confint(pairs(emm.binmodel.cond)) confint(pairs(emm.binmodel.posi)) summary(binmodel.PCT) #===Additional analyses based on the reviewer comments: including the interaction term between condition and position #using GLMM wit interaction term binmodel.PCTint <- glmer(pros ~ condition+position+condition:position+(1|actor),data = binom, family = binomial(link="logit"),control=contr) binmodel.PCTint.red1 <- glmer(pros ~ condition+position+(1|actor),data = binom, family = binomial(link="logit"),control=contr) anova(binmodel.PCTint,binmodel.PCTint.red1) drop1(binmodel.PCTint, test='Chisq') Anova(binmodel.PCTint, test='Chisq') #post-hoc emm_cond_pos3 <- emmeans(binmodel.PCTint, ~position*condition, type="response") contrast(emm_cond_pos3, "pairwise", by="condition", adjust="tukey") confint(binmodel.PCTint, by = "condition") emm1.1 = emmeans(binmodel.PCTint, specs = pairwise ~ condition:position, type = "response", adjust = "none") emm1.1 emm1.1$emmeans %>% as.data.frame() emm2 = emmeans(binmodel.test3.3, specs = pairwise ~ position|condition, type = "response") emm2 emm2$contrasts %>% as.data.frame() #Using wilcoxon signed rank: #1) does the poportion of prosocial choices in test sessions differ from control sessions testcontrpros=read.table(file="testcontrproslong.txt", header=T, sep="\t") wilcox.pros.testcontr<-wilcox.test(prop~condition,data=testcontrpros,paired=TRUE,exact=FALSE) wilcox.pros.testcontr Zstatwilcox.pros.testcontr<-qnorm(wilcox.pros.testcontr$p.value/2) Zstatwilcox.pros.testcontr abs(Zstatwilcox.pros.testcontr)/sqrt(7) prostestcontrdata <- ddply(testcontrpros, c("condition"), summarise, N = length(prop), median = median(prop), mean = mean(prop), sd = sd(prop), se = sd / sqrt(7)) prostestcontrdata #2) In the test sessions, does the proportion of prosocial choices and selfish choices differ? testdata=read.table(file="testpasswilcox.txt", header=T, sep="\t") wilcox.test(prop~choice,data=testdata,paired=TRUE,exact=FALSE) Zstatwilcox.test.test <-qnorm(wilcox.test.test $p.value/2) Zstatwilcox.test.test abs(Zstatwilcox.test.test )/sqrt(7) #7 individuen/paren tdata <- ddply(testdata, c("choice"), summarise, N = length(prop), median = median(prop), mean = mean(prop), sd = sd(prop), se = sd / sqrt(7)) tdata #3) In the control sessions, does the proportion of prosocial choices and selfish choices differ? controldata=read.table(file="contrpasswilcox.txt", header=T, sep="\t") wilcox.test.control<-wilcox.test(prop~choice,data=controldata,paired=TRUE,exact=FALSE) wilcox.test.control Zstatwilcox.test.control <-qnorm(wilcox.test.control$p.value/2) Zstatwilcox.test.control abs(Zstatwilcox.test.control)/sqrt(7) #7 individuen/paren cdata <- ddply(controldata, c("choice"), summarise, N = length(prop), median = median(prop), mean = mean(prop), sd = sd(prop), se = sd / sqrt(7)) #---------------------------------------------------- #Binomial GLMM onlt including first '20' test trials & all control trials binomfirst20=read.table(file="binomfirst20datapasscriterion.txt", header=T, sep="\t") attach(binomfirst20) binmodelfirst20<- glmer(pros ~ condition+position+condition:position+(1|actor),data = binomfirst20, family = binomial(link="logit"),control=contr) vif(binmodelfirst20) drop1(binmodelfirst20, test='Chisq') Anova(binmodelfirst20) emm_cond_posfirst20 <- emmeans(binmodelfirst20.test1, ~position*condition, type="response") emm_cond_posfirst20 contrast(emm_cond_posfirst20, "pairwise", by="condition", adjust="tukey") #Using wilcoxon signed rank: #1) does the poportion of prosocial choices in first 20 test trials differ from control trials test20contrpros=read.table(file="testfirst20vscontr.txt", header=T, sep="\t") wilcox.pros.test20contr<-wilcox.test(prop~condition,data=test20contrpros,paired=TRUE,exact=FALSE) wilcox.pros.test20contr Zstatwilcox.pros.test20contr<-qnorm(wilcox.pros.test20contr$p.value/2) Zstatwilcox.pros.test20contr abs(Zstatwilcox.pros.test20contr)/sqrt(7) prostest20 <- ddply(test20contrpros, c("condition"), summarise, N = length(prop), median = median(prop), mean = mean(prop), sd = sd(prop), se = sd / sqrt(7)) #Individual analyses (but see comment of reviewer 2) library(multcomp) djpos <- subset(binom, actor=='DJ') bspos <- subset(binom, actor=='BS') hbpos <- subset(binom, actor=='HB') kgpos <- subset(binom, actor=='KG') kkpos <- subset(binom, actor=='KK') mzpos <- subset(binom, actor=='MZ') nypos <- subset(binom, actor=='NY') djmodel <- glm(pros ~ condition+position, data = djpos, family = binomial(link="logit")) drop1(djmodel, test='Chisq') Anova(djmodel) emm.djmodel.cond<-emmeans(djmodel, "condition",type="response") emm.djmodel.pos<-emmeans(djmodel, "position",type="response") pairs(emm.djmodel.cond) pairs(emm.djmodel.pos) confint(pairs(emm.djmodel.cond)) confint(pairs(emm.djmodel.pos)) summary(djmodel) hbmodel <- glm(pros ~ condition+position, data = hbpos, family = binomial(link="logit")) drop1(hbmodel, test='Chisq') Anova(hbmodel) emm.hbmodel.cond<-emmeans(hbmodel, "condition",type="response") emm.hbmodel.pos<-emmeans(hbmodel, "position",type="response") pairs(emm.hbmodel.cond) pairs(emm.hbmodel.pos) confint(pairs(emm.hbmodel.cond)) confint(pairs(emm.hbmodel.pos)) summary(hbmodel) bsmodel <- glm(pros ~ condition+position, data = bspos, family = binomial(link="logit")) drop1(bsmodel, test='Chisq') # significant Anova(bsmodel) emm.bsmodel.cond<-emmeans(bsmodel, "condition",type="response") emm.bsmodel.pos<-emmeans(bsmodel, "position",type="response") pairs(emm.bsmodel.cond) pairs(emm.bsmodel.pos) confint(pairs(emm.bsmodel.cond)) confint(pairs(emm.bsmodel.pos)) summary(bsmodel) kgmodel <- glm(pros ~ condition+position, data = kgpos, family = binomial(link="logit")) drop1(kgmodel, test='Chisq') Anova(kgmodel) emm.kgmodel.cond<-emmeans(kgmodel, "condition",type="response") emm.kgmodel.pos<-emmeans(kgmodel, "position",type="response") pairs(emm.kgmodel.cond) pairs(emm.kgmodel.pos) confint(pairs(emm.kgmodel.cond)) confint(pairs(emm.kgmodel.pos)) summary(kgmodel) kkmodel <- glm(pros ~ condition+position, data = kkpos, family = binomial(link="logit")) drop1(kkmodel, test='Chisq') Anova(kkmodel) emm.kkmodel.cond<-emmeans(kkmodel, "condition",type="response") emm.kkmodel.pos<-emmeans(kkmodel, "position",type="response") pairs(emm.kkmodel.cond) pairs(emm.kkmodel.pos) confint(pairs(emm.kkmodel.cond)) confint(pairs(emm.kkmodel.pos)) summary(kkmodel) mzmodel <- glm(pros ~ condition+position, data = mzpos, family = binomial(link="logit")) drop1(mzmodel, test='Chisq') Anova(mzmodel) emm.mzmodel.cond<-emmeans(mzmodel, "condition",type="response") emm.mzmodel.pos<-emmeans(mzmodel, "position",type="response") pairs(emm.mzmodel.cond) pairs(emm.mzmodel.pos) confint(pairs(emm.mzmodel.cond)) confint(pairs(emm.mzmodel.pos)) summary(mzmodel) nymodel <- glm(pros ~ condition+position, data = nypos, family = binomial(link="logit")) drop1(nymodel, test='Chisq') Anova(nymodel) emm.nymodel.cond<-emmeans(nymodel, "condition",type="response") emm.nymodel.pos<-emmeans(nymodel, "position",type="response") pairs(emm.nymodel.cond) confint(pairs(emm.nymodel.cond)) pairs(emm.nymodel.pos) confint(pairs(emm.nymodel.pos)) summary(nymodel) #--------------------------------------------------- #Group service paradigm #--------------------------------------------------- #Import dataset totalpull=read.table(file="totalpulls.txt", header=T, sep="\t") #Does proportion of pulls differ between test and control sessions? wilcox.test.proptotal<-wilcox.test(prop~type,data=totalpull,paired=TRUE,exact=FALSE) wilcox.test.proptotal Zstatwilcox.test.proptotal <-qnorm(wilcox.test.proptotal$p.value/2) Zstatwilcox.test.proptotal abs(Zstatwilcox.test.proptotal)/sqrt(7) groupservicetotal <- ddply(totalpull,c("type"), summarise,N=length(prop),median = median(prop),sd= sd(prop),se= sd / sqrt(7)) #7 individuen groupservicetotal