library(dplyr) library(tidyr) library(lme4) library(lsmeans) library(multcomp)
sessionInfo()
This analysis will look at the revenue by watershed which accounts for the amount of crop area that has been converted to prairie. The assumed prairie cost per acre is
prairie_cost_per_acre = 95
This analysis will calculate the mean difference in revenue for both corn and soybean comparing the all crop (control) to the prairie treatment.
df <- 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(crop_prop = 1-prairie_pct/100, revenue = crop_prop * revenue - prairie_cost_per_acre * (1-crop_prop), any_prairie = ifelse(treatment == "0% prairie", "No", "Yes")) df %>% 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
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
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
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.