Arthropod associations with residual dry matter and foundations shrubs in California desert ecosystems

Pitfall traps were deployed along an aridity gradient in California to measure associations of ground-dwelling arthropods with foundation shrubs, and to assess the ability of residual dry matter (RDM) as an indicator of arthropod community structure.

URL of repository with cleaning scripts and data: https://github.com/jennabraun/rdm.gradient

env <- read.csv("clean_data/cov.csv")
comm <- read.csv("clean_data/comm.csv")
rdm <- read.csv("clean_data/veg_covariates.csv")
rdm <- filter(rdm, Microsite != "larrea")
eph <- filter(env)
sum(eph$abun)
## [1] 5820
str(eph)
## 'data.frame':    135 obs. of  31 variables:
##  $ X             : chr  "C1 ephedra 10" "C1 ephedra 14" "C1 ephedra 18" "C1 ephedra 2" ...
##  $ Region        : chr  "Cuyama" "Cuyama" "Cuyama" "Cuyama" ...
##  $ site          : chr  "C1" "C1" "C1" "C1" ...
##  $ ID            : int  10 14 18 2 26 30 6 10 18 2 ...
##  $ Microsite     : chr  "ephedra" "ephedra" "ephedra" "ephedra" ...
##  $ RDM           : num  5.076 2.938 0.294 2.561 1.674 ...
##  $ x             : int  295 425 375 410 380 425 225 NA NA NA ...
##  $ y             : int  245 400 365 405 280 375 195 NA NA NA ...
##  $ z             : int  170 195 120 210 115 190 160 NA NA NA ...
##  $ bur.drip      : int  2 3 1 1 2 1 1 1 0 1 ...
##  $ bur.buf       : int  0 0 0 0 0 0 0 NA NA NA ...
##  $ max.veg.height: int  49 32 20 31 16 NA 32 20 7 19 ...
##  $ rdm.cov       : int  85 50 7 80 20 0 45 40 5 20 ...
##  $ green.cov     : int  0 0 0 0 0 0 0 5 0 0 ...
##  $ non.veg.cov   : int  15 20 91 5 75 0 55 55 95 80 ...
##  $ largewood.cov : int  NA 30 NA 15 5 0 NA NA NA NA ...
##  $ lrgrocks.cov  : logi  NA NA NA NA NA NA ...
##  $ AT            : chr  "Y" NA NA NA ...
##  $ PF            : chr  "Y" "Y" "Y" "Y" ...
##  $ PF.out        : chr  "june.28" "june.28" "june.28" "june.28" ...
##  $ PF.in         : chr  "july.1" "july.1" "july.1" "july.1" ...
##  $ RDM.date      : chr  "july.1" "july.1" "july.1" "july.1" ...
##  $ comments      : chr  NA NA NA NA ...
##  $ uniID         : chr  "C1 ephedra 10" "C1 ephedra 14" "C1 ephedra 18" "C1 ephedra 2" ...
##  $ abun          : int  15 32 53 72 25 45 5 8 22 32 ...
##  $ H             : num  1.08 1.08 0.83 1.36 1.79 ...
##  $ Simpson       : num  0.578 0.43 0.335 0.646 0.749 ...
##  $ Species       : int  4 9 8 8 10 14 4 6 8 5 ...
##  $ Even          : num  0.777 0.493 0.399 0.653 0.776 ...
##  $ abun.nos      : int  15 31 51 71 25 44 5 8 22 31 ...
##  $ Species.nos   : int  4 8 6 7 10 13 4 6 8 4 ...
eph$Microsite <- as.factor(eph$Microsite)
rdm$Microsite <- as.factor(rdm$Microsite)
eph$Microsite <- relevel(eph$Microsite, ref = "open")
rdm$Microsite <- relevel(rdm$Microsite, ref = "open")

se <- function(x) sd(x, na.rm = TRUE)/sqrt(length(x)) ## SE

rii <- read.csv("clean_data/rii.csv")
indrii <- read.csv("clean_data/rii_individual.csv")

effect <- dplyr::select(rii, site, ESI,Max,arid)

eph <- left_join(eph, effect, by = "site")
eph <- distinct(eph)

rdm <- left_join(rdm, effect, by = "site")
rdm <- distinct(rdm)

#env <- left_join(env, effect, by = "site")
#env <- distinct(env)

eph.abun <- filter(eph, abun < 350)
eph.abun <- distinct(eph.abun)

eph %>% group_by(Region) %>% summarise(sum = sum(abun))
## # A tibble: 3 x 2
##   Region    sum
##   <chr>   <int>
## 1 Cuyama   1616
## 2 Mojave   1622
## 3 Panoche  2582

EDA

The sites are sorted by aridity - from most arid to least arid.

Residual dry matter

ggplot(rdm, aes(Microsite, RDM)) + geom_boxplot() + facet_grid(~reorder(site, arid))+ theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

Does shrub size influence RDM?

ggplot(rdm, aes(x, RDM)) + geom_smooth() + facet_grid(~reorder(site, arid)) + xlab("Shrub Width")+ theme_classic()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 272 rows containing non-finite values (stat_smooth).

ggplot(rdm, aes(z, RDM)) + geom_smooth() + facet_grid(~reorder(site, arid)) + xlab("Shrub Height")+ theme_classic()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 272 rows containing non-finite values (stat_smooth).

Insect morphospecies richness

ggplot(eph, aes(Species, fill = Microsite)) + geom_density(alpha = 0.7) + facet_grid(~reorder(site, arid))+ theme_classic()

ggplot(eph, aes(RDM, Species)) + geom_smooth(method = lm) + geom_point() + theme_classic()
## `geom_smooth()` using formula 'y ~ x'

ggplot(eph, aes(Microsite, Species)) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggplot(eph, aes(RDM, Species, fill = Microsite)) + geom_smooth(method = lm) + geom_point() + theme_classic()
## `geom_smooth()` using formula 'y ~ x'

ggplot(eph, aes(RDM, Species, fill = Microsite)) + geom_smooth(method = lm) + geom_point() + theme_classic() + facet_grid(~Region)
## `geom_smooth()` using formula 'y ~ x'

Insect abundance

There is one sample with 1200 insects in it, almost all ants that really throws everything off.

eph %>% filter(abun < 350) %>% ggplot(aes(abun, fill = Microsite)) + geom_density(alpha = 0.7) + facet_grid(~reorder(site, arid))+ theme_classic()

ggplot(filter(eph, abun <350), aes(Microsite, abun)) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggplot(filter(eph, abun <350), aes(Microsite, abun)) + geom_boxplot() + facet_grid(~Region) + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggplot(filter(eph, abun < 350), aes(RDM, abun)) + geom_smooth(method = lm) + geom_point() + theme_classic()
## `geom_smooth()` using formula 'y ~ x'

ggplot(filter(eph, abun < 350), aes(RDM, abun, fill = Microsite)) + geom_smooth(method = lm) + geom_point() + theme_classic()
## `geom_smooth()` using formula 'y ~ x'

ggplot(filter(eph, abun < 350), aes(RDM, abun, fill = Microsite)) + geom_smooth(method = lm) + geom_point() + theme_classic() + facet_grid(~Region)
## `geom_smooth()` using formula 'y ~ x'

Site level RII

Residual dry matter

ggplot(rii, aes(arid, average, color = response)) + geom_errorbar(aes(arid, ymin = average - error, ymax = average + error)) + geom_point() + geom_jitter() 

ggplot(rii, aes(arid, average, color = response)) + geom_smooth(method = lm)
## `geom_smooth()` using formula 'y ~ x'

Individual RII (rdm)

ggplot(indrii, aes(site, rii)) + geom_boxplot()

ggplot(indrii, aes(x, rii)) + geom_point() + geom_smooth(method = lm) + xlab("Shrub Width")
## `geom_smooth()` using formula 'y ~ x'

Counts and summaries

Simple counts

arth <- read.csv("clean_data/arth_long.csv")
arth.eph <- filter(arth, Microsite != "larrea")

sum(arth.eph$Quantity)
## [1] 5820
ants <- filter(arth.eph, Family == "Formicidae") 
sum(ants$Quantity)
## [1] 3632
ants %>% count(highest.rtu)
##                  highest.rtu  n
## 1    Aphaenogaster megommata 20
## 2           Camponotus black  2
## 3   Camponotus semitestaceus  1
## 4      Crematogaster depilis  7
## 5         Dorymyrmex insanus 30
## 6           Forelius mccooki  3
## 7         Forelius pruinosus 28
## 8              Messor andrei  7
## 9           Messor pergandei 33
## 10    Myrmecocystus kennedyi 26
## 11    Myrmecocystus lugubris  1
## 12   Myrmecocystus mexicanus  1
## 13   Myrmecocystus testaceus  6
## 14          Pheidole barbata  1
## 15       Pheidole desertorum  4
## 16           Pheidole hyatti 60
## 17 Pogonomyrmex californicus 26
## 18          Solenopsis aurea  1
## 19         Solenopsis xyloni 39
## 20          Tapinoma sessile  1
## 21        Temnothorax andrei  7
mean(eph.abun$abun)
## [1] 31.7594
sum(eph.abun$abun)
## [1] 4224
sd(eph.abun$abun)
## [1] 39.15764
order <- arth.eph %>% group_by(Order) %>% summarise(sum = sum(Quantity))
sum(order$sum)
## [1] 5820
beetles <- filter(arth.eph, Order == "Coleoptera")
beetles %>% count(highest.rtu)
##                   highest.rtu  n
## 1  Aethecerinus latecinctus\n  1
## 2                    Alaephus  7
## 3                   Aleochara  3
## 4        Anepsius delicatulus 25
## 5              Apsena rufipes  1
## 6                     Attalus  1
## 7                  Auchmobius 11
## 8                  Blapstinus  4
## 9         Calathus ruficollis  1
## 10     Carpophilus hemipterus  3
## 11         Cerenopus concoler  4
## 12                   Conibius  1
## 13         Cononotus sericans 19
## 14      Cryptoglossa muricata 12
## 15           Edrotes rotundus  1
## 16                  Eleodes 2  3
## 17             Eleodes armata 23
## 18           Eleodes dentipes 10
## 19        Eleodes longicollis  8
## 20                 Histeridae  3
## 21       Horistonotus simplex 10
## 22                  Hymenorus  1
## 23              Hyperaspidius  1
## 24                 Metoponium  2
## 25         Niptus ventriculus 13
## 26        Nocibiotes caudatus  1
## 27      Notibius puncticollis  1
## 28      Octinodes amplicollis  1
## 29               Phyllotreta   1
## 30               Pterostichus  1
## 31                     Ptinus  1
## 32    Tetragonoderus pallidus  2
## 33        Thinoxenus squalens  7
## 34                 Triorophus 13
## 35                 Trogoderma  1
## 36          Typhaea stercorea  7
fam <- arth.eph %>% group_by(Family) %>% summarise(sum = sum(Quantity))

#how many plots contain green veg

#eph %>% count(green.cov == 0)

Models

Drivers of arthropod abundance

eph.abun$x[is.na(eph.abun$x == TRUE)] <- 0
eph$x[is.na(eph$x == TRUE)] <- 0
m1 <- glmmTMB(abun ~ Microsite + RDM + rdm.cov +  arid + ESI + (1|site), family = "poisson", data = eph.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#summary(m1)
#looks overdispersed
check_overdispersion(m1)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## # Overdispersion test
## 
##        dispersion ratio =   22.526
##   Pearson's Chi-Squared = 2838.309
##                 p-value =  < 0.001
## Overdispersion detected.
m1.nb <- glmmTMB(abun ~ Microsite + RDM + rdm.cov + arid + ESI +(1|site), family = "nbinom2", data = eph.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m1.nb)
##  Family: nbinom2  ( log )
## Formula:          abun ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph.abun
## 
##      AIC      BIC   logLik deviance df.resid 
##   1159.0   1182.2   -571.5   1143.0      125 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1865   0.4319  
## Number of obs: 133, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.75 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.004877   0.367399   8.179 2.87e-16 ***
## Micrositeephedra  0.586142   0.173345   3.381 0.000721 ***
## RDM               0.035951   0.034697   1.036 0.300136    
## rdm.cov          -0.008758   0.004015  -2.181 0.029162 *  
## arid              0.031663   0.023284   1.360 0.173881    
## ESI              -2.180906   2.266085  -0.962 0.335843    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.nb1 <- glmmTMB(abun ~ Microsite + RDM + rdm.cov +  arid + ESI + (1|site), family = "nbinom1", data = eph.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
AIC(m1, m1.nb, m1.nb1)
##        df      AIC
## m1      7 3114.700
## m1.nb   8 1159.042
## m1.nb1  8 1184.955
#nbinom2 best family

#test each rdm interaction

m2 <- glmmTMB(abun ~ Microsite * RDM + rdm.cov  + arid + ESI + (1|site), family = "nbinom2", data = eph.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m2)
##  Family: nbinom2  ( log )
## Formula:          abun ~ Microsite * RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph.abun
## 
##      AIC      BIC   logLik deviance df.resid 
##   1159.8   1185.9   -570.9   1141.8      124 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1887   0.4344  
## Number of obs: 133, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.76 
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.111731   0.382577   8.134 4.17e-16 ***
## Micrositeephedra      0.415532   0.231399   1.796   0.0725 .  
## RDM                   0.014737   0.038613   0.382   0.7027    
## rdm.cov              -0.008278   0.003995  -2.072   0.0383 *  
## arid                  0.031002   0.023382   1.326   0.1849    
## ESI                  -2.623553   2.305880  -1.138   0.2552    
## Micrositeephedra:RDM  0.055679   0.050754   1.097   0.2726    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- glmmTMB(abun ~ Microsite * rdm.cov + RDM + arid + ESI + (1|site), family = "nbinom2", data = eph.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m3)
##  Family: nbinom2  ( log )
## Formula:          abun ~ Microsite * rdm.cov + RDM + arid + ESI + (1 | site)
## Data: eph.abun
## 
##      AIC      BIC   logLik deviance df.resid 
##   1158.4   1184.5   -570.2   1140.4      124 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1924   0.4387  
## Number of obs: 133, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.78 
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               3.150702   0.383592   8.214  < 2e-16 ***
## Micrositeephedra          0.306962   0.240772   1.275  0.20234    
## rdm.cov                  -0.014116   0.005046  -2.797  0.00516 ** 
## RDM                       0.047844   0.034837   1.373  0.16964    
## arid                      0.029886   0.023566   1.268  0.20474    
## ESI                      -2.470902   2.293169  -1.078  0.28125    
## Micrositeephedra:rdm.cov  0.008749   0.005362   1.632  0.10274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(m1.nb, m2, m3)
##       df      AIC
## m1.nb  8 1159.042
## m2     9 1159.842
## m3     9 1158.444
#interactions are non-significant and do not improve the model. Dropping.
anova(m1.nb, m2)
## Data: eph.abun
## Models:
## m1.nb: abun ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
## m2: abun ~ Microsite * RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m1.nb  8 1159.0 1182.2 -571.52   1143.0                         
## m2     9 1159.8 1185.8 -570.92   1141.8 1.2002      1     0.2733
anova(m1.nb, m3)
## Data: eph.abun
## Models:
## m1.nb: abun ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
## m3: abun ~ Microsite * rdm.cov + RDM + arid + ESI + (1 | site), zi=~0, disp=~1
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m1.nb  8 1159.0 1182.2 -571.52   1143.0                         
## m3     9 1158.4 1184.5 -570.22   1140.4 2.5974      1      0.107
check_collinearity(m1.nb)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Term  VIF Increased SE Tolerance
##  Microsite 1.46         1.21      0.69
##        RDM 1.73         1.31      0.58
##    rdm.cov 2.33         1.53      0.43
##       arid 1.25         1.12      0.80
##        ESI 1.26         1.12      0.79
summary(m1.nb)
##  Family: nbinom2  ( log )
## Formula:          abun ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph.abun
## 
##      AIC      BIC   logLik deviance df.resid 
##   1159.0   1182.2   -571.5   1143.0      125 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1865   0.4319  
## Number of obs: 133, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.75 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.004877   0.367399   8.179 2.87e-16 ***
## Micrositeephedra  0.586142   0.173345   3.381 0.000721 ***
## RDM               0.035951   0.034697   1.036 0.300136    
## rdm.cov          -0.008758   0.004015  -2.181 0.029162 *  
## arid              0.031663   0.023284   1.360 0.173881    
## ESI              -2.180906   2.266085  -0.962 0.335843    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(m1.nb, type = 2))
Chisq Df Pr(>Chisq)
Microsite 11.4335866 1 0.0007213
RDM 1.0735891 1 0.3001362
rdm.cov 4.7580001 1 0.0291624
arid 1.8491650 1 0.1738805
ESI 0.9262357 1 0.3358430

Insect morphospecies richness

m4 <- glmmTMB(Species ~ Microsite + RDM + rdm.cov +  arid + ESI + (1|site), family = "nbinom2", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
summary(m4)
##  Family: nbinom2  ( log )
## Formula:          Species ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##       NA       NA       NA       NA      127 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03885  0.1971  
## Number of obs: 135, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.04e+07 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       2.0086656  0.1679635  11.959  < 2e-16 ***
## Micrositeephedra  0.1774510  0.0749860   2.366  0.01796 *  
## RDM               0.0044239  0.0166705   0.265  0.79072    
## rdm.cov          -0.0003282  0.0018611  -0.176  0.86001    
## arid              0.0121462  0.0106088   1.145  0.25224    
## ESI              -2.9892639  1.0635163  -2.811  0.00494 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#nb models don't coverge -> likely suggests overdispersion not an issue but should add a formal check

eph <- filter(eph, Microsite != "larrea")
m4.p <- glmmTMB(Species ~ Microsite+RDM + rdm.cov +  arid + ESI +(1|site), family = "poisson", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m4.p)
##  Family: poisson  ( log )
## Formula:          Species ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##    653.0    673.3   -319.5    639.0      128 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03885  0.1971  
## Number of obs: 135, groups:  site, 9
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       2.0086635  0.1679633  11.959  < 2e-16 ***
## Micrositeephedra  0.1774498  0.0749859   2.366  0.01796 *  
## RDM               0.0044244  0.0166705   0.265  0.79070    
## rdm.cov          -0.0003282  0.0018611  -0.176  0.86000    
## arid              0.0121464  0.0106088   1.145  0.25223    
## ESI              -2.9892683  1.0635148  -2.811  0.00494 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_overdispersion(m4.p)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## # Overdispersion test
## 
##        dispersion ratio =   0.878
##   Pearson's Chi-Squared = 112.362
##                 p-value =   0.836
## No overdispersion detected.
#wow wow, I can actually use a poisson model 

#check each interaction

m5 <- glmmTMB(Species ~ Microsite * RDM + rdm.cov +  arid + ESI + (1|site), family = "poisson", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m5)
##  Family: poisson  ( log )
## Formula:          Species ~ Microsite * RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##    655.0    678.2   -319.5    639.0      127 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03886  0.1971  
## Number of obs: 135, groups:  site, 9
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.0079247  0.1728937  11.614  < 2e-16 ***
## Micrositeephedra      0.1786359  0.0995794   1.794  0.07283 .  
## RDM                   0.0045879  0.0189565   0.242  0.80876    
## rdm.cov              -0.0003303  0.0018644  -0.177  0.85938    
## arid                  0.0121466  0.0106104   1.145  0.25230    
## ESI                  -2.9858511  1.0802178  -2.764  0.00571 ** 
## Micrositeephedra:RDM -0.0004482  0.0247744  -0.018  0.98557    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- glmmTMB(Species ~ Microsite * rdm.cov + RDM +  arid + ESI + (1|site), family = "poisson", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
summary(m6)
##  Family: poisson  ( log )
## Formula:          Species ~ Microsite * rdm.cov + RDM + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##    654.2    677.4   -319.1    638.2      127 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.0375   0.1936  
## Number of obs: 135, groups:  site, 9
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.043432   0.170035  12.018  < 2e-16 ***
## Micrositeephedra          0.110769   0.106358   1.041  0.29765    
## rdm.cov                  -0.001821   0.002534  -0.719  0.47238    
## RDM                       0.007735   0.017132   0.452  0.65161    
## arid                      0.012136   0.010459   1.160  0.24589    
## ESI                      -3.086352   1.055792  -2.923  0.00346 ** 
## Micrositeephedra:rdm.cov  0.002327   0.002633   0.884  0.37692    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m7 <- glmmTMB(Species ~ Microsite * ESI+ rdm.cov + RDM +  arid + (1|site), family = "poisson", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
AIC(m4.p, m5, m6)
##      df      AIC
## m4.p  7 652.9706
## m5    8 654.9703
## m6    8 654.1833
anova(m4.p, m5)
## Data: eph
## Models:
## m4.p: Species ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
## m5: Species ~ Microsite * RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
##      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## m4.p  7 652.97 673.31 -319.49   638.97                        
## m5    8 654.97 678.21 -319.49   638.97 3e-04      1     0.9856
anova(m4.p, m6)
## Data: eph
## Models:
## m4.p: Species ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site), zi=~0, disp=~1
## m6: Species ~ Microsite * rdm.cov + RDM + arid + ESI + (1 | site), zi=~0, disp=~1
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m4.p  7 652.97 673.31 -319.49   638.97                         
## m6    8 654.18 677.43 -319.09   638.18 0.7873      1     0.3749
#no interactions 

check_collinearity(m4.p)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##       Term  VIF Increased SE Tolerance
##  Microsite 1.31         1.14      0.76
##        RDM 1.65         1.29      0.60
##    rdm.cov 2.06         1.44      0.49
##       arid 1.25         1.12      0.80
##        ESI 1.27         1.13      0.79
summary(m4.p)
##  Family: poisson  ( log )
## Formula:          Species ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##    653.0    673.3   -319.5    639.0      128 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03885  0.1971  
## Number of obs: 135, groups:  site, 9
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       2.0086635  0.1679633  11.959  < 2e-16 ***
## Micrositeephedra  0.1774498  0.0749859   2.366  0.01796 *  
## RDM               0.0044244  0.0166705   0.265  0.79070    
## rdm.cov          -0.0003282  0.0018611  -0.176  0.86000    
## arid              0.0121464  0.0106088   1.145  0.25223    
## ESI              -2.9892683  1.0635148  -2.811  0.00494 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(m4.p, type = 2))
Chisq Df Pr(>Chisq)
Microsite 5.6000444 1 0.0179600
RDM 0.0704384 1 0.7906996
rdm.cov 0.0311081 1 0.8599993
arid 1.3108917 1 0.2522333
ESI 7.9002850 1 0.0049427 
## Richness no singletons

m1nos <- glmmTMB(Species.nos ~ Microsite+RDM + rdm.cov +  arid + ESI +(1|site), family = "poisson", data = eph)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m1nos)
##  Family: poisson  ( log )
## Formula:          
## Species.nos ~ Microsite + RDM + rdm.cov + arid + ESI + (1 | site)
## Data: eph
## 
##      AIC      BIC   logLik deviance df.resid 
##    645.9    666.2   -316.0    631.9      128 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03739  0.1934  
## Number of obs: 135, groups:  site, 9
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       1.9783008  0.1664840  11.883  < 2e-16 ***
## Micrositeephedra  0.1538798  0.0770390   1.997  0.04578 *  
## RDM               0.0061796  0.0169556   0.364  0.71552    
## rdm.cov           0.0001445  0.0019122   0.076  0.93978    
## arid              0.0092079  0.0105255   0.875  0.38167    
## ESI              -2.9138193  1.0551737  -2.761  0.00575 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(m1nos, type = 2))
Chisq Df Pr(>Chisq)
Microsite 3.9897181 1 0.0457787
RDM 0.1328291 1 0.7155166
rdm.cov 0.0057068 1 0.9397825
arid 0.7653145 1 0.3816705
ESI 7.6256576 1 0.0057544 
check_overdispersion(m1nos)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## # Overdispersion test
## 
##        dispersion ratio =   0.870
##   Pearson's Chi-Squared = 111.371
##                 p-value =   0.852
## No overdispersion detected.

Influence of shrub size

#redo models
#shrub microsite only
shr <- filter(eph, Microsite == "ephedra")

m1.shr <- glmmTMB(Species ~ x + z + (1|site), family = "poisson", data = shr)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
summary(m1.shr)
##  Family: poisson  ( log )
## Formula:          Species ~ x + z + (1 | site)
## Data: shr
## 
##      AIC      BIC   logLik deviance df.resid 
##    368.1    377.1   -180.1    360.1       66 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03835  0.1958  
## Number of obs: 70, groups:  site, 9
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.0329051  0.2034844   9.990   <2e-16 ***
## x           -0.0004386  0.0006006  -0.730    0.465    
## z            0.0008097  0.0016605   0.488    0.626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(m1.shr, type = 2))
Chisq Df Pr(>Chisq)
x 0.5332531 1 0.4652424
z 0.2377630 1 0.6258265 
shr.abun <- filter(shr, abun < 350)
m2.shr <- glmmTMB(abun ~ x + z + (1|site), family = "nbinom2", data = shr.abun)
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you

## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
summary(m2.shr)
##  Family: nbinom2  ( log )
## Formula:          abun ~ x + z + (1 | site)
## Data: shr.abun
## 
##      AIC      BIC   logLik deviance df.resid 
##    634.4    645.7   -312.2    624.4       65 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.4685   0.6844  
## Number of obs: 70, groups:  site, 9
## 
## Overdispersion parameter for nbinom2 family (): 1.83 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.384052   0.557074   6.075 1.24e-09 ***
## x           -0.000739   0.001433  -0.516    0.606    
## z            0.001646   0.003736   0.441    0.659    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(m2.shr, type = 2))
Chisq Df Pr(>Chisq)
x 0.2659704 1 0.6060477
z 0.1942839 1 0.6593749

RII

Data wrangling and calculations for RII are in RII.R

Bootstrap the mean and 95% CI for RII rdm biomass and rdm cover. These calculations are made at the individual level (n = 30 for each of nine sites) and the RII represents all nine sites

shapiro.test(indrii$rii)
## 
##  Shapiro-Wilk normality test
## 
## data:  indrii$rii
## W = 0.90461, p-value = 4.143e-12
t.test(indrii$rii)
## 
##  One Sample t-test
## 
## data:  indrii$rii
## t = 10.369, df = 271, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2831132 0.4158190
## sample estimates:
## mean of x 
## 0.3494661
#not normal though

library(boot)
bootmean <- function(d, i) mean(d[i])

#rdm biomass

RDM.rii <- boot(indrii$rii, bootmean, R=10000, stype ="i")
ci.RDM.rii <- boot.ci(RDM.rii)
## Warning in boot.ci(RDM.rii): bootstrap variances needed for studentized
## intervals
bca.RDM <- ci.RDM.rii$bca

#rdm cover

rii.cover <- read.csv("clean_data/rii_individual_rdmcover.csv")
str(rii.cover)
## 'data.frame':    271 obs. of  30 variables:
##  $ X.1           : int  1 3 5 7 9 11 13 15 17 19 ...
##  $ rii           : num  0.5238 -0.0526 -0.027 -0.1111 -0.0323 ...
##  $ X             : int  1 3 5 7 9 11 13 15 17 19 ...
##  $ Region        : chr  "Panoche" "Panoche" "Panoche" "Panoche" ...
##  $ site          : chr  "PAN3" "PAN3" "PAN3" "PAN3" ...
##  $ ID            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Microsite     : chr  "ephedra" "ephedra" "ephedra" "ephedra" ...
##  $ RDM           : num  9.96 5.79 3.41 1.8 4.26 ...
##  $ x             : int  355 430 520 370 455 275 275 310 135 355 ...
##  $ y             : int  275 420 455 325 375 255 250 245 115 345 ...
##  $ z             : int  160 160 170 120 140 130 95 110 110 130 ...
##  $ bur.drip      : int  3 5 5 3 1 1 3 1 2 0 ...
##  $ bur.buf       : int  4 1 1 2 2 1 2 4 0 1 ...
##  $ max.veg.height: int  59 42 37 21 43 40 43 62 28 44 ...
##  $ rdm.cov       : int  80 90 90 60 90 85 80 80 95 90 ...
##  $ green.cov     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ non.veg.cov   : int  20 10 10 40 10 15 20 20 0 10 ...
##  $ largewood.cov : int  0 0 0 0 0 0 0 0 5 0 ...
##  $ lrgrocks.cov  : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ AT            : chr  NA NA NA NA ...
##  $ PF            : chr  NA NA NA NA ...
##  $ PF.out        : chr  NA NA NA NA ...
##  $ PF.in         : chr  NA NA NA NA ...
##  $ RDM.date      : chr  "june.24.morning" "june.24.morning" "june.24.morning" "june.24.morning" ...
##  $ comments      : chr  NA NA NA NA ...
##  $ ESI           : num  0.222 0.222 0.222 0.222 0.222 ...
##  $ Max           : num  35.3 35.3 35.3 35.3 35.3 35.3 35.3 35.3 35.3 35.3 ...
##  $ arid          : num  15.5 15.5 15.5 15.5 15.5 ...
##  $ label         : chr  "Semi-arid 6" "Semi-arid 6" "Semi-arid 6" "Semi-arid 6" ...
##  $ uniID         : chr  "PAN3 1" "PAN3 2" "PAN3 3" "PAN3 4" ...
mean(rii.cover$rii)
## [1] 0.4568816
cover.rii <- boot(rii.cover$rii, bootmean, R=1000, stype ="i")
ci.cover.rii <- boot.ci(cover.rii)
## Warning in boot.ci(cover.rii): bootstrap variances needed for studentized
## intervals
bca.cover <- ci.cover.rii$bca
t.test(rii.cover$rii)
## 
##  One Sample t-test
## 
## data:  rii.cover$rii
## t = 18.162, df = 270, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.4073540 0.5064093
## sample estimates:
## mean of x 
## 0.4568816
site.rii <- read.csv("clean_data/rii.csv")
#want to make a plot of site level rii vs average rii

site.rii <- dplyr::select(site.rii, 3:6, 13, 14, ESI, arid, Max)

#for each site, calculate the mean RDM biomass and cover for the whole site and also just under shrubs

#need to fix this
site.rii <- rdm %>% dplyr::select(site, RDM, rdm.cov) %>% group_by(site) %>% summarise(mean.rdm.biomass = mean(RDM), mean.rdm.cov = mean(rdm.cov)) %>% right_join(., site.rii, by = "site")

site.rii <- rdm %>% filter(Microsite == "ephedra") %>% dplyr::select(site, RDM, rdm.cov) %>% group_by(site) %>% summarise(shr.rdm.biomass = mean(RDM), shr.rdm.cov = mean(rdm.cov)) %>% right_join(., site.rii, by = "site")

Correlations

#relative effects
rel.rii <- pivot_wider(site.rii, site, names_from = response, values_from = average)

cor.test(rel.rii$richness, rel.rii$rdm)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$richness and rel.rii$rdm
## t = -0.40531, df = 7, p-value = 0.6973
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7410266  0.5700194
## sample estimates:
##        cor 
## -0.1514268
cor.test(rel.rii$richness, rel.rii$cov)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$richness and rel.rii$cov
## t = -0.84102, df = 7, p-value = 0.4281
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8050869  0.4521495
## sample estimates:
##        cor 
## -0.3029396
cor.test(rel.rii$abun, rel.rii$rdm)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$abun and rel.rii$rdm
## t = 0.41989, df = 7, p-value = 0.6872
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5663322  0.7434717
## sample estimates:
##       cor 
## 0.1567424
cor.test(rel.rii$abun, rel.rii$cov)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$abun and rel.rii$cov
## t = -0.043779, df = 7, p-value = 0.9663
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6732688  0.6547714
## sample estimates:
##         cor 
## -0.01654473
cor.test(rel.rii$rdm, rel.rii$cov)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$rdm and rel.rii$cov
## t = 5.1266, df = 7, p-value = 0.001359
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5478209 0.9764779
## sample estimates:
##       cor 
## 0.8886384
cor.test(rel.rii$abun, rel.rii$richness)
## 
##  Pearson's product-moment correlation
## 
## data:  rel.rii$abun and rel.rii$richness
## t = 0.66533, df = 7, p-value = 0.5271
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5014641  0.7814343
## sample estimates:
##       cor 
## 0.2438766
#environmental: aridity and stress
site.rii %>% filter(response == "richness") %>% cor.test(~ESI + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and average
## t = 0.81902, df = 7, p-value = 0.4398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4584381  0.8022778
## sample estimates:
##       cor 
## 0.2957173
m <- site.rii %>% filter(response == "richness") %>% lm(average ~ ESI,.)
summary(m)
## 
## Call:
## lm(formula = average ~ ESI, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16821 -0.04948 -0.00767  0.05944  0.13097 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.04775    0.06155   0.776    0.463
## ESI          0.34680    0.42343   0.819    0.440
## 
## Residual standard error: 0.101 on 7 degrees of freedom
## Multiple R-squared:  0.08745,    Adjusted R-squared:  -0.04292 
## F-statistic: 0.6708 on 1 and 7 DF,  p-value: 0.4398
site.rii %>% filter(response == "richness") %>% cor.test(~arid + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and average
## t = -0.50959, df = 7, p-value = 0.626
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7580385  0.5432228
## sample estimates:
##        cor 
## -0.1891309
site.rii %>% filter(response == "richness") %>% cor.test(~Max + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  Max and average
## t = 0.64203, df = 7, p-value = 0.5413
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5078350  0.7780838
## sample estimates:
##       cor 
## 0.2358205
site.rii %>% filter(response == "richness") %>% cor.test(~shr.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.biomass and average
## t = -0.17915, df = 7, p-value = 0.8629
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7002612  0.6245870
## sample estimates:
##         cor 
## -0.06755779
site.rii %>% filter(response == "richness") %>% cor.test(~shr.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.cov and average
## t = 0.1886, df = 7, p-value = 0.8558
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6224096  0.7020721
## sample estimates:
##        cor 
## 0.07110297
site.rii %>% filter(response == "richness") %>% cor.test(~mean.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.biomass and average
## t = 0.34978, df = 7, p-value = 0.7368
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.583880  0.731513
## sample estimates:
##      cor 
## 0.131064
site.rii %>% filter(response == "richness") %>% cor.test(~mean.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.cov and average
## t = 0.43952, df = 6, p-value = 0.6757
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6031214  0.7837451
## sample estimates:
##       cor 
## 0.1766125
site.rii %>% filter(response == "abun") %>% cor.test(~ESI + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and average
## t = -0.67725, df = 7, p-value = 0.52
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7831296  0.4981866
## sample estimates:
##        cor 
## -0.2479814
site.rii %>% filter(response == "abun") %>% cor.test(~arid + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and average
## t = -0.9018, df = 7, p-value = 0.3971
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8126300  0.4346206
## sample estimates:
##        cor 
## -0.3226235
site.rii %>% filter(response == "abun") %>% cor.test(~Max + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  Max and average
## t = 0.79579, df = 7, p-value = 0.4523
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4650481  0.7992638
## sample estimates:
##       cor 
## 0.2880319
site.rii %>% filter(response == "abun") %>% cor.test(~shr.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.biomass and average
## t = -0.80372, df = 7, p-value = 0.448
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8002979  0.4627961
## sample estimates:
##        cor 
## -0.2906612
site.rii %>% filter(response == "abun") %>% cor.test(~shr.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.cov and average
## t = -0.79074, df = 7, p-value = 0.455
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7986025  0.4664798
## sample estimates:
##        cor 
## -0.2863543
site.rii %>% filter(response == "abun") %>% cor.test(~mean.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.biomass and average
## t = -0.76906, df = 7, p-value = 0.467
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7957383  0.4726050
## sample estimates:
##        cor 
## -0.2791249
site.rii %>% filter(response == "abun") %>% cor.test(~mean.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.cov and average
## t = -0.5613, df = 6, p-value = 0.5949
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8018284  0.5712209
## sample estimates:
##        cor 
## -0.2233596
site.rii %>% filter(response == "rdm") %>% cor.test(~ESI + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and average
## t = -0.43658, df = 7, p-value = 0.6756
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7462433  0.5620882
## sample estimates:
##        cor 
## -0.1628094
site.rii %>% filter(response == "rdm") %>% cor.test(~arid + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and average
## t = -2.0059, df = 7, p-value = 0.08488
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9051156  0.1001476
## sample estimates:
##        cor 
## -0.6041567
m <- site.rii %>% filter(response == "rdm") %>% lm(average~arid, data =.)
summary(m)
## 
## Call:
## lm(formula = average ~ arid, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.312373 -0.100267  0.000974  0.153395  0.257657 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.598366   0.148208   4.037  0.00495 **
## arid        -0.018860   0.009402  -2.006  0.08488 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2173 on 7 degrees of freedom
## Multiple R-squared:  0.365,  Adjusted R-squared:  0.2743 
## F-statistic: 4.024 on 1 and 7 DF,  p-value: 0.08488
site.rii %>% filter(response == "rdm") %>% cor.test(~Max + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  Max and average
## t = 1.9479, df = 7, p-value = 0.09245
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1174963  0.9018910
## sample estimates:
##       cor 
## 0.5928897
site.rii %>% filter(response == "rdm") %>% cor.test(~shr.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.biomass and average
## t = -0.27996, df = 7, p-value = 0.7876
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7190967  0.6008855
## sample estimates:
##       cor 
## -0.105229
site.rii %>% filter(response == "rdm") %>% cor.test(~shr.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.cov and average
## t = -0.6024, df = 7, p-value = 0.5659
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7722651  0.5185747
## sample estimates:
##        cor 
## -0.2220047
site.rii %>% filter(response == "rdm") %>% cor.test(~mean.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.biomass and average
## t = -1.6674, df = 7, p-value = 0.1394
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8842032  0.2027436
## sample estimates:
##        cor 
## -0.5331672
site.rii %>% filter(response == "rdm") %>% cor.test(~mean.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.cov and average
## t = -1.1814, df = 6, p-value = 0.2821
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8721181  0.3894789
## sample estimates:
##       cor 
## -0.434424
site.rii %>% filter(response == "cov") %>% cor.test(~ESI + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and average
## t = -0.60923, df = 7, p-value = 0.5616
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7732783  0.5167334
## sample estimates:
##        cor 
## -0.2243947
site.rii %>% filter(response == "cov") %>% cor.test(~arid + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and average
## t = -1.2092, df = 7, p-value = 0.2658
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8461993  0.3431783
## sample estimates:
##        cor 
## -0.4156825
site.rii %>% filter(response == "cov") %>% cor.test(~Max + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  Max and average
## t = 1.1267, df = 7, p-value = 0.297
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3681049  0.8378976
## sample estimates:
##       cor 
## 0.3917983
site.rii %>% filter(response == "cov") %>% cor.test(~shr.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.biomass and average
## t = -0.39447, df = 7, p-value = 0.705
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7391945  0.5727479
## sample estimates:
##        cor 
## -0.1474661
site.rii %>% filter(response == "cov") %>% cor.test(~shr.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  shr.rdm.cov and average
## t = -0.74966, df = 7, p-value = 0.4779
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7931385  0.4780605
## sample estimates:
##        cor 
## -0.2726133
site.rii %>% filter(response == "cov") %>% cor.test(~mean.rdm.biomass + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.biomass and average
## t = -1.9238, df = 7, p-value = 0.09579
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9005083  0.1247493
## sample estimates:
##        cor 
## -0.5880954
site.rii %>% filter(response == "cov") %>% cor.test(~mean.rdm.cov + average, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  mean.rdm.cov and average
## t = -1.6455, df = 6, p-value = 0.151
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9062106  0.2422235
## sample estimates:
##        cor 
## -0.5576304
site.rii %>% filter(response == "rdm") %>% cor.test(~ESI+ mean.rdm.biomass, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and mean.rdm.biomass
## t = 3.7477, df = 7, p-value = 0.007189
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3340583 0.9601396
## sample estimates:
##      cor 
## 0.816938
m2 <- site.rii %>% filter(response == "rdm") %>% lm(mean.rdm.biomass ~ ESI, data =.) 
summary(m2)
## 
## Call:
## lm(formula = mean.rdm.biomass ~ ESI, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72274 -0.49571  0.07935  0.54021  0.71961 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   1.0053     0.3687   2.727  0.02947 * 
## ESI           9.5061     2.5365   3.748  0.00719 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6052 on 7 degrees of freedom
## Multiple R-squared:  0.6674, Adjusted R-squared:  0.6199 
## F-statistic: 14.05 on 1 and 7 DF,  p-value: 0.007189
site.rii %>% filter(response == "rdm") %>% cor.test(~arid+ mean.rdm.biomass, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and mean.rdm.biomass
## t = 2.3414, df = 7, p-value = 0.05174
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.002506808  0.921336281
## sample estimates:
##       cor 
## 0.6627182
site.rii %>% filter(response == "rdm") %>% cor.test(~ESI+ mean.rdm.cov, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and mean.rdm.cov
## t = 4.5229, df = 6, p-value = 0.004005
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4591782 0.9779963
## sample estimates:
##       cor 
## 0.8793269
site.rii %>% filter(response == "rdm") %>% cor.test(~arid+ mean.rdm.cov, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and mean.rdm.cov
## t = 2.1027, df = 6, p-value = 0.08017
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0985471  0.9294286
## sample estimates:
##       cor 
## 0.6513584
site.rii %>% filter(response == "rdm") %>% cor.test(~ESI+ shr.rdm.biomass, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and shr.rdm.biomass
## t = 4.4989, df = 7, p-value = 0.002802
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4628129 0.9705211
## sample estimates:
##       cor 
## 0.8619899
site.rii %>% filter(response == "rdm") %>% cor.test(~arid+ shr.rdm.biomass, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and shr.rdm.biomass
## t = 1.806, df = 7, p-value = 0.1139
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1603852  0.8933998
## sample estimates:
##       cor 
## 0.5637886
site.rii %>% filter(response == "rdm") %>% cor.test(~ESI+ shr.rdm.cov, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  ESI and shr.rdm.cov
## t = 5.7814, df = 7, p-value = 0.0006763
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6189772 0.9810072
## sample estimates:
##       cor 
## 0.9093052
site.rii %>% filter(response == "rdm") %>% cor.test(~arid+ shr.rdm.cov, data =.)
## 
##  Pearson's product-moment correlation
## 
## data:  arid and shr.rdm.cov
## t = 2.3075, df = 7, p-value = 0.05439
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01212914  0.91986891
## sample estimates:
##       cor 
## 0.6572872
ggplot(filter(rii,response == "abun"), aes(reorder(site, arid), arid)) + geom_bar(stat="identity")

Figures

Figure 1

#richness
lbls <- read.csv("clean_data/sitelabels.csv")
eph <- left_join(eph, lbls, by = "site")

eph$label <- as.factor(eph$label)
eph$label <- factor(eph$label, levels = c("Semi-arid 1", "Semi-arid 2", "Semi-arid 3", "Semi-arid 4", "Semi-arid 5", "Semi-arid 6", "Arid 7", "Arid 8", "Arid 9"))

c <- ggplot(eph, aes(label, Species, fill = Microsite)) + geom_boxplot() + ylab("Arthropod Richness") + theme(legend.position = "bottom", legend.box = "horizontal") + scale_fill_manual(values=c("white", "gray"))
l1 <- get_legend(c)

rdm <- left_join(rdm, lbls, by = "site")

rdm$label <- factor(rdm$label, levels = c("Semi-arid 1", "Semi-arid 2", "Semi-arid 3", "Semi-arid 4", "Semi-arid 5", "Semi-arid 6", "Arid 7", "Arid 8", "Arid 9"))

rdm$label <- as.factor(rdm$label)
eph$label <- as.factor(eph$label)

a <- ggplot(rdm, aes(label, RDM, fill = Microsite)) +  
 stat_boxplot(geom ='errorbar')+ geom_boxplot() + ylab("RDM Biomass (g)") + theme(legend.position = "none") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + stat_summary(fun=mean, colour="black", geom="point", shape=18, size=3, position=position_dodge(.8)) + scale_fill_manual(values=c("white", "gray")) + theme(legend.box = "horizontal", legend.position = c(0.9, 0.9), legend.title = element_blank())
a

b <- ggplot(rdm, aes(label, rdm.cov, fill = Microsite))  +  
 stat_boxplot(geom ='errorbar')+ geom_boxplot() + ylab("RDM Cover (%)") + theme(legend.position = "none")+ xlab("")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + stat_summary(fun=mean, colour="black", geom="point", shape=18, size=3, position=position_dodge(.8))+ scale_fill_manual(values=c("white", "gray"))


c <- ggplot(eph, aes(label, Species, fill = Microsite)) +  
 stat_boxplot(geom ='errorbar') + geom_boxplot() + ylab("Arthropod Richness") + theme(legend.position = "none")+ xlab("")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + stat_summary(fun=mean, colour="black", geom="point", shape=18, size=3, position=position_dodge(.8))+ scale_fill_manual(values=c("white", "gray"))

d <- ggplot(filter(eph, abun < 350), aes(label, abun, fill = Microsite)) +  
 stat_boxplot(geom ='errorbar') + geom_boxplot() + ylab("Arthropod Abundance (captures/trap)") + theme(legend.position = "none") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + xlab("Site") + stat_summary(fun=mean, colour="black", geom="point", shape=18, size=3, position=position_dodge(.8))+ scale_fill_manual(values=c("white", "gray"))


plots <- list(a,b,c,d)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing non-finite values (stat_summary).
maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 1))

p
## TableGrob (4 x 1) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
## 3 3 (3-3,1-1) arrange gtable[layout]
## 4 4 (4-4,1-1) arrange gtable[layout]
#grid.arrange(p, l1, ncol = 1,  heights = c(1, 0.1))

Figure 2 RII

full.rii <- read.csv("clean_data/fullrii.csv")
full.rii <- mutate(full.rii, ci = 1.96*error)

full.rii$response <- factor(full.rii$metric, levels = c("abun", "richness", "cov", "rdm"))
full.rii$site <- "Overall"
test <- bind_rows(rii, full.rii)

rii <- mutate(rii, ci = error*1.96)

test$response <- factor(test$response, levels = c("rdm", "cov", "richness", "abun"))

test <- left_join(test, lbls, by = "site")
test$label[is.na(test$label)] <- "Pooled Sites"

test$label <- factor(test$label, levels = c("Pooled Sites", "Semi-arid 1", "Semi-arid 2", "Semi-arid 3", "Semi-arid 4", "Semi-arid 5", "Semi-arid 6", "Arid 7", "Arid 8", "Arid 9"))
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
group.colors <- c("#000000", hue_pal()(9))
names( group.colors ) <- LETTERS[1:4]; group.colors
##         A         B         C         D      <NA>      <NA>      <NA>      <NA> 
## "#000000" "#F8766D" "#D39200" "#93AA00" "#00BA38" "#00C19F" "#00B9E3" "#619CFF" 
##      <NA>      <NA> 
## "#DB72FB" "#FF61C3"
rii.plot <- ggplot(test, aes(x=response, y=average)) + geom_point(aes(response, average,  color = label), position = position_dodge(width = 0.7),  size=4)+ geom_errorbar(aes(ymin=average-error, ymax=average + error, color = label), width=.1, position= position_dodge(width = 0.7))  + ylim(-1, 1) + theme(axis.text.x=element_text(angle=90, vjust=.5)) + ylab("RII") + xlab("") + geom_hline(aes(yintercept=0))  + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + scale_x_discrete(labels = c("RDM Biomass (g)", "RDM Cover (%)", "Arthropod Richness", "Arthropod Abundance \n (captures/trap)")) +  scale_colour_manual( values= group.colors, name = "Sites")+ theme(text = element_text(size = 25))

rii.plot 

#geom_errorbar(aes(ymin=average-ci, ymax=average + ci), width=.1, position= position_dodge(width = 0.5)) + geom_point(aes(metric, average), position=position_dodge(width=0.5), colour="black", fill = "black", pch=21, size=4) + ylim(1, -1) + theme(axis.text.x=element_text(angle=90, vjust=.5)) + ylab("RII") + xlab("") + geom_hline(aes(yintercept=0))  + coord_flip()

Figure 3 Scatterplots/correlation coefficients

fit1<- lm(richness ~ rdm, rel.rii)
a <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0, y = 0.23),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5)))+ ylab("Mean RII Arthropod Richness") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))


fit1 <- lm(richness ~ cov, rel.rii)
b <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.25, y = 0.23),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5)))  + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))
b
## `geom_smooth()` using formula 'y ~ x'

fit1 <- lm(abun ~ rdm, rel.rii)
c <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = -0, y = 0.65),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + xlab("Mean RII RDM Biomass") + ylab("Mean RII Arthropod Abundance") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))
c
## `geom_smooth()` using formula 'y ~ x'

fit1 <- lm(abun ~ cov, rel.rii)
d <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.24, y = 0.65),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + xlab("Mean RII RDM cover") + ylab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))

d
## `geom_smooth()` using formula 'y ~ x'

grid.arrange(a, b, c, d, ncol = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Figure 4 - Aridty and ESI vs RII

#so many scatters, let's do a 3 by 2


fit1 <- filter(site.rii, response == "rdm") %>% lm(mean.rdm.biomass ~ arid,.)
a <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 3.5),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("Mean RDM Biomass (g)") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))
a
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "rdm") %>% lm(mean.rdm.biomass ~ ESI,.)
b <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.03, y = 3.5),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
b
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "rdm") %>% lm(mean.rdm.cov ~ arid,.)
c <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 60),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("Mean RDM Cover (%)") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
c
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "rdm") %>% lm(mean.rdm.cov ~ ESI,.)
d <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.03, y = 60),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
d
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "rdm") %>% lm(average ~ arid,.)
e <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 1),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + ylab("Mean RII RDM Biomass")+ theme(text = element_text(size = 18))
e
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "rdm") %>% lm(average ~ ESI,.)
f <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.03, y = 0.78),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + theme(text = element_text(size = 18))
f
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "cov") %>% lm(average ~ arid,.)
g <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 0.75),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("Mean RII Cover") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
g
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "cov") %>% lm(average ~ ESI,.)
h <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.04, y = 0.7),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
h
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "richness") %>% lm(average ~ arid,.)
i <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 0.24),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("Mean RII Richness") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))

i
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "richness") %>% lm(average ~ ESI,.)
j <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.03, y = 0.3),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
j
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "abun") %>% lm(average ~ arid,.)
k <- ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 4, y = 0.8),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + xlab("Aridity") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm")) + ylab("Mean RII Abundance")+ theme(text = element_text(size = 18))
k
## `geom_smooth()` using formula 'y ~ x'

fit1 <- filter(site.rii, response == "abun") %>% lm(average ~ ESI,.)
l <-  ggplot(fit1$model, aes_string(x = names(fit1$model)[2], y = names(fit1$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  geom_label(aes(x = 0.04, y = 0.8),vjust=1, hjust = 0, 
             label = paste("Adj r2 = ",signif(summary(fit1)$adj.r.squared, 5),
 " \nP =",signif(summary(fit1)$coef[2,4], 5))) + ylab("") + xlab("Evaporative Stress Index") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))  + theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))+ theme(text = element_text(size = 18))
l
## `geom_smooth()` using formula 'y ~ x'

#grid.arrange(a, b, c, d, e, f, g, h, i, j, k, l, ncol = 2)


plots <- list(a,b,c,d,e,f,g,h,i,j,k,l)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 2))

p
## TableGrob (6 x 2) "arrange": 12 grobs
##     z     cells    name           grob
## 1   1 (1-1,1-1) arrange gtable[layout]
## 2   2 (1-1,2-2) arrange gtable[layout]
## 3   3 (2-2,1-1) arrange gtable[layout]
## 4   4 (2-2,2-2) arrange gtable[layout]
## 5   5 (3-3,1-1) arrange gtable[layout]
## 6   6 (3-3,2-2) arrange gtable[layout]
## 7   7 (4-4,1-1) arrange gtable[layout]
## 8   8 (4-4,2-2) arrange gtable[layout]
## 9   9 (5-5,1-1) arrange gtable[layout]
## 10 10 (5-5,2-2) arrange gtable[layout]
## 11 11 (6-6,1-1) arrange gtable[layout]
## 12 12 (6-6,2-2) arrange gtable[layout]
full.rii <- mutate(full.rii, upperCI = average + ci, lowerci = average - ci)

Appendix figures

ggplot(eph, aes(RDM, fill = Microsite)) + geom_density(alpha = 0.7) + theme_classic() + xlab("RDM (g)")

# Conclusions

  • Ephedra signficantly influences insect abundance and richness
  • Ephedra supports greater plant biomass