######Data and analyses for Carter, Wilkinson, Page (2017) Food-sharing vampire bats are more nepotistic under conditions of perceived risk. Behavioral Ecology. ##### Gerald Carter, January 5, 2017 ##### contact: gerry@socialbat.org, socialbat.org #install packages install.packages("lme4") install.packages("car") install.packages("broman") ########################### # Is nepotism stronger during food sharing in trapped condition? ########################### #Paired t-test of nepotism scores #importing data (nepotism scores, see methods) as text scores= "recipient normal trapped bella 0.459881424 0.5 cerce 0.387671233 0.5 cerceson 0.474158606 NA countess 0.195930095 NA flapper 0 NA gremlin NA 0.43800813 houdini NA 0.039563863 lily 0.41432031 NA lilyson 0.459432799 0.5 lucy 0.314011946 0.5 lucyson 0.420031881 0.45385514 marce NA 0.195 mya 0.071176069 NA pixus 0.445157152 0.5 snow NA 0.045696154 sonofzaria NA 0.5 streblick NA 0 thecount NA 0 vampirella 0.259464227 NA veronica 0.12483254 NA whimsey 0.218829457 NA zaria 0.432926829 0.5" scores= read.table(textConnection(scores), header=TRUE) #get differences scores$diffs=scores$normal-scores$trapped #t-test t.test(scores$diffs,mu=0) #results: t = -3.6632, df = 6, p-value = 0.01054 ########################### #Permutation test #load packages library(lme4) library(car) # import all possible trial donations as text (see methods) # includes possible donations from all bats that fed that recipient during the experiment temp="obs recipient donor trt trial donation total.donations donation.rate prop.donation kinship dyad 1 bella countess normal 68 11 1418 5.5 0.007757405 0 bellacountess 2 bella countess trapped 68 0 0 0 NA 0 bellacountess 3 bella gelfling trapped 56 0 179 0 0 0.1772 bellagelfling 4 bella gelfling normal 56 75 1113 37.5 0.067385445 0.1772 bellagelfling 5 bella lily normal 68 1286 1418 643 0.906911142 0.5 bellalily 6 bella lily trapped 68 0 0 0 NA 0.5 bellalily 7 bella lily normal 56 996 1113 498 0.894878706 0.5 bellalily 8 bella lily trapped 56 179 179 179 1 0.5 bellalily 9 bella mya trapped 68 0 0 0 NA 0 bellamya 10 bella mya normal 68 61 1418 30.5 0.043018336 0 bellamya 11 bella pixus trapped 56 0 179 0 0 0.5 bellapixus 12 bella pixus normal 56 42 1113 21 0.037735849 0.5 bellapixus 13 bella veronica normal 68 36 1418 18 0.02538787 0 bellaveronica 14 bella veronica trapped 68 0 0 0 NA 0 bellaveronica 15 cerce cerceson normal 52 539 730 269.5 0.738356164 0.5 cercecerceson 16 cerce cerceson trapped 52 9 9 9 1 0.5 cercecerceson 17 cerce lily normal 52 54 730 27 0.073972603 0.25 cercelily 18 cerce lily trapped 52 0 9 0 0 0.25 cercelily 19 cerce veronica normal 52 153 730 76.5 0.209589041 0 cerceveronica 20 cerce veronica trapped 52 0 9 0 0 0 cerceveronica 21 cerceson cerce trapped 64 0 0 0 NA 0.5 cercesoncerce 22 cerceson cerce normal 64 648 703 324 0.921763869 0.5 cercesoncerce 23 cerceson gomez normal 64 55 703 27.5 0.078236131 0.1697 cercesongomez 24 cerceson gomez trapped 64 0 0 0 NA 0.1697 cercesongomez 25 countess flapper normal 54 72 1175 36 0.061276596 0.5 countessflapper 26 countess flapper trapped 65 0 0 0 NA 0.5 countessflapper 27 countess flapper trapped 54 0 0 0 NA 0.5 countessflapper 28 countess flapper normal 65 205 1174 102.5 0.174616695 0.5 countessflapper 29 countess lily normal 65 158 1174 79 0.134582624 0 countesslily 30 countess lily trapped 65 0 0 0 NA 0 countesslily 31 countess lucy trapped 65 0 0 0 NA 0.2319 countesslucy 32 countess lucy normal 65 109 1174 54.5 0.092844974 0.2319 countesslucy 33 countess mina trapped 54 0 0 0 NA 0.0062 countessmina 34 countess mina normal 54 27 1175 13.5 0.022978723 0.0062 countessmina 35 countess mya normal 65 352 1174 176 0.299829642 0.0795 countessmya 36 countess mya trapped 54 0 0 0 NA 0.0795 countessmya 37 countess mya normal 54 147 1175 73.5 0.125106383 0.0795 countessmya 38 countess mya trapped 65 0 0 0 NA 0.0795 countessmya 39 countess pixus normal 54 165 1175 82.5 0.140425532 0 countesspixus 40 countess pixus trapped 54 0 0 0 NA 0 countesspixus 41 countess thecount normal 54 229 1175 114.5 0.194893617 0.3793 countessthecount 42 countess thecount trapped 54 0 0 0 NA 0.3793 countessthecount 43 countess veronica normal 65 169 1174 84.5 0.1439523 0.4081 countessveronica 44 countess veronica normal 54 247 1175 123.5 0.210212766 0.4081 countessveronica 45 countess veronica trapped 54 0 0 0 NA 0.4081 countessveronica 46 countess veronica trapped 65 0 0 0 NA 0.4081 countessveronica 47 countess zaria normal 54 37 1175 18.5 0.031489362 0 countesszaria 48 countess zaria trapped 54 0 0 0 NA 0 countesszaria 49 flapper gomez trapped 60 0 0 0 NA 0 flappergomez 50 flapper gomez normal 60 1068 1068 534 1 0 flappergomez 51 lily bella normal 70 784 933 392 0.840300107 0.5 lilybella 52 lily bella normal 59 1200 1590 600 0.754716981 0.5 lilybella 53 lily bella trapped 59 0 0 0 NA 0.5 lilybella 54 lily bella trapped 70 0 0 0 NA 0.5 lilybella 55 lily gelfling normal 59 192 1590 96 0.120754717 0.5 lilygelfling 56 lily gelfling trapped 59 0 0 0 NA 0.5 lilygelfling 57 lily pixus trapped 59 0 0 0 NA 0.25 lilypixus 58 lily pixus normal 59 198 1590 99 0.124528302 0.25 lilypixus 59 lilyson bella trapped 196 0 218 0 0 0.325 lilysonbella 60 lilyson bella normal 196 188 811 94 0.231812577 0.325 lilysonbella 61 lilyson lily normal 196 623 811 311.5 0.768187423 0.5 lilysonlily 62 lilyson lily trapped 196 218 218 218 1 0.5 lilysonlily 63 lucy bella normal 71 47 2164 23.5 0.021719039 0 lucybella 64 lucy bella trapped 71 0 0 0 NA 0 lucybella 65 lucy bella normal 57 12 2171 6 0.005527407 0 lucybella 66 lucy bella trapped 57 0 564 0 0 0 lucybella 67 lucy mina normal 57 39 2171 19.5 0.017964072 0 lucymina 68 lucy mina trapped 57 0 564 0 0 0 lucymina 69 lucy mya trapped 71 0 0 0 NA 0 lucymya 70 lucy mya normal 71 473 2164 236.5 0.21857671 0 lucymya 71 lucy mya trapped 57 0 564 0 0 0 lucymya 72 lucy mya normal 57 742 2171 371 0.341777982 0 lucymya 73 lucy pixus trapped 57 0 564 0 0 0 lucypixus 74 lucy pixus normal 57 214 2171 107 0.098572087 0 lucypixus 75 lucy veronica trapped 57 0 564 0 0 0.2752 lucyveronica 76 lucy veronica normal 57 103 2171 51.5 0.047443574 0.2752 lucyveronica 77 lucy zaria trapped 57 564 564 564 1 0.5 lucyzaria 78 lucy zaria trapped 71 0 0 0 NA 0.5 lucyzaria 79 lucy zaria normal 71 1604 2164 802 0.741219963 0.5 lucyzaria 80 lucy zaria normal 57 1061 2171 530.5 0.488714878 0.5 lucyzaria 81 lucyson bella normal 198 301 1882 150.5 0.159936238 0 lucysonbella 82 lucyson bella trapped 198 158 1712 158 0.09228972 0 lucysonbella 83 lucyson lucy normal 198 1581 1882 790.5 0.840063762 0.5 lucysonlucy 84 lucyson lucy trapped 198 1554 1712 1554 0.90771028 0.5 lucysonlucy 85 mya bella trapped 53 0 0 0 NA 0 myabella 86 mya bella normal 53 240 1847 120 0.129940444 0 myabella 87 mya bella normal 63 271 1145 135.5 0.236681223 0 myabella 88 mya bella trapped 63 0 0 0 NA 0 myabella 89 mya lily normal 63 808 1145 404 0.705676856 0.0174 myalily 90 mya lily trapped 63 0 0 0 NA 0.0174 myalily 91 mya lily normal 53 435 1847 217.5 0.235517055 0.0174 myalily 92 mya lily trapped 53 0 0 0 NA 0.0174 myalily 93 mya lucy normal 53 9 1847 4.5 0.004872767 0 myalucy 94 mya lucy trapped 53 0 0 0 NA 0 myalucy 95 mya mina trapped 63 0 0 0 NA 0.2234 myamina 96 mya mina normal 63 51 1145 25.5 0.044541485 0.2234 myamina 97 mya pixus normal 53 114 1847 57 0.061721711 0.1887 myapixus 98 mya pixus trapped 53 0 0 0 NA 0.1887 myapixus 99 mya streblick trapped 53 0 0 0 NA 0.5 myastreblick 100 mya streblick normal 53 192 1847 96 0.103952355 0.5 myastreblick 101 mya thecount trapped 53 0 0 0 NA 0.2547 myathecount 102 mya thecount normal 53 380 1847 190 0.205739036 0.2547 myathecount 103 mya veronica trapped 63 0 0 0 NA 0 myaveronica 104 mya veronica normal 53 426 1847 213 0.230644288 0 myaveronica 105 mya veronica trapped 53 0 0 0 NA 0 myaveronica 106 mya veronica normal 63 15 1145 7.5 0.013100437 0 myaveronica 107 pixus bella trapped 62 0 8 0 0 0.5 pixusbella 108 pixus bella normal 62 1388 1559 694 0.890314304 0.5 pixusbella 109 pixus lily trapped 62 8 8 8 1 0.25 pixuslily 110 pixus lily normal 62 0 1159 0 0 0.25 pixuslily 111 pixus lucy normal 62 171 1559 85.5 0.109685696 0 pixuslucy 112 pixus lucy trapped 62 0 8 0 0 0 pixuslucy 113 vampirella bella trapped 69 0 0 0 NA 0.1686 vampirellabella 114 vampirella bella normal 69 219 647 109.5 0.338485317 0.1686 vampirellabella 115 vampirella cerce normal 55 10 592 5 0.016891892 0.3623 vampirellacerce 116 vampirella cerce trapped 55 0 0 0 NA 0.3623 vampirellacerce 117 vampirella countess normal 55 128 592 64 0.216216216 0.0904 vampirellacountess 118 vampirella countess trapped 55 0 0 0 NA 0.0904 vampirellacountess 119 vampirella mya trapped 69 0 0 0 NA 0.1892 vampirellamya 120 vampirella mya normal 69 209 647 104.5 0.323029366 0.1892 vampirellamya 121 vampirella thecount normal 55 176 592 88 0.297297297 0.2114 vampirellathecount 122 vampirella thecount trapped 55 0 0 0 NA 0.2114 vampirellathecount 123 vampirella veronica trapped 55 0 0 0 NA 0.0669 vampirellaveronica 124 vampirella veronica trapped 69 0 0 0 NA 0.0669 vampirellaveronica 125 vampirella veronica normal 69 48 647 24 0.074188563 0.0669 vampirellaveronica 126 vampirella veronica normal 55 10 592 5 0.016891892 0.0669 vampirellaveronica 127 vampirella whimsey trapped 55 0 0 0 NA 0.5 vampirellawhimsey 128 vampirella whimsey normal 55 206 592 103 0.347972973 0.5 vampirellawhimsey 129 vampirella whimsey trapped 69 0 0 0 NA 0.5 vampirellawhimsey 130 vampirella whimsey normal 69 171 647 85.5 0.264296754 0.5 vampirellawhimsey 131 veronica bella trapped 61 0 0 0 NA 0 veronicabella 132 veronica bella normal 61 185 882 92.5 0.209750567 0 veronicabella 133 veronica countess trapped 61 0 0 0 NA 0.4081 veronicacountess 134 veronica countess normal 61 186 882 93 0.210884354 0.4081 veronicacountess 135 veronica mya trapped 61 0 0 0 NA 0 veronicamya 136 veronica mya normal 61 389 882 194.5 0.441043084 0 veronicamya 137 veronica pixus trapped 61 0 0 0 NA 0.4441 veronicapixus 138 veronica pixus normal 61 77 882 38.5 0.087301587 0.4441 veronicapixus 139 whimsey gomez normal 67 83 129 41.5 0.643410853 0.063 whimseygomez 140 whimsey gomez trapped 67 0 0 0 NA 0.063 whimseygomez 141 whimsey vampirella trapped 67 0 0 0 NA 0.5 whimseyvampirella 142 whimsey vampirella normal 67 46 129 23 0.356589147 0.5 whimseyvampirella 143 zaria lucy trapped 58 302 302 302 1 0.5 zarialucy 144 zaria lucy normal 58 720 984 360 0.731707317 0.5 zarialucy 145 zaria lucy trapped 66 36 36 36 1 0.5 zarialucy 146 zaria lucy normal 66 1615 1615 807.5 1 0.5 zarialucy" data= read.table(textConnection(temp), header=TRUE) #model predicts the proportion of total food from a donor as influenced by: #kinship [0-0.5], treatment=trt [trapped vs free],and their interaction #while controlling for random effects: trial and dyad #data=possible donations #data2=actual donations #extract actual donations data2=subset(data, donation>0) data2$trial=as.factor(data2$trial) #proportions are logit-transformed and have rounded minimum value subtracted to prevent 0 denominators data2$response=logit(data2$prop.donation-0.004) #fit model using observed data fit= lmer(response~kinship*trt+(1|trial/dyad), data=data2) summary(fit) #get effect size observed= fixef(fit) #create data frame of fixed effects obs=data.frame(t(as.matrix(observed))) colnames(obs)= c("intercept", "kinship","trt", "interaction") #set number of randomizations reps=5000 # create initial random data # randomly draw kinship from all possible donors # randomly assign treatment from all possible donations rand.data= transform (data, kinship=sample(kinship), trt=sample(trt)) #create empty matrix to fill with randomized fixed effects coefs.rand <- matrix(NA,nrow=reps,ncol=length(observed)) #fill matrix with results of test under null hypothesis (randomized kinship and trt) set.seed(1) for(i in 1:reps){ rand.data= transform (rand.data, kinship=sample(kinship),trt=sample(trt)) rand.data2=subset(rand.data, donation>0) rand.data2$response=logit(rand.data2$prop.donation-0.004) rand.fit= lmer(response~kinship*trt+(1|trial/dyad), data=rand.data2) rand.observed=fixef(rand.fit) coefs.rand[i,1]= rand.observed[1] coefs.rand[i,2]= rand.observed[2] coefs.rand[i,3]= rand.observed[3] coefs.rand[i,4]= rand.observed[4] } colnames(coefs.rand)= c("intercept", "kinship","trt", "interaction") expected= data.frame(coefs.rand) #print observed fixed effects obs #calculate one-tailed p-values p.kinship= sum(expected$kinship>obs$kinship)/reps #assumes positive obs value p.trt= sum (expected$trt>obs$trt)/reps #assuming positive obs value p.interaction= sum (expected$interaction>obs$interaction)/reps #assumes positive values p.values=cbind(p.kinship, p.trt, p.interaction) p.values #results #observed # intercept kinship trt(trapped) interaction #1 -2.492463 4.793516 1.775636 7.200634 # p.kinship p.trt p.interaction #[1,] 0.0148 0.0322 0.0086 #kin give a higher proportion of food received overall #kin give an even higher proportion of food received in the trapped condition ##################################### #compare total food received in trapped and normal cages ##################################### #import data as text tempdata="subject trap free1 free2 total.free bella 179 933 180 1113 bella 0 1357 61 1418 bellasson 8 1559 0 1559 cerce 9 289 441 730 cerceson 0 504 199 703 countess 0 643 532 1175 countess 0 839 335 1174 flapper 0 1068 0 1068 lily 0 948 642 1590 lily 0 916 17 933 lilyson 218 690 121 811 lucy 564 1734 437 2171 lucy 0 1357 807 2164 lucyson 1712 1882 0 1882 mya 0 744 1103 1847 mya 0 286 859 1145 vampirella 0 354 238 592 vampirella 0 476 171 647 veronica 0 682 200 882 whimsey 0 129 0 129 zaria 302 768 216 984 zaria 36 703 912 1615" totals= read.table(textConnection(tempdata), header=TRUE) #take the mean for multiple observations from one subject subjects=aggregate(cbind(trap, free1, free2)~subject,data=totals,FUN=mean, na.rm=F) sdiff=subjects$trap-subjects$free1 #permuted paired t-test library(broman) paired.perm.test(sdiff) #confidence intervals for making figure 1 fig1="subject treatment share log.share bella atrap 179 5.192956851 bella free1 933 6.839476438 bella free2 180 5.198497031 bella atrap 0 0 bella free1 1357 7.213768308 bella free2 61 4.127134385 bellasson atrap 8 2.197224577 bellasson free1 1559 7.3524411 bellasson free2 0 0 cerce atrap 9 2.302585093 cerce free1 289 5.669880923 cerce free2 441 6.091309882 cerceson atrap 0 0 cerceson free1 504 6.224558429 cerceson free2 199 5.298317367 countess atrap 0 0 countess free1 643 6.467698726 countess free2 532 6.278521424 countess atrap 0 0 countess free1 839 6.733401892 countess free2 335 5.81711116 flapper atrap 0 0 flapper free1 1068 6.974478911 flapper free2 0 0 lily atrap 0 0 lily free1 948 6.855408799 lily free2 642 6.466144724 lily atrap 0 0 lily free1 916 6.821107472 lily free2 17 2.890371758 lilyson atrap 218 5.38907173 lilyson free1 690 6.538139824 lilyson free2 121 4.804021045 lucy atrap 564 6.336825731 lucy free1 1734 7.458762692 lucy free2 437 6.08221891 lucy atrap 0 0 lucy free1 1357 7.213768308 lucy free2 807 6.694562059 lucyson atrap 1712 7.446001498 lucyson free1 1882 7.540621529 lucyson free2 0 0 mya atrap 0 0 mya free1 744 6.613384218 mya free2 1103 7.006695227 mya atrap 0 0 mya free1 286 5.659482216 mya free2 859 6.756932389 vampirella atrap 0 0 vampirella free1 354 5.872117789 vampirella free2 238 5.476463552 vampirella atrap 0 0 vampirella free1 476 6.167516491 vampirella free2 171 5.147494477 veronica atrap 0 0 veronica free1 682 6.52649486 veronica free2 200 5.303304908 whimsey atrap 0 0 whimsey free1 129 4.86753445 whimsey free2 0 0 zaria atrap 302 5.713732806 zaria free1 768 6.64509097 zaria free2 216 5.379897354 zaria atrap 36 3.610917913 zaria free1 703 6.556778356 zaria free2 912 6.816735881" fig1= read.table(textConnection(fig1), header=TRUE) #bootstrap mean within groups (trapped hour, free hour 1, free hour 2) library(boot) mean.w=function(x,w) sum(x*w) #95% CI for trapped group atrap=fig1[which(fig1$treatment=="atrap"),] trt.boot=boot(data=atrap$log.share,statistic=mean.w,R=1000, stype="w") boot.ci(trt.boot, type=c("bca")) mean(atrap$log.share) #95% CI for first free hour free1=fig1[which(fig1$treatment=="free1"),] trt.boot=boot(data=free1$log.share,statistic=mean.w,R=1000, stype="w") boot.ci(trt.boot, type=c("bca")) mean(free1$log.share) #95% CI for second free hour free2=fig1[which(fig1$treatment=="free2"),] trt.boot=boot(data=free2$log.share,statistic=mean.w,R=1000, stype="w") boot.ci(trt.boot, type=c("bca")) mean(free2$log.share) #########end