LBH coor singing analysis oct-15
Marcelo Araya Salas
10/24/2015
Simulated data
Prob of overlap simulated data
## Prob for overlapping bouts## [1] 0.6975813## SD Prob for overlapping bouts## [1] 0.04527745## Prob for alternating bouts## [1] 0.3636657## SD Prob for alternating bouts## [1] 0.04916876## Prob for non-coordinated bouts## [1] 0.4732574## SDProb for non-coordinated bouts## [1] 0.06136546## classification performance## pvalclass
## <0.025 >0.025 & <0.975 <0.975
## no.ovlp 477 523 0
## notsync 22 950 28
## ovlp 0 16 984## overall classification performance## [1] 0.8036667## [[1]]
##
## [[2]]
##
## [[3]]
Descriptive stats
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 04 00 00 00 00 00 00 0d 3b 01 00 02 00 02 00 00 00 05 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 04 00 00 00 00 00 00 0d 3b 01 00 02 00 02 00 00 00 05 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 04 00 00 00 00 00 00 0d 3b 01 00 02 00 02 00 00 00 05 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 04 00 00 00 00 00 00 0d 3b 01 00 02 00 02 00 00 00 05 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 04 00 00 00 00 00 00 0d 3b 01 00 02 00 02 00 00 00 05 00## number of events## [1] 185## Start and end date## [1] "2012-02-01 UTC" "2012-05-19 UTC"## Number of recordings## [1] 9668## Time periods## [1] "0530" "0600" "0630" "0700" "0730" "0800" "0830" "0900" "0930" "1000"
## [11] "1030" "1100" "1130" "1200" "1230" "1300" "1330" "1400" "1430" "1500"
## [21] "1530" "1600" "1630" "1700" "1730" "0500"## Number of time periods## [1] 26## Number of days recorded## [1] 79## Number of days recorded per lek## Group.1 x
## 1 CCE 12
## 2 CCL 12
## 3 LOC 6
## 4 SUR 49## average distance between mics## [1] 14.93317## Number of days analyzed per lek## Group.1 x
## 1 CCE 5
## 2 CCL 3
## 3 LOC 2
## 4 SUR 7## Mean and sd perch height## [1] 2.150768## [1] 1.130179## Number of pairs## [1] 19## Number of individuals in pairs## [1] 22## NUmber of same individual i different pairs##
## -1 100 103 104 106 108 109 111 117 143 148 149 24 32 54 60 72 81
## 1 2 1 1 2 3 1 2 2 1 3 1 1 4 2 3 1 1
## 86 88 90 98
## 3 1 1 1## Mean number of times the same individual was in different pairs## [1] 1.727273## Number of events with only one male known## [1] 35## Number of males in events with only one mae knwow## [1] 21## Number of identified males in all events## [1] 42## mean duration of events## [1] 66.99738## duration range of events## [1] 5.103 438.772## mean number of songs## [1] 220.3351## range number of songs## [1] 24 1315## counts of events per pattern##
## No.ovlp No.pattern Ovlp
## 43 104 38## counts of events per pattern for known dyads##
## No.ovlp No.pattern Ovlp
## 29 65 27## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00## [1] 88 -1 153 146 36 95 149 147 151 83 148 86 143 32 152 109 111
## [18] 117 108 106 100 104 103 65 98 72 97 81 43 24 90 53 9 54
## [35] 60## proportion of dyads with perches at both distance categories## close far
## 0.3455357 0.8782738## min proportion of dyads with perches at both distance categories## close far
## 0.03571429 0.57142857## max proportion of dyads with perches at both distance categories## close far
## 0.5714286 1.0000000## mean distance between territory centroids## [1] 39.18379## sd distance between territory centroids## [1] 26.46792## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 02 00 00 00 00 00 00 0d 3b 00 00 00 00 52 01 00 00 15 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 01 00 00 00 00 00 00 0d 3b 01 00 00 00 ac 01 00 00 25 00## number of possible pairs## [1] 251## number of observed singing pairs## [1] 19## proportion## [1] 0.07569721Real data
## cutoff for distance## 0% 25% 50% 75% 100%
## 0.00000 13.39981 21.50203 26.96696 75.57736##
## 103//104 108//109 108//111 108//117 111//117 32//148 32//54 32//60
## Close 1 14 1 1 0 1 0 0
## Far 1 0 0 0 2 0 3 4
##
## 32//86 54//60 60//81 98//106
## Close 0 0 0 0
## Far 2 1 8 5
## Singing pattern vs song type sharing ALL DATA
contextual factors effect ALL DATA
cat("sample size")## sample sizenrow(stres)## [1] 61Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+ Time+songtype+distcat, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time+songtype+distcat, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time+songtype, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+songtype, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+distcat, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~songtype+distcat, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~songtype, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~distcat, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur, random=~1 |lek , data = stres, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~1, random=~1 |lek , data = stres, method = "ML")
# Cand.models[[1]]<-glm(sing.perf~Time+songtype+distcat, data = stres)
# Cand.models[[2]]<-glm(sing.perf~Time+songtype, data = stres)
# Cand.models[[3]]<-glm(sing.perf~dur+songtype+distcat, data = stres)
# Cand.models[[4]]<-glm(sing.perf~dur+Time, data = stres)
# Cand.models[[5]]<-glm(sing.perf~dur, data = stres)
# Cand.models[[6]]<-glm(sing.perf~dur+songtype, data = stres)
# Cand.models[[7]]<-glm(sing.perf~Time, data = stres)
# Cand.models[[8]]<-glm(sing.perf~Time+distcat, data = stres)
# Cand.models[[9]]<-glm(sing.perf~songtype+distcat, data = stres)
# Cand.models[[10]]<-glm(sing.perf~songtype, data = stres)
# Cand.models[[11]]<-glm(sing.perf~distcat, data = stres)
# Cand.models[[12]]<-glm(sing.perf~1, data = stres)
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~sing.perf","sing.perf~",Modnames)
aictabs<-list()
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt
## sing.perf~dur + distcat 5 36.51 0.00 0.58 0.58
## sing.perf~distcat 4 38.64 2.13 0.20 0.78
## sing.perf~songtype + distcat 5 40.01 3.50 0.10 0.88
## sing.perf~dur + Time + songtype + distcat 7 40.45 3.94 0.08 0.96
## sing.perf~Time + songtype + distcat 6 41.71 5.21 0.04 1.00
## sing.perf~dur + songtype 5 63.71 27.20 0.00 1.00
## sing.perf~dur 4 64.96 28.45 0.00 1.00
## sing.perf~1 3 65.91 29.40 0.00 1.00
## sing.perf~songtype 4 66.48 29.98 0.00 1.00
## sing.perf~Time 4 66.68 30.17 0.00 1.00
## sing.perf~Time + songtype 5 67.44 30.93 0.00 1.00
## LL
## sing.perf~dur + distcat -12.71
## sing.perf~distcat -14.96
## sing.perf~songtype + distcat -14.46
## sing.perf~dur + Time + songtype + distcat -12.17
## sing.perf~Time + songtype + distcat -14.08
## sing.perf~dur + songtype -26.31
## sing.perf~dur -28.12
## sing.perf~1 -29.74
## sing.perf~songtype -28.88
## sing.perf~Time -28.98
## sing.perf~Time + songtype -28.17cs<-list()
cs[[length(cs)+1]]<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## sing.perf~dur + distcat 5 36.51 0.00 0.58
## sing.perf~distcat 4 38.64 2.13 0.20
## sing.perf~songtype + distcat 5 40.01 3.50 0.10
## sing.perf~dur + Time + songtype + distcat 7 40.45 3.94 0.08
##
## Model probabilities sum to 0.96confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
#effect size for song type
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Different"), distcat = c("Far", "Far"),Time = c(0,0), dur = c(1,1)), second.ord = T)
b<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Different"), distcat = c("Close", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
es<-list()
es[[length(es)+1]]<-data.frame(data="All data",factor="coord perf~song type",effect_size=mean(a$Mod.avg.eff,b$Mod.avg.eff), lowerCI=mean(a$Lower.CL,b$Lower.CL), upperCI=mean(a$Upper.CL, b$Upper.CL))
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 All data coord perf~song type 0.01283926 -0.2020502 0.2277288#effect size for distance
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Same"), distcat = c("Far", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
b<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Different", "Different"), distcat = c("Far", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
es[[length(es)+1]]<-data.frame(data="All data",factor="coord perf~distance",effect_size=mean(a$Mod.avg.eff,b$Mod.avg.eff), lowerCI=mean(a$Lower.CL,b$Lower.CL), upperCI=mean(a$Upper.CL, b$Upper.CL))
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 All data coord perf~distance 0.5563044 0.3214828 0.7911261d<-modavg(parm = "Time", cand.set = confmodels, modnames = confmodnames)
es[[length(es)+1]]<-data.frame(data="All data",factor="coord perf~time",effect_size=d$Mod.avg.beta,lowerCI= d$Lower.CL , upperCI=d$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 All data coord perf~time 0.01307057 -0.01367228 0.03981342# modavg(parm = "dur", cand.set = confmodels, modnames = confmodnames)
dur.models<-confmodels[grep("dur", confmodnames, invert = F)]
d<-summary(dur.models[[1]])$coefficients$fixed[2]
d1<-confint(dur.models[[1]], parm="dur", level=0.95)[[1]][1:2]
es[[length(es)+1]]<-data.frame(data="All data",factor="coord perf~dur",effect_size=d,lowerCI= d1[1], upperCI=d1[2])
es[[length(es)]]## data factor effect_size lowerCI upperCI
## dur All data coord perf~dur -0.001524101 -0.003028919 -1.928248e-05contextual factors effect known dyads
stresKP <- droplevels(stres[grep("//", stres$pair.x),])
cat("bout sample size")## bout sample sizenrow(stresKP)## [1] 34cat("pair sample size")## pair sample sizelength(unique(stresKP$pair.x))## [1] 12Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+ Time+songtype+distcat, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time+songtype+distcat, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time+songtype, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+songtype, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur+distcat, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~songtype+distcat, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~Time, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~songtype, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~distcat, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~dur, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(sing.perf~1, random=~1 |lek/pair.x , data = stresKP, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~sing.perf","sing.perf~",Modnames)
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt
## sing.perf~distcat 5 2.87 0.00 0.43 0.43
## sing.perf~songtype + distcat 6 3.72 0.85 0.28 0.71
## sing.perf~dur + distcat 6 4.28 1.41 0.21 0.92
## sing.perf~Time + songtype + distcat 7 6.88 4.01 0.06 0.98
## sing.perf~dur + Time + songtype + distcat 8 9.50 6.63 0.02 0.99
## sing.perf~1 4 12.70 9.83 0.00 1.00
## sing.perf~dur 5 14.17 11.30 0.00 1.00
## sing.perf~songtype 5 15.37 12.50 0.00 1.00
## sing.perf~Time 5 15.47 12.60 0.00 1.00
## sing.perf~dur + songtype 6 16.75 13.88 0.00 1.00
## sing.perf~Time + songtype 6 18.33 15.47 0.00 1.00
## LL
## sing.perf~distcat 4.64
## sing.perf~songtype + distcat 5.69
## sing.perf~dur + distcat 5.42
## sing.perf~Time + songtype + distcat 5.71
## sing.perf~dur + Time + songtype + distcat 6.13
## sing.perf~1 -1.66
## sing.perf~dur -1.01
## sing.perf~songtype -1.61
## sing.perf~Time -1.66
## sing.perf~dur + songtype -0.82
## sing.perf~Time + songtype -1.61cs[[length(cs)+1]]<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## sing.perf~distcat 5 2.87 0.00 0.43
## sing.perf~songtype + distcat 6 3.72 0.85 0.28
## sing.perf~dur + distcat 6 4.28 1.41 0.21
## sing.perf~Time + songtype + distcat 7 6.88 4.01 0.06
##
## Model probabilities sum to 0.98confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
#effect size for song type
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Different"), distcat = c("Far", "Far"),Time = c(0,0), dur = c(1,1)), second.ord = T)
b<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Different"), distcat = c("Close", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
es[[length(es)+1]]<-data.frame(data="Known dyads",factor="coord perf~song type",effect_size=mean(a$Mod.avg.eff,b$Mod.avg.eff), lowerCI=mean(a$Lower.CL,b$Lower.CL), upperCI=mean(a$Upper.CL, b$Upper.CL))
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 Known dyads coord perf~song type 0.05160043 -0.1819734 0.2851743#effect size for distance
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Same", "Same"), distcat = c("Far", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
b<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(songtype = c("Different", "Different"), distcat = c("Far", "Close"),Time = c(0,0), dur = c(1,1)), second.ord = T)
es[[length(es)+1]]<-data.frame(data="Known dyads",factor="coord perf~distance",effect_size=mean(a$Mod.avg.eff,b$Mod.avg.eff), lowerCI=mean(a$Lower.CL,b$Lower.CL), upperCI=mean(a$Upper.CL, b$Upper.CL))
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 Known dyads coord perf~distance 0.5081319 0.2714378 0.7448259d<-modavg(parm = "Time", cand.set = confmodels, modnames = confmodnames)
es[[length(es)+1]]<-data.frame(data="Known dyads",factor="coord perf~time",effect_size=d$Mod.avg.beta,lowerCI= d$Lower.CL , upperCI=d$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI upperCI
## 1 Known dyads coord perf~time 0.002373733 -0.02205481 0.02680228# modavg(parm = "dur", cand.set = confmodels, modnames = confmodnames)
dur.models<-confmodels[grep("dur", confmodnames, invert = F)]
d<-summary(dur.models[[1]])$coefficients$fixed[2]
d1<-confint(dur.models[[1]], parm="dur", level=0.95)[[1]][1:2]
es[[length(es)+1]]<-data.frame(data="Known dyads",factor="coord perf~dur",effect_size=d,lowerCI= d1[1], upperCI=d1[2])
es[[length(es)]]## data factor effect_size lowerCI upperCI
## dur Known dyads coord perf~dur -0.000837914 -0.002136698 0.0004608705cat("pattern by distance ALL DATA")## pattern by distance ALL DATAtable(stres$distcat, stres$pattern)[-2,]##
## No.ovlp No.pattern Ovlp
## Close 9 7 0
## Far 3 31 11cat("pattern by distance known dyads")## pattern by distance known dyadstable(stresKP$distcat, stresKP$pattern)##
## No.ovlp No.pattern Ovlp
## Close 4 4 0
## Far 2 19 5cat("distance by song type sharing ALL DATA")## distance by song type sharing ALL DATAtable(stres$distcat, stres$songtype)[-2,]##
## Different Same
## Close 0 16
## Far 20 25cat("distance by song type sharing known dyads")## distance by song type sharing known dyadstable(stresKP$distcat, stresKP$songtype)##
## Different Same
## Close 0 8
## Far 16 10cat("pattern by song type sharing known dyads")## pattern by song type sharing known dyadstable(stres$pattern, stres$songtype)##
## Different Same
## No.ovlp 2 10
## No.pattern 14 24
## Ovlp 4 7cat("pattern by song type sharing known dyads")## pattern by song type sharing known dyadstable(stresKP$pattern, stresKP$songtype)[]##
## Different Same
## No.ovlp 2 4
## No.pattern 11 12
## Ovlp 3 2ssresults <- droplevels(ssresults)
a <- table(ssresults$pair, ssresults$pattern)
cat("proportion of dyads showing a coordination pattern")## proportion of dyads showing a coordination pattern1- length(which(apply(a[,-2],1,sum)/apply(a,1,sum) == 0))/nrow(a)## [1] 0.7368421cat("proportion of dyads showing alternation")## proportion of dyads showing alternationlength(which(a[,1] > 0))/nrow(a)## [1] 0.4736842cat("proportion of dyads showing overlapping")## proportion of dyads showing overlappinglength(which(a[,3] > 0))/nrow(a)## [1] 0.5263158b <- a[apply(a,1,sum)>2,]
cat("pattern by song type sharing known dyads")## pattern by song type sharing known dyads1- length(which(apply(b[,-2],1,sum)/apply(b,1,sum) == 0))/nrow(b)## [1] 0.9230769cat("mean bouts per dyad")## mean bouts per dyadmean(table(ssresults$pair))## [1] 6.368421cat("sd bouts per dyad")## sd bouts per dyadsd(table(ssresults$pair))## [1] 6.112015# altchi <- a[,2:1]
# altchi[,1] <-apply(a,1,sum) * 0.025
#
# chisq.test(altchi)
#
#
# ovlchi <- a[,3:2]
# ovlchi[,1] <-apply(a,1,sum) * 0.025
#
# chisq.test(ovlchi)
bothchi <- a[,3:2]
bothchi[,1] <- apply(a[,-2],1,sum)
cat("chi square for occurrence of coordination in dyads")## chi square for occurrence of coordination in dyadschisq.test(bothchi, simulate.p.value = F) ##
## Pearson's Chi-squared test
##
## data: bothchi
## X-squared = 42.246, df = 18, p-value = 0.001022cat("proportions of songs overlapped")## proportions of songs overlappede <- sapply(1:nrow(results),function(x)
{results$obs.ovlps[x]/results$n[x]*2
})
tapply(e, INDEX = results$pattern, mean)## No.ovlp No.pattern Ovlp
## 0.1381852 0.4112937 0.5732274cat("proportion of unidentified birds")## proportion of unidentified birdsd <- table(results$class)[-1]
(sum(d[1:2]) + (sum(d[3:5])*2))/(nrow(results)*2)## [1] 0.2486486cat("proportion of bouts with at least 1 unidentified bird")## proportion of bouts with at least 1 unidentified bird1 - 121/185## [1] 0.3459459Gap duration
gap duration for ALL DATA
me<-aggregate(gaps[,grep("gap|indivsolo",colnames(gaps))],by=list(gaps$sound.files,gaps$pattern),mean)
me<-me[!duplicated(me),]
me1<-data.frame(gaptype="solo gap",me[,c(1:4)])
me2<-data.frame(gaptype="duet gap",me[,c(1:2,5,4)])
colnames(me1)<-c("gaptype", "sound.files", "pattern", "gap", "indiv")
colnames(me2)<-c("gaptype", "sound.files", "pattern", "gap", "indiv")
me<-rbind(me1, me2)
me<-merge(me, results[,grep("sing.event|pair",colnames(results))], by.x= "sound.files", by.y= "sing.event")
me<-merge(me,songtype1[,grep("songtype|pair", colnames(songtype1))], by= "pair")
me<-me[!duplicated(me),]
cat("sample size")## sample sizenrow(me)## [1] 132Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype * pattern, random=~1 | sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~pattern, random=~1 | sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype , random=~ 1 | sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~1 , random=~1 | sound.files, data = me, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~gap","gap duration~",Modnames)
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## gap duration~gaptype * pattern 8 -384.47 0.00 1 1 200.82
## gap duration~gaptype 4 -366.55 17.92 0 1 187.43
## gap duration~pattern 5 -328.55 55.92 0 1 169.51
## gap duration~1 3 -323.42 61.05 0 1 164.80cs[[length(cs)+1]]<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## gap duration~gaptype * pattern 8 -384.47 0 1
##
## Model probabilities sum to 1confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~gaptype+NO.OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.dur~gaptype+NO.OVLP 0.1044533 0.06514806
## upperCI
## Mod.avg.diff 0.1437586a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.pattern", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~gaptype+NO.PATT",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.dur~gaptype+NO.PATT 0.0399524 0.01147333
## upperCI
## Mod.avg.diff 0.06843146a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("Ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~gaptype+OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.dur~gaptype+OVLP 0.04299138 -0.003259029
## upperCI
## Mod.avg.diff 0.0892418a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("Ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~Ovlp.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~Ovlp.VS.No.ovlp+SOLOGAP -0.01310553
## lowerCI upperCI
## Mod.avg.diff -0.0560241 0.02981303a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~No.pattern.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~No.pattern.VS.No.ovlp+SOLOGAP -0.01682756
## lowerCI upperCI
## Mod.avg.diff -0.0511493 0.01749417a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~No.pattern.VS.Ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~No.pattern.VS.Ovlp+SOLOGAP -0.003722032
## lowerCI upperCI
## Mod.avg.diff -0.04212879 0.03468472a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~Ovlp.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~Ovlp.VS.No.ovlp+DUETGAP 0.0745675
## lowerCI upperCI
## Mod.avg.diff 0.03164893 0.1174861a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~No.pattern.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~No.pattern.VS.No.ovlp+DUETGAP 0.08132852
## lowerCI upperCI
## Mod.avg.diff 0.04700678 0.1156503a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="All data",factor="gap.dur~No.pattern.VS.Ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.dur~No.pattern.VS.Ovlp+DUETGAP -0.00676102
## lowerCI upperCI
## Mod.avg.diff -0.04516778 0.03164574gap duration for known dyads
# by pattern for Known individuals
me<-me[!is.na(me$indiv),]
cat("sample size number of bouts")## sample size number of boutsnrow(me)## [1] 102cat("sample size number of dyads")## sample size number of dyadslength(unique(me$indiv))## [1] 12Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype * pattern, random=~1 | indiv/sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~pattern, random=~1 | indiv/sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype , random=~1 | indiv/sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~1 , random=~1 | indiv/sound.files, data = me, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~gap","gap duration~",Modnames)
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## gap duration~gaptype * pattern 9 -296.13 0.00 0.98 0.98 158.04
## gap duration~gaptype 5 -288.24 7.89 0.02 1.00 149.43
## gap duration~pattern 6 -254.60 41.53 0.00 1.00 133.74
## gap duration~1 4 -251.08 45.05 0.00 1.00 129.75cs[[length(cs)+1]]<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## gap duration~gaptype * pattern 9 -296.13 0 0.98
##
## Model probabilities sum to 0.98confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~gaptype+NO.OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.dur~gaptype+NO.OVLP 0.09268589 0.04163014
## upperCI
## Mod.avg.diff 0.1437416a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.pattern", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~gaptype+NO.PATT",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.dur~gaptype+NO.PATT 0.04677359 0.008824618
## upperCI
## Mod.avg.diff 0.08472257a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("Ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~gaptype+OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.dur~gaptype+OVLP 0.03751271 -0.01742934
## upperCI
## Mod.avg.diff 0.09245476a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("Ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~Ovlp.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~Ovlp.VS.No.ovlp+SOLOGAP -0.02746772
## lowerCI upperCI
## Mod.avg.diff -0.08050223 0.02556679a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~No.pattern.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~No.pattern.VS.No.ovlp+SOLOGAP -0.03334329
## lowerCI upperCI
## Mod.avg.diff -0.07832559 0.01163901a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~Ovlp.VS.No.pattern+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~Ovlp.VS.No.pattern+SOLOGAP -0.005875568
## lowerCI upperCI
## Mod.avg.diff -0.05309184 0.04134071a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~Ovlp.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~Ovlp.VS.No.ovlp+DUETGAP 0.0826409
## lowerCI upperCI
## Mod.avg.diff 0.02960639 0.1356754a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~No.pattern.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~No.pattern.VS.No.ovlp+DUETGAP 0.07925559
## lowerCI upperCI
## Mod.avg.diff 0.03427329 0.1242379a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.dur~Ovlp.VS.No.pattern+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.dur~Ovlp.VS.No.pattern+DUETGAP 0.003385315
## lowerCI upperCI
## Mod.avg.diff -0.04383096 0.05060159variation in gap duration for ALL DATA
cv.gaps <- readRDS("~/Dropbox/LBH data/Vocalizations/Coordinated singing/cv.gaps.RDS")
me.cv <- reshape(cv.gaps[,grep("sound.files|gap|pattern|indivsolo", colnames(cv.gaps))],idvar = c("sound.files"),ids = cv.gaps$sound.files, times=names(cv.gaps)[c(2,4)],timevar = "gaptype", varying = list(names(cv.gaps)[c(2,4)]), direction = "long", v.names = "gap")
me.cv$gaptype <- as.factor(me.cv$gaptype)
cat("sample size")## sample sizenrow(me.cv)## [1] 134Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype * pattern, random=~1 | sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~pattern, random=~1 | sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype , random=~ 1 | sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~1 , random=~1 | sound.files, data = me.cv, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~gap","gap variation~",Modnames)
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## gap variation~gaptype 4 -381.35 0.00 0.9 0.9 194.83
## gap variation~gaptype * pattern 8 -376.92 4.43 0.1 1.0 197.04
## gap variation~1 3 -358.34 23.02 0.0 1.0 182.26
## gap variation~pattern 5 -354.06 27.29 0.0 1.0 182.27cs[[length(cs)+1]] <- confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## gap variation~gaptype 4 -381.35 0.00 0.9
## gap variation~gaptype * pattern 8 -376.92 4.43 0.1
##
## Model probabilities sum to 1confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("No.ovlp", "No.ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~gaptype+NO.OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.cv~gaptype+NO.OVLP 0.04472304 0.0180233
## upperCI
## Mod.avg.diff 0.07142279a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("No.pattern", "No.pattern")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~gaptype+NO.PATT",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.cv~gaptype+NO.PATT 0.04773163 0.02733494
## upperCI
## Mod.avg.diff 0.06812833a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("Ovlp", "Ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~gaptype+OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff All data gap.cv~gaptype+OVLP 0.04970973 0.02263569
## upperCI
## Mod.avg.diff 0.07678376a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("Ovlp", "No.ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~Ovlp.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~Ovlp.VS.No.ovlp+SOLOGAP -0.002666117
## lowerCI upperCI
## Mod.avg.diff -0.03012683 0.02479459a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("No.pattern", "No.ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~No.pattern.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~No.pattern.VS.No.ovlp+SOLOGAP -0.001586391
## lowerCI upperCI
## Mod.avg.diff -0.02465744 0.02148466a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~No.pattern.VS.Ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~No.pattern.VS.Ovlp+SOLOGAP 0.001079726
## lowerCI upperCI
## Mod.avg.diff -0.02151804 0.02367749a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.ovlp", "Ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~Ovlp.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~Ovlp.VS.No.ovlp+DUETGAP -0.002320566
## lowerCI upperCI
## Mod.avg.diff -0.02865142 0.02401028a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.ovlp", "No.pattern")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~No.pattern.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~No.pattern.VS.No.ovlp+DUETGAP -0.001422199
## lowerCI upperCI
## Mod.avg.diff -0.02411251 0.02126811a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]] <- data.frame(data="All data",factor="gap.cv~No.pattern.VS.Ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff All data gap.cv~No.pattern.VS.Ovlp+DUETGAP -0.000898367
## lowerCI upperCI
## Mod.avg.diff -0.02321441 0.02141767gap duration for known dyads
# by pattern for Known individuals
me.cv<-me.cv[!is.na(me.cv$indivsolo),]
names(me.cv)[grep("indivsolo", names(me.cv))] <- "indiv"
cat("sample size number of bouts")## sample size number of boutsnrow(me.cv)## [1] 102cat("sample size number of individuals")## sample size number of individualslength(unique(me.cv$indiv))## [1] 12Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype * pattern, random=~1 | indiv/sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~pattern, random=~1 | indiv/sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype , random=~1 | indiv/sound.files, data = me.cv, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~1 , random=~1 | indiv/sound.files, data = me.cv, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
Modnames<-gsub("~gap","gap variation~",Modnames)
aictabs[[length(aictabs)+1]]<-aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
aictabs[[length(aictabs)]]##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## gap variation~gaptype 5 -300.26 0.00 0.93 0.93 155.44
## gap variation~gaptype * pattern 9 -295.03 5.23 0.07 1.00 157.49
## gap variation~1 4 -283.37 16.89 0.00 1.00 145.89
## gap variation~pattern 6 -279.23 21.03 0.00 1.00 146.06cs[[length(cs)+1]]<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs[[length(cs)]]##
## Confidence set for the best model
##
## Method: raw sum of model probabilities
##
## 95% confidence set:
## K AICc Delta_AICc AICcWt
## gap variation~gaptype 5 -300.26 0.00 0.93
## gap variation~gaptype * pattern 9 -295.03 5.23 0.07
##
## Model probabilities sum to 1confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("No.ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~gaptype+NO.OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.cv~gaptype+NO.OVLP 0.04004719 0.01454456
## upperCI
## Mod.avg.diff 0.06554983a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("No.pattern", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~gaptype+NO.PATT",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.cv~gaptype+NO.PATT 0.04150475 0.019708
## upperCI
## Mod.avg.diff 0.0633015a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "sologap"), pattern = c("Ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~gaptype+OVLP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size lowerCI
## Mod.avg.diff known dyads gap.cv~gaptype+OVLP 0.04331794 0.01637223
## upperCI
## Mod.avg.diff 0.07026365a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("Ovlp", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~Ovlp.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.cv~Ovlp.VS.No.ovlp+SOLOGAP -0.002128939
## lowerCI upperCI
## Mod.avg.diff -0.03014032 0.02588245a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("No.pattern", "No.ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~No.pattern.VS.No.ovlp+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.cv~No.pattern.VS.No.ovlp+SOLOGAP -0.001319688
## lowerCI upperCI
## Mod.avg.diff -0.02575237 0.023113a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("sologap", "sologap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~Ovlp.VS.No.pattern+SOLOGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.cv~Ovlp.VS.No.pattern+SOLOGAP 0.0008092507
## lowerCI upperCI
## Mod.avg.diff -0.02265158 0.02427008a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.ovlp", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="cgap.cv~Ovlp.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads cgap.cv~Ovlp.VS.No.ovlp+DUETGAP -0.001141808
## lowerCI upperCI
## Mod.avg.diff -0.02593717 0.02365355a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.ovlp", "No.pattern")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~No.pattern.VS.No.ovlp+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor
## Mod.avg.diff known dyads gap.cv~No.pattern.VS.No.ovlp+DUETGAP
## effect_size lowerCI upperCI
## Mod.avg.diff -0.0001378696 -0.02264006 0.02236432a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duetgap", "duetgap"), pattern = c("No.pattern", "Ovlp")))
es[[length(es)+1]]<-data.frame(data="known dyads",factor="gap.cv~Ovlp.VS.No.pattern+DUETGAP",effect_size=a$Mod.avg.eff, lowerCI= a$Lower.CL, upperCI=a$Upper.CL)
es[[length(es)]]## data factor effect_size
## Mod.avg.diff known dyads gap.cv~Ovlp.VS.No.pattern+DUETGAP -0.001003939
## lowerCI upperCI
## Mod.avg.diff -0.02485721 0.02284933##save results
es <- lapply(es, function(x)
{colnames(x)<-c("Data set", "Fixed effect (predictor)", "Effect size", "Lower 95 CI", "Upper 95% CI")
return(x)}
)
esT<-do.call("rbind", es)
esT$`Fixed effect (predictor)`<- sapply(1:nrow(esT), function(x){
if(any(all(c(esT$`Lower 95 CI`[x], esT$`Upper 95% CI`[x]) > 0), all(c(esT$`Lower 95 CI`[x], esT$`Upper 95% CI`[x]) < 0)))
return(paste(esT$`Fixed effect (predictor)`[x], "*",sep = "")) else
return(as.character(esT$`Fixed effect (predictor)`[x]))})
colnames(esT)<-c("Data set", "Fixed effect (predictor)", "Effect size", "Lower 95 CI", "Upper 95% CI")
write.csv(esT, "~/Dropbox/LBH data/Vocalizations/Coordinated singing/effect size results.csv", row.names = F)
round.esT <- esT
round.esT[,3:5] <- apply(esT[,3:5],2,round, 3)
library(grid)
library(gridExtra)
g <- tableGrob(round.esT, rows = NULL)
grid.newpage()
grid.draw(g)
names(aictabs) <-rep(c("All data", "known dyads"), 3)
aicT<-list()
for(i in 1:length(aictabs))
{
aicT[[length(aicT)+1]]<-rbind(data.frame(data=names(aictabs)[i],aictabs[[i]]))
}
aicT<-do.call("rbind",aicT)
write.csv(aicT, "~/Dropbox/LBH data/Vocalizations/Coordinated singing/All models.csv", row.names = F)
# csT<-list()
# for(i in 1:length(cs))
# {
# csT[[length(csT)+1]]<-rbind(data.frame(data=names(cs)[i],cs[[i]]$table))
# }
#
# csT<-do.call("rbind",csT)
#
# write.csv(esT, "~/Dropbox/LBH data/Vocalizations/Coordinated singing/95 prob models.csv", row.names = F)gap duration for known dyads controlling for pair
# by pattern for Known individuals
me<-me[grep("//",me$pair),]
cat("sample size")
nrow(me)
Cand.models <- list( )
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype * pattern, random=~1 | pair/sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~pattern, random=~1 | pair/sound.files, data = me, method = "ML")
Cand.models[[length(Cand.models)+1]]<-lme(gap~gaptype , random=~1 | pair/sound.files, data = me, method = "ML")
Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
# Modnames<-gsub("sing.perf","",Modnames)
aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
cs<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
cs
confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
a<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.ovlp", "No.ovlp")))
a
b<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("No.pattern", "No.pattern")))
b
c<-modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "solo gap"), pattern = c("Ovlp", "Ovlp")))
c
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("Ovlp", "No.ovlp")))
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "No.ovlp")))
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("solo gap", "solo gap"), pattern = c("No.pattern", "Ovlp")))
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "Ovlp")))
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.ovlp", "No.pattern")))
modavgEffect(cand.set = confmodels, modnames = confmodnames, newdata = data.frame(gaptype = c("duet gap", "duet gap"), pattern = c("No.pattern", "Ovlp")))Graphs on effect of duration on test performance on real data
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Graphs testing coordination performance along events (on real data)
Graphs testing coordination performance along events (on simulated data)
AOresults<-droplevels(ssresults[ssresults$pattern!="No.pattern",])
a<-a1<-table(ssresults$pair, ssresults$pattern)[,c(1,3)]
a[a>1]<-1
cat("pairs with two both patterns")## pairs with two both patternsa1[apply(a,1,sum)>1,]##
## No.ovlp Ovlp
## 103//104 1 9
## 108//111 3 1
## 111//117 3 6
## 143//148 1 2
## -1//90 3 2nrow(a1[apply(a,1,sum)>1,])## [1] 5cat("pairs with single pattern")## pairs with single patterna1[apply(a,1,sum)==1,]##
## No.ovlp Ovlp
## 108//109 13 0
## 108//117 0 3
## 24//88 2 0
## 32//54 0 1
## 32//60 0 1
## 54//60 0 1
## 60//81 2 0
## 72//100 0 1
## 86//149 1 0nrow(a1[apply(a,1,sum)==1,])## [1] 9#chi square
a<-a1<-table(ssresults$pair, ssresults$pattern)[,c(1,3)]
a1<-a1[apply(a,1,sum)>0,]
# library(MASS)
# a1[a1>1]<-1
# chisq.test(a1)
# spagg<-aggregate(AOresults[,grep("logdist|dur|Time", colnames(AOresults))], list(AOresults$pattern, AOresults$pair), mean)
# colnames(spagg)[1:2]<-c("Pattern", "Pair")
#
#
# pattern <- as.integer(ifelse(spagg$Pattern=="Ovlp",1,0))
# dist<-spagg$Distance
# pair<-spagg$Pair
#
# Cand.models <- list( )
#
# Cand.models[[length(Cand.models)+1]]<-glmer(Pattern ~ logdist + Time + (1 | Pair), family = binomial, data = spagg)
# Cand.models[[length(Cand.models)+1]]<-glmer(Pattern ~ Time + (1 | Pair), family = binomial, data = spagg)
# Cand.models[[length(Cand.models)+1]]<-glmer(Pattern ~ logdist +(1 | Pair), family = binomial, data = spagg)
# Cand.models[[length(Cand.models)+1]]<-glmer(Pattern ~ 1 + (1 | Pair), family = binomial, data = spagg)
# Cand.models[[length(Cand.models)+1]]<-glmer(Pattern ~ dur +(1 | Pair), family = binomial, data = spagg)
#
# Modnames<-sapply(Cand.models, function(x) paste(as.formula(x), collapse = ""))
#
# cat("coor pattern vs distance repeated measures by pair")
#
# aictab(cand.set = Cand.models, modnames = Modnames, second.ord = T)
#
# cs<-confset(cand.set = Cand.models, modnames = Modnames, second.ord = T, method = "raw")
#
# cs
#
# confmodels<-Cand.models[Modnames %in% cs[[length(cs)]]$table[,1]]
#
# confmodnames<-Modnames[Modnames %in% cs[[length(cs)]]$table[,1]]
#
# modavg(parm = "logdist", cand.set = confmodels, modnames = confmodnames)meanSNR<-aggregate(SNRdata[,grep("SNR", colnames(SNRdata))], list(SNRdata$coord.sing.event, SNRdata$sound.files), mean)
colnames(meanSNR)<-c("coord.sing.event", "channel", "SNR")
meanSNR2<-merge(meanSNR, results[,grep("coord.sing.event|dur|sing.perf", colnames(results))], by= "coord.sing.event")
meanSNR2<-meanSNR2[meanSNR2$SNR>2,]
a<-table(meanSNR2$coord.sing.event)
a<-a[a>1]
meanSNR2<-meanSNR2[meanSNR2$coord.sing.event %in% names(a),]
meanSNR2<-meanSNR2[order(meanSNR2$sing.perf),]
# View(meanSNR2[,c(1,3,4,5)])
#best alternating ones
meanSNR2[1:30,c(1,3,4,5)]
#best overlapping ones
meanSNR2[nrow(meanSNR2)-30:30,c(1,3,4)]
#ovlp 391 from 1:02 to 1:09
#alternation 390 from 43 to 52
library(seewave)
library(tuneR)
setwd("/media/m/Eutoxeres/LBH song recordings/Coordinated singing/cleaned/split wavs")
wavs<-list.files(pattern = ".wav$",ignore.case = T)
wavleft<-grep("left", wavs,value = T)
wavright<-grep("right", wavs,value = T)library(RColorBrewer)
tiff(filename = "/home/m/Dropbox/LBH data/Vocalizations/Coordinated singing/graphs/spectrograms.tiff",res = 200,
width = 2000, height = 1000)
par(mfrow=c(4,1))
# col1<-colorRampPalette(brewer.pal(9,"Blues"))(10)[9]
# col2<-colorRampPalette(brewer.pal(9,"Blues"))(10)[5]
col1<-"#147777"
col2<-"#A13737"
ovlp<-90
#first event
flim<-c(2.5,9.8)
tlim<-c(62.5, 72)
event<-391
par(mar=c(0,4.2,4,0.6))
spectro(readWave(grep(event,wavleft,value = T)),f = 44100, ovlp = ovlp, wl = 300, flim = flim, tlim = tlim,palette = reverse.gray.colors.2, scale = F, grid = F, axisX = F, flab = "", tlab = "",trel = F)
start <- alldat$start[alldat$sound.files == grep(event,wavleft,value = T)]-tlim[1]
end <- alldat$end[alldat$sound.files == grep(event,wavleft,value = T)]-tlim[1]
start <- start[start > 0 & start < tlim[2]]
end <- end[end > 0 & end < tlim[2]]
invisible(sapply(1:length(start),FUN = function(x)
{
rect(xleft = start[x], ybottom = flim[1], xright = end[x], ytop = flim[1]+0.5, col = col1, border = col1)
})
)
text("Frequency (kHz)", x = -0.04*(tlim[2]-tlim[1]) , y = 1.95, srt = 90, xpd=T, cex=1.5)
text("A", x = 0.1 , y = 9, xpd=T, cex=1.5, font = 2)
#second channel
par(mar=c(4,4.2, 0, 0.6))
spectro(readWave(grep(event,wavright,value = T)),f = 44100, ovlp = ovlp, wl = 300, flim = flim, tlim = tlim,palette = reverse.gray.colors.2, scale = F, grid = F, trel = F, flab= "", tlab = "")
start <- alldat$start[alldat$sound.files == grep(event,wavright,value = T)]-tlim[1]
end <- alldat$end[alldat$sound.files == grep(event,wavright,value = T)]-tlim[1]
start <- start[start > 0 & start < tlim[2]]
end <- end[end > 0 & end < tlim[2]]
invisible(sapply(1:length(start),FUN = function(x)
{
rect(xleft = start[x], ybottom = flim[2]-0.5, xright = end[x], ytop = flim[2], col = col2, border = col2)
})
)
text("Frequency (kHz)", x = -0.04*(tlim[2]-tlim[1]), y = flim[2] - flim[1]+1.95, srt = 90, xpd=T, cex=1.5)
text("Time (s)", x = (tlim[2]-tlim[1])/2, y = (flim[2]-flim[1]/1.2)-(flim[2]-flim[1]), xpd=T, cex=1.5)
#second event
tlim<-c(43, 43 + tlim[2] - tlim[1])
event<-390
par(mar=c(0,4.2,4,0.6))
spectro(readWave(grep(event,wavleft,value = T)),f = 44100, ovlp = ovlp, wl = 300, flim = flim, tlim = tlim,palette = reverse.gray.colors.2, scale = F, grid = F, axisX = F, flab = "", tlab = "",trel = F)
start <- alldat$start[alldat$sound.files == grep(event,wavleft,value = T)]-tlim[1]
end <- alldat$end[alldat$sound.files == grep(event,wavleft,value = T)]-tlim[1]
start <- start[start > 0 & start < tlim[2]]
end <- end[end > 0 & end < tlim[2]]
invisible(sapply(1:length(start),FUN = function(x)
{
rect(xleft = start[x], ybottom = flim[1], xright = end[x], ytop = flim[1]+0.5, col = col1, border = col1)
})
)
text("Frequency (kHz)", x = -0.04*(tlim[2]-tlim[1]) , y = 1.95, srt = 90, xpd=T, cex=1.5)
text("B", x = 0.1 , y = 9, xpd=T, cex=1.5, font = 2)
#second channel
par(mar=c(4,4.2, 0, 0.6))
spectro(readWave(grep(event,wavright,value = T)),f = 44100, ovlp = ovlp, wl = 300, flim = flim, tlim = tlim,palette = reverse.gray.colors.2, scale = F, grid = F, trel = F, flab= "", tlab = "")
start <- alldat$start[alldat$sound.files == grep(event,wavright,value = T)]-tlim[1]
end <- alldat$end[alldat$sound.files == grep(event,wavright,value = T)]-tlim[1]
# start <- start[start > 0 & start < tlim[2]]
# end <- end[end > 0 & end < tlim[2]]
invisible(sapply(1:length(start),FUN = function(x)
{
rect(xleft = start[x], ybottom = flim[2]-0.5, xright = end[x], ytop = flim[2], col = col2, border = col2)
})
)
text("Frequency (kHz)", x = -0.04*(tlim[2]-tlim[1]), y = flim[2] - flim[1]+1.95, srt = 90, xpd=T, cex=1.5)
text("Time (s)", x = (tlim[2]-tlim[1])/2, y = (flim[2]-flim[1]/1.2)-(flim[2]-flim[1]), xpd=T, cex=1.5)
dev.off()require(gridExtra)
tiff(filename = "/home/m/Dropbox/LBH data/Vocalizations/Coordinated singing/graphs/multipanel fig.tiff",res = 300,
width = 4000, height = 2400)
ggs<-grid.arrange(gg1.col,gg9.col,gg22.col,gg23.col, ncol=2)
dev.off()
jpeg(filename = "/home/m/Dropbox/LBH data/Vocalizations/Coordinated singing/graphs/multipanel fig.jpeg",res = 300,
width = 4000, height = 2400)
ggs<-grid.arrange(gg1.col,gg9.col,gg22.col,gg23.col, ncol=2)
dev.off()