vicia_ndvi.R

# library("googledrive")
library("googlesheets4")
library("tidyverse")
library("ggpubr")

# drive_auth(use_oob = TRUE, cache = TRUE)
# 
mysheets <- googledrive::drive_find(type = "spreadsheet")
mysheets
# reuse token to Sheet Auth. 
# gs4_auth(token = drive_token())
# sheets_read(vicia)
# vicia <- Sheets$id[1]
# gs4_browse(vicia)

dat <- sheets_read("1n1dMV7CbXkdodh1t_QcbFC0L1w_e4ph70EbQthJYpSE") %>% 
  rename(protect_1 = Tra1, protect_2 = Tra2, lote_test=Tot) %>% 
  pivot_longer(-(Date:Days), names_to = "franja", values_to = "ndvi") %>%
  mutate(trat=factor(ifelse(str_detect(franja,"protect"),"Protegido", "Testigo")))

str(dat)

formula <- y ~ x

dat  %>% 
  ggplot(aes(Days, ndvi, col=trat)) +
  geom_point()+
  geom_smooth(se = FALSE, method = lm)+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),
    formula = formula, label.y.npc = "middle") +  
  cowplot::theme_minimal_grid()+
  labs(x = "Fecha desde aplicación de fungicida", 
       title = "Progreso de ndvi desde aplicación de tratamiento\nComparación de pendientes", 
       subtitle = "Con fungicida -0,002 vs Testigo -0,0047 (P=0,0042)", 
       col = "Tratamiento")

library(emmeans)
mod <- lm(ndvi ~  Days*trat, data=dat)
# Obtain slopes for each trat ...
( mod.emt <- emtrends(mod, "trat", var = "Days") )
pairs(mod.emt)