Cost of facilitation: invasive grasses limit recruitment of benefactor shrubs.

ecoblender

alex.filazzola

Abstract

Shrubs facilitate the abundance and productive of annual plants in desert ecosystems.However, these shrub microhabitats favour plant species with competitive life histories. Ephedra californica is a dominant shrub in the San Joaquin desert that has been identified as a facilitator. Here, we explore the factors that limit the Ephedra californica recruitment into the San Joaquin desert including substrate, water availability, and herbivory. We also explore the role of the invasive grass Bromus madritensis on limiting establishment of E. californica. Vegetation surveys were conducted in the field during 2013 and collected seed was then used to conduct two greenhouse trials. The first explore germination and establishment techniques for E. california and the optimal substrate. The second examined E. californica establishment in present of the invasive B. madritensis responding to different water levels and herbivory. These results can have implications for land managers in the San Joaquin Valley to maintain native shrub biodiversity in the region.

---

## load libraries
library(tidyverse)
# library(OIsurv)
library(emmeans)
library(ggthemes)
library(ape)

## load data
substrate <-read.csv("data/Ephedra.substrate.csv")
recruit <-read.csv("data/Ephedra.recruitment.csv")
landscape <-read.csv("data/ephedra.landscape.csv")

##load functions
## inverse hyperbolic sine transformation 
ihs <- function(x) {
    y <- log(x + sqrt(x ^ 2 + 1))
    return(y)
}
se <- function(x, ...) sd(x, na.rm=T)/sqrt(length(na.omit(x)))
source("functions.r")

Ephedra characteristics at the landscape

avg.size <- landscape %>% summarize(h=mean(H), d1=mean(D1), area=mean(Area), density=mean(Shrub.density), rdm=mean(RDM.2013),h.se=se(H), d1.se=se(D1), area.se=se(Area), density.se=se(Shrub.density), rdm.se=se(RDM.2013))
avg.size

##          h       d1     area     density      rdm       h.se      d1.se
## 1 1.332925 3.446412 9.321902 0.004318859 7.042223 0.01506635 0.05137079
##     area.se   density.se    rdm.se
## 1 0.2708882 6.272347e-05 0.2109119

Test of limiting factors - brome

recruit[is.na(recruit)] <- 0

## detailed names for treatments
recruit[,"Lvl"] <- as.character(recruit[,"Lvl"]) ## convert variables to characters
treats <- data.frame(Lvl=c("none","low","high","med","once","twice"), treat=c("control","low water","high water","partial shade","clipped once","clipped twice"), stringsAsFactors =F)
recruit <- merge(recruit, treats, by="Lvl")
recruit[recruit$Treatment=="shade" & recruit$Lvl == "high","treat"] <- "full shade"
recruit[,"treat"] <- as.factor(recruit[,"treat"]) ## convert back to factor

## Ephedra establishment & survival
recruit[,"emergence"] <- ifelse(recruit$ephedra.emergence>0, 1,0)
recruit[,"survival"] <- ifelse(recruit$ephedra.end==recruit$emergence,1,0)
recruit[recruit$emergence==0,"survival"] <- 0
names(recruit)[16] <- "biomass" ## simplify biomass name for Ephedra

## water
water <- subset(recruit, Treatment == "control" | Treatment=="water")
water.long <- water %>% gather(measure, value, c(19,20)) ## collapse different ephedra measures

## Plot all metrics as an average of brome density
water.germ <- water.long %>% group_by(Density, treat, measure) %>% summarize(value=mean(value, na.rm=T))

plot1 <-  ggplot(water.germ, aes(x=Density, y=value,color=treat)) +ylab("E. californica response")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + geom_smooth( method="lm", formula= y~x, se=FALSE, lwd=2) + facet_grid(measure~"water", scales= "free_y") + scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))

## shade
shade <- subset(recruit, Treatment == "control" | Treatment=="shade")
shade.long <- shade %>% gather(measure, value, c(19,20)) ## collapse different ephedra measures

plot2 <- ggplot(shade.long, aes(x=Density, y=value, color=treat)) +ylab("E. californica response")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + geom_smooth( method="lm", formula= y~x, se=FALSE, lwd=2) + facet_grid(measure~"shade", scales= "free_y") + scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))


## clipping
clip <- subset(recruit, Treatment == "control" | Treatment=="clipped")
clip.long <- clip %>% gather(measure, value, c(19,20)) ## collapse different ephedra measures

plot3 <- ggplot(clip.long , aes(x=Density, y=value, color=treat)) +ylab("E. californica response")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + geom_smooth( method="lm", formula= y~x, se=T, lwd=2) + facet_grid(measure~"clipping", scales= "free_y") + scale_color_manual(values=c("#E69F00", "#56B4E9","#999999"))

# library(gridExtra)
# grid.arrange(plot1, plot2, plot3, ncol=3)
plot1

plot2

plot3

## models for survival and emergence
library(emmeans)
recruit$treat <- factor(recruit$treat, levels(recruit$treat)[c(1:2,7:3)]) ## reorder control to last treatment

## Compare Emergence
m1 <- glm(ephedra.begin ~ treat * brome.begin+ Density, data=recruit, family="binomial")
car::Anova(m1, type="3")

## Analysis of Deviance Table (Type III tests)
## 
## Response: ephedra.begin
##                   LR Chisq Df Pr(>Chisq)  
## treat              16.0833  6    0.01331 *
## brome.begin         2.0576  1    0.15145  
## Density             3.3518  1    0.06713 .
## treat:brome.begin   9.1364  6    0.16605  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m1, specs = trt.vs.ctrlk ~treat, var="survival")

## $emmeans
##  treat          emmean    SE  df asymp.LCL asymp.UCL
##  clipped once  -0.0993 0.202 Inf   -0.4945     0.296
##  clipped twice -0.2949 0.204 Inf   -0.6948     0.105
##  partial shade -1.1095 0.263 Inf   -1.6248    -0.594
##  low water     -0.1017 0.201 Inf   -0.4962     0.293
##  high water     0.5016 0.209 Inf    0.0913     0.912
##  full shade    -0.5494 0.210 Inf   -0.9617    -0.137
##  control        0.1927 0.202 Inf   -0.2028     0.588
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df z.ratio p.value
##  clipped once - control    -0.292 0.285 Inf -1.023  0.7613 
##  clipped twice - control   -0.488 0.287 Inf -1.701  0.3408 
##  partial shade - control   -1.302 0.332 Inf -3.925  0.0005 
##  low water - control       -0.294 0.285 Inf -1.034  0.7553 
##  high water - control       0.309 0.291 Inf  1.061  0.7394 
##  full shade - control      -0.742 0.291 Inf -2.550  0.0544 
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: dunnettx method for 6 tests

out1 <- data.frame(emmeans(m1, specs = trt.vs.ctrlk ~treat, var="ephedra.begin")$contrast)

pR2 = 1 - m1$deviance / m1$null.deviance

## Compare Survival
m2 <- glm(survival ~ brome.begin * treat + Density, data=recruit, family="binomial")
car::Anova(m2, type="3")

## Analysis of Deviance Table (Type III tests)
## 
## Response: survival
##                   LR Chisq Df Pr(>Chisq)   
## brome.begin         2.5040  1   0.113559   
## treat              17.6628  6   0.007133 **
## Density             6.2313  1   0.012551 * 
## brome.begin:treat   3.3961  6   0.757743   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pR2 = 1 - m2$deviance / m2$null.deviance

emmeans(m2, specs = trt.vs.ctrlk ~treat, var="survival")

## $emmeans
##  treat         emmean    SE  df asymp.LCL asymp.UCL
##  clipped once  -0.933 0.225 Inf     -1.37    -0.493
##  clipped twice -0.901 0.225 Inf     -1.34    -0.461
##  partial shade -2.277 0.382 Inf     -3.03    -1.528
##  low water     -1.044 0.230 Inf     -1.49    -0.593
##  high water    -1.015 0.234 Inf     -1.47    -0.556
##  full shade    -2.279 0.340 Inf     -2.95    -1.612
##  control       -0.829 0.218 Inf     -1.26    -0.402
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df z.ratio p.value
##  clipped once - control   -0.1042 0.313 Inf -0.333  0.9893 
##  clipped twice - control  -0.0725 0.312 Inf -0.233  0.9963 
##  partial shade - control  -1.4478 0.440 Inf -3.289  0.0056 
##  low water - control      -0.2147 0.316 Inf -0.680  0.9178 
##  high water - control     -0.1859 0.320 Inf -0.581  0.9470 
##  full shade - control     -1.4498 0.403 Inf -3.601  0.0018 
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: dunnettx method for 6 tests

out2 <- data.frame(emmeans(m2, specs = trt.vs.ctrlk ~treat, var="survival")$contrast)

## Compare growth
m3 <- lm(sqrt(Ephedra.above) ~ Density + treat * brome.begin, data= subset(recruit, Ephedra.above >0))
car::Anova(m3, type="3")

## Anova Table (Type III tests)
## 
## Response: sqrt(Ephedra.above)
##                    Sum Sq  Df  F value    Pr(>F)    
## (Intercept)       0.34992   1 140.1699 < 2.2e-16 ***
## Density           0.00010   1   0.0411  0.839610    
## treat             0.04588   6   3.0633  0.006916 ** 
## brome.begin       0.00029   1   0.1171  0.732583    
## treat:brome.begin 0.01540   6   1.0284  0.408043    
## Residuals         0.48430 194                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

out3 <- data.frame(emmeans(m3, specs = trt.vs.ctrlk ~treat, var="Ephedra.above")$contrast)


m4 <- lm(sqrt(Ephedra.below) ~ Density + treat*brome.begin, data= subset(recruit, Ephedra.below >0))
car::Anova(m4, type="3")

## Anova Table (Type III tests)
## 
## Response: sqrt(Ephedra.below)
##                     Sum Sq  Df  F value Pr(>F)    
## (Intercept)       0.067103   1 124.2331 <2e-16 ***
## Density           0.000001   1   0.0015 0.9693    
## treat             0.007784   6   2.4020 0.0293 *  
## brome.begin       0.000096   1   0.1769 0.6745    
## treat:brome.begin 0.003843   6   1.1858 0.3156    
## Residuals         0.101006 187                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

out4 <- data.frame(emmeans(m4, specs = trt.vs.ctrlk ~treat, var="Ephedra.below")$contrast)

## effect of treatment on brome
m5 <- lm(sqrt(brome.biomass) ~ treat, data=subset(recruit, brome.biomass >0 ))
anova(m5)

## Analysis of Variance Table
## 
## Response: sqrt(brome.biomass)
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## treat       6 25.304  4.2173  66.969 < 2.2e-16 ***
## Residuals 530 33.376  0.0630                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

out5 <- data.frame(emmeans(m5, specs = trt.vs.ctrlk ~treat, var="brome.biomass")$contrast)


m6 <- glm(brome.prop.begin/100 ~ treat, data=recruit, family="binomial")
anova(m6, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: brome.prop.begin/100
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                    699     460.48         
## treat  6   6.6555       693     453.82   0.3539

emmeans(m6, pairwise~treat, var="brome.biomass")

## $emmeans
##  treat          emmean    SE  df asymp.LCL asymp.UCL
##  clipped once  -0.0921 0.200 Inf   -0.4845     0.300
##  clipped twice  0.1342 0.200 Inf   -0.2587     0.527
##  partial shade -0.2554 0.202 Inf   -0.6506     0.140
##  low water      0.1926 0.201 Inf   -0.2012     0.586
##  high water    -0.0861 0.200 Inf   -0.4784     0.306
##  full shade     0.3743 0.204 Inf   -0.0246     0.773
##  control        0.1181 0.200 Inf   -0.2745     0.511
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                      estimate    SE  df z.ratio p.value
##  clipped once - clipped twice  -0.22627 0.283 Inf -0.799  0.9852 
##  clipped once - partial shade   0.16331 0.284 Inf  0.575  0.9975 
##  clipped once - low water      -0.28466 0.284 Inf -1.004  0.9534 
##  clipped once - high water     -0.00601 0.283 Inf -0.021  1.0000 
##  clipped once - full shade     -0.46637 0.285 Inf -1.634  0.6604 
##  clipped once - control        -0.21020 0.283 Inf -0.742  0.9899 
##  clipped twice - partial shade  0.38958 0.284 Inf  1.370  0.8178 
##  clipped twice - low water     -0.05839 0.284 Inf -0.206  1.0000 
##  clipped twice - high water     0.22025 0.283 Inf  0.777  0.9871 
##  clipped twice - full shade    -0.24011 0.286 Inf -0.841  0.9807 
##  clipped twice - control        0.01606 0.283 Inf  0.057  1.0000 
##  partial shade - low water     -0.44797 0.285 Inf -1.574  0.6993 
##  partial shade - high water    -0.16933 0.284 Inf -0.596  0.9969 
##  partial shade - full shade    -0.62969 0.286 Inf -2.198  0.2964 
##  partial shade - control       -0.37352 0.284 Inf -1.314  0.8457 
##  low water - high water         0.27865 0.284 Inf  0.982  0.9579 
##  low water - full shade        -0.18172 0.286 Inf -0.635  0.9957 
##  low water - control            0.07446 0.284 Inf  0.262  1.0000 
##  high water - full shade       -0.46036 0.285 Inf -1.613  0.6741 
##  high water - control          -0.20419 0.283 Inf -0.721  0.9914 
##  full shade - control           0.25617 0.286 Inf  0.897  0.9731 
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 7 estimates

outall <- bind_rows(out1,out2,out3,out4, out5)
outall[,"comparison"] <- c(rep("emergence", nrow(out1)), rep("survival", nrow(out2)), rep("ABG", nrow(out3)),rep("BWG", nrow(out4)), rep("bromeBio", nrow(out5)))


outall[,"treatment"] <- sub("\\ -.*", "", outall$contrast) ## simplify treatment names
outall[,"treatment"] <- factor(outall[,"treatment"], levels=c("clipped once","clipped twice", "low water", "high water", "partial shade", "full shade"))

## Plot EMMEANS for emergence and survival
ggplot(outall %>% filter(comparison %in% c("emergence", "survival")), aes(x=treatment, y=estimate)) + geom_point(size=3) + geom_errorbar(aes(ymin=estimate-SE*1.96, ymax=estimate+SE*1.96), width=0) + facet_grid(comparison~.) + theme_Publication() + ylab("Estimated marginal means") + geom_hline(yintercept = 0, lwd=1, lty=2, color="#00000070")

## Plot EMMEANS for ABG and BWG biomass with brome
outall[,"biomass"]<- ifelse(outall[,"comparison"] == "ABG", "E. californica - AGB", 
                               ifelse(outall[,"comparison"] == "BWG", "E. californica - BGB", "B. madritensis"))
ggplot(outall %>% filter(comparison %in% c("ABG", "BWG","bromeBio")), aes(x=treatment, y=estimate)) + geom_point(size=3)  + geom_errorbar(aes(ymin=estimate-SE*1.96, ymax=estimate+SE*1.96), width=0) + facet_grid(biomass~.,  scales= "free_y") + theme_Publication() + ylab("Estimated marginal means")+ geom_hline(yintercept = 0, lwd=1, lty=2, color="#00000070")

#### plot densities of brome

## calculate means for points on figure
avgRates <- recruit %>% group_by(Density) %>% summarize(emerge=mean(ephedra.begin), emerge.se=se(ephedra.begin), surv=mean(survival), surv.se=se(survival))


m1 <- glm(ephedra.begin ~ Density* treat, data=recruit, family="binomial")
car::Anova(m1, type="2")

## Analysis of Deviance Table (Type II tests)
## 
## Response: ephedra.begin
##               LR Chisq Df Pr(>Chisq)    
## Density         3.7238  1    0.05364 .  
## treat          30.5072  6  3.147e-05 ***
## Density:treat   6.4480  6    0.37492    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m1, specs = trt.vs.ctrlk ~treat, var="survival")

## $emmeans
##  treat         emmean    SE  df asymp.LCL asymp.UCL
##  clipped once  -0.120 0.201 Inf   -0.5135    0.2727
##  clipped twice -0.286 0.203 Inf   -0.6849    0.1125
##  partial shade -1.061 0.249 Inf   -1.5487   -0.5740
##  low water     -0.080 0.200 Inf   -0.4723    0.3123
##  high water     0.493 0.207 Inf    0.0876    0.8992
##  full shade    -0.490 0.206 Inf   -0.8936   -0.0858
##  control        0.201 0.201 Inf   -0.1933    0.5951
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df z.ratio p.value
##  clipped once - control    -0.321 0.284 Inf -1.131  0.6973 
##  clipped twice - control   -0.487 0.286 Inf -1.703  0.3396 
##  partial shade - control   -1.262 0.320 Inf -3.947  0.0005 
##  low water - control       -0.281 0.284 Inf -0.990  0.7798 
##  high water - control       0.292 0.289 Inf  1.013  0.7669 
##  full shade - control      -0.691 0.288 Inf -2.398  0.0803 
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: dunnettx method for 6 tests

out1 <- data.frame(emmeans(m1, specs = trt.vs.ctrlk ~treat, var="ephedra.begin")$contrast)

m1 <- glm(ephedra.begin ~ Density, data=recruit, family="binomial")

ggplot() + geom_smooth( data = recruit, aes(x=Density, y= survival), method="glm", method.args = list(family = "binomial"), col="black", lwd=2, fullrange=T) + geom_smooth( data = recruit, aes(x=Density, y= ephedra.begin), method="glm", method.args = list(family = "binomial"), col="black", lwd=2, fullrange=T) + ylab(expression("Proportion of "*paste(italic("E. californica"))*" individuals")) + theme_Publication() + 
  annotate("text", x = 15, y = 0.26, label = "survival probability", size=5) +
  annotate("text", x = 15, y = 0.5, label = "emergence probability", size=5) + 
  xlab(expression(paste(italic("B. madritensis "))*"density")) + xlim(-1,21) +
  geom_point(data=avgRates, aes(x=Density, y=emerge), size=2) + geom_errorbar(data=avgRates, aes(x=Density, ymin=emerge-emerge.se, ymax=emerge+emerge.se), width=0)+
  geom_point(data=avgRates, aes(x=Density, y=surv), size=2) + geom_errorbar(data=avgRates, aes(x=Density, ymin=surv-surv.se, ymax=surv+surv.se), width=0)

Patterns of Brome Emergence

## Relationship of seeds to emerged individuals
plot1 <- ggplot(recruit %>% filter(Density != 0), aes(x=Density, y= brome.begin)) + geom_jitter(height=0, width=0.5) + xlim(1,21) + ylim(0,20) + geom_smooth(method="glm", method.args = list(family = "quasipoisson"), fullrange=T, color="Grey50") + theme_Publication() + xlab(expression("Number of "*paste(italic("B. madritensis"))*" seeds")) + ylab("Number of emerged individuals")

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

plot2 <- ggplot(recruit %>% filter(Density != 0), aes(as.factor(Density), y=brome.prop.begin))  + geom_boxplot(fill="Grey70") + theme_Publication() + xlab(expression("Number of "*paste(italic("B. madritensis"))*" seeds")) + ylab(expression("Proportion of emerged"*paste(italic("B. madritensis"))*" individuals")) 

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

gridExtra::grid.arrange(plot1, plot2, ncol=2)

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Relationship of seeds to biomass
ggplot(recruit %>% filter(Density != 0), aes(x=Density, y= log(brome.biomass))) + geom_jitter(height=0.5, width=0.5) + xlim(1,21)  + geom_smooth(method="glm", method.args = list(family = "quasipoisson"), fullrange=T) + theme_Publication() + xlab(expression("Number of "*paste(italic("B. madritensis"))*" seeds")) + ylab("Number of emerged individuals")

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

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

## Warning: Computation failed in `stat_smooth()`:
## negative values not allowed for the 'quasiPoisson' family

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

m1 <- glm(brome.begin ~ Density * treat, data=subset(recruit, Density !=0) , family="quasipoisson")
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: quasipoisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                            559    1638.56              
## Density        1  1075.25       558     563.31 < 2.2e-16 ***
## treat          6    20.93       552     542.38 0.0003357 ***
## Density:treat  6     0.67       546     541.71 0.9920988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Average  emergence among treatments
emergePatterns <- recruit %>% group_by(treat) %>% summarize(propEmerge = mean(ephedra.begin),Emerge.se = se(ephedra.begin), propSurv = mean(survival), Surv.se=se(survival), bromeEmerge=mean(brome.prop.begin/100), brome.se=se(brome.prop.begin/100))
write.csv(emergePatterns, "emergePatterns.csv", row.names=F)

## Average emergence across everything
recruit %>% summarize(propEmerge = mean(ephedra.begin),Emerge.se = se(emergence), propSurv = mean(survival), Surv.se=se(survival), bromeEmerge=mean(brome.prop.begin/100), brome.se=se(brome.prop.begin/100))

##   propEmerge  Emerge.se  propSurv    Surv.se bromeEmerge   brome.se
## 1  0.4585714 0.01650094 0.2328571 0.01598617   0.5136429 0.01325073

Shrub-plant facilitation results

## load plant comp data
plant <- read.csv("data//ephedra.plantcomp.csv")

library(lmerTest)

## community measures
plant[,"abd"] <- apply(plant[,9:48], 1, function(x) sum(x, na.rm=T)) ## abundance 
plant[,"rich"] <- apply(plant[,9:48], 1, function(x) sum(x>0, na.rm=T))

## percent brome
plant[is.na(plant)] <- 0
mean(plant[,"B.rubens"]/plant[,"abd"])*100

## [1] 69.22444

## fit a linear model
m1 <- glmer.nb(abd ~ microsite + (1|year) + (1|site), data=plant)
summary(m1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(9.0118)  ( log )
## Formula: abd ~ microsite + (1 | year) + (1 | site)
##    Data: plant
## 
##      AIC      BIC   logLik deviance df.resid 
##   4883.1   4903.7  -2436.6   4873.1      445 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4761 -0.6159 -0.0603  0.5463  3.2187 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  year   (Intercept) 0.003222 0.05677 
##  site   (Intercept) 0.078784 0.28069 
## Number of obs: 450, groups:  year, 3; site, 3
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     4.91598    0.16702  29.433   <2e-16 ***
## micrositeshrub  0.31954    0.03261   9.798   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## microstshrb -0.098

car::Anova(m1, test="Chisq")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: abd
##            Chisq Df Pr(>Chisq)    
## microsite 95.992  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2 <- glmer.nb(rich ~ microsite + (1|year) + (1|site), data=plant)
summary(m2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(329887.7)  ( log )
## Formula: rich ~ microsite + (1 | year) + (1 | site)
##    Data: plant
## 
##      AIC      BIC   logLik deviance df.resid 
##   1399.3   1419.9   -694.7   1389.3      445 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.10029 -0.48707  0.01885  0.34541  1.89031 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  year   (Intercept) 0.0008052 0.02838 
##  site   (Intercept) 0.0374271 0.19346 
## Number of obs: 450, groups:  year, 3; site, 3
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.19133    0.11531  10.331   <2e-16 ***
## micrositeshrub -0.57803    0.05999  -9.636   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## microstshrb -0.177

car::Anova(m2, test="Chisq")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: rich
##            Chisq Df Pr(>Chisq)    
## microsite 92.855  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### calculate averages and variance
plantavg <- plant %>% group_by(microsite) %>% summarize(ABD=mean(abd), RICH=mean(rich),ABD.se=se(abd)*1.96,  RICH.se=se(rich)*1.96 )

##abundance
abd <- plantavg[,c(1:2,4)]
abd[,2:3] <- abd[,2:3]/0.25

## plot
plot1 <- ggplot(abd, aes(x=microsite, y=ABD, fill=microsite)) + geom_bar(stat="identity", position=position_dodge(width=0.9), color="black") + ylab(expression("plant density (m"^2*")")) + theme_Publication() + geom_errorbar(aes(x=microsite, ymin=ABD-ABD.se, ymax=ABD+ABD.se), position=position_dodge(width=0.9), width=0) + xlab("") + scale_fill_brewer()

## richness
rich <- plantavg[,c(1,3,5)]


## plot
plot2 <- ggplot(rich, aes(x=microsite, y=RICH, fill=microsite)) + geom_bar(stat="identity", position=position_dodge(width=0.9), color="black") + ylab("species richness") + theme_Publication() + geom_errorbar(aes(x=microsite, ymin=RICH-RICH.se, ymax=RICH+RICH.se), position=position_dodge(width=0.9), width=0) + xlab("") + scale_fill_brewer()

require(gridExtra)

grid.arrange(plot1, plot2, nrow=2)

### compare just brome
plant[,"perBrome"] <- plant[,"B.rubens"]/plant[,"abd"]*100

m1 <- glmer(B.rubens/abd ~ microsite + (1|site) + (1|year), data=plant, family="binomial")
summary(m1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: B.rubens/abd ~ microsite + (1 | site) + (1 | year)
##    Data: plant
## 
##      AIC      BIC   logLik deviance df.resid 
##    305.3    321.8   -148.7    297.3      446 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3433 -0.1943  0.0857  0.2643  1.2188 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 2.6488   1.6275  
##  year   (Intercept) 0.6417   0.8011  
## Number of obs: 450, groups:  site, 3; year, 3
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.2411     1.0640  -0.227    0.821    
## micrositeshrub   4.2433     0.4690   9.048   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## microstshrb -0.099

car::Anova(m1, test="Chisq")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: B.rubens/abd
##            Chisq Df Pr(>Chisq)    
## microsite 81.858  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2 <- glmer.nb(B.rubens ~ microsite + (1|site) + (1|year), data=plant)
summary(m2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(3.1316)  ( log )
## Formula: B.rubens ~ microsite + (1 | site) + (1 | year)
##    Data: plant
## 
##      AIC      BIC   logLik deviance df.resid 
##   4921.8   4942.3  -2455.9   4911.8      445 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6697 -0.6630 -0.1346  0.5268  2.8941 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.26562  0.5154  
##  year   (Intercept) 0.03112  0.1764  
## Number of obs: 450, groups:  site, 3; year, 3
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     4.13618    0.31718   13.04   <2e-16 ***
## micrositeshrub  0.99767    0.05745   17.37   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## microstshrb -0.093

car::Anova(m2, test="Chisq")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: B.rubens
##            Chisq Df Pr(>Chisq)    
## microsite 301.57  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## simplify dataframe to just brome and not-brome
bromeOther <- plant[,c("year","microsite","abd","B.rubens")]
bromeOther$B.rubens <- bromeOther$B.rubens/0.25 ## convert to individuals per m2
bromeOther <- bromeOther %>% gather(plants, density, 3:4)

avgBrome <- bromeOther %>% group_by(year, microsite, plants) %>% summarize(dens=mean(density), error=se(density))


ggplot(avgBrome %>% filter(plants != "abd"), aes(x=microsite, y=dens, fill=microsite)) + geom_bar(stat="identity", position=position_dodge(width=0.9), color="black") +  ylab(expression(paste(italic("B. madritensis"))*"density (m"^2*")")) + theme_Publication() + geom_errorbar(aes(x=microsite, ymin=dens-error, ymax=dens+error), position=position_dodge(width=0.9), width=0) + xlab("") + scale_fill_brewer() + facet_grid(~year)

Adjusted null hypothesis based on germ rate

library(car)

### Determine emergence and survival rate of Ephera pre-trial
substrate <-read.csv("data/Ephedra.substrate.csv")

## determine survival rate
survRate <- substrate %>% filter(date == "2/25/2015") %>%  ## last survey instance
   group_by(Micro, Sand) %>% summarize(propSurv = mean(survival))

## determine emergence rate
emergeRate <- substrate %>% group_by(Micro, Sand, Rep) %>% summarize(max.emerge = max(emergence, na.rm=T)) %>%  ## number of emerged individuals
  mutate(emerge.occ = ifelse(max.emerge>0, 1, 0)) %>%  ## yes/no emergence
  group_by(Micro, Sand) %>% summarize(prop.emerge = mean(emerge.occ))


### Specify new Null

## Compare Emergence
m1 <- glm(ephedra.begin ~ Density * treat, data=recruit, family="binomial")
car::Anova(m1, type="2")

## Analysis of Deviance Table (Type II tests)
## 
## Response: ephedra.begin
##               LR Chisq Df Pr(>Chisq)    
## Density         3.7238  1    0.05364 .  
## treat          30.5072  6  3.147e-05 ***
## Density:treat   6.4480  6    0.37492    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m1, specs = trt.vs.ctrlk ~treat, var="survival")

## $emmeans
##  treat         emmean    SE  df asymp.LCL asymp.UCL
##  clipped once  -0.120 0.201 Inf   -0.5135    0.2727
##  clipped twice -0.286 0.203 Inf   -0.6849    0.1125
##  partial shade -1.061 0.249 Inf   -1.5487   -0.5740
##  low water     -0.080 0.200 Inf   -0.4723    0.3123
##  high water     0.493 0.207 Inf    0.0876    0.8992
##  full shade    -0.490 0.206 Inf   -0.8936   -0.0858
##  control        0.201 0.201 Inf   -0.1933    0.5951
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df z.ratio p.value
##  clipped once - control    -0.321 0.284 Inf -1.131  0.6973 
##  clipped twice - control   -0.487 0.286 Inf -1.703  0.3396 
##  partial shade - control   -1.262 0.320 Inf -3.947  0.0005 
##  low water - control       -0.281 0.284 Inf -0.990  0.7798 
##  high water - control       0.292 0.289 Inf  1.013  0.7669 
##  full shade - control      -0.691 0.288 Inf -2.398  0.0803 
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: dunnettx method for 6 tests

out1 <- data.frame(emmeans(m1, specs = trt.vs.ctrlk ~treat, var="ephedra.begin")$contrast)


m1.lht <- lht(m1, c("(Intercept) = 0.68"))
m1.lht

## Linear hypothesis test
## 
## Hypothesis:
## (Intercept) = 0.68
## 
## Model 1: restricted model
## Model 2: ephedra.begin ~ Density * treat
## 
##   Res.Df Df  Chisq Pr(>Chisq)  
## 1    687                       
## 2    686  1 6.0349    0.01403 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1