library(dplyr) library(tidyr) library(lme4) library(lsmeans) library(multcomp)
sessionInfo()
This vignette compares crop revenue in the crop area of each watershed, i.e. it does not account for the fact that part of the watershed is now prairie in the treatment watersheds. Thus, the aim is to investigate the effect of prairie on the crop.
This analysis will calculate the mean difference in revenue for both corn and soybean comparing the all crop (control) to the prairie treatment.
STRIPSTyndall::revenue %>% left_join(STRIPSMeta::watersheds) %>% left_join(STRIPSMeta::crop_seed_info) %>% filter(year > 2007) %>% # for PNAS paper we are not including 2007 (establishment) mutate(any_prairie = ifelse(treatment == "0% prairie", "No", "Yes")) %>% group_by(block, any_prairie, crop_species) %>% summarize(mean = mean(revenue)) %>% spread(any_prairie, mean) %>% mutate(Diff = No - Yes) %>% group_by(crop_species) %>% summarize(`No prairie - Any prairie` = mean(Diff, na.rm=TRUE))
Corn analysis
df <- STRIPSTyndall::revenue %>% left_join(STRIPSMeta::watersheds) %>% left_join(STRIPSMeta::crop_seed_info) corn <- df %>% filter(crop_species == "corn") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + (1|watershed) + (1|year), data=corn) # summary(m) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Soybean analysis
soybean <- df %>% filter(crop_species == "soybean") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + (1|watershed) + (1|year), data=soybean) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Corn analysis
df <- STRIPSTyndall::revenue %>% left_join(STRIPSMeta::watersheds) %>% left_join(STRIPSMeta::crop_seed_info) corn <- df %>% filter(crop_species == "corn") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + size_ha + (1|watershed) + (1|year), data=corn) # summary(m) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Soybean analysis
soybean <- df %>% filter(crop_species == "soybean") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + size_ha + (1|watershed) + (1|year), data=soybean) summary(m) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Corn analysis
df <- STRIPSTyndall::revenue %>% left_join(STRIPSMeta::watersheds) %>% left_join(STRIPSMeta::crop_seed_info) corn <- df %>% filter(crop_species == "corn") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + log(size_ha) + (1|watershed) + (1|year), data=corn) # summary(m) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Soybean analysis
soybean <- df %>% filter(crop_species == "soybean") %>% dplyr::select(revenue, treatment, block, size_ha, watershed, year) m <- lmer(revenue ~ treatment + block + log(size_ha) + (1|watershed) + (1|year), data=soybean) summary(m) m %>% lsmeans("treatment") %>% contrast("trt.vs.ctrl", ref=2:4) %>% confint
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.