Neighbours shapes pollination network topology

The front page of the repository containing data and additional code can be accessed here: [https://github.com/jennabraun/foundation-pollination]

Build network from visitation dataset

#dataset with aggregated visits (i.e. number of visits per rep, no pollinator ID) with joined covariates & calculated densities
#contains reps with zero visits
#wrangling all done using visits.R

data <- read.csv("Data/Output/visitation_cleaned.csv")

sum(data$Quantity)
## [1] 635
#long version with pollinator ID, list of visits only for building network 
visits <- read.csv("Data/visitation_data.csv")
visits$Quantity <- as.numeric(as.character(visits$Quantity))
#some non-pollinator interactions were also recorded but designated with quantity = 0
visits <- dplyr::filter(visits, Quantity>0)
visits$Site <- "site"
visits$ID <- gsub(" ","", visits$ID)
visits$ID <- gsub('\"',"", visits$ID)
visits$Species <- gsub(" ","", visits$Species)
visits$Species <- gsub('\"',"", visits$Species)
sum(visits$Quantity)
## [1] 635
visits$uniID <- paste(visits$Species, visits$WP)

Functions

#this function creates an adjencancy matrix from the long visit only dataset
#uniID - plant ID, can use to create species/individual/RTU/functional groups networks from dataset
#y is bottom - uniID or plant Species
#z is top - pollinator RTU or functional group
net <- function(x, y, z){
  y <- enquo(y)
  z <- enquo(z)
  fnet <- dplyr::select(x, !!y, !!z, Site, Quantity)
  fnet$Quantity <- as.numeric(fnet$Quantity)
  print(sum(fnet$Quantity))
  long.ag <- fnet %>% group_by(!!y, !!z) %>% summarise(Quantity = sum(Quantity)) 
  print(sum(long.ag$Quantity))
  wide <- spread(long.ag, !!z, Quantity)
  wide[is.na(wide)] <- 0
  wide <- as.data.frame(wide)
  rownames(wide) = wide[,1 ] # the first row will be the header
  rownames(wide) <- wide[,1]
  wide <- select(wide, -!!y)
}

#calculates z-scores
net.zscore = function(obsval, nullval) {
  (obsval - mean(nullval))/sd(nullval)  
}

#creates random network using vegan with fixed marginal totals. Creates 1000 and puts in one big dataframe
ranindices <- function(x, index){
  y <- permatfull(x, fixedmar = "both", times = 1000)
  y <- y$perm
  y_indices <- map(y, networklevel,index = index)
  ran <-  do.call("rbind", y_indices)
  ran <- as.data.frame(ran)
}




#calculates the modularity of 1000 random networks with fixed marginal totals of input and puts in a big dataframe. Bipartite modularity 
ranModules <- function(x){
  y <- permatfull(x, fixedmar = "both", times = 1000)
  y <- y$perm
  z <- list()
  null.mod <- map(y, computeModules)
  for (i in 1:1000){
    ind <- null.mod[[i]]@likelihood
    z <- c(z, ind)
  } 
  do.call(rbind, lapply(z, as.data.frame))
}

Calculations

Centrality measures

#make an individual network using functional groups
#ind1 <- net(visits, uniID, fun.grp)

ind1 <- net(visits, uniID, ID)
## [1] 635
## `summarise()` has grouped output by 'uniID'. You can override using the `.groups` argument.
## [1] 635
#output dataset for repo
write.csv(ind1, "Data/individual_network_bipartite.csv")

sum(ind1)
## [1] 635
#project into a unipartite network
one1 <- as.one.mode(ind1, fill = 0, project = "lower", weight = "TRUE")
write.csv(one1, "Data/individual_plant_onemode.csv")
#calculate centrality measures
one1 <- graph_from_adjacency_matrix(one1)


cls1 <- closeness(one1)
## Warning in closeness(one1): At centrality.c:2784 :closeness centrality is not
## well-defined for disconnected graphs
ei1 <- eigen_centrality(one1)
bt1 <- betweenness(one1)
dg1 <- degree(one1)

#cbind and join to main covariate dataset
dat1 <- as.data.frame(cbind(cls1, ei1$vector, bt1, dg1))
dat1$uniID <- row.names(dat1)
all <- left_join(dat1, data, by = "uniID")

#betweenness is zero inflated so a binary variable
all <- mutate(all, btbin = ifelse(bt1 == 0, 0, 1))




#filter out the over zero part to model separately
bt <- filter(all, bt1 > 0)

Effective partners

#calculate species level indices
sp <- specieslevel(ind1, index = "effective partners")
pl <- sp[2]
pl <- as.data.frame(pl)
pl$uniID <- row.names(pl)
all <- left_join(pl, all, by = "uniID")

Proportion of conspecific/heterospecific interactors

#unweighted network
un <- as.one.mode(ind1, fill = 0, project = "lower", weight = "FALSE")

#need to subset
un <- as.data.frame(un)
un$uniID <- row.names(un)
dat <- left_join(all, un, by = "uniID") 
dat <- select(dat, -(18:28))
#calculate number of conspecific interactors
dat <- mutate(dat, cons = ifelse(Species == "PP", rowSums(select (dat, starts_with("PP"))), +
  ifelse(Species == "SC", rowSums(select(dat, starts_with("SC"))),+    ifelse(Species == "HH", rowSums(select(dat, starts_with("HH"))),  +      ifelse(Species == "LT",rowSums(select(dat, starts_with("LT"))), +       ifelse(Species == "AS", rowSums(select(dat, starts_with("AS"))), +      ifelse(Species == "EC", rowSums(select(dat, starts_with("EC"))), +      ifelse(Species == "BW", rowSums(select(dat, starts_with("BW"))), +        ifelse(Species == "EL", rowSums(select(dat, starts_with("EL"))), +  ifelse(Species == "SM", rowSums(select(dat, starts_with("SM"))),
  + ifelse(Species == "SD", rowSums(select(dat, starts_with("SD"))), 0)))))))))))
#calculate heterospecific interactors

dat <- mutate(dat, hets = ifelse(Species == "PP", rowSums(select(dat, 37:188, 207:265)), +
   ifelse(Species == "SC", rowSums(select(dat, 256:265, 37:206)),+    
 ifelse(Species == "HH", rowSums(select(dat, 37:131, 136:265)),  +      ifelse(Species == "LT",rowSums(select(dat, 37:135, 189:265)), +       ifelse(Species == "AS", rowSums(select(dat, 78:265)), +          ifelse(Species == "EC", rowSums(select(dat, 37:96, 131:265)), +      ifelse(Species == "BW", rowSums(select(dat, 37:77, 96:265)), +        ifelse(Species == "EL", rowSums(select(dat, 37:130, 132:265)), +  ifelse(Species == "SM", rowSums(select(dat, 37:261)),
 + ifelse(Species == "SD", rowSums(select(dat, 37:255, 262:265)), 0)))))))))))


#standardize by abundances
dat <- select(dat, -(37:265))
abun <- dat %>% group_by(Species) %>% count()
dat <- left_join(dat, abun, by = "Species")
dat <- mutate(dat, hetn = 231 - n)
dat <- mutate(dat, st.con = cons/n)
dat <- mutate(dat, st.het = hets/hetn)
dat <- select(dat, uniID, 37:42)
all <- left_join(all, dat, by = "uniID")

Contribution to nestedness

#this step is slow & can be skipped w/o issues. load from csv but code is here.  takes 30 min 
# ne <- nestedcontribution(ind1)
# ne.low <- ne[[2]]
# ne.top <- ne[[1]]
# ne.low$uniID <- row.names(ne.low)
# write.csv(ne.low, "data/Output/nelow.csv")
# write.csv(ne.top, "data/Output/netop.csv")
ne.low <- read.csv("data/Output/nelow.csv")
ne.low <- select(ne.low, -X)
all <- left_join(all, ne.low, by = "uniID")

Unipartitite modularity

one1 <- as.one.mode(ind1, fill = 0, project = "lower", weighted = TRUE)
mods <- netcarto(one1)
n1data <- mods[[1]]
n1data <- rename(n1data, uniID = name)
all <- left_join(all, n1data, by = "uniID")
all$module <- as.factor(all$module)

These sections take a long time to run. Therefore, the output is stored in .csv but the code is commented out here.

Individual-based network

#make 1000 random networks, project to one-mode, compute modularity
#null.mod <- permatfull(ind1, fixedmar = "both", times = 1000)
#null.mod <- null.mod$perm
#null.one <- map(null.mod, as.one.mode, fill = 0, project = "lower", weighted = TRUE)
#null.module <- map(null.one, netcarto)

#t <- list()  
#for (i in 1:1000){
#    ind <- null.module[[i]][2]
#   t <- c(t, ind)
#  } 
#t <- unlist(t, use.names=FALSE)
#t <- as.data.frame(t)
#write.csv(t, "Data/Randomization_Outputs/randommoduleoutput_unipartite.csv")
t <- read.csv("Data/Randomization_Outputs/randommoduleoutput_unipartite.csv")
mods[[2]]
## [1] 0.5774919
net.zscore(mods[[2]], t$t)
## [1] 41.44413
mean(t$t)
## [1] 0.1564183
#print(null.module[[3]][2])

EDA

labels = c("Acamptopappus sphaerocephalus", "Eriogonum fasciculatum", "Ericameria cooperi", "Ericameria linearifolia", "Echinocereus engelmannii", "Larrea tridentata", "Opuntia basilaris", "Cylindropuntia echinocarpa", "Salvia dorri", "Scutellaria mexicana")

ggplot(data, aes(N.flowers, Quantity)) + geom_point(aes(color = Species)) + geom_smooth(method = lm, se = FALSE, color = "black") + scale_color_discrete("", labels = labels) + ylab("Visitation Rate") + xlab("Floral display size")
## `geom_smooth()` using formula 'y ~ x'

ggplot(data, aes(N.flowers, Height)) + geom_point(aes(color = Species)) + geom_smooth(method = lm, se = FALSE, color = "black") + scale_color_discrete("", labels = labels) + xlab("Floral display size")
## `geom_smooth()` using formula 'y ~ x'

ggplot(data, aes(shrub.density, Quantity)) + geom_point(aes(color = Species)) + geom_smooth(method = lm, se = FALSE, color = "black") + scale_color_discrete("", labels = labels) + ylab("Visitation rate") + xlab("Neighbourhood shrub density")
## `geom_smooth()` using formula 'y ~ x'

pl <- select(all, Species, shrub.density, dg1, st.con, st.het)
pl <- gather(pl,type, prop, 4:5)

ggplot(pl, aes(dg1, prop, color = type)) + geom_point()  +xlab("Degree Centrality") + ylab("Proportion of Population")

ggplot(pl, aes(shrub.density, prop, color = type)) + geom_point() + geom_smooth(method = lm) +xlab("Shrub Density") + ylab("Proportion of Population")
## `geom_smooth()` using formula 'y ~ x'

c <- select(all, 3:6, Quantity)
c <- rename(c, closeness = cls1, eigancentrality = V2, betweenness = bt1, degree = dg1, visitation.rate = Quantity)
M <- cor(c)
corrplot(M, method = "number") # Display the correlation coefficient

Modelling

Visitation rates

m1 <- MASS::glm.nb(Quantity ~density + N.flowers + imp.den + Species, data = data)       
m1p <- glm(Quantity ~density + N.flowers + imp.den + Species, family = "poisson", data = data)    
check_overdispersion(m1p)
## # Overdispersion test
## 
##        dispersion ratio =   2.421
##   Pearson's Chi-Squared = 922.346
##                 p-value = < 0.001
## Overdispersion detected.
AIC(m1, m1p)
##     df      AIC
## m1  14 1327.841
## m1p 13 1500.967
#negative binomial #2 is the best due to overdispersion

#correlation test
cor.test(data$shrub.density, data$cactus.density)
## 
##  Pearson's product-moment correlation
## 
## data:  data$shrub.density and data$cactus.density
## t = -3.2772, df = 392, p-value = 0.001142
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.25793733 -0.06556399
## sample estimates:
##        cor 
## -0.1633025
m1 <- MASS::glm.nb(Quantity ~density * N.flowers + day + Species, data = data)
m2 <- MASS::glm.nb(Quantity ~ density + N.flowers + day + Species,  data = data)       
m3 <- MASS::glm.nb(Quantity ~ density + N.flowers * day + Species, data = data) 
m4 <-MASS::glm.nb(Quantity ~shrub.density * N.flowers + day + Species, data = data)
m5 <- MASS::glm.nb(Quantity ~ shrub.density + N.flowers + day + Species,data = data)    
m6 <- MASS::glm.nb(Quantity ~ shrub.density + N.flowers * day + Species,  data = data)

m6.2 <- MASS::glm.nb(Quantity ~ shrub.density + Height+day + N.flowers * day + Species,  data = data)
summary(m6.2)
## 
## Call:
## MASS::glm.nb(formula = Quantity ~ shrub.density + Height + day + 
##     N.flowers * day + Species, data = data, init.theta = 1.281937207, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9349  -1.1398  -0.4666   0.3391   2.6458  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.3193294  0.2755295  -4.788 1.68e-06 ***
## shrub.density  0.0808752  0.0357399   2.263 0.023643 *  
## Height         0.0049071  0.0024917   1.969 0.048911 *  
## day            0.0454833  0.0177519   2.562 0.010402 *  
## N.flowers      0.0119291  0.0020781   5.740 9.44e-09 ***
## SpeciesBW      1.0168923  0.2639745   3.852 0.000117 ***
## SpeciesEC     -0.1579723  0.2678134  -0.590 0.555285    
## SpeciesEL     -1.3408917  1.1653249  -1.151 0.249873    
## SpeciesHH      1.1434189  0.5534674   2.066 0.038836 *  
## SpeciesLT     -0.3254614  0.3953542  -0.823 0.410386    
## SpeciesPP      1.3085603  0.2968011   4.409 1.04e-05 ***
## SpeciesSC      1.0944596  0.2064465   5.301 1.15e-07 ***
## SpeciesSD      0.1699269  0.4341951   0.391 0.695531    
## SpeciesSM     -0.2573513  0.4491798  -0.573 0.566688    
## day:N.flowers -0.0005851  0.0001533  -3.818 0.000135 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.2819) family taken to be 1)
## 
##     Null deviance: 500.16  on 393  degrees of freedom
## Residual deviance: 406.11  on 379  degrees of freedom
## AIC: 1316.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.282 
##           Std. Err.:  0.195 
## 
##  2 x log-likelihood:  -1284.081
check_collinearity(m6.2)
## Warning: Model has interaction terms. VIFs might be inflated. You may check
##   multicollinearity among predictors of a model without interaction terms.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##  shrub.density 1.19         1.09      0.84
##            day 2.64         1.62      0.38
##      N.flowers 3.93         1.98      0.25
##  day:N.flowers 2.23         1.49      0.45
## 
## Moderate Correlation
## 
##    Term  VIF Increased SE Tolerance
##  Height 9.12         3.02      0.11
## 
## High Correlation
## 
##     Term   VIF Increased SE Tolerance
##  Species 27.97         5.29      0.04
#including height increases colinearity

m6.3 <- MASS::glm.nb(Quantity ~ shrub.density + day + N.flowers * day + Species,  data = data)
AIC(m6.2, m6.3)
##      df      AIC
## m6.2 16 1316.081
## m6.3 15 1317.507
#including height increases colinearity without improving fit

m7 <- MASS::glm.nb(Quantity ~ cactus.density + N.flowers * day + Species, data = data)
m8 <- MASS::glm.nb(Quantity ~ con.density + N.flowers * day+ Species, data = data)

m9 <- MASS::glm.nb(Quantity ~ shrub.density + N.flowers * day + Species, data = data)
m10 <- MASS::glm.nb(Quantity ~ N.flowers * day + Species,  data = data)
m11 <- MASS::glm.nb(Quantity ~ density + N.flowers * day + Species,  data = data)

#m9.imp <- MASS::glm.nb(Quantity ~ shrub.density + N.flowers * imp.den + Species, data = data)
#AIC(m9.imp)


knitr::kable(anova(m9, m10, test = "Chisq"))
Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
N.flowers * day + Species 1.223626 381 -1292.470 NA NA NA
shrub.density + N.flowers * day + Species 1.248591 380 -1287.507 1 vs 2 1 4.962725 0.0258994 
knitr::kable(anova(m10, m11, test = "Chisq"))
Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
N.flowers * day + Species 1.223626 381 -1292.470 NA NA NA
density + N.flowers * day + Species 1.243864 380 -1288.315 1 vs 2 1 4.154188 0.0415317 
#adding shrub density significantly improves model fit more than density

knitr::kable(AIC(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11))
df AIC
m1 15 1327.017
m2 14 1327.794
m3 15 1318.315
m4 15 1326.006
m5 14 1326.763
m6 15 1317.507
m7 15 1321.756
m8 15 1318.881
m9 15 1317.507
m10 14 1320.470
m11 15 1318.315 
#anova(m6, mnull)
library(tab)
## Loading required package: knitr
summary(m6)
## 
## Call:
## MASS::glm.nb(formula = Quantity ~ shrub.density + N.flowers * 
##     day + Species, data = data, init.theta = 1.248590687, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9552  -1.1597  -0.4494   0.3880   2.5977  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.1099982  0.2493269  -4.452 8.51e-06 ***
## shrub.density  0.0832383  0.0359451   2.316 0.020574 *  
## N.flowers      0.0129446  0.0020718   6.248 4.16e-10 ***
## day            0.0446048  0.0177417   2.514 0.011933 *  
## SpeciesBW      1.0181016  0.2656378   3.833 0.000127 ***
## SpeciesEC     -0.2646716  0.2682290  -0.987 0.323771    
## SpeciesEL     -1.3147551  1.1813828  -1.113 0.265754    
## SpeciesHH      1.0816165  0.5552322   1.948 0.051410 .  
## SpeciesLT      0.2841706  0.2186900   1.299 0.193799    
## SpeciesPP      1.2219151  0.2936316   4.161 3.16e-05 ***
## SpeciesSC      1.1636821  0.2054516   5.664 1.48e-08 ***
## SpeciesSD      0.2213435  0.4345732   0.509 0.610517    
## SpeciesSM     -0.1294727  0.4452333  -0.291 0.771206    
## N.flowers:day -0.0005800  0.0001531  -3.789 0.000151 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.2486) family taken to be 1)
## 
##     Null deviance: 494.61  on 393  degrees of freedom
## Residual deviance: 405.17  on 380  degrees of freedom
## AIC: 1317.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.249 
##           Std. Err.:  0.187 
## 
##  2 x log-likelihood:  -1287.507
tabglm(m6)
## Waiting for profiling to be done...
Variable Beta (SE) 95% CI P
Intercept -1.11 (0.25) (-1.63, -0.61) <0.001
shrub.density 0.08 (0.04) (0.01, 0.16) 0.02
N.flowers 0.01 (0.00) (0.01, 0.02) <0.001
day 0.04 (0.02) (0.01, 0.08) 0.01
Species
     AS (ref) – – –
     BW 1.02 (0.27) (0.50, 1.54) <0.001
     EC -0.26 (0.27) (-0.80, 0.27) 0.32
     EL -1.31 (1.18) (-4.31, 0.61) 0.27
     HH 1.08 (0.56) (-0.02, 2.23) 0.05
     LT 0.28 (0.22) (-0.15, 0.72) 0.19
     PP 1.22 (0.29) (0.64, 1.81) <0.001
     SC 1.16 (0.21) (0.76, 1.57) <0.001
     SD 0.22 (0.43) (-0.63, 1.06) 0.61
     SM -0.13 (0.45) (-1.03, 0.73) 0.77
N.flowers by day -0.00 (0.00) (-0.00, -0.00) <0.001 
knitr::kable(car::Anova(m6, type = 3))
LR Chisq Df Pr(>Chisq)
shrub.density 4.980705 1 0.0256315
N.flowers 31.690178 1 0.0000000
day 6.119279 1 0.0133715
Species 56.752376 9 0.0000000
N.flowers:day 11.399443 1 0.0007347 
check_collinearity(m6)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##  shrub.density 1.18         1.09      0.84
##      N.flowers 3.74         1.93      0.27
##            day 2.31         1.52      0.43
##        Species 3.49         1.87      0.29
##  N.flowers:day 2.34         1.53      0.43
#outputting "data" file to .csv for publication

Centrality - degree

#degree centrality

#choosing between error families
ggplot(all, aes(dg1)) + geom_density()

m2 <- MASS::glm.nb(dg1 ~ N.flowers+ Quantity + shrub.density + day + Species, data = all)
m3 <- glm(dg1 ~ N.flowers+ Quantity + shrub.density + day + Species, family = "poisson", data = all)
check_overdispersion(m3)
## # Overdispersion test
## 
##        dispersion ratio =   22.668
##   Pearson's Chi-Squared = 4919.051
##                 p-value =  < 0.001
## Overdispersion detected.
AIC(m2, m3)
##    df      AIC
## m2 15 2368.756
## m3 14 6834.353
#nbinom2
summary(m2)
## 
## Call:
## MASS::glm.nb(formula = dg1 ~ N.flowers + Quantity + shrub.density + 
##     day + Species, data = all, init.theta = 2.171691063, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4171  -0.7540  -0.0473   0.4566   2.6632  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.046e+00  1.625e-01  18.750  < 2e-16 ***
## N.flowers      7.434e-04  8.197e-04   0.907  0.36444    
## Quantity       2.383e-01  2.116e-02  11.264  < 2e-16 ***
## shrub.density  4.515e-03  2.797e-02   0.161  0.87175    
## day           -1.336e-02  9.628e-03  -1.387  0.16534    
## SpeciesBW      1.893e-01  2.087e-01   0.907  0.36447    
## SpeciesEC      1.683e-01  1.824e-01   0.922  0.35630    
## SpeciesEL     -3.466e+01  3.803e+06   0.000  0.99999    
## SpeciesHH      1.096e+00  3.633e-01   3.017  0.00256 ** 
## SpeciesLT      7.743e-01  1.601e-01   4.836 1.32e-06 ***
## SpeciesPP      1.069e+00  2.087e-01   5.122 3.02e-07 ***
## SpeciesSC      1.193e+00  1.546e-01   7.716 1.20e-14 ***
## SpeciesSD     -1.758e-01  3.126e-01  -0.562  0.57384    
## SpeciesSM     -4.214e-01  4.192e-01  -1.005  0.31484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1717) family taken to be 1)
## 
##     Null deviance: 511.16  on 230  degrees of freedom
## Residual deviance: 260.67  on 217  degrees of freedom
## AIC: 2368.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.172 
##           Std. Err.:  0.211 
## 
##  2 x log-likelihood:  -2338.756
#check if interaction should be included
m7 <- MASS::glm.nb(dg1 ~ Quantity+ shrub.density + day  * N.flowers + Species, data = all)
AIC(m2, m7)
##    df      AIC
## m2 15 2368.756
## m7 16 2369.264
#mnull <- glmmTMB(dg1 ~  (1|Species), family = "nbinom2", data = all)



#model without interaction is best
summary(m2)
## 
## Call:
## MASS::glm.nb(formula = dg1 ~ N.flowers + Quantity + shrub.density + 
##     day + Species, data = all, init.theta = 2.171691063, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4171  -0.7540  -0.0473   0.4566   2.6632  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.046e+00  1.625e-01  18.750  < 2e-16 ***
## N.flowers      7.434e-04  8.197e-04   0.907  0.36444    
## Quantity       2.383e-01  2.116e-02  11.264  < 2e-16 ***
## shrub.density  4.515e-03  2.797e-02   0.161  0.87175    
## day           -1.336e-02  9.628e-03  -1.387  0.16534    
## SpeciesBW      1.893e-01  2.087e-01   0.907  0.36447    
## SpeciesEC      1.683e-01  1.824e-01   0.922  0.35630    
## SpeciesEL     -3.466e+01  3.803e+06   0.000  0.99999    
## SpeciesHH      1.096e+00  3.633e-01   3.017  0.00256 ** 
## SpeciesLT      7.743e-01  1.601e-01   4.836 1.32e-06 ***
## SpeciesPP      1.069e+00  2.087e-01   5.122 3.02e-07 ***
## SpeciesSC      1.193e+00  1.546e-01   7.716 1.20e-14 ***
## SpeciesSD     -1.758e-01  3.126e-01  -0.562  0.57384    
## SpeciesSM     -4.214e-01  4.192e-01  -1.005  0.31484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1717) family taken to be 1)
## 
##     Null deviance: 511.16  on 230  degrees of freedom
## Residual deviance: 260.67  on 217  degrees of freedom
## AIC: 2368.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.172 
##           Std. Err.:  0.211 
## 
##  2 x log-likelihood:  -2338.756
knitr::kable(car::Anova(m2, type = 2))
LR Chisq Df Pr(>Chisq)
N.flowers 0.6942087 1 0.4047365
Quantity 106.5245963 1 0.0000000
shrub.density 0.0244659 1 0.8757054
day 1.8975463 1 0.1683532
Species 111.1358216 9 0.0000000 
check_collinearity(m2)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##      N.flowers 1.00         1.00      1.00
##       Quantity 1.00         1.00      1.00
##  shrub.density 1.00         1.00      1.00
##            day 1.00         1.00      1.00
##        Species 1.00         1.00      1.00

Betweenness

#betweeness
ggplot(all, aes(bt1)) + geom_density()

b2 <-glm(btbin ~ shrub.density + Quantity + day + N.flowers + Species, family = "binomial", data = all)

summary(b2)
## 
## Call:
## glm(formula = btbin ~ shrub.density + Quantity + day + N.flowers + 
##     Species, family = "binomial", data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0337  -0.7101  -0.4394   0.8343   1.9811  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.613e+00  6.285e-01  -4.158 3.21e-05 ***
## shrub.density  2.420e-01  9.893e-02   2.447   0.0144 *  
## Quantity       7.245e-01  1.208e-01   6.000 1.97e-09 ***
## day            7.569e-03  3.448e-02   0.220   0.8262    
## N.flowers      3.205e-03  3.112e-03   1.030   0.3030    
## SpeciesBW      2.762e-01  7.119e-01   0.388   0.6980    
## SpeciesEC     -5.248e-01  6.479e-01  -0.810   0.4179    
## SpeciesEL     -1.555e+01  2.400e+03  -0.006   0.9948    
## SpeciesHH     -1.588e+01  1.162e+03  -0.014   0.9891    
## SpeciesLT     -5.433e-01  5.464e-01  -0.994   0.3201    
## SpeciesPP      1.448e-01  6.981e-01   0.207   0.8356    
## SpeciesSC     -7.934e-01  5.584e-01  -1.421   0.1554    
## SpeciesSD     -2.509e+00  2.177e+00  -1.152   0.2493    
## SpeciesSM     -5.662e-01  1.413e+00  -0.401   0.6886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 314.28  on 230  degrees of freedom
## Residual deviance: 227.08  on 217  degrees of freedom
## AIC: 255.08
## 
## Number of Fisher Scoring iterations: 15
car::Anova(b2, type = 2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: btbin
##               LR Chisq Df Pr(>Chisq)    
## shrub.density    6.170  1    0.01299 *  
## Quantity        57.434  1  3.496e-14 ***
## day              0.048  1    0.82623    
## N.flowers        1.087  1    0.29705    
## Species          9.155  9    0.42309    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#bnull <- glmmTMB(btbin ~ 1 + (1|Species), family = "binomial", data = all)
#anova(b2, bnull)

summary(b2)
## 
## Call:
## glm(formula = btbin ~ shrub.density + Quantity + day + N.flowers + 
##     Species, family = "binomial", data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0337  -0.7101  -0.4394   0.8343   1.9811  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.613e+00  6.285e-01  -4.158 3.21e-05 ***
## shrub.density  2.420e-01  9.893e-02   2.447   0.0144 *  
## Quantity       7.245e-01  1.208e-01   6.000 1.97e-09 ***
## day            7.569e-03  3.448e-02   0.220   0.8262    
## N.flowers      3.205e-03  3.112e-03   1.030   0.3030    
## SpeciesBW      2.762e-01  7.119e-01   0.388   0.6980    
## SpeciesEC     -5.248e-01  6.479e-01  -0.810   0.4179    
## SpeciesEL     -1.555e+01  2.400e+03  -0.006   0.9948    
## SpeciesHH     -1.588e+01  1.162e+03  -0.014   0.9891    
## SpeciesLT     -5.433e-01  5.464e-01  -0.994   0.3201    
## SpeciesPP      1.448e-01  6.981e-01   0.207   0.8356    
## SpeciesSC     -7.934e-01  5.584e-01  -1.421   0.1554    
## SpeciesSD     -2.509e+00  2.177e+00  -1.152   0.2493    
## SpeciesSM     -5.662e-01  1.413e+00  -0.401   0.6886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 314.28  on 230  degrees of freedom
## Residual deviance: 227.08  on 217  degrees of freedom
## AIC: 255.08
## 
## Number of Fisher Scoring iterations: 15
knitr::kable(car::Anova(b2, type = 2))
LR Chisq Df Pr(>Chisq)
shrub.density 6.1700927 1 0.0129928
Quantity 57.4337122 1 0.0000000
day 0.0481998 1 0.8262260
N.flowers 1.0873921 1 0.2970498
Species 9.1549728 9 0.4230937 
check_collinearity(b2)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##  shrub.density 1.32         1.15      0.76
##       Quantity 1.12         1.06      0.89
##            day 1.41         1.19      0.71
##      N.flowers 2.34         1.53      0.43
##        Species 4.46         2.11      0.22
#gaussian part



b4 <- glm(log(bt1) ~ Quantity + shrub.density + day + N.flowers + Species, family = "gaussian", data = bt)
shapiro.test(resid(b4))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(b4)
## W = 0.94969, p-value = 0.0009764
#b4null <- glmmTMB(log(bt1) ~ (1|Species), family = "gaussian", data = bt)
#anova(b4, b4null)
check_collinearity(b4)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##       Quantity 1.16         1.07      0.87
##  shrub.density 2.03         1.42      0.49
##            day 1.79         1.34      0.56
##      N.flowers 3.09         1.76      0.32
## 
## High Correlation
## 
##     Term   VIF Increased SE Tolerance
##  Species 11.10         3.33      0.09
summary(b4)
## 
## Call:
## glm(formula = log(bt1) ~ Quantity + shrub.density + day + N.flowers + 
##     Species, family = "gaussian", data = bt)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.7773  -0.6321   0.0009   1.0717   3.2165  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.5907910  0.6132236   7.486 6.05e-11 ***
## Quantity       0.2281380  0.0665368   3.429 0.000937 ***
## shrub.density -0.1500114  0.0983848  -1.525 0.131036    
## day            0.0275796  0.0344578   0.800 0.425718    
## N.flowers      0.0009064  0.0028772   0.315 0.753505    
## SpeciesBW     -0.3196672  0.6355436  -0.503 0.616277    
## SpeciesEC      0.6342991  0.7403700   0.857 0.394002    
## SpeciesLT     -0.2892603  0.5765147  -0.502 0.617148    
## SpeciesPP      1.4177174  0.7185191   1.973 0.051733 .  
## SpeciesSC     -0.1385590  0.5573058  -0.249 0.804252    
## SpeciesSD      2.0075691  1.6953356   1.184 0.239646    
## SpeciesSM      1.7220659  1.2741964   1.351 0.180124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.310566)
## 
##     Null deviance: 246.63  on 96  degrees of freedom
## Residual deviance: 196.40  on 85  degrees of freedom
## AIC: 369.7
## 
## Number of Fisher Scoring iterations: 2
knitr::kable(car::Anova(b4, type = 2))
LR Chisq Df Pr(>Chisq)
Quantity 11.7563110 1 0.0006064
shrub.density 2.3248375 1 0.1273235
day 0.6406186 1 0.4234869
N.flowers 0.0992488 1 0.7527329
Species 8.9063914 7 0.2594468

Eigancentrality

ggplot(all, aes(sqrt(V2))) + geom_density()

v2 <- glm(sqrt(V2) ~ Quantity+ shrub.density + N.flowers + day +Species, family = "gaussian", data = all)

shapiro.test(resid(v2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(v2)
## W = 0.88547, p-value = 3.222e-12
hist(resid(v2))

#histogram looks fine

check_collinearity(v2)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##       Quantity 1.14         1.07      0.88
##  shrub.density 1.40         1.18      0.71
##      N.flowers 4.04         2.01      0.25
##            day 1.41         1.19      0.71
## 
## Moderate Correlation
## 
##     Term  VIF Increased SE Tolerance
##  Species 7.31         2.70      0.14
summary(v2)
## 
## Call:
## glm(formula = sqrt(V2) ~ Quantity + shrub.density + N.flowers + 
##     day + Species, family = "gaussian", data = all)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.51283  -0.04446   0.00596   0.05249   0.52419  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.412e-02  3.248e-02   0.435  0.66411    
## Quantity       3.185e-02  4.302e-03   7.403 2.90e-12 ***
## shrub.density -3.320e-04  5.626e-03  -0.059  0.95299    
## N.flowers     -2.751e-05  1.651e-04  -0.167  0.86783    
## day           -2.184e-03  1.935e-03  -1.129  0.26025    
## SpeciesBW     -3.934e-02  4.168e-02  -0.944  0.34631    
## SpeciesEC      5.718e-03  3.633e-02   0.157  0.87509    
## SpeciesEL     -3.772e-02  1.425e-01  -0.265  0.79147    
## SpeciesHH      4.325e-01  7.356e-02   5.881 1.53e-08 ***
## SpeciesLT      1.019e-01  3.207e-02   3.176  0.00171 ** 
## SpeciesPP      3.949e-01  4.209e-02   9.381  < 2e-16 ***
## SpeciesSC      5.166e-01  3.103e-02  16.649  < 2e-16 ***
## SpeciesSD      6.902e-02  6.165e-02   1.120  0.26415    
## SpeciesSM     -5.679e-02  8.318e-02  -0.683  0.49554    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.01941603)
## 
##     Null deviance: 16.4538  on 230  degrees of freedom
## Residual deviance:  4.2133  on 217  degrees of freedom
## AIC: -239.42
## 
## Number of Fisher Scoring iterations: 2
knitr::kable(car::Anova(v2, type = 2))
LR Chisq Df Pr(>Chisq)
Quantity 54.8044655 1 0.0000000
shrub.density 0.0034828 1 0.9529400
N.flowers 0.0277617 1 0.8676702
day 1.2740897 1 0.2590009
Species 435.7077823 9 0.0000000 
ggplot(all, aes(Species, V2)) + geom_boxplot() + ylab("Eigancentrality")

#maybe just use the number of sampled plants as an offset?? #conspecifics

#all <- rename(all, cons = cons.x, hets = hets.x, n = n.x)
all <- mutate(all, uncon = n - cons)
y <- cbind(all$cons, all$uncon)

m4 <- glm(y ~ shrub.density+ dg1 + N.flowers + day +Species, family = "binomial"(link = "logit"), data = all)
summary(m4)
## 
## Call:
## glm(formula = y ~ shrub.density + dg1 + N.flowers + day + Species, 
##     family = binomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.0982  -1.4347  -0.1204   1.1166   3.8393  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.572e+00  9.816e-02 -26.206  < 2e-16 ***
## shrub.density  5.741e-02  1.572e-02   3.652  0.00026 ***
## dg1            1.350e-02  4.324e-04  31.227  < 2e-16 ***
## N.flowers      7.147e-04  4.303e-04   1.661  0.09670 .  
## day           -1.334e-03  5.224e-03  -0.255  0.79840    
## SpeciesBW      1.374e+00  1.373e-01  10.001  < 2e-16 ***
## SpeciesEC      2.741e-01  1.155e-01   2.372  0.01767 *  
## SpeciesEL     -1.119e+01  5.354e+02  -0.021  0.98333    
## SpeciesHH      7.691e-01  5.227e-01   1.471  0.14118    
## SpeciesLT      1.074e+00  9.132e-02  11.756  < 2e-16 ***
## SpeciesPP      7.343e-01  1.443e-01   5.088 3.63e-07 ***
## SpeciesSC      1.005e+00  9.573e-02  10.494  < 2e-16 ***
## SpeciesSD     -3.774e-01  5.567e-01  -0.678  0.49785    
## SpeciesSM      2.837e+00  6.000e-01   4.729 2.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3285.35  on 230  degrees of freedom
## Residual deviance:  830.03  on 217  degrees of freedom
## AIC: 1623.6
## 
## Number of Fisher Scoring iterations: 12
#interact_plot(m4, dg1, Quantity)


#check_overdispersion(m4)
#performance use to work on binomial too, residual deviance still suggests over though
#overdispersed


cor.test(all$dg1, all$Quantity)
## 
##  Pearson's product-moment correlation
## 
## data:  all$dg1 and all$Quantity
## t = 13.237, df = 229, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5784854 0.7257965
## sample estimates:
##       cor 
## 0.6584006
cor.test(all$bt1, all$dg1)
## 
##  Pearson's product-moment correlation
## 
## data:  all$bt1 and all$dg1
## t = 7.1098, df = 229, p-value = 1.467e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3133588 0.5254719
## sample estimates:
##       cor 
## 0.4252367
m4.2 <- glm(y ~ shrub.density+ dg1+ N.flowers + day + Species + Quantity, family = "quasibinomial"(link = "logit"), data = all)
summary(m4.2)
## 
## Call:
## glm(formula = y ~ shrub.density + dg1 + N.flowers + day + Species + 
##     Quantity, family = quasibinomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.8592  -1.3015  -0.0016   1.1509   3.8882  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -2.485e+00  1.747e-01 -14.225  < 2e-16 ***
## shrub.density  5.948e-02  2.777e-02   2.142  0.03329 *  
## dg1            1.767e-02  1.190e-03  14.856  < 2e-16 ***
## N.flowers      1.202e-03  7.672e-04   1.567  0.11854    
## day            2.264e-03  9.344e-03   0.242  0.80879    
## SpeciesBW      1.567e+00  2.461e-01   6.366 1.15e-09 ***
## SpeciesEC      2.589e-01  2.036e-01   1.272  0.20479    
## SpeciesEL     -1.119e+01  9.507e+02  -0.012  0.99062    
## SpeciesHH      5.854e-01  9.289e-01   0.630  0.52921    
## SpeciesLT      9.114e-01  1.654e-01   5.509 1.02e-07 ***
## SpeciesPP      5.726e-01  2.629e-01   2.178  0.03050 *  
## SpeciesSC      7.652e-01  1.790e-01   4.274 2.88e-05 ***
## SpeciesSD     -2.871e-01  9.783e-01  -0.294  0.76941    
## SpeciesSM      3.069e+00  1.078e+00   2.848  0.00483 ** 
## Quantity      -1.576e-01  3.253e-02  -4.843 2.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 3.152627)
## 
##     Null deviance: 3285.35  on 230  degrees of freedom
## Residual deviance:  756.49  on 216  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 12
AIC(m4, m4.2)
##      df      AIC
## m4   14 1623.644
## m4.2 15       NA
knitr::kable(car::Anova(m4.2, type = 2))
LR Chisq Df Pr(>Chisq)
shrub.density 4.5836027 1 0.0322793
dg1 261.3418684 1 0.0000000
N.flowers 2.4576295 1 0.1169552
day 0.0586953 1 0.8085701
Species 76.8114949 9 0.0000000
Quantity 23.3249572 1 0.0000014 
#check_collinearity(m4)
#check_overdispersion(m4)

summary(m4)
## 
## Call:
## glm(formula = y ~ shrub.density + dg1 + N.flowers + day + Species, 
##     family = binomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.0982  -1.4347  -0.1204   1.1166   3.8393  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.572e+00  9.816e-02 -26.206  < 2e-16 ***
## shrub.density  5.741e-02  1.572e-02   3.652  0.00026 ***
## dg1            1.350e-02  4.324e-04  31.227  < 2e-16 ***
## N.flowers      7.147e-04  4.303e-04   1.661  0.09670 .  
## day           -1.334e-03  5.224e-03  -0.255  0.79840    
## SpeciesBW      1.374e+00  1.373e-01  10.001  < 2e-16 ***
## SpeciesEC      2.741e-01  1.155e-01   2.372  0.01767 *  
## SpeciesEL     -1.119e+01  5.354e+02  -0.021  0.98333    
## SpeciesHH      7.691e-01  5.227e-01   1.471  0.14118    
## SpeciesLT      1.074e+00  9.132e-02  11.756  < 2e-16 ***
## SpeciesPP      7.343e-01  1.443e-01   5.088 3.63e-07 ***
## SpeciesSC      1.005e+00  9.573e-02  10.494  < 2e-16 ***
## SpeciesSD     -3.774e-01  5.567e-01  -0.678  0.49785    
## SpeciesSM      2.837e+00  6.000e-01   4.729 2.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3285.35  on 230  degrees of freedom
## Residual deviance:  830.03  on 217  degrees of freedom
## AIC: 1623.6
## 
## Number of Fisher Scoring iterations: 12
m4 <- glm(y ~ shrub.density+ dg1 + N.flowers + day + Species, family = "quasibinomial"(link = "logit"), data = all)
summary(m4)
## 
## Call:
## glm(formula = y ~ shrub.density + dg1 + N.flowers + day + Species, 
##     family = quasibinomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.0982  -1.4347  -0.1204   1.1166   3.8393  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -2.572e+00  1.807e-01 -14.238  < 2e-16 ***
## shrub.density  5.741e-02  2.893e-02   1.984   0.0485 *  
## dg1            1.350e-02  7.959e-04  16.966  < 2e-16 ***
## N.flowers      7.147e-04  7.920e-04   0.902   0.3678    
## day           -1.334e-03  9.615e-03  -0.139   0.8898    
## SpeciesBW      1.374e+00  2.528e-01   5.434 1.48e-07 ***
## SpeciesEC      2.741e-01  2.127e-01   1.289   0.1988    
## SpeciesEL     -1.119e+01  9.855e+02  -0.011   0.9910    
## SpeciesHH      7.691e-01  9.621e-01   0.799   0.4249    
## SpeciesLT      1.074e+00  1.681e-01   6.387 1.01e-09 ***
## SpeciesPP      7.343e-01  2.657e-01   2.764   0.0062 ** 
## SpeciesSC      1.005e+00  1.762e-01   5.701 3.85e-08 ***
## SpeciesSD     -3.774e-01  1.025e+00  -0.368   0.7130    
## SpeciesSM      2.837e+00  1.104e+00   2.569   0.0109 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 3.387697)
## 
##     Null deviance: 3285.35  on 230  degrees of freedom
## Residual deviance:  830.03  on 217  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 12
knitr::kable(car::Anova(m4, type = 2))
LR Chisq Df Pr(>Chisq)
shrub.density 3.9350041 1 0.0472910
dg1 346.3087067 1 0.0000000
N.flowers 0.8145845 1 0.3667683
day 0.0192591 1 0.8896263
Species 82.6356289 9 0.0000000 
car::vif(m4)
##                   GVIF Df GVIF^(1/(2*Df))
## shrub.density 1.120002  1        1.058301
## dg1           1.247408  1        1.116874
## N.flowers     2.364213  1        1.537600
## day           1.340436  1        1.157772
## Species       3.710468  9        1.075561
cor.test(all$shrub.density, all$dg1)
## 
##  Pearson's product-moment correlation
## 
## data:  all$shrub.density and all$dg1
## t = -1.3662, df = 229, p-value = 0.1732
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.21648182  0.03962084
## sample estimates:
##         cor 
## -0.08991668
#heterospecific access
all <- mutate(all, unhet = hetn - hets)
z <- cbind(all$hets, all$unhet)

m1 <- glm(z ~ day + shrub.density+ dg1+ N.flowers+ Species, family = "binomial"(link = "logit"), data = all)
summary(m1)
## 
## Call:
## glm(formula = z ~ day + shrub.density + dg1 + N.flowers + Species, 
##     family = binomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.4554  -1.6463  -0.5594   0.9963   7.3405  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.288e+00  7.605e-02 -43.234  < 2e-16 ***
## day           -1.020e-02  4.394e-03  -2.322   0.0202 *  
## shrub.density  6.003e-03  1.180e-02   0.509   0.6108    
## dg1            8.974e-03  3.241e-04  27.688  < 2e-16 ***
## N.flowers      6.493e-04  3.970e-04   1.636   0.1019    
## SpeciesBW      1.635e-01  8.992e-02   1.818   0.0691 .  
## SpeciesEC     -6.125e-02  8.682e-02  -0.705   0.4805    
## SpeciesEL     -1.592e+01  5.721e+02  -0.028   0.9778    
## SpeciesHH      8.057e-01  1.114e-01   7.235 4.66e-13 ***
## SpeciesLT     -1.055e+00  8.880e-02 -11.879  < 2e-16 ***
## SpeciesPP      6.096e-01  8.323e-02   7.325 2.39e-13 ***
## SpeciesSC     -5.096e-01  8.506e-02  -5.991 2.08e-09 ***
## SpeciesSD      1.819e-01  1.285e-01   1.416   0.1569    
## SpeciesSM     -2.064e-01  2.040e-01  -1.012   0.3115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2475.71  on 230  degrees of freedom
## Residual deviance:  907.19  on 217  degrees of freedom
## AIC: 1803.5
## 
## Number of Fisher Scoring iterations: 13
#check_overdispersion(m1)
#also overdisperse


m5.2 <- glm(z ~ dg1 + shrub.density + N.flowers + day +Species, family = "quasibinomial"(link = "logit"), data = all)
summary(m5.2)
## 
## Call:
## glm(formula = z ~ dg1 + shrub.density + N.flowers + day + Species, 
##     family = quasibinomial(link = "logit"), data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.4554  -1.6463  -0.5594   0.9963   7.3405  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -3.288e+00  1.566e-01 -20.997  < 2e-16 ***
## dg1            8.974e-03  6.674e-04  13.447  < 2e-16 ***
## shrub.density  6.003e-03  2.429e-02   0.247 0.805021    
## N.flowers      6.493e-04  8.174e-04   0.794 0.427857    
## day           -1.020e-02  9.047e-03  -1.128 0.260697    
## SpeciesBW      1.635e-01  1.851e-01   0.883 0.378268    
## SpeciesEC     -6.125e-02  1.788e-01  -0.343 0.732213    
## SpeciesEL     -1.592e+01  1.178e+03  -0.014 0.989233    
## SpeciesHH      8.057e-01  2.293e-01   3.514 0.000538 ***
## SpeciesLT     -1.055e+00  1.829e-01  -5.769 2.73e-08 ***
## SpeciesPP      6.096e-01  1.714e-01   3.558 0.000460 ***
## SpeciesSC     -5.096e-01  1.751e-01  -2.910 0.003993 ** 
## SpeciesSD      1.819e-01  2.645e-01   0.687 0.492526    
## SpeciesSM     -2.064e-01  4.200e-01  -0.491 0.623583    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 4.23966)
## 
##     Null deviance: 2475.71  on 230  degrees of freedom
## Residual deviance:  907.19  on 217  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 13
knitr::kable(car::Anova(m5.2, type = 2))
LR Chisq Df Pr(>Chisq)
dg1 183.3044968 1 0.0000000
shrub.density 0.0609745 1 0.8049623
N.flowers 0.6250378 1 0.4291814
day 1.2714973 1 0.2594860
Species 140.1316632 9 0.0000000 
cor.test(all$dg1, all$Quantity)
## 
##  Pearson's product-moment correlation
## 
## data:  all$dg1 and all$Quantity
## t = 13.237, df = 229, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5784854 0.7257965
## sample estimates:
##       cor 
## 0.6584006

Effective partners

ggplot(all, aes(effective.partners)) + geom_density()

m1 <- glm(effective.partners ~ Quantity + shrub.density + N.flowers + day +Species, family = "gaussian", data = all)
hist(resid(m1))

m2 <- glm(effective.partners ~ Quantity + shrub.density + N.flowers + day + Species, family = Gamma(link = "log"), data = all)


AIC(m1,m2)
##    df      AIC
## m1 15 449.4442
## m2 15 361.9216
shapiro.test(all$effective.partners)
## 
##  Shapiro-Wilk normality test
## 
## data:  all$effective.partners
## W = 0.75919, p-value < 2.2e-16
summary(m2)
## 
## Call:
## glm(formula = effective.partners ~ Quantity + shrub.density + 
##     N.flowers + day + Species, family = Gamma(link = "log"), 
##     data = all)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9621  -0.2571  -0.1184   0.2371   0.8432  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.0328035  0.0831648   0.394  0.69364    
## Quantity       0.1082462  0.0110170   9.825  < 2e-16 ***
## shrub.density  0.0380488  0.0144060   2.641  0.00886 ** 
## N.flowers      0.0010411  0.0004228   2.462  0.01458 *  
## day           -0.0013285  0.0049537  -0.268  0.78882    
## SpeciesBW      0.1939648  0.1067310   1.817  0.07055 .  
## SpeciesEC     -0.1080376  0.0930282  -1.161  0.24678    
## SpeciesEL     -0.3362941  0.3648966  -0.922  0.35775    
## SpeciesHH     -0.2853827  0.1883492  -1.515  0.13118    
## SpeciesLT     -0.1157679  0.0821091  -1.410  0.15999    
## SpeciesPP      0.0661002  0.1077858   0.613  0.54035    
## SpeciesSC     -0.0661505  0.0794453  -0.833  0.40595    
## SpeciesSD      0.0277277  0.1578620   0.176  0.86074    
## SpeciesSM     -0.0571744  0.2129991  -0.268  0.78863    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1273064)
## 
##     Null deviance: 46.398  on 230  degrees of freedom
## Residual deviance: 26.784  on 217  degrees of freedom
## AIC: 361.92
## 
## Number of Fisher Scoring iterations: 7
knitr::kable(car::Anova(m2, type = 2))
LR Chisq Df Pr(>Chisq)
Quantity 82.7193092 1 0.0000000
shrub.density 6.7162601 1 0.0095538
N.flowers 5.6694639 1 0.0172627
day 0.0736917 1 0.7860357
Species 13.3983665 9 0.1453934 
check_collinearity(m2)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Term  VIF Increased SE Tolerance
##       Quantity 1.14         1.07      0.88
##  shrub.density 1.40         1.18      0.71
##      N.flowers 4.04         2.01      0.25
##            day 1.41         1.19      0.71
## 
## Moderate Correlation
## 
##     Term  VIF Increased SE Tolerance
##  Species 7.31         2.70      0.14
ggplot(all, aes(Species, effective.partners)) + geom_boxplot()

Nestedness

ggplot(all, aes(nestedcontribution)) + geom_density()

shapiro.test(all$nestedcontribution)
## 
##  Shapiro-Wilk normality test
## 
## data:  all$nestedcontribution
## W = 0.9416, p-value = 5.516e-08
n1 <- glm(nestedcontribution ~ N.flowers+day  + Quantity  + shrub.density + 
            Species, family = "gaussian", data = all)
shapiro.test(resid(n1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(n1)
## W = 0.9807, p-value = 0.003047
car::vif(n1)
##                   GVIF Df GVIF^(1/(2*Df))
## N.flowers     2.332391  1        1.527217
## day           1.401935  1        1.184033
## Quantity      1.116855  1        1.056814
## shrub.density 1.379016  1        1.174315
## Species       4.348473  9        1.085083
summary(n1)
## 
## Call:
## glm(formula = nestedcontribution ~ N.flowers + day + Quantity + 
##     shrub.density + Species, family = "gaussian", data = all)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.54230  -0.34911  -0.04695   0.36838   2.56875  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.4224514  0.1514005  -2.790  0.00573 ** 
## N.flowers      0.0006303  0.0007697   0.819  0.41369    
## day           -0.0132097  0.0090182  -1.465  0.14443    
## Quantity       0.1294838  0.0200563   6.456 6.93e-10 ***
## shrub.density -0.0026897  0.0262259  -0.103  0.91841    
## SpeciesBW      0.0204989  0.1943025   0.106  0.91608    
## SpeciesEC      0.1241651  0.1693567   0.733  0.46425    
## SpeciesEL     -0.4867736  0.6642896  -0.733  0.46449    
## SpeciesHH      0.8599533  0.3428873   2.508  0.01287 *  
## SpeciesLT      0.6377295  0.1494786   4.266 2.97e-05 ***
## SpeciesPP      1.0211541  0.1962227   5.204 4.51e-07 ***
## SpeciesSC      1.2668395  0.1446291   8.759 5.67e-16 ***
## SpeciesSD     -0.0216937  0.2873858  -0.075  0.93990    
## SpeciesSM     -0.4178677  0.3877622  -1.078  0.28239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.4219155)
## 
##     Null deviance: 176.291  on 230  degrees of freedom
## Residual deviance:  91.556  on 217  degrees of freedom
## AIC: 471.77
## 
## Number of Fisher Scoring iterations: 2
knitr::kable(car::Anova(n1, type = 2))
LR Chisq Df Pr(>Chisq)
N.flowers 0.6707428 1 0.4127928
day 2.1456086 1 0.1429783
Quantity 41.6799550 1 0.0000000
shrub.density 0.0105185 1 0.9183126
Species 127.2782871 9 0.0000000

#Modules and modularity modelling

library(nnet)

m1 <- multinom(module ~ Species + shrub.density  + Quantity+ N.flowers + day, data = all)
## # weights:  56 (39 variable)
## initial  value 313.302526 
## iter  10 value 133.730573
## iter  20 value 67.058444
## iter  30 value 64.100438
## iter  40 value 63.732325
## iter  50 value 63.629500
## iter  60 value 63.613698
## iter  70 value 63.609926
## iter  80 value 63.609800
## iter  80 value 63.609800
## final  value 63.609800 
## converged
summary(m1)
## Warning in sqrt(diag(vc)): NaNs produced
## Call:
## multinom(formula = module ~ Species + shrub.density + Quantity + 
##     N.flowers + day, data = all)
## 
## Coefficients:
##   (Intercept)  SpeciesBW SpeciesEC  SpeciesHH SpeciesLT SpeciesPP  SpeciesSC
## 1    134.7520 -194.87193 -167.6419  -55.22961 -162.4881  137.9230 -129.65998
## 2    139.3741  191.75534 -167.2072 -116.86158 -169.8534  135.6282   71.46729
## 3    135.1930  -97.68008 -369.7410  168.39503 -166.1612  140.6470   77.12401
##   SpeciesSD  SpeciesSM shrub.density Quantity     N.flowers       day
## 1  126.3942   2.574872   -0.01461747 28.72668 -7.952179e-06 0.2782226
## 2  121.4919  14.462814   -0.03611962 28.34664  3.521932e-03 0.1996707
## 3  125.1843 -26.381954   -0.04694358 28.69329  1.691513e-03 0.1885051
## 
## Std. Errors:
##   (Intercept)     SpeciesBW    SpeciesEC    SpeciesHH SpeciesLT SpeciesPP
## 1    1.384115  1.556517e-11 1.973154e+00 1.189784e-12  1.108174 1.0485915
## 2    1.330712           NaN 1.582369e+00 1.778241e-15  1.310817 0.8199640
## 3    1.556226 3.212889e-125 2.104059e-87 3.175894e-96  1.637366 0.9625024
##   SpeciesSC SpeciesSD    SpeciesSM shrub.density Quantity  N.flowers       day
## 1       NaN 0.9376496 3.555986e-06     0.7195968 1.034941 0.01412017 0.2312982
## 2 0.6544695 0.9522779 3.555986e-06     0.7212355 1.045156 0.01323238 0.2266176
## 3 0.6544695 1.1071535 1.099133e-18     0.7382535 1.042893 0.01775685 0.2363944
## 
## Residual Deviance: 127.2196 
## AIC: 205.2196
options(scipen=999)
knitr::kable(car::Anova(m1, type = 2))
LR Chisq Df Pr(>Chisq)
Species 231.7233592 24 0.0000000
shrub.density 0.0151803 3 0.9995048
Quantity 8.7665960 3 0.0325606
N.flowers 0.1942784 3 0.9785076
day 2.5649047 3 0.4636755

Network comparisons

The random network parts take a really long time to run, therefore are output as csv and read back in. Code is commented out.

These are all bipartite networks

#calculate network level indices for the bipartite networks. 

#make a species network
sp <- net(visits, Species, ID)
## [1] 635
## `summarise()` has grouped output by 'Species'. You can override using the `.groups` argument.
## [1] 635
write.csv(sp, "Data/species_network.csv")
n <- networklevel(sp)
n <- as.data.frame(n)

#output networks repo
write.csv(ind1, "Data/individal_network.csv")


#use existing individual network
#i <- as.data.frame(networklevel(ind1))
#write.csv(i, "data/Output/ind_network_indices.csv")
i <- read.csv("data/Output/ind_network_indices.csv")

sp <- net(visits, Species, ID)
## [1] 635
## `summarise()` has grouped output by 'Species'. You can override using the `.groups` argument.
## [1] 635
sp.m <- computeModules(sp)

i.m <- computeModules(ind1)

mod <- i.m@moduleWeb

#run the random modules function on the species network
#sp.ran.mod <- ranModules(sp)
#write.csv(sp.ran.mod, "data/Randomization_outputs/randommodule_sp_bi.csv")

sp.ran.mod <- read.csv("data/Randomization_outputs/randommodule_sp_bi.csv")
mean(sp.ran.mod$`X..`)
## [1] 0.1494235
net.zscore(sp.m@likelihood, sp.ran.mod$`X..`)
## [1] 46.87052
#run the random modules function on the individual network, takes ~72 hours
#ind.ran.mod <- ranModules(ind1)
#write.csv(ind.ran.mod, "Data/Randomization_outputs/randommoduleoutput_bipartite.csv")

ind.ran.mod <- read.csv("Data/Randomization_outputs/randommoduleoutput_bipartite.csv")
ind.ran.mod$t <- ind.ran.mod$X..i..
net.zscore(i.m@likelihood, ind.ran.mod$t)
## [1] 32.16754
mean(ind.ran.mod$t)
## [1] 0.4286197
max(ind.ran.mod$t)
## [1] 0.4564127
#individual networks

ind.WC <- ranindices(ind1, "weighted connectance")
ind.WNODF <-ranindices(ind1, "weighted NODF")
ind.H2 <- ranindices(ind1, "H2")

net.zscore(i[10,2], ind.WNODF$`weighted NODF`)
## [1] -4.691895
mean(ind.WNODF$`weighted NODF`)
## [1] 3.699326
max(ind.WNODF$`weighted NODF`)
## [1] 4.661721
min(ind.WNODF$`weighted NODF`)
## [1] 2.699631
net.zscore(i[14,2], ind.WC$`weighted connectance`)
## [1] -21.28123
mean(ind.WC$`weighted connectance`)
## [1] 0.07369418
max(ind.WC$`weighted connectance`)
## [1] 0.07918955
min(ind.WC$`weighted connectance`)
## [1] 0.06792919
net.zscore(i[19,2], ind.H2$H2)
## [1] 35.59199
mean(ind.H2$H2)
## [1] 0.1320496
max(ind.H2$H2)
## [1] 0.167811
#species networks

sp.WNODF <- ranindices(sp, "weighted NODF")
net.zscore(n[10,], sp.WNODF$`weighted NODF`)
## [1] -11.14848
mean(sp.WNODF$`weighted NODF`)
## [1] 35.30176
max(sp.WNODF$`weighted NODF`)
## [1] 41.54066
sp.WC <- ranindices(sp, "weighted connectance")
n[14,]
## [1] 0.06720467
net.zscore(n[14,], sp.WC$`weighted connectance`)
## [1] -34.16967
mean(sp.WC$`weighted connectance`)
## [1] 0.1545774
max(sp.WC$`weighted connectance`)
## [1] 0.1632196
sp.H2 <- ranindices(sp, "H2")
n[19,]
## [1] 0.6150591
net.zscore(n[19,], sp.H2$H2)
## [1] 63.0344
mean(sp.H2$H2)
## [1] 0.08900217
#violindles


labels = c("Acamptopappus sphaerocephalus", "Eriogonum fasciculatum", "Ericameria cooperi", "Echinocereus engelmannii", "Larrea tridentata", "Opuntia basilaris", "Cylindropuntia echinocarpa", "Salvia dorri", "Scutellaria mexicana")

ggplot(filter(all, module != "NA"), aes(module, Quantity)) + geom_violin() + stat_summary(fun.y=mean, geom="point", shape=23, size=2) + geom_jitter(aes(color = Species, shape = Species), position=position_jitter(0.2)) + ylab("Visitation Rate") + xlab("Module") + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + scale_shape_manual(values=c(3, 16, 17, 3, 16, 17, 3, 16, 17), labels = labels) +scale_color_discrete(name = "Species", labels = labels)
## Warning: `fun.y` is deprecated. Use `fun` instead.

labels = c("Acamptopappus sphaerocephalus", "Eriogonum fasciculatum", "Ericameria cooperi", "Ericameria linearifolia", "Echinocereus engelmannii", "Larrea tridentata", "Opuntia basilaris", "Cylindropuntia echinocarpa", "Salvia dorri", "Scutellaria mexicana")

ggplot(all, aes(module, N.flowers)) + geom_boxplot()

ggplot(all, aes(module, shrub.density)) + geom_boxplot()

#centrality species comparisons
#make a smaller dataframe of means and se for figure 1
library(ggthemes)
se <- function(x) {sd(x)/sqrt(length(x))} ## SE
central.means <- all %>% group_by(Species) %>% summarise(dg = mean(dg1), visits = mean(Quantity), eig = mean(V2), bt = mean(bt1)) %>% gather(var, mean, 2:5)
central.se <- all %>% group_by(Species) %>% summarise(dg = se(dg1), visits = se(Quantity), eig = se(V2), bt = se(bt1)) %>% gather(var, se, 2:5)
cent <- cbind(central.means, central.se)
cent <- cent[,c(1,2,3,6)]
cent$var <- factor(cent$var, levels = c("visits", "dg", "eig", "bt", labels = c("Visitation Rates", "Degree centrality", "Eigenvector Centrality", "Betweenness")))

ggplot(cent, aes(Species, mean)) + geom_point(stat = "identity", aes(fill = Species)) +  geom_errorbar(aes(Species, ymin = mean - se, ymax = mean + se)) + theme_classic() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + facet_wrap(.~var, scales = "free")  

visit_plot <- filter(cent, var == "visits")
dg_plot <- filter(cent, var == "dg")
eig_plot <- filter(cent, var == "eig")
bt_plot <- filter(cent, var == "bt")

a <- ggplot(visit_plot, aes(Species, mean)) + geom_point() +  geom_errorbar(aes(Species, ymin = mean - se, ymax = mean + se)) + theme_classic() + ylab("Mean Visitation Rate") + xlab("") + theme(axis.text.x = element_blank(),panel.border = element_rect(colour = "black", fill=NA, size=1), aspect.ratio = 1) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))


b <- ggplot(dg_plot, aes(Species, mean)) + geom_point() +  geom_errorbar(aes(Species, ymin = mean - se, ymax = mean + se)) + theme_classic() + ylab("Mean Degree Centrality") + xlab("") + theme(axis.text.x = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), aspect.ratio = 1) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))

c <- ggplot(eig_plot, aes(Species, mean, fill = Species)) + geom_point() +  geom_errorbar(aes(Species, ymin = mean - se, ymax = mean + se)) + theme_classic() + ylab("Mean Eigenvector Centrality") + scale_x_discrete(labels = labels) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x = element_text(angle = 90, hjust = 0.95, vjust = .1), legend.position = "none", aspect.ratio = 1) + xlab("") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))

d <- ggplot(bt_plot, aes(Species, mean, fill = Species)) + geom_point() +  geom_errorbar(aes(Species, ymin = mean - se, ymax = mean + se)) + theme_classic() + ylab("Mean Betweenness Centrality") + scale_x_discrete(labels = labels) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x = element_text(angle = 90, hjust = 0.95, vjust = .1), legend.position = "none") + xlab("") + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))


library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
fig1 <- grid.arrange(a, b, c, d, ncol = 2, heights = c(1, 1.4), widths = c(1, 1.4))

#Figure 2


e <- ggplot(pl, aes(dg1, prop, color = type)) + geom_point()  +xlab("Degree Centrality")+ geom_smooth(method = lm) + ylab("Proportion of Population") + scale_color_discrete(name = "", labels = c("Conspecific", "Heterospecific")) + theme_classic() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + theme(legend.position = "none", aspect.ratio = 1)

f <- ggplot(pl, aes(shrub.density, prop, color = type)) + geom_point() + geom_smooth(method = lm) +xlab("Shrub Density") + ylab("Proportion of Population") + scale_color_discrete(name = "", labels = c("Conspecific", "Heterospecific")) + theme_classic() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), legend.position = "bottom", aspect.ratio = 1) 
e
## `geom_smooth()` using formula 'y ~ x'

grid.arrange(e, f, ncol =1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Figure 1

labels <- select(all, uniID, Species)
labels <- labels[match(row.names(ind1), labels$uniID),]
all.equal(labels$uniID, row.names(ind1))
all.equal(labels$uniID, V(one1)$name)
#ok, these are the same. let's overwrite the network labels with these simplified species labels

V(one1)$name <- labels$Species
labels$Species
V(one1)$Species <- labels$Species
pal1 <- heat.colors(10, alpha=1)

library(colorspace)
library(RColorBrewer)
V(one1)$Species <- as.factor(V(one1)$Species)
levels(V(one1)$Species) <- 1:10
comps <- as.numeric(V(one1)$Species)
comps
colbar <- brewer.pal(n = 11, name = "Paired")
#(max(comps)+1)
V(one1)$color <- colbar[comps]

#V(one1)$color <- pal1[V(one1)$Species]

#making a figure for the unipartite
plot(one1, rescale = TRUE, edge.arrow.size=0, edge.curved=0, vertex.shape= "circle",vertex.size = 2, vertex.frame.color="#555555", edge.color = "black", layout = layout.fruchterman.reingold, margin = 0, vertex.label = NA)

#legend(x=-1.5, y=-1.1, pch=21,
#col="#777777", pt.cex=2, cex=.8, bty="n", ncol=1)
unique(sp)
plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("center",fill = colbar, legend=c("Acamptopappus sphaerocephalus", "Eriogonum fasciculatum", "Ericameria cooperi", "Ericameria linearfolia", "Echinocereus engelmannii", "Larrea tridentata", "Opuntia basilaris", "Cylindropuntia echinocarpa", "Salvia dorri", "Scutellaria mexicana"))




plot(l)
sp <- V(one1)$Species
V(one1)$Species
#vertex.color="red", vertex.size = 3,
all.equal(labels$uniID, V(one1)$name)
vertex.label = NA
V(one1)$name



#let's try the bipartite network now
plot(sp)





plotweb(sp)
visweb(ind1)


plotweb(sp, method = "cca", text.rot = 90, labsize =1.5, ybig = 0.7, low.y = 0.7, high.y = 0.98, plot.axes = FALSE, y.width.low = 0.05, y.width.high = 0.05, col.low = colbar, high.lablength = 30, low.lablength = 0)


plotweb(ind1, method = "cca", text.rot = 90, labsize =1.5, ybig = 0.7, low.y = 0.7, high.y = 0.98, plot.axes = FALSE, y.width.low = 0.05, y.width.high = 0.05, col.low = colbar[comps], bor.col.low = colbar[comps], high.lablength = 00, low.lablength = 0)

colbar
colbar[comps]

levels(V(one1)$Species) <- 1:10
comps <- as.numeric(V(one1)$Species)
comps
colbar <- brewer.pal(n = 11, name = "Paired")
#(max(comps)+1)
V(one1)$color <- colbar[comps]