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()