In ISU-STRIPS/STRIPSTyndall: STRIPS Tyndall Economic data

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.

Simple analysis

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))

Using contrasts

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

Using contrasts including watershed size

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

Using contrasts including logarithm of watershed size

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