Set working directory, load and subset all the data tables:
setwd("~/Grasshopper Sparrows/Results and Analysis/Aggregation Manuscript")
#results of cooperative defense trials
behavior_all<-read.csv ("cooperative_defense_trials.csv")
BLRS<-subset(behavior_all, Predator %in% c("BLRS"))
BHCO<-subset(behavior_all, Predator %in% c("BHCO"))
HOSP<-subset(behavior_all, Predator %in% c("HOSP"))
threat<-subset(behavior_all, Predator %in% c("BLRS", "BHCO"))
#all nests, with partial predation, brood parasitism, and density columns
predation<-read.csv("nests_predation_and_parasitism.csv")
#observed nests, with cooperative care and density columns
coop<-read.csv("nests_cooperative_care.csv")
#all nests, with EPP column and density column
nests<-read.csv("nests_withEPP.csv")
#density of all males from GIS rasters
alldensity<-read.csv("density_all_males.csv")
#density of all males classified as "aggregated" (density values in upper quartile)
agg_all<-read.csv("density_aggregated_only.csv")
#density of all males classified as "related" (r values above 0.25)
related_all<-read.csv("density_related_only.csv")
#relatedness values of all cuckolded males
EPP_withrelatedness_all<-read.csv("density_cuckolded_only.csv")
#relatedness values of cuckolded males and the likely genetic fathers
EPP_dads<-read.csv("density_genetic_father_relatedness.csv")
#density of all males, without null values
agg_all_noNulls<-read.csv("density_all_males_noNulls.csv")
#RMark nest success file, 2014 and 2015 nests
GRSP = read.csv("nest_success_model.CSV")
GRSPuntransformed = read.csv("nest_success_model.CSV")
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
library(gridExtra)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked _by_ '.GlobalEnv':
##
## coop
library(forcats)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Compare densities over time periods:
density<-aov(density~time.of.year, data=alldensity)
summary(density)
## Df Sum Sq Mean Sq F value Pr(>F)
## time.of.year 8 37.71 4.713 45.91 <2e-16 ***
## Residuals 1061 108.93 0.103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc <- TukeyHSD(x=density, conf.level=0.95)
posthoc
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = density ~ time.of.year, data = alldensity)
##
## $time.of.year
## diff lwr upr p adj
## early2014-early2013 0.48630795 0.32791404 0.644701863 0.0000000
## early2015-early2013 0.03077493 -0.12376503 0.185314891 0.9995081
## late2013-early2013 -0.16469153 -0.33256269 0.003179623 0.0593554
## late2014-early2013 -0.08614144 -0.25745898 0.085176104 0.8247698
## late2015-early2013 0.07374853 -0.08149836 0.228995433 0.8664300
## mid 2015-early2013 0.03764046 -0.11548629 0.190767215 0.9977393
## mid2013-early2013 -0.18531980 -0.34510019 -0.025539421 0.0098626
## mid2014-early2013 -0.01678431 -0.17656469 0.142996076 0.9999963
## early2015-early2014 -0.45553302 -0.57457449 -0.336491555 0.0000000
## late2013-early2014 -0.65099949 -0.78690017 -0.515098799 0.0000000
## late2014-early2014 -0.57244939 -0.71258494 -0.432313842 0.0000000
## late2015-early2014 -0.41255942 -0.53251721 -0.292601629 0.0000000
## mid 2015-early2014 -0.44866749 -0.56586849 -0.331466492 0.0000000
## mid2013-early2014 -0.67162776 -0.79739756 -0.545857957 0.0000000
## mid2014-early2014 -0.50309226 -0.62886206 -0.377322460 0.0000000
## late2013-early2015 -0.19546646 -0.32685507 -0.064077860 0.0001462
## late2014-early2015 -0.11691637 -0.25268065 0.018847914 0.1573901
## late2015-early2015 0.04297360 -0.07184730 0.157794507 0.9638531
## mid 2015-early2015 0.00686553 -0.10507214 0.118803205 0.9999999
## mid2013-early2015 -0.21609474 -0.33697489 -0.095214583 0.0000012
## mid2014-early2015 -0.04755924 -0.16843939 0.073320914 0.9515039
## late2014-late2013 0.07855010 -0.07221484 0.229315034 0.7942044
## late2015-late2013 0.23844007 0.10622068 0.370659453 0.0000009
## mid 2015-late2013 0.20233199 0.07260856 0.332055428 0.0000502
## mid2013-late2013 -0.02062827 -0.15814241 0.116885864 0.9999417
## mid2014-late2013 0.14790722 0.01039309 0.285421361 0.0241132
## late2015-late2014 0.15988997 0.02332152 0.296458419 0.0087304
## mid 2015-late2014 0.12378190 -0.01037154 0.257935332 0.0977728
## mid2013-late2014 -0.09917837 -0.24087915 0.042522417 0.4225598
## mid2014-late2014 0.06935713 -0.07234366 0.211057914 0.8455837
## mid 2015-late2015 -0.03610807 -0.14901974 0.076803591 0.9864568
## mid2013-late2015 -0.25906834 -0.38085098 -0.137285696 0.0000000
## mid2014-late2015 -0.09053284 -0.21231549 0.031249801 0.3364857
## mid2013-mid 2015 -0.22296027 -0.34202838 -0.103892155 0.0000003
## mid2014-mid 2015 -0.05442477 -0.17349288 0.064643342 0.8899380
## mid2014-mid2013 0.16853550 0.04102399 0.296047002 0.0014166
Densities by time of year:
alldensity$time.of.year = fct_relevel(alldensity$time.of.year, "early2013", "mid2013", "late2013", "early2014", "mid2014", "late2014", "early2015", "mid 2015", "late2015")
ggplot(data=alldensity, aes(time.of.year, density)) +
geom_boxplot(fill="darkgray")+
xlab("Time of Year")+
ylab("Density at each territory (males/ha)")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
Histogram of densities:
fig1<-ggplot(alldensity, aes(x=density)) +
geom_histogram(color="black", fill="gray")+
xlab("Density of breeding males \n (territories/ha)")+
ylab("Number of centroids")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure1.jpeg", width = 8, height = 5, units = 'in', res = 600)
fig1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#dev.off()
Hypothesis 1, Prediction i: are aggregated birds more likely to approach threat models? Welch’s t-tests comparing aggregated males’ and unaggregated males’ closest approach distance to snakes, cowbirds, and House Sparrow controls:
t.test(Closest.Approach.Distance..m.~Aggregated., data=BLRS) #p=0.3508, t=1.0965, df=3.0941
##
## Welch Two Sample t-test
##
## data: Closest.Approach.Distance..m. by Aggregated.
## t = 1.0965, df = 3.0941, p-value = 0.3508
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -14.58612 30.32898
## sample estimates:
## mean in group N mean in group Y
## 14.800000 6.928571
t.test(Closest.Approach.Distance..m.~Aggregated., data=BHCO) #p=0.7274, t=0.37584, df=3.7313
##
## Welch Two Sample t-test
##
## data: Closest.Approach.Distance..m. by Aggregated.
## t = 0.37584, df = 3.7313, p-value = 0.7274
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.77008 23.15342
## sample estimates:
## mean in group N mean in group Y
## 12.36667 9.67500
t.test(Closest.Approach.Distance..m.~Aggregated., data=HOSP) #p=0.9503, t=0.067665, df=3
##
## Welch Two Sample t-test
##
## data: Closest.Approach.Distance..m. by Aggregated.
## t = 0.067665, df = 3, p-value = 0.9503
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.904821 7.204821
## sample estimates:
## mean in group N mean in group Y
## 5.00 4.85
Hypothesis 1, Prediction ii: are aggregated birds’ approach latencies lower? Welch’s t-tests comparing aggregated males’ and unaggregated males’ approach time to snakes, cowbirds, and House Sparrow controls:
t.test(ApproachTime~Aggregated., data=BLRS) #p=0.3826, t=-0.91926, df=8.7676
##
## Welch Two Sample t-test
##
## data: ApproachTime by Aggregated.
## t = -0.91926, df = 8.7676, p-value = 0.3826
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -12.185079 5.163672
## sample estimates:
## mean in group N mean in group Y
## 5.191667 8.702370
t.test(ApproachTime~Aggregated., data=BHCO) #p=0.943, t=-0.075237, df=4.861
##
## Welch Two Sample t-test
##
## data: ApproachTime by Aggregated.
## t = -0.075237, df = 4.861, p-value = 0.943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -15.81022 14.91856
## sample estimates:
## mean in group N mean in group Y
## 8.233333 8.679165
t.test(ApproachTime~Aggregated., data=HOSP) #p=0.5541, t=-0.64493, df=3.9978
##
## Welch Two Sample t-test
##
## data: ApproachTime by Aggregated.
## t = -0.64493, df = 3.9978, p-value = 0.5541
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -8.732682 5.441015
## sample estimates:
## mean in group N mean in group Y
## 1.583333 3.229166
Graphing the results of the first two predictions:
p1<- ggplot(behavior_all, aes(x=Predator, y=Closest.Approach.Distance..m., fill=Aggregated.)) +
geom_boxplot()+
scale_x_discrete(labels=c("BHCO" = "Cowbird", "BLRS" = "Snake",
"HOSP" = "Control"))+
xlab("")+
ylab("Closest approach (m)")+
scale_fill_manual(name = "Aggregated?", values = c("gray", "black")
, labels = c("N" = "No", "Y" = "Yes"))+
annotate("text", x = 0.6, y = 38, label = "A", size=12)+
theme(
legend.position=c(.85, .85),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
p2<- ggplot(behavior_all, aes(x=Predator, y=ApproachTime, fill=Aggregated.)) +
geom_boxplot()+
scale_x_discrete(labels=c("BHCO" = "Cowbird", "BLRS" = "Snake",
"HOSP" = "Control"))+
xlab("")+
ylab("Time to closest approach (min)")+
scale_fill_manual(name = "Aggregated?", values = c("gray", "black")
, labels = c("N" = "No", "Y" = "Yes"))+
annotate("text", x = 0.6, y = 30, label = "B", size=12)+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position=c(.85, .85),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure2.jpeg", width = 8, height = 5, units = 'in', res = 600)
grid.arrange(p1, p2, ncol = 2)
#dev.off()
Hypothesis 1, Prediction iii: are aggregated birds more likely to make alarm vocalizations? Chi-squared analysis:
tbl = table(threat$Chipping, threat$Aggregated.)
tbl # the contingency table
##
## N Y
## No 6 7
## Yes 1 4
chisq.test(tbl) #x2=0.23017, df=1, p=0.6314
## Warning in chisq.test(tbl): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tbl
## X-squared = 0.23017, df = 1, p-value = 0.6314
Hypothesis 1, Prediction iv: are aggregated birds less likely to be depredated or parasitized? Welch’s t-test comparing densities (log200) of partially depredated or parasitized nests to successful nests:
t.test(log200~Predation., predation)
##
## Welch Two Sample t-test
##
## data: log200 by Predation.
## t = -1.6886, df = 114.03, p-value = 0.09403
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.32432159 0.02584033
## sample estimates:
## mean in group FALSE mean in group TRUE
## -0.4266253 -0.2773846
t.test(log200~Parasitism, predation)
##
## Welch Two Sample t-test
##
## data: log200 by Parasitism
## t = 2.2185, df = 204.19, p-value = 0.02762
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.0185572 0.3150097
## sample estimates:
## mean in group Parasitized mean in group Unparasitized
## -0.2222689 -0.3890523
Graphing the results of prediction iv:
p1<- ggplot(predation, aes(Predation., log200, group=Predation.))+
geom_boxplot(fill="lightgray")+
annotate("text", x = 0.6, y = 0.5, label = "A", size=12)+
scale_x_discrete(labels=c("FALSE" = "Successful", "TRUE" = "Depredated"))+
xlab("")+
ylab("Territory density (log-territories/ha)")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
p2<- ggplot(predation, aes(Parasitism, log200, group=Parasitism))+
geom_boxplot(fill="lightgray")+
annotate("text", x = 0.6, y = 0.5, label = "B", size=12)+
xlab("")+
ylab("Territory density (log-territories/ha)")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure3.jpeg", width = 8, height = 5, units = 'in', res = 600)
grid.arrange(p1, p2, ncol = 2)
#dev.off()
Hypothesis 2, Prediction v: do aggregated nests experience higher rates of cooperative care? Welch’s t-test to compare the densities of nests with and without cooperative care:
t.test(log200~CooperativeCare, data=coop)
##
## Welch Two Sample t-test
##
## data: log200 by CooperativeCare
## t = 0.04429, df = 25.679, p-value = 0.965
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4313205 0.4503051
## sample estimates:
## mean in group no mean in group yes
## -0.3522496 -0.3617419
Graph results of prediction v:
fig4<- ggplot(coop, aes(CooperativeCare, log200, group=CooperativeCare))+
geom_boxplot(fill="lightgray")+
xlab("Extra-pair defense")+
ylab("Density (log of males/ha)")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure4.jpeg", width = 5, height = 5, units = 'in', res = 600)
fig4
#dev.off()
Hypothesis 3, Prediction vi: are nests with extra-pair young located in areas of higher densities? ANOVA to compare densities of nests with different numbers of EPY:
an_density_numEPP<-aov(log200~as.factor(NumberEPP), data=nests)
summary(an_density_numEPP) #df=2, 30, f=2.054, p=0.146
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(NumberEPP) 2 1.759 0.8796 2.054 0.146
## Residuals 30 12.848 0.4283
## 189 observations deleted due to missingness
Graph of prediction vi:
fig5<-ggplot(nests, aes(NumberEPP, log200, group=NumberEPP))+
geom_boxplot(fill="lightgray")+
xlab("Number of extra-pair nestlings")+
ylab("Density (log of males/ha)")+
scale_x_continuous(breaks=c(0,1,2))+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure5.jpeg", width = 5, height = 5, units = 'in', res = 600)
fig5
## Warning: Removed 189 rows containing missing values (stat_boxplot).
#dev.off()
Hypothesis 4, Prediction viii: is the average distance between related males’ territories lower than the distance between unrelated males’ territories? Fit a linear model comparing relatedness to distance with all males:
fit <- lm(Queller...Goodnight..1989~Distance, data=agg_all)
summary(fit) #r2=0.0002789
##
## Call:
## lm(formula = Queller...Goodnight..1989 ~ Distance, data = agg_all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80089 -0.14352 -0.02026 0.12431 1.00726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.010e-03 1.804e-03 -0.560 0.5756
## Distance -1.535e-06 5.012e-07 -3.062 0.0022 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1953 on 33608 degrees of freedom
## (348 observations deleted due to missingness)
## Multiple R-squared: 0.0002789, Adjusted R-squared: 0.0002492
## F-statistic: 9.376 on 1 and 33608 DF, p-value: 0.0022
Graph of relatedness to distance:
fig6<-ggplot(data=agg_all, aes(Queller...Goodnight..1989, Distance)) +
geom_point(colour="darkgray", size=1)+
xlab("Relationship Coefficient")+
ylab("Distance apart (m)")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure6.jpeg", width = 7, height = 5, units = 'in', res = 600)
fig6
## Warning: Removed 348 rows containing missing values (geom_point).
#dev.off()
Compare the relatedness of aggregated and unaggregated individuals:
t.test(Queller...Goodnight..1989~agg_all, data=agg_all)
##
## Welch Two Sample t-test
##
## data: Queller...Goodnight..1989 by agg_all
## t = -7.0117, df = 30797, p-value = 2.404e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.01937081 -0.01090695
## sample estimates:
## mean in group no mean in group yes
## -0.012076508 0.003062372
Compare the distance between related and unrelated individuals:
t.test(Distance~related, data=related_all)
##
## Welch Two Sample t-test
##
## data: Distance by related
## t = 3.5855, df = 4524.7, p-value = 0.00034
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 59.86089 204.30240
## sample estimates:
## mean in group no mean in group yes
## 2916.698 2784.617
Graph the last two results:
p1<-ggplot(agg_all, aes(agg_all,Queller...Goodnight..1989))+
geom_boxplot(fill="lightgray")+
xlab("Aggregated?")+
ylab("Relationship Coefficient")+
annotate("text", x = 0.6, y = 1, label = "A", size=12)+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
p2<- ggplot(related_all, aes(related, Distance))+
geom_boxplot(fill="lightgray")+
xlab("Related?")+
ylab("Distance between (m)")+
annotate("text", x = 0.6, y = 11000, label = "B", size=12)+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("S1.jpeg", width = 8, height = 5, units = 'in', res = 600)
grid.arrange(p1, p2, ncol = 2)
## Warning: Removed 348 rows containing non-finite values (stat_boxplot).
#dev.off()
Hypothesis 4, Prediction ix: are genetic fathers related to cuckolded males?
Compare relatedness values of males that were cuckolded to those that were not:
t.test(Queller...Goodnight..1989~EPP, data=EPP_withrelatedness_all)
##
## Welch Two Sample t-test
##
## data: Queller...Goodnight..1989 by EPP
## t = -8.1076, df = 1701.8, p-value = 9.769e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.07311396 -0.04462985
## sample estimates:
## mean in group no mean in group yes
## -0.06908837 -0.01021646
Graph it:
s2<-ggplot(EPP_withrelatedness_all, aes(EPP,Queller...Goodnight..1989))+
geom_boxplot(fill="lightgray")+
xlab("Extra-pair nestlings?")+
ylab("Relationship Coefficient")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("S2.jpeg", width = 5, height = 5, units = 'in', res = 600)
s2
#dev.off()
Compare the distribution of all relatedness values (between all males) to the distribution of relatedness values generated from comparing cuckolded males to the true genetic fathers of offspring:
original<-(agg_all_noNulls$Queller...Goodnight..1989)
resample<-(EPP_dads$Queller...Goodnight..1989)
ks.test(original, resample, alternative = c("two.sided"), exact = NULL)
## Warning in ks.test(original, resample, alternative = c("two.sided"), exact
## = NULL): p-value will be approximate in the presence of ties
##
## Two-sample Kolmogorov-Smirnov test
##
## data: original and resample
## D = 0.19236, p-value = 0.2173
## alternative hypothesis: two-sided
Graph it:
p1<- ggplot(EPP_dads, aes(x=Queller...Goodnight..1989)) +
geom_histogram(color="black", fill="gray", bins=10)+
xlab("Relatedness of social and genetic fathers")+
ylab(" \n Number of pairs")+
coord_cartesian(xlim = c(-0.5, 1))+
scale_y_continuous(breaks=c(0,2,4,6,8,10))+
annotate("text", x = -0.5, y = 10, label = "A", size=12)+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
p2<- ggplot(agg_all_noNulls, aes(x=Queller...Goodnight..1989)) +
geom_histogram(color="black", fill="gray")+
xlab("Relatedness between all males")+
ylab("Number of pairs")+
coord_cartesian(xlim = c(-0.5, 1))+
annotate("text", x = -0.5, y = 4000, label = "B", size=12)+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("S3.jpeg", width = 8, height = 5, units = 'in', res = 600)
grid.arrange(p1, p2, ncol = 1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#dev.off()
Nest success
Set up rmark file and look for errors
# look at first and last occasions
# needs to be 1 for dataset
(min = min(GRSP$FirstFound))
# this gives maximum count for number of occasions
(occ = max(GRSP$LastChecked))
# check input for errors
GRSP$Diff = GRSP$LastChecked-GRSP$LastPresent
GRSP$Error="No"
GRSP$Error[(GRSP$LastPresent-GRSP$FirstFound)<0] <- "Yes"
GRSP$Error[(GRSP$LastChecked-GRSP$LastPresent)<0] <- "Yes"
GRSP$Error[GRSP$Fate==1 & GRSP$Diff<=0] <- "Yes"
GRSP$Error[GRSP$Fate==0 & GRSP$Diff!=0] <- "Yes"
(Error <- subset(GRSP, Error=="Yes"))
# drop errors if present
GRSP <- GRSP[which(GRSP$Error=='No'), ]
Treat factors as factors and scale continuous variables:
GRSP$Year= factor(GRSP$Year)
GRSP$Treatment=factor(GRSP$Treatment)
GRSP$Edge= scale(GRSP$Edge)
GRSP$Road= scale(GRSP$Road)
GRSP$Tree= scale(GRSP$Tree)
GRSP$Burn= scale(GRSP$Burn)
GRSP$Density=scale(GRSP$Density)
GRSP$AbsForb=scale(GRSP$AbsForb)
GRSP$Height=scale(GRSP$Height)
Run a set of models to compare temporal effects on nest success:
run.GRSP.temporal =function()
{ # Daily survival rate (DSR) is constant
t1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~1)))
# DSR follows a linear trend through time
t2=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Time)))
# DSR follows a quadratic trend through time
t3=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Time+I(Time^2))))
# DSR is different for every day of exposure period
t4=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~time)))
# Year
t5=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Year)), groups=c("Year"))
#Year+Date
t6=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Year+Time)), groups=c("Year"))
#Year*Date
t7=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Year*Time)), groups=c("Year"))
return(collect.models())
}
GRSP.results.temporal=run.GRSP.temporal() # This runs the models above and takes a minute or 2
##
## Note: only 51 parameters counted of 102 specified parameters
## AICc and parameter count have been adjusted upward
Examine table of model results:
GRSP.results.temporal # print model selection table to screen
## model npar AICc DeltaAICc weight Deviance
## 5 S(~Year) 2 801.7144 0.000000 0.47322044 797.7059
## 6 S(~Year + Time) 3 802.7273 1.012895 0.28517810 796.7102
## 7 S(~Year * Time) 4 804.4095 2.695127 0.12297707 796.3811
## 3 S(~Time + I(Time^2)) 3 805.4198 3.705395 0.07420740 799.4027
## 1 S(~1) 1 807.2571 5.542676 0.02961329 805.2542
## 2 S(~Time) 2 808.6438 6.929370 0.01480369 804.6353
## 4 S(~time) 102 918.7080 116.993627 0.00000000 698.6438
#options(width=100) # set page width to 100 characters
#sink("results.table.txt") # capture screen output to file
#print.marklist(GRSP.results.temporal) # send output to file
#sink() # return output to screen
#system("notepad results.table.txt",invisible=FALSE)
run a model set looking at structural landscape variables:
run.GRSP.environmental =function()
{
# Daily survival rate (DSR) is constant
t1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~1)))
#Edge
E1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge)))
# Tree
E2=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree)))
# Forb
E3=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb)))
# Height
E4=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height)))
# Forb+Height
E5=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height)))
# Forb*Height
E6=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*Height)))
#Edge+AbsForb
E7=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+AbsForb)))
#Edge*AbsForb
E8=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*AbsForb)))
#Edge+Height
E9=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+Height)))
#Edge*Height
E10=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Height)))
#Tree *Edge
E11=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Edge)))
#Tree +Edge
E12=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Edge)))
#Tree+AbsForb
E13=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb)))
#Tree*AbsForb
E14=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb)))
#Tree+Height
E15=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height)))
#Tree*Height
E16=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height)))
#Edge+Height+AbsForb
E17=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+Height+AbsForb)))
#Edge*Height+AbsForb
E18=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+Height*AbsForb)))
#Edge*Height*AbsForb
E19=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Height*AbsForb)))
#Edge+AbsForb+Height
E20=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+AbsForb+Height)))
#Edge*AbsForb+Height
E21=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*AbsForb+Height)))
#Edge+AbsForb*Height
E22=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+AbsForb*Height)))
#Tree+AbsForb
E23=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb)))
#Tree*AbsForb
E24=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb)))
#Tree+Height
E25=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height)))
#Tree*Height
E26=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height)))
#Tree+Height+AbsForb
E27=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height+AbsForb)))
#Tree*Height+AbsForb
E28=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height+AbsForb)))
#Tree+Height*AbsForb
E29=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height*AbsForb)))
#Tree*AbsForb+Height*Edge
E30=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb+Height*Edge)))
#Tree*Height*AbsForb
E31=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+Height)))
#Tree+AbsForb+Edge
E32=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+Edge)))
#Tree*AbsForb+Edge
E33=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb+Edge)))
#Tree+Height+Edge
E34=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height+Edge)))
#Tree*Height+Edge
E35=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height+Edge)))
#Tree+Height+AbsForb+Edge
E36=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height+AbsForb+Edge)))
#Tree*Height+AbsForb+Edge
E37=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height+AbsForb+Edge)))
#Tree+Height*AbsForb+Edge
E38=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height*AbsForb+Edge)))
#Tree*Height*AbsForb+Edge
E39=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height*AbsForb+Edge)))
#Tree+AbsForb+Height*Edge
E40=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+Height*Edge)))
#Tree*AbsForb+Height+Edge
E41=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb+Height+Edge)))
#Tree+AbsForb*Edge
E42=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*Edge)))
#Tree*AbsForb*Edge
E43=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*AbsForb*Edge)))
#Tree+Height*Edge
E44=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height*Edge)))
#Tree*Height*Edge
E45=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height*Edge)))
#Tree+Height+AbsForb*Edge
E46=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height+AbsForb*Edge)))
#Tree*Height+AbsForb*Edge
E47=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height+AbsForb*Edge)))
#Tree+Height*AbsForb*Edge
E48=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Height*AbsForb*Edge)))
#Tree*Height*AbsForb*Edge
E49=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Height*AbsForb*Edge)))
return(collect.models())
}
GRSP.results.environmental=run.GRSP.environmental() # This runs the models above and takes a minute or 2
show model table:
GRSP.results.environmental # print model selection table to screen
## model npar AICc DeltaAICc
## 23 S(~AbsForb) 2 807.2502 0.000000000
## 50 S(~1) 1 807.2571 0.006936166
## 12 S(~Tree) 2 807.6950 0.444800000
## 5 S(~Tree + AbsForb) 3 807.8060 0.555874842
## 16 S(~Tree + AbsForb) 3 807.8060 0.555874842
## 45 S(~AbsForb + Height) 3 809.0104 1.760264842
## 34 S(~Height) 2 809.1086 1.858420000
## 1 S(~Edge) 2 809.2121 1.961910000
## 47 S(~Edge + AbsForb) 3 809.2499 1.999774842
## 7 S(~Tree + Height) 3 809.3479 2.097714842
## 18 S(~Tree + Height) 3 809.3479 2.097714842
## 20 S(~Tree + Height + AbsForb) 4 809.3585 2.108316775
## 25 S(~Tree + AbsForb + Height) 4 809.3585 2.108316775
## 26 S(~Tree + AbsForb + Edge) 4 809.5072 2.257066775
## 4 S(~Tree + Edge) 3 809.6477 2.397554842
## 6 S(~Tree * AbsForb) 4 809.7846 2.534436775
## 17 S(~Tree * AbsForb) 4 809.7846 2.534436775
## 3 S(~Tree * Edge) 4 810.3551 3.104996775
## 48 S(~Edge * AbsForb) 4 810.5632 3.313056775
## 37 S(~Tree + AbsForb * Edge) 5 810.8221 3.571951899
## 46 S(~AbsForb * Height) 4 810.8556 3.605486775
## 8 S(~Tree * Height) 4 810.8872 3.637006775
## 19 S(~Tree * Height) 4 810.8872 3.637006775
## 30 S(~Tree + Height + AbsForb + Edge) 5 810.9853 3.735161899
## 9 S(~Edge + Height + AbsForb) 4 811.0101 3.759966775
## 13 S(~Edge + AbsForb + Height) 4 811.0101 3.759966775
## 22 S(~Tree + Height * AbsForb) 5 811.0513 3.801161899
## 49 S(~Edge + Height) 3 811.0685 3.818354842
## 21 S(~Tree * Height + AbsForb) 5 811.1146 3.864451899
## 28 S(~Tree + Height + Edge) 4 811.2774 4.027236775
## 27 S(~Tree * AbsForb + Edge) 5 811.4638 4.213641899
## 35 S(~Tree + AbsForb + Height * Edge) 6 811.6379 4.387756333
## 2 S(~Edge * Height) 4 811.6432 4.393096775
## 39 S(~Tree + Height * Edge) 5 811.8726 4.622431899
## 41 S(~Tree + Height + AbsForb * Edge) 6 812.1779 4.927756333
## 14 S(~Edge * AbsForb + Height) 5 812.2501 4.999991899
## 32 S(~Tree + Height * AbsForb + Edge) 6 812.6397 5.389546333
## 31 S(~Tree * Height + AbsForb + Edge) 6 812.6756 5.425436333
## 29 S(~Tree * Height + Edge) 5 812.7637 5.513591899
## 10 S(~Edge + Height * AbsForb) 5 812.8572 5.607021899
## 15 S(~Edge + AbsForb * Height) 5 812.8572 5.607021899
## 36 S(~Tree * AbsForb + Height + Edge) 6 812.8784 5.628256333
## 24 S(~Tree * AbsForb + Height * Edge) 7 813.5177 6.267526211
## 42 S(~Tree * Height + AbsForb * Edge) 7 813.8788 6.628606211
## 38 S(~Tree * AbsForb * Edge) 8 814.6348 7.384687686
## 11 S(~Edge * Height * AbsForb) 8 816.5129 9.262717686
## 43 S(~Tree + Height * AbsForb * Edge) 9 816.5663 9.316176930
## 40 S(~Tree * Height * Edge) 8 816.7952 9.545047686
## 33 S(~Tree * Height * AbsForb + Edge) 9 817.9903 10.740196930
## 44 S(~Tree * Height * AbsForb * Edge) 16 827.2026 19.952471175
## weight Deviance
## 23 8.567615e-02 803.2416
## 50 8.537953e-02 805.2542
## 12 6.859190e-02 803.6864
## 5 6.488634e-02 801.7890
## 16 6.488634e-02 801.7890
## 45 3.553230e-02 802.9934
## 34 3.383056e-02 805.1001
## 1 3.212452e-02 805.2035
## 47 3.152204e-02 803.2329
## 7 3.001559e-02 803.3308
## 18 3.001559e-02 803.3308
## 20 2.985690e-02 801.3300
## 25 2.985690e-02 801.3300
## 26 2.771687e-02 801.4788
## 4 2.583673e-02 803.6306
## 6 2.412759e-02 801.7561
## 17 2.412759e-02 801.7561
## 3 1.813921e-02 802.3267
## 48 1.634703e-02 802.5348
## 37 1.436218e-02 800.7794
## 46 1.412337e-02 802.8272
## 8 1.390253e-02 802.8587
## 19 1.390253e-02 802.8587
## 30 1.323670e-02 800.9426
## 9 1.307355e-02 802.9817
## 13 1.307355e-02 802.9817
## 22 1.280702e-02 801.0086
## 49 1.269740e-02 805.0515
## 21 1.240809e-02 801.0719
## 28 1.143817e-02 803.2489
## 27 1.042028e-02 801.4211
## 35 9.551482e-03 799.5781
## 2 9.526011e-03 803.6148
## 39 8.493987e-03 801.8299
## 41 7.291405e-03 800.1181
## 14 7.032755e-03 802.2074
## 32 5.788084e-03 800.5799
## 31 5.685143e-03 800.6158
## 29 5.439997e-03 802.7210
## 10 5.191712e-03 802.8145
## 15 5.191712e-03 802.8145
## 36 5.136882e-03 800.8186
## 24 3.731504e-03 799.4379
## 42 3.115131e-03 799.7989
## 38 2.134496e-03 798.5321
## 11 8.346153e-04 800.4102
## 43 8.126019e-04 798.4379
## 40 7.247349e-04 800.6925
## 33 3.987088e-04 799.8619
## 44 3.983234e-06 794.8124
#options(width=100) # set page width to 100 characters
#sink("results.table.txt") # capture screen output to file
#print.marklist(GRSP.results.environmental) # send output to file
#sink() # return output to screen
#system("notepad results.table.txt",invisible=FALSE)
Run a set of models looking at the effect of density:
run.GRSP.density =function()
{
# Daily survival rate (DSR) is constant
d1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~1)))
# Year
d2=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200)))
# Year
d3=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200+I(log200^2))))
return(collect.models())
}
GRSP.results.density=run.GRSP.density() # This runs the models above and takes a minute or 2
##
## Note: only 2 parameters counted of 3 specified parameters
## AICc and parameter count have been adjusted upward
Show model output:
GRSP.results.density # print model selection table to screen
## model npar AICc DeltaAICc weight Deviance
## 2 S(~log200) 2 807.2367 0.00000000 0.4244155 803.2282
## 1 S(~1) 1 807.2571 0.02036617 0.4201156 805.2542
## 3 S(~log200 + I(log200^2)) 3 809.2453 2.00853484 0.1554689 803.2282
#options(width=100) # set page width to 100 characters
#sink("results.table.txt") # capture screen output to file
#print.marklist(GRSP.results.density) # send output to file
#sink() # return output to screen
#system("notepad results.table.txt",invisible=FALSE)
Show the beta predictions for the density model:
GRSP.results.density$d2$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.1225429 0.0985619 1.9293615 2.3157243
## S:log200 -0.2427605 0.1760618 -0.5878416 0.1023207
Run the full combined model set:
run.GRSP.combined =function()
{
# Daily survival rate (DSR) is constant
t1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~1)))
# Year
c1=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Year)), groups=c("Year"))
#Year+Date
c2=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Year+Time)), groups=c("Year"))
#Daily Survival rate (DSR) varies by density (200)
c3=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200)))
#Edge
c4=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge)))
# Tree
c5=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree)))
# Forb
c6=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb)))
# Height
c7=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height)))
# Forb+Height
c8=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height)))
#Tree+AbsForb
c9=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb)))
#Edge+log200
c10=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+log200)))
# Tree+log200
c11=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+log200)))
# Forb+log200
c12=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+log200)))
# Height+log200
c13=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+log200)))
# Forb+Height+log200
c14=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+log200)))
#Tree+AbsForb+log200
c15=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+log200)))
#Edge*log200
c16=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*log200)))
# Tree*log200
c17=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*log200)))
# Forb*log200
c18=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*log200)))
# Height*log200
c19=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*log200)))
# Forb+Height*log200
c20=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*log200)))
# Forb+Height*log200
c21=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+AbsForb*log200)))
#Tree+AbsForb*log200
c22=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*log200)))
#Tree+AbsForb*log200
c23=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Tree*log200)))
#Daily Survival rate (DSR) varies by density (200) +Year
c24=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200+Year)), groups=c("Year"))
#Edge+Year
c25=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+Year)), groups=c("Year"))
# Tree+Year
c26=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Year)), groups=c("Year"))
# Forb+Year
c27=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Year)), groups=c("Year"))
# Height+Year
c28=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+Year)), groups=c("Year"))
# Forb+Height+Year
c29=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+Year)), groups=c("Year"))
#Tree+AbsForb+Year
c30=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+Year)), groups=c("Year"))
#Edge+log200+Year
c31=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+log200+Year)), groups=c("Year"))
# Tree+log200+Year
c32=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+log200+Year)), groups=c("Year"))
# Forb+log200+Year
c33=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+log200+Year)), groups=c("Year"))
# Height+log200+Year
c34=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+log200+Year)), groups=c("Year"))
# Forb+Height+log200+Year
c35=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+log200+Year)), groups=c("Year"))
#Tree+AbsForb+log200+Year
c36=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+log200+Year)), groups=c("Year"))
#Edge*log200+Year
c37=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*log200+Year)), groups=c("Year"))
# Tree*log200+Year
c38=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*log200+Year)), groups=c("Year"))
# Forb*log200+Year
c39=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*log200+Year)), groups=c("Year"))
# Height*log200+Year
c40=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*log200+Year)), groups=c("Year"))
# Forb+Height*log200+Year
c41=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*log200+Year)), groups=c("Year"))
# Forb+Height*log200+Year
c42=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+AbsForb*log200+Year)), groups=c("Year"))
#Tree+AbsForb*log200+Year
c43=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*log200+Year)), groups=c("Year"))
#Tree+AbsForb*log200+Year
c44=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~ AbsForb+Tree*log200+Year)), groups=c("Year"))
#Daily Survival rate (DSR) varies by density (200) *Year
c45=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200*Year)), groups=c("Year"))
#Edge*Year
c46=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Year)), groups=c("Year"))
# Tree*Year
c47=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Year)), groups=c("Year"))
# Forb*Year
c48=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*Year)), groups=c("Year"))
# Height*Year
c49=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*Year)), groups=c("Year"))
#Forb+Height*Year
c50=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*Year)), groups=c("Year"))
#Tree+AbsForb*Year
c51=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*Year)), groups=c("Year"))
#Edge+log200*Year
c52=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+log200*Year)), groups=c("Year"))
# Tree+log200*Year
c53=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+log200*Year)), groups=c("Year"))
# Forb+log200*Year
c54=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+log200*Year)), groups=c("Year"))
# Height+log200*Year
c55=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+log200*Year)), groups=c("Year"))
# Forb+Height+log200*Year
c56=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+log200*Year)), groups=c("Year"))
#Tree+AbsForb+log200*Year
c57=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+log200*Year)), groups=c("Year"))
#Edge*log200*Year
c58=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*log200*Year)), groups=c("Year"))
# Tree*log200*Year
c59=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*log200*Year)), groups=c("Year"))
# Forb*log200*Year
c60=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*log200*Year)), groups=c("Year"))
# Height*log200*Year
c61=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*log200*Year)), groups=c("Year"))
# Forb+Height*log200*Year
c62=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*log200*Year)), groups=c("Year"))
# Forb+Height*log200*Year
c63=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+AbsForb*log200*Year)), groups=c("Year"))
#Tree+AbsForb*log200*Year
c64=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*log200*Year)), groups=c("Year"))
#Tree+AbsForb*log200*Year
c65=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~ AbsForb+Tree*log200*Year)), groups=c("Year"))
#Daily Survival rate (DSR) varies by density (200) +Year+Time
c66=mark(GRSP,nocc=occ,model ="Nest",model.parameters=list(S=list(formula=~log200+Year+Time)), groups=c("Year"))
#Edge+Year+Time
c67=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+Year+Time)), groups=c("Year"))
# Tree+Year+Time
c68=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+Year+Time)), groups=c("Year"))
# Forb+Year+Time
c69=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Year+Time)), groups=c("Year"))
# Height+Year
c70=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+Year+Time)), groups=c("Year"))
# Forb+Height+Year+Time
c71=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+Year+Time)), groups=c("Year"))
#Tree+AbsForb+Year+Time
c72=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+Year+Time)), groups=c("Year"))
#Edge+log200+Year+Time
c73=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+log200+Year+Time)), groups=c("Year"))
# Tree+log200+Year+Time
c74=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+log200+Year+Time)), groups=c("Year"))
# Forb+log200+Year+Time
c75=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+log200+Year+Time)), groups=c("Year"))
# Height+log200+Year+Time
c76=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+log200+Year+Time)), groups=c("Year"))
# Forb+Height+log200+Year+Time
c77=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+log200+Year+Time)), groups=c("Year"))
#Tree+AbsForb+log200+Year+Time
c78=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+log200+Year+Time)), groups=c("Year"))
#Edge*log200+Year+Time
c79=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*log200+Year+Time)), groups=c("Year"))
# Tree*log200+Year+Time
c80=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*log200+Year+Time)), groups=c("Year"))
# Forb*log200+Year+Time
c81=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*log200+Year+Time)), groups=c("Year"))
# Height*log200+Year+Time
c82=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*log200+Year+Time)), groups=c("Year"))
# Forb+Height*log200+Year +Time
c83=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*log200+Year+Time)), groups=c("Year"))
# Forb+Height*log200+Year+Time
c84=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+AbsForb*log200+Year+Time)), groups=c("Year"))
#Tree+AbsForb*log200+Year+Time
c85=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*log200+Year+Time)), groups=c("Year"))
#Tree+AbsForb*log200+Year+Time
c86=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~ AbsForb+Tree*log200+Year+Time)), groups=c("Year"))
#Daily Survival rate (DSR) varies by density (200) *Year+Time
c87=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~log200*Year+Time)), groups=c("Year"))
#Edge*Year+Time
c88=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Year+Time)), groups=c("Year"))
# Tree*Year+Time
c89=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Year+Time)), groups=c("Year"))
# Forb*Year+Time
c90=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*Year+Time)), groups=c("Year"))
# Height*Year+Time
c91=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*Year+Time)), groups=c("Year"))
# Forb+Height*Year+Time
c92=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*Year+Time)), groups=c("Year"))
#Tree+AbsForb*Year+Time
c93=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*Year+Time)), groups=c("Year"))
#Edge*Year+log200
c94=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Year+log200)), groups=c("Year"))
# Tree*Year+log200
c95=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Year+log200)), groups=c("Year"))
# Forb*Year+log200
c96=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*Year+log200)), groups=c("Year"))
# Height*Year+log200
c97=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*Year+log200)), groups=c("Year"))
# Forb+Height*Year+log200
c98=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*Year+log200)), groups=c("Year"))
#Tree+AbsForb*Year+log200
c99=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*Year+log200)), groups=c("Year"))
#Edge+log200*Year+Time
c100=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge+log200*Year+Time)), groups=c("Year"))
# Tree+log200*Year+Time
c101=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+log200*Year+Time)), groups=c("Year"))
# Forb+log200*Year+Time
c102=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+log200*Year+Time)), groups=c("Year"))
# Height+log200*Year+Time
c103=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+log200*Year+Time)), groups=c("Year"))
# Forb+Height+log200*Year+Time
c104=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height+log200*Year+Time)), groups=c("Year"))
#Tree+AbsForb+log200*Year+Time
c105=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb+log200*Year+Time)), groups=c("Year"))
#Edge*log200*Year+Time
c106=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*log200*Year+Time)), groups=c("Year"))
# Tree*log200*Year+Time
c107=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*log200*Year+Time)), groups=c("Year"))
# Forb*log200*Year+Time
c108=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*log200*Year+Time)), groups=c("Year"))
# Height*log200*Year+Time
c109=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*log200*Year+Time)), groups=c("Year"))
# Forb+Height*log200*Year +Time
c119=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*log200*Year+Time)), groups=c("Year"))
# Forb+Height*log200*Year+Time
c110=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height+AbsForb*log200*Year+Time)), groups=c("Year"))
#Tree+AbsForb*log200*Year+Time
c111=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*log200*Year+Time)), groups=c("Year"))
#Tree+AbsForb*log200*Year+Time
c112=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~ AbsForb+Tree*log200*Year+Time)), groups=c("Year"))
#Edge*Year+log200+Time
c113=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Edge*Year+log200+Time)), groups=c("Year"))
# Tree*Year+log200+Time
c114=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree*Year+log200+Time)), groups=c("Year"))
# Forb*Year+log200+Time
c115=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb*Year+log200+Time)), groups=c("Year"))
# Height*Year+log200+Time
c116=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Height*Year+log200+Time)), groups=c("Year"))
# Forb+Height*Year+log200+Time
c117=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~AbsForb+Height*Year+log200+Time)), groups=c("Year"))
#Tree+AbsForb*Year+log200+Time
c118=mark(GRSP,nocc=occ,model="Nest",model.parameters=list(S=list(formula=~Tree+AbsForb*Year+log200+Time)), groups=c("Year"))
return(collect.models())
}
GRSP.results.combined=run.GRSP.combined() # This runs the models above and takes a minute or 2
Examine the full model table:
GRSP.results.combined # print model selection table to screen
## model npar AICc DeltaAICc
## 60 S(~log200 * Year) 4 796.1157 0.0000000
## 74 S(~Edge * log200 * Year) 8 796.9797 0.8639809
## 71 S(~Height + log200 * Year) 5 797.6031 1.4873351
## 70 S(~AbsForb + log200 * Year) 5 797.6153 1.4995351
## 106 S(~log200 * Year + Time) 5 797.9412 1.8254951
## 68 S(~Edge + log200 * Year) 5 798.0226 1.9068251
## 69 S(~Tree + log200 * Year) 5 798.0509 1.9351351
## 51 S(~Edge * log200 + Year) 5 798.7682 2.6524451
## 9 S(~Edge * log200 * Year + Time) 9 798.9515 2.8357202
## 72 S(~AbsForb + Height + log200 * Year) 6 799.0538 2.9380896
## 6 S(~Height + log200 * Year + Time) 6 799.2475 3.1317696
## 77 S(~AbsForb * log200 * Year) 8 799.3886 3.2728709
## 5 S(~AbsForb + log200 * Year + Time) 6 799.4223 3.3065996
## 73 S(~Tree + AbsForb + log200 * Year) 6 799.5585 3.4427096
## 3 S(~Edge + log200 * Year + Time) 6 799.8578 3.7420396
## 4 S(~Tree + log200 * Year + Time) 6 799.8693 3.7535496
## 97 S(~Edge * log200 + Year + Time) 6 800.3471 4.2313296
## 63 S(~AbsForb * Year) 4 800.4143 4.2985400
## 7 S(~AbsForb + Height + log200 * Year + Time) 7 800.6492 4.5334294
## 80 S(~Height + AbsForb * log200 * Year) 9 800.7221 4.6063402
## 78 S(~Height * log200 * Year) 8 800.8098 4.6940609
## 11 S(~AbsForb * log200 * Year + Time) 9 801.3233 5.2075402
## 81 S(~Tree + AbsForb * log200 * Year) 9 801.3559 5.2401602
## 8 S(~Tree + AbsForb + log200 * Year + Time) 7 801.3585 5.2427394
## 1 S(~Year) 2 801.7144 5.5986632
## 110 S(~AbsForb * Year + Time) 5 801.7514 5.6356151
## 67 S(~Tree + AbsForb * Year) 5 801.7714 5.6556551
## 116 S(~AbsForb * Year + log200) 5 801.9825 5.8667151
## 79 S(~AbsForb + Height * log200 * Year) 9 802.4838 6.3680202
## 14 S(~Height + AbsForb * log200 * Year + Time) 10 802.4962 6.3804434
## 75 S(~Tree * log200 * Year) 8 802.6181 6.5023409
## 12 S(~Height * log200 * Year + Time) 9 802.6380 6.5222602
## 32 S(~Year + Time) 3 802.7273 6.6115581
## 39 S(~Tree + Year) 3 802.7635 6.6477981
## 37 S(~log200 + Year) 3 802.9159 6.8001581
## 113 S(~Tree + AbsForb * Year + Time) 6 803.1088 6.9930196
## 52 S(~Tree * log200 + Year) 5 803.1372 7.0214051
## 38 S(~Edge + Year) 3 803.2637 7.1479981
## 15 S(~Tree + AbsForb * log200 * Year + Time) 10 803.2863 7.1705634
## 119 S(~Tree + AbsForb * Year + log200) 6 803.4454 7.3296196
## 41 S(~Height + Year) 3 803.4471 7.3313581
## 19 S(~AbsForb * Year + log200 + Time) 6 803.4657 7.3499396
## 40 S(~AbsForb + Year) 3 803.5476 7.4318481
## 85 S(~Tree + Year + Time) 4 803.7895 7.6737300
## 61 S(~Edge * Year) 4 803.7962 7.6804800
## 82 S(~AbsForb + Tree * log200 * Year) 9 803.9691 7.8533902
## 64 S(~Height * Year) 4 804.0800 7.9643000
## 46 S(~Tree + log200 + Year) 4 804.1193 8.0035700
## 88 S(~Height + Year + Time) 4 804.1455 8.0297500
## 99 S(~Tree * log200 + Year + Time) 6 804.1644 8.0486796
## 83 S(~log200 + Year + Time) 4 804.1742 8.0585000
## 23 S(~AbsForb + Height * log200 * Year + Time) 10 804.2737 8.1579834
## 84 S(~Edge + Year + Time) 4 804.3430 8.2272300
## 10 S(~Tree * log200 * Year + Time) 9 804.3648 8.2491002
## 48 S(~Height + log200 + Year) 4 804.4899 8.3741200
## 45 S(~Edge + log200 + Year) 4 804.5037 8.3879700
## 86 S(~AbsForb + Year + Time) 4 804.5512 8.4354200
## 44 S(~Tree + AbsForb + Year) 4 804.5856 8.4698100
## 114 S(~Edge * Year + log200) 5 804.7416 8.6258351
## 62 S(~Tree * Year) 4 804.7669 8.6511200
## 47 S(~AbsForb + log200 + Year) 4 804.8557 8.7399600
## 107 S(~Edge * Year + Time) 5 804.8986 8.7828551
## 22 S(~Tree + AbsForb * Year + log200 + Time) 7 804.9162 8.8004894
## 59 S(~AbsForb + Tree * log200 + Year) 6 804.9724 8.8566796
## 42 S(~AbsForb + Height + Year) 4 805.2564 9.1406500
## 111 S(~Height * Year + Time) 5 805.2600 9.1442151
## 117 S(~Height * Year + log200) 5 805.2758 9.1600551
## 92 S(~Tree + log200 + Year + Time) 5 805.3753 9.2595651
## 94 S(~Height + log200 + Year + Time) 5 805.4401 9.3243951
## 90 S(~Tree + AbsForb + Year + Time) 5 805.6047 9.4889351
## 28 S(~Edge * log200) 4 805.6451 9.5293400
## 16 S(~AbsForb + Tree * log200 * Year + Time) 10 805.6807 9.5649034
## 108 S(~Tree * Year + Time) 5 805.7972 9.6814551
## 91 S(~Edge + log200 + Year + Time) 5 805.8121 9.6963051
## 89 S(~AbsForb + Height + Year + Time) 5 805.9179 9.8021351
## 29 S(~Tree * log200) 4 805.9292 9.8134600
## 66 S(~AbsForb + Height * Year) 5 805.9320 9.8162651
## 105 S(~AbsForb + Tree * log200 + Year + Time) 7 805.9574 9.8416094
## 50 S(~Tree + AbsForb + log200 + Year) 5 806.0469 9.9311751
## 93 S(~AbsForb + log200 + Year + Time) 5 806.0963 9.9805051
## 17 S(~Edge * Year + log200 + Time) 6 806.1211 10.0053796
## 115 S(~Tree * Year + log200) 5 806.1260 10.0102451
## 55 S(~Height * log200 + Year) 5 806.3861 10.2703551
## 49 S(~AbsForb + Height + log200 + Year) 5 806.4199 10.3041351
## 36 S(~AbsForb + Tree * log200) 5 806.5306 10.4148051
## 20 S(~Height * Year + log200 + Time) 6 806.6338 10.5180696
## 53 S(~AbsForb * log200 + Year) 5 806.8153 10.6995751
## 112 S(~AbsForb + Height * Year + Time) 6 807.0749 10.9591496
## 118 S(~AbsForb + Height * Year + log200) 6 807.2275 11.1117496
## 43 S(~log200) 2 807.2367 11.1209732
## 76 S(~AbsForb) 2 807.2502 11.1344032
## 120 S(~1) 1 807.2571 11.1413394
## 96 S(~Tree + AbsForb + log200 + Year + Time) 6 807.2833 11.1675496
## 101 S(~Height * log200 + Year + Time) 6 807.3284 11.2126696
## 95 S(~AbsForb + Height + log200 + Year + Time) 6 807.3339 11.2181796
## 18 S(~Tree * Year + log200 + Time) 6 807.3470 11.2312596
## 65 S(~Tree) 2 807.6950 11.5792032
## 109 S(~Tree + AbsForb) 3 807.8060 11.6902781
## 24 S(~AbsForb + log200) 3 807.9741 11.8583281
## 58 S(~Tree + AbsForb * log200 + Year) 6 808.0006 11.8848196
## 100 S(~AbsForb * log200 + Year + Time) 6 808.0524 11.9366396
## 13 S(~Tree + log200) 3 808.1134 11.9976481
## 56 S(~AbsForb + Height * log200 + Year) 6 808.3196 12.2038796
## 57 S(~Height + AbsForb * log200 + Year) 6 808.3906 12.2748896
## 21 S(~AbsForb + Height * Year + log200 + Time) 7 808.5543 12.4385794
## 25 S(~Height + log200) 3 808.8522 12.7364481
## 27 S(~Tree + AbsForb + log200) 4 808.8611 12.7453800
## 98 S(~AbsForb + Height) 3 809.0104 12.8946681
## 87 S(~Height) 2 809.1086 12.9928232
## 2 S(~Edge + log200) 3 809.2029 13.0871481
## 54 S(~Edge) 2 809.2121 13.0963132
## 102 S(~AbsForb + Height * log200 + Year + Time) 7 809.2246 13.1088694
## 104 S(~Tree + AbsForb * log200 + Year + Time) 7 809.2320 13.1162694
## 103 S(~Height + AbsForb * log200 + Year + Time) 7 809.3025 13.1867994
## 26 S(~AbsForb + Height + log200) 4 809.5347 13.4189600
## 30 S(~AbsForb * log200) 4 809.9347 13.8189900
## 31 S(~Height * log200) 4 810.5434 14.4276600
## 35 S(~Tree + AbsForb * log200) 5 810.8076 14.6918851
## 33 S(~AbsForb + Height * log200) 5 811.2639 15.1481051
## 34 S(~Height + AbsForb * log200) 5 811.5086 15.3928751
## weight Deviance
## 60 1.308461e-01 788.0873
## 74 8.494733e-02 780.8770
## 71 6.219995e-02 787.5604
## 70 6.182169e-02 787.5726
## 106 5.252421e-02 787.8985
## 68 5.043116e-02 787.9799
## 69 4.972234e-02 788.0082
## 51 3.473679e-02 788.7255
## 9 3.169509e-02 780.8230
## 72 3.011360e-02 786.9940
## 6 2.733416e-02 787.1877
## 77 2.547217e-02 783.2859
## 5 2.504620e-02 787.3625
## 73 2.339839e-02 787.4986
## 3 2.014593e-02 787.7980
## 4 2.003032e-02 787.8095
## 97 1.577391e-02 788.2872
## 63 1.525263e-02 792.3858
## 7 1.356248e-02 786.5693
## 80 1.307696e-02 782.5936
## 78 1.251579e-02 784.7071
## 11 9.681838e-03 783.1948
## 81 9.525208e-03 783.2274
## 8 9.512932e-03 787.2787
## 1 7.962079e-03 797.7059
## 110 7.816323e-03 791.7087
## 67 7.738395e-03 791.7287
## 116 6.963375e-03 791.9398
## 79 5.419544e-03 784.3553
## 14 5.385985e-03 782.3391
## 75 5.067519e-03 786.5154
## 12 5.017299e-03 784.5095
## 32 4.798209e-03 796.7102
## 39 4.712049e-03 796.7465
## 37 4.366417e-03 796.8989
## 113 3.965025e-03 791.0489
## 52 3.909147e-03 793.0945
## 38 3.669380e-03 797.2467
## 15 3.628213e-03 783.1292
## 119 3.350843e-03 791.3855
## 41 3.347932e-03 797.4301
## 19 3.316971e-03 791.4059
## 40 3.183871e-03 797.5305
## 85 2.821184e-03 795.7610
## 61 2.811679e-03 795.7678
## 82 2.578807e-03 785.8407
## 64 2.439692e-03 796.0516
## 46 2.392256e-03 796.0909
## 88 2.361145e-03 796.1170
## 99 2.338903e-03 792.1046
## 83 2.327447e-03 796.1458
## 23 2.214508e-03 784.1166
## 84 2.139146e-03 796.3145
## 10 2.115882e-03 786.2364
## 48 1.987667e-03 796.4614
## 45 1.973950e-03 796.4753
## 86 1.927670e-03 796.5227
## 44 1.894807e-03 796.5571
## 114 1.752607e-03 794.6989
## 62 1.730589e-03 796.7384
## 47 1.655399e-03 796.8273
## 107 1.620272e-03 794.8559
## 22 1.606049e-03 790.8364
## 59 1.561555e-03 792.9126
## 42 1.354858e-03 797.2279
## 111 1.352445e-03 795.2173
## 117 1.341776e-03 795.2331
## 92 1.276650e-03 795.3326
## 94 1.235931e-03 795.3974
## 90 1.138321e-03 795.5620
## 28 1.115555e-03 797.6166
## 16 1.095894e-03 785.5235
## 108 1.033855e-03 795.7545
## 91 1.026207e-03 795.7694
## 89 9.733167e-04 795.8752
## 29 9.678209e-04 797.9008
## 66 9.664645e-04 795.8893
## 105 9.542946e-04 791.8775
## 50 9.125013e-04 796.0042
## 93 8.902698e-04 796.0535
## 17 8.792659e-04 794.0613
## 115 8.771294e-04 796.0833
## 55 7.701610e-04 796.3434
## 49 7.572622e-04 796.3772
## 36 7.164973e-04 796.4878
## 20 6.804418e-04 794.5740
## 53 6.214090e-04 796.7726
## 112 5.457726e-04 795.0151
## 118 5.056791e-04 795.1677
## 43 5.033524e-04 803.2282
## 76 4.999837e-04 803.2416
## 120 4.982527e-04 805.2542
## 96 4.917657e-04 795.2235
## 101 4.807956e-04 795.2686
## 95 4.794729e-04 795.2741
## 18 4.763474e-04 795.2872
## 65 4.002845e-04 803.6864
## 109 3.786598e-04 801.7890
## 24 3.481429e-04 801.9570
## 58 3.435619e-04 795.9407
## 100 3.347745e-04 795.9926
## 13 3.247167e-04 802.0963
## 56 2.929018e-04 796.2598
## 57 2.826847e-04 796.3308
## 21 2.604699e-04 794.4745
## 25 2.244276e-04 802.8351
## 27 2.234275e-04 800.8327
## 98 2.073572e-04 802.9934
## 87 1.974263e-04 805.1001
## 2 1.883314e-04 803.1858
## 54 1.874703e-04 805.2035
## 102 1.862970e-04 795.1448
## 104 1.856090e-04 795.1522
## 103 1.791776e-04 795.2227
## 26 1.595404e-04 801.5063
## 30 1.306187e-04 801.9063
## 31 9.634612e-05 802.5150
## 35 8.442255e-05 800.7649
## 33 6.720345e-05 801.2211
## 34 5.946212e-05 801.4659
#options(width=100) # set page width to 100 characters
#sink("results.table.xls") # capture screen output to file
#print.marklist(GRSP.results.combined) # send output to file
#sink() # return output to screen
#system("notepad results.table.xls",invisible=FALSE)
Look at beta values for each of the top models
density*year:
GRSP.results.combined$c45$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2407644 0.1369698 1.9723036 2.5092252
## S:log200 -0.5016157 0.2333160 -0.9589150 -0.0443163
## S:Year2 -0.1737687 0.1997983 -0.5653734 0.2178359
## S:log200:Year2 1.1838319 0.3892436 0.4209145 1.9467493
density*year+edge:
GRSP.results.combined$c52$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2528458 0.1421624 1.9742075 2.5314840
## S:Edge -0.0305306 0.0928234 -0.2124645 0.1514033
## S:log200 -0.4972879 0.2339640 -0.9558572 -0.0387185
## S:Year2 -0.1967190 0.2119150 -0.6120724 0.2186344
## S:log200:Year2 1.1702727 0.3917505 0.4024416 1.9381037
density*year+height:
GRSP.results.combined$c55$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2365880 0.1372988 1.9674824 2.5056937
## S:Height 0.0676104 0.0936210 -0.1158868 0.2511076
## S:log200 -0.5277409 0.2387284 -0.9956487 -0.0598332
## S:Year2 -0.1719394 0.2002235 -0.5643775 0.2204986
## S:log200:Year2 1.1987828 0.3924315 0.4296172 1.9679485
density*year+forb:
GRSP.results.combined$c54$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2223259 0.1390573 1.9497736 2.4948783
## S:AbsForb -0.0665414 0.0917125 -0.2462979 0.1132152
## S:log200 -0.4874209 0.2332851 -0.9446597 -0.0301822
## S:Year2 -0.1182230 0.2147037 -0.5390423 0.3025963
## S:log200:Year2 1.2233253 0.3918750 0.4552503 1.9914003
density*year+date:
GRSP.results.combined$c87$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.1816934 0.1921052 1.8051672 2.5582197
## S:log200 -0.4767372 0.2376937 -0.9426169 -0.0108575
## S:Year2 -0.1897134 0.2030376 -0.5876670 0.2082402
## S:Time 0.0016480 0.0037979 -0.0057959 0.0090919
## S:log200:Year2 1.1540924 0.3937546 0.3823334 1.9258513
density*year+tree:
GRSP.results.combined$c53$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2407515 0.1369413 1.9723465 2.5091565
## S:Tree -0.0269063 0.0955351 -0.2141551 0.1603425
## S:log200 -0.4901846 0.2364342 -0.9535957 -0.0267735
## S:Year2 -0.1753146 0.1998792 -0.5670778 0.2164486
## S:log200:Year2 1.1588328 0.3991985 0.3764038 1.9412618
densityyearedge:
GRSP.results.combined$c58$results$beta
## estimate se lcl ucl
## S:(Intercept) 2.2710491 0.1467543 1.9834106 2.5586875
## S:Edge -0.0944049 0.1365072 -0.3619591 0.1731493
## S:log200 -0.5274318 0.2306500 -0.9795057 -0.0753578
## S:Year2 -0.2986059 0.2140374 -0.7181192 0.1209075
## S:Edge:log200 -0.5720414 0.2863925 -1.1333707 -0.0107122
## S:Edge:Year2 -0.2025645 0.2261790 -0.6458754 0.2407464
## S:log200:Year2 1.0856287 0.4122138 0.2776896 1.8935679
## S:Edge:log200:Year2 0.2123654 0.4331761 -0.6366597 1.0613906
Graph the density*year interaction: First, generate predictions based on the model:
#create a sequence of values within the range of the density values
min.density200=min(GRSP$log200)
max.density200=max(GRSP$log200)
avg.log200=mean(GRSP$log200)
log200.values=seq(from=min.density200, to=max.density200, length=50)
#predict nest success based on different densities for a given year
pred.density1=covariate.predictions(GRSP.results.combined$c45, data=data.frame(log200=log200.values, Time=50), indices=c(1))
#put those predictions into a vector
pred.density.1=pred.density1$estimates
#make your densities "real" instead of logs
density.values= (10^log200.values)+0.01
#make a data frame to graph
plot.density.1=data.frame(pred.density.1$estimate, density.values, pred.density.1$lcl, pred.density.1$ucl, pred.density.1$se)
#do the thing for year 2
pred.density2=covariate.predictions(GRSP.results.combined$c45, data=data.frame(log200=log200.values, Time=50), indices=c(103))
pred.density.2=pred.density2$estimates
plot.density.2=data.frame(pred.density.2$estimate, density.values, pred.density.2$lcl, pred.density.2$ucl, pred.density.2$se)
then graph the predictions for the density*year interaction:
#plot the relationship between density and year
densityplotyr1=ggplot()+geom_line(data=plot.density.1, aes(x=density.values, y=pred.density.1$estimate, color="2014"), size=1)+
xlab("Male Densities (territories/ha)")+ylab("Estimated Daily Nest Survival ± 95% CL")+
geom_line(data=plot.density.2, aes(x=density.values, y=pred.density.2$estimate, colour="2015"), size=1, linetype=2)+
geom_ribbon (aes(x=density.values, ymin=pred.density.2$lcl, ymax=pred.density.2$ucl, alpha=0.2))+
xlim(0,2.9)+geom_ribbon(aes(x=density.values, ymin=pred.density.1$lcl, ymax=pred.density.1$ucl), alpha=0.2)+
guides(alpha=FALSE, colour=guide_legend(override.aes=list(linetype=c(1,2), fill=c(NA,NA))))+
scale_color_manual(name = "Year",
values = c("2014" = "red", "2015" = "blue"))+
theme(
legend.key=element_rect(colour = "transparent", fill="white"),
legend.position = c(0.8, 0.2),
legend.text=element_text(size=15),
legend.title=element_text(size=15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("Figure7.jpeg", width = 8, height = 5, units = 'in', res = 600)
densityplotyr1
#dev.off()
Graph the density-by-edge-by-year interaction First, generate a prediction based on the edge values, by year:
#create a sequence of values within the range of the edge values
min.edge=min(GRSP$Edge)
max.edge=max(GRSP$Edge)
avg.edge=mean(GRSPuntransformed$Edge)
std.edge=sd(GRSPuntransformed$Edge)
edge.values=seq(from=min.edge, to=max.edge, length=50)
edge.values_transformed= (edge.values*std.edge)+avg.edge
#predict nest success based on different edge values for a given year
pred.edge1=covariate.predictions(GRSP.results.combined$c58, data=data.frame(Edge=edge.values, Time=50), indices=c(1))
#put those predictions into a vector
pred.edge.1=pred.edge1$estimates
#make a data frame to graph
plot.edge.1=data.frame(pred.edge.1$estimate, edge.values_transformed, pred.edge.1$lcl, pred.edge.1$ucl, pred.edge.1$se)
#do the thing for year 2
pred.edge2=covariate.predictions(GRSP.results.combined$c58, data=data.frame(Edge=edge.values, Time=50), indices=c(103))
pred.edge.2=pred.edge2$estimates
plot.edge.2=data.frame(pred.edge.2$estimate, edge.values_transformed, pred.edge.2$lcl, pred.edge.2$ucl, pred.edge.2$se)
Then graph the year*edge interaction
#plot the relationship between edge and year
edgeplotyr1=ggplot()+geom_line(data=plot.edge.1, aes(x=edge.values_transformed, y=pred.edge.1$estimate, color="2014"), size=1)+
xlab("Distance to nearest edge (m)")+ylab("Estimated Daily Nest Survival ± 95% CL")+
geom_line(data=plot.edge.2, aes(x=edge.values_transformed, y=pred.edge.2$estimate, colour="2015"), size=1, linetype=2)+
geom_ribbon (aes(x=edge.values_transformed, ymin=pred.edge.2$lcl, ymax=pred.edge.2$ucl, alpha=0.2))+
xlim(0,450)+geom_ribbon(aes(x=edge.values_transformed, ymin=pred.edge.1$lcl, ymax=pred.edge.1$ucl), alpha=0.2)+
guides(alpha=FALSE, colour=guide_legend(override.aes=list(linetype=c(1,2), fill=c(NA,NA))))+
scale_color_manual(name = "Year",
values = c("2014" = "red", "2015" = "blue"))+
theme(
legend.key=element_rect(colour = "transparent", fill="white"),
legend.position = c(0.3, 0.2),
legend.text=element_text(size=15),
legend.title=element_text(size=15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("s4.jpeg", width = 8, height = 5, units = 'in', res = 600)
edgeplotyr1
#dev.off()
Then make predictions based on density values, holding edge constant:
#predict nest success based on different density values for a given edge variable
pred.edge3=covariate.predictions(GRSP.results.combined$c58, data=data.frame(Edge=min.edge, log200=log200.values), indices=c(1))
#put those predictions into a vector
pred.edge.3=pred.edge3$estimates
#make a data frame to graph
plot.edge.3=data.frame(pred.edge.3$estimate, density.values, pred.edge.3$lcl, pred.edge.3$ucl, pred.edge.3$se)
#do the thing for year 2
pred.edge4=covariate.predictions(GRSP.results.combined$c58, data=data.frame(Edge=max.edge, log200=log200.values), indices=c(1))
pred.edge.4=pred.edge4$estimates
plot.edge.4=data.frame(pred.edge.4$estimate, density.values, pred.edge.4$lcl, pred.edge.4$ucl, pred.edge.4$se)
Then graph the density*edge interaction:
#plot the relationship between edge and density, holding edge constant
edgeplot1=ggplot()+geom_line(data=plot.edge.3, aes(x=density.values, y=pred.edge.3$estimate, color="Min (0.03m)"), size=1)+
xlab("Density (territories/ha)")+ylab("Estimated Daily Nest Survival ± 95% CL")+
geom_line(data=plot.edge.4, aes(x=density.values, y=pred.edge.4$estimate, colour="Max (448m)"), size=1, linetype=2)+
geom_ribbon (aes(x=density.values, ymin=pred.edge.4$lcl, ymax=pred.edge.4$ucl, alpha=0.2))+
xlim(0,2.9)+geom_ribbon(aes(x=density.values, ymin=pred.edge.3$lcl, ymax=pred.edge.3$ucl), alpha=0.2)+
guides(alpha=FALSE, colour=guide_legend(override.aes=list(linetype=c(2,1), fill=c(NA,NA))))+
scale_color_manual(name = "Distance to edge",
values = c("Min (0.03m)" = "red", "Max (448m)" = "blue"))+
theme(
legend.key=element_rect(colour = "transparent", fill="white"),
legend.position = c(0.3, 0.2),
legend.text=element_text(size=15),
legend.title=element_text(size=15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(colour = 'black', size = 15),
axis.text.y = element_text(colour = 'black', size = 15),
axis.line = element_line(colour = "black", size = 1),
axis.title.x = element_text(colour = "black", size=15),
axis.title.y = element_text(colour = "black", size=15),
panel.background = element_rect(fill="white", colour="white"),
plot.background = element_rect(fill = "white", colour="white"))
#jpeg("s5.jpeg", width = 8, height = 5, units = 'in', res = 600)
edgeplot1
#dev.off()