Shrub Density influences on Animals

The following is a review of shrubs that have been reported as benefactor species on local animal communities. Web of Science resources were queried using the following key terms.

"facil* density* shrub”, “density shrub* facilitate” and “density shrub* animal* facil*"

Questions:

Is there evidence that shrub studies that reported facilitation of other plants or of animals influence measures of animal community composition?
Is there are a relationship between animal abundance or richness and shrub or tree density? Is there a difference between the measures of animal community composition in studies that reported facilitation or not?

Hypothesis:

Higher density of foundation species that have reported facilitation support higher densities or more rich animal communities.

Predictions:

1. Reported benefactor plants better support animal communities.
   
2. Higher density of benefactor plants predict higher measures of animal community structure or composition.

Good Figures (June 8th 2021)

ggplot(data = data_e1, aes(yi)) +
  geom_histogram(bins = 40) + labs(x = "Incidence Rate for Animal Abundance", y = "Frequency") + theme_classic()

ggplot(summary_data, aes(Study.ID, abundance_rate)) +
  geom_point(size =2) +
  geom_errorbar(aes(ymin = abundance - var_abundance, ymax = abundance + var_abundance), size= 0.7, width=0.1, position = position_dodge(width = 0.5)) +
  coord_flip() +
  geom_hline(yintercept = 1, colour="grey", linetype = "longdash", size = 2) + 
  labs(x = "Study", y = "Incident Rate of Abundance") + theme_classic()

Meta Models

library(broom)
#x = density
#y = abundance
#With density
mod1.1 <- rma(yi = yi, sei = vi, method = "ML", test = "knha", control=list(stepadj=0.5), data = data_e1)

## Warning in rma(yi = yi, sei = vi, method = "ML", test = "knha", control =
## list(stepadj = 0.5), : Ratio of largest to smallest sampling variance extremely
## large. May not be able to obtain stable results.

summary(mod1.1)

## 
## Random-Effects Model (k = 113; tau^2 estimator: ML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -281.0118  1239.2260   566.0235   571.4783   566.1326   
## 
## tau^2 (estimated amount of total heterogeneity): 7.6846 (SE = 1.0334)
## tau (square root of estimated tau^2 value):      2.7721
## I^2 (total heterogeneity / total variability):   100.00%
## H^2 (total variability / sampling variability):  25834509573.18
## 
## Test for Heterogeneity:
## Q(df = 112) = 1049029.0227, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval    pval   ci.lb   ci.ub 
##   1.9008  0.2743  6.9287  <.0001  1.3572  2.4444  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod1.2 <- rma(yi = yi, sei = vi, method = "ML", mods = ~macrohabitat-1, test = "knha", control=list(stepadj=0.5), data = data_e1)

## Warning in rma(yi = yi, sei = vi, method = "ML", mods = ~macrohabitat - : Ratio
## of largest to smallest sampling variance extremely large. May not be able to
## obtain stable results.

summary(mod1.2)

## 
## Mixed-Effects Model (k = 113; tau^2 estimator: ML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -280.7982  1238.7988   567.5963   575.7785   567.8165   
## 
## tau^2 (estimated amount of residual heterogeneity):     7.6299 (SE = 1.0261)
## tau (square root of estimated tau^2 value):             2.7622
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   24746718787.32
## 
## Test for Residual Heterogeneity:
## QE(df = 111) = 775272.6960, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 111) = 24.0660, p-val < .0001
## 
## Model Results:
## 
##                        estimate      se    tval    pval   ci.lb   ci.ub 
## macrohabitatdesert       2.1345  0.4659  4.5815  <.0001  1.2113  3.0577  *** 
## macrohabitatgrassland    1.7750  0.3407  5.2097  <.0001  1.0999  2.4501  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.1 <- tidy(mod1.1) 
m1.1

## # A tibble: 1 x 6
##   term    type    estimate std.error statistic  p.value
##   <chr>   <chr>      <dbl>     <dbl>     <dbl>    <dbl>
## 1 overall summary     1.90     0.274      6.93 2.84e-10

plot(mod1.1)

###Weighted regression with all data

gmodel <- lm(standardized.density ~ abundance, data=data_all, weights = n_days)
gmodel

## 
## Call:
## lm(formula = standardized.density ~ abundance, data = data_all, 
##     weights = n_days)
## 
## Coefficients:
## (Intercept)    abundance  
##     4851.00        43.49

summary(gmodel)

## 
## Call:
## lm(formula = standardized.density ~ abundance, data = data_all, 
##     weights = n_days)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -197129  -31209  -16262  -14951  351038 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4851.00     990.55   4.897 3.33e-06 ***
## abundance      43.49      25.08   1.734   0.0856 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66830 on 111 degrees of freedom
## Multiple R-squared:  0.02639,    Adjusted R-squared:  0.01762 
## F-statistic: 3.008 on 1 and 111 DF,  p-value: 0.08562

Fig1 with Line

ggplot(data = data_e1, aes(abundance, standardized.density)) +
  geom_point() +
  geom_smooth(method = "lm", mapping = aes(weight = n_days), se = F) + labs(x = "Abundance", y = "Density") + theme_classic()

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

### Fig1 With Curve

fig1 <- ggplot(data = data_e1, aes(abundance, standardized.density, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Abundance", y = "Density") + theme_classic() + facet_wrap(~macrohabitat)
fig1

color_points <- ggplot(data = data_e1, aes(abundance, standardized.density, color = macrohabitat)) +
  geom_point() + scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue"))
color_points + theme(legend.title = element_blank()) + labs(x = "Abundance", y = "Density") +theme_classic()

model1 <- lm(standardized.density ~ abundance, data=data_all)
summary(model1)

## 
## Call:
## lm(formula = standardized.density ~ abundance, data = data_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -742.1  -682.4  -563.8  -480.3 12273.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 742.2721   254.0139   2.922  0.00421 **
## abundance    -0.4242     4.1289  -0.103  0.91835   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2193 on 111 degrees of freedom
## Multiple R-squared:  9.509e-05,  Adjusted R-squared:  -0.008913 
## F-statistic: 0.01056 on 1 and 111 DF,  p-value: 0.9184

plot(fitted(model1), resid(model1), xlab='Fitted Values', ylab='Residuals') + abline(0,0)

## integer(0)

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(model1)

## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 0.012607, df = 1, p-value = 0.9106

wt <- 1 / lm(abs(model1$residuals) ~ model1$fitted.values)$fitted.values^2
wls_model <- lm(standardized.density ~ abundance, data = data_all, weights=wt)
summary(wls_model)

## 
## Call:
## lm(formula = standardized.density ~ abundance, data = data_all, 
##     weights = wt)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7073 -0.6890 -0.5400 -0.4571 11.8927 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 705.8472   250.4437   2.818  0.00572 **
## abundance     0.6203     4.4588   0.139  0.88961   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.125 on 111 degrees of freedom
## Multiple R-squared:  0.0001743,  Adjusted R-squared:  -0.008833 
## F-statistic: 0.01935 on 1 and 111 DF,  p-value: 0.8896

data.frame(y = rstandard(model1),
           x = model1$fitted.values) %>%
  ggplot(aes(x = x, y = y)) + 
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(title = "Standardized Residuals vs Fitted Values Plot")

ggplot(data = data_all, aes(y = standardized.density, x = abundance)) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE,
              color = "black", 
              size = 0.5, 
              linetype = "dashed") +
  geom_smooth(method = lm, se = FALSE, 
              aes(weight = wt),
              color = "red", 
              size = 0.5,
              linetype = "dashed") +
  labs(title = "Scatterplot of Density ~ abundance")

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

###Trying to do regressions with summary data

cmodel <- lm(density ~ abundance_rate, data=summary_data, weights = )
cmodel

## 
## Call:
## lm(formula = density ~ abundance_rate, data = summary_data)
## 
## Coefficients:
##    (Intercept)  abundance_rate  
##           2460           -1015

summary(cmodel)

## 
## Call:
## lm(formula = density ~ abundance_rate, data = summary_data)
## 
## Residuals:
##       1       2       3       4       5       6 
## -1308.0  5552.8   650.8  -262.7 -2183.8 -2449.0 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)        2460       1808   1.361    0.245
## abundance_rate    -1016       1171  -0.867    0.435
## 
## Residual standard error: 3309 on 4 degrees of freedom
## Multiple R-squared:  0.1582, Adjusted R-squared:  -0.05221 
## F-statistic: 0.7519 on 1 and 4 DF,  p-value: 0.4348

ggplot(data = summary_data, aes(abundance, density)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

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

model <- lm(density ~ abundance, data=data_e1)
summary(model)

## 
## Call:
## lm(formula = density ~ abundance, data = data_e1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1836.7 -1469.3  -141.3   669.9  3776.7 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1839.576    181.629  10.128  < 2e-16 ***
## abundance     -9.481      2.952  -3.211  0.00173 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1568 on 111 degrees of freedom
## Multiple R-squared:  0.08502,    Adjusted R-squared:  0.07677 
## F-statistic: 10.31 on 1 and 111 DF,  p-value: 0.001728

plot(fitted(model), resid(model), xlab='Fitted Values', ylab='Residuals') + abline(0,0)

## integer(0)

library(lmtest)
bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1.6003, df = 1, p-value = 0.2059

wt <- 1 / lm(abs(model$residuals) ~ model$fitted.values)$fitted.values^2
wls_model <- lm(density ~ abundance, data = data_e1, weights=wt)
summary(wls_model)

## 
## Call:
## lm(formula = density ~ abundance, data = data_e1, weights = wt)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4650 -1.2501 -0.0978  0.5592  3.3273 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1817.007    182.458   9.958  < 2e-16 ***
## abundance     -8.885      2.464  -3.606 0.000467 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.318 on 111 degrees of freedom
## Multiple R-squared:  0.1049, Adjusted R-squared:  0.09682 
## F-statistic: 13.01 on 1 and 111 DF,  p-value: 0.0004666

data.frame(y = rstandard(model),
           x = model$fitted.values) %>%
  ggplot(aes(x = x, y = y)) + 
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(title = "Standardized Residuals vs Fitted Values Plot")

#Sampling effort
ggplot(data_e1, aes(n_days, yi)) +
  geom_point(aes(color = macrohabitat)) + 
  labs(x = "total days sampled per study", y = "abundance capture rates")

MyData <- read.csv("MyData.csv")

###Forest Plot for paper
figure3 <- ggplot(data_all, aes(as.factor(rep), abundance_rate, color = macrohabitat)) +
  geom_point(size =2) +
  geom_errorbar(aes(ymin = abundance_rate - var_abundancerate, ymax = abundance_rate + var_abundancerate), size= 0.4, width=0.1, position = position_dodge(width = 0.5)) +
  coord_flip() +
  geom_hline(yintercept = 1, colour="grey", linetype = "longdash", size = 0.5) + 
  labs(x = "Rep", y = "Incident Rate of Abundance") + theme_classic() + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),legend.title = element_blank(), axis.text = element_text(size = 10)) + labs(x = "Observations") + scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + facet_wrap(~macrohabitat)
  
figure3

gd <- data_all %>% group_by(macrohabitat) %>%
  summarise(abundance_rate = mean(abundance_rate), rep = rep/rep)

## `summarise()` has grouped output by 'macrohabitat'. You can override using the `.groups` argument.

grassland <- gd[-c(1:39),]
desert <- gd[-c(40:113),]

desert.forest <- ggplot(subset(data_all, macrohabitat %in% c("desert")), aes(as.factor(rep), abundance_rate, color = macrohabitat)) + 
  geom_point(size =2) +
  geom_errorbar(aes(ymin = abundance_rate - var_abundancerate, ymax = abundance_rate + var_abundancerate), size= 0.4, width=0.1, position = position_dodge(width = 0.5)) +
  coord_flip() + facet_wrap(~microsite, scales = "free") +
  geom_hline(yintercept = 1, colour="grey", linetype = "longdash", size = 1.5) + 
  labs(x = "Rep", y = "Incident Rate of Abundance") + theme_classic() + 
  theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank(),legend.title = element_blank(), axis.text = element_text(size = 10)) + labs(x = "Observations") + scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + facet_wrap(~macrohabitat, scales = "free") + geom_point(data = desert, size = 5, shape = "triangle")+scale_y_continuous(breaks=c(0,2,4,6,8))
desert.forest 

grassland.forest <- ggplot(subset(data_all, macrohabitat %in% c("grassland")), aes(as.factor(rep), abundance_rate, color = macrohabitat)) +
  geom_point() +
  geom_errorbar(aes(ymin = abundance_rate - var_abundancerate, ymax = abundance_rate + var_abundancerate), size= 0.4, width=0.1, position = position_dodge(width = 0.5)) +
  coord_flip() + facet_wrap(~microsite, scales = "free") +
  geom_hline(yintercept = 1, colour="grey", linetype = "longdash", size = 1.5) + 
  labs(x = "Rep", y = "Incident Rate of Abundance") + theme_classic() + 
  theme(axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), legend.position = "none", legend.title = element_blank(), axis.text = element_text(size = 10)) + labs(x = "Observations") + scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + facet_wrap(~macrohabitat, scales = "free") + geom_point(data = grassland, size = 5, shape = "triangle")
grassland.forest

###Figure is good! :D
forestfig <- ggarrange(desert.forest, grassland.forest, nrow = 1, align = "v", heights = 10)
forestfig

### Has abundance on x axis so does not work!!!
grassland <- ggplot(subset(data_e1, macrohabitat %in% c("grassland")), aes(abundance, standardized.density, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Animal Abundance", y = "Shrub Density") + theme_classic() + scale_color_manual(values = c("grassland" = "dark blue")) +theme(legend.position = "none") 
grassland

### Has abundance on x axis so does not work!!!
desert <- ggplot(subset(data_e1, macrohabitat %in% c("desert")), aes(abundance, standardized.density, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Animal Abundance", y = " Shrub Density") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none") + ylim(0,15)
desert

## Warning: Removed 26 rows containing non-finite values (stat_smooth).

## Warning: Removed 26 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_smooth).

scatterplot <- ggarrange(desert, grassland, nrow = 1, align = "v", heights = 10)

## Warning: Removed 26 rows containing non-finite values (stat_smooth).

## Warning: Removed 26 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_smooth).

scatterplot

##Try with natural log (Do not use this!)

data_1$natlog <- log(data_1$abundance)
log.grassland <- ggplot(subset(data_1, macrohabitat %in% c("grassland")), aes(standardized.density, natlog,  color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = "Log Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) +theme(legend.position = "none")

log.desert <- ggplot(subset(data_1, macrohabitat %in% c("desert")), aes(standardized.density, natlog, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = "Log Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")

logplot <- ggarrange(log.desert, log.grassland, nrow = 1, align = "v", heights = 10)
logplot

data_1$natlog.dense <- log(data_1$standardized.density)
log.grassland.dense <- ggplot(subset(data_1, macrohabitat %in% c("grassland")), aes(natlog.dense, abundance,  color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Log Shrub Density", y = "Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) +theme(legend.position = "none")

log.desert.dense <- ggplot(subset(data_1, macrohabitat %in% c("desert")), aes(natlog.dense, abundance, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Log Shrub Density", y = "Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")
logplot.dense <- ggarrange(log.desert.dense, log.grassland.dense, nrow = 1, align = "v", heights = 10)

## Warning: Removed 4 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

logplot.dense

Alex_paper <- read.csv("Carrizo_telemetry_with_shrub_density.csv")
Alex_paper_mean <- aggregate(Alex_paper[, 10], list(Alex_paper$year), mean)

Alex_paper_count <- Alex_paper %>%
 count(year, lizard)

library(dplyr)
Alex_paper_count <- Alex_paper_count %>% 
  group_by(year) %>% 
  summarise(Sum = sum(n))
colnames(Alex_paper_mean) <- c("year", "Average Density")

Alex_final <- merge(Alex_paper_mean, Alex_paper_count)

Test <- Alex_paper %>%
 count(lizard, year)

###Figures have density on x axis
desertnew <- ggplot(subset(data_e1, macrohabitat %in% c("desert")), aes(standardized.density, yi, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = " Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")
desertnew

grasslandnew <- ggplot(subset(data_e1, macrohabitat %in% c("grassland")), aes(standardized.density, yi, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = " Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")
grasslandnew

scatterplotnew <- ggarrange(desertnew, grasslandnew, nrow = 1, align = "v", heights = 10)
scatterplotnew

log.facet <- ggplot(data_1, aes(standardized.density, natlog,  color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = "Log Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) +theme(legend.position = "none")
log.facet

plot <- ggplot(data_e1, aes(standardized.density, yi, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x +I(x^2), mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density", y = "Animal Abundance") + theme_classic() + facet_wrap(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")
plot

###Main figure!
library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

yi.plot <- ggplot(data_e1, aes(standardized.density, yi, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density (per km^2)", y = "Animal Abundance") + theme_classic() + facet_grid(~macrohabitat, scales = "free")+ scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none") +  scale_y_continuous(breaks= pretty_breaks())
yi.plot

####Testing main data with no zeros

exampleDF <- data_e1 %>%
  group_split(zeroes = standardized.density == 0 | yi == 0)

nozeros <- ggplot(exampleDF[[1]], aes(x = standardized.density, y = yi)) + 
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F)  + 
  geom_point() +
  geom_point(data = exampleDF[[2]], aes(x = standardized.density, y = yi), color = "red") +
  xlab("Shrub Density") + 
  ylab("Abundance")
nozeros

yi.plot <- ggplot(data_e1, aes(standardized.density, yi, color = macrohabitat)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F) + labs(x = "Shrub Density (per km^2)", y = "Animal Abundance") + theme_classic() + scale_color_manual(values = c("desert" = "red", "grassland" = "dark blue")) + theme(legend.position = "none")
yi.plot

data_new <- data_e1 %>% filter(plant.species !="Unknown")

Caragana.K <- ggplot(subset(data_new, plant.species %in% c("Caragana korshinskii")), aes( standardized.density, yi, color = plant.species)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F) + labs(x = "", y = "Animal Abundance") + theme_classic() + facet_wrap(~plant.species) + scale_color_manual(values = c("Caragana korshinskii" = "dark blue")) + theme(legend.position = "none") +  scale_y_continuous(breaks= pretty_breaks())
Caragana.K

Caragana.M <- ggplot(subset(data_new, plant.species %in% c("Caragana microphylla")), aes( standardized.density, yi, color = plant.species)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F) + labs(x = "", y = " ") + theme_classic() + facet_wrap(~plant.species) + scale_color_manual(values = c("Caragana microphylla" = "red")) + theme(legend.position = "none") + ylim(0,0.05) + scale_y_continuous(breaks = pretty_breaks()) 

## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Caragana.M

Ephedra<- ggplot(subset(data_new, plant.species %in% c("Ephedra californica")), aes( standardized.density, yi, color = plant.species)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y~x, mapping = aes(weight = n_days), se = F) + labs(x = "", y = "") + theme_classic() + facet_wrap(~plant.species)  + scale_color_manual(values = c("Ephedra californica" = "red")) + theme(legend.position = "none") + scale_y_continuous(breaks = pretty_breaks()) 
Ephedra

fig3 <- ggarrange(Caragana.K, Caragana.M, Ephedra, ncol = 2, nrow = 2)
fig3 <- annotate_figure(fig3, bottom = text_grob("Shrub Density (per km^2)", 
               color = "black", face = "bold", size = 10))
fig3

###By shrub species
mod1.3 <- rma(yi = yi, sei = vi, method = "ML", mods = ~plant.species-1, test = "knha", control=list(stepadj=0.5), data = data_new)

## Warning in rma(yi = yi, sei = vi, method = "ML", mods = ~plant.species - : Ratio
## of largest to smallest sampling variance extremely large. May not be able to
## obtain stable results.

summary(mod1.3)

## 
## Mixed-Effects Model (k = 112; tau^2 estimator: ML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -275.6087  1206.8080   559.2175   570.0914   559.5913   
## 
## tau^2 (estimated amount of residual heterogeneity):     7.1971 (SE = 0.9729)
## tau (square root of estimated tau^2 value):             2.6827
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   23543254453.90
## 
## Test for Residual Heterogeneity:
## QE(df = 109) = 166085.2924, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 109) = 18.6144, p-val < .0001
## 
## Model Results:
## 
##                                    estimate      se    tval    pval    ci.lb 
## plant.speciesCaragana korshinskii    1.7927  0.3363  5.3303  <.0001   1.1261 
## plant.speciesCaragana microphylla    0.0130  1.0070  0.0129  0.9897  -1.9829 
## plant.speciesEphedra californica     2.6828  0.5122  5.2375  <.0001   1.6676 
##                                     ci.ub 
## plant.speciesCaragana korshinskii  2.4593  *** 
## plant.speciesCaragana microphylla  2.0089      
## plant.speciesEphedra californica   3.6980  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretations

1. Density matters in both grasslands and deserts
2. The net effect of increasing density is positive and it’s significantly positive (1.92)