Get first half of synthesis vignette working with full data

Author: KristinaRiemer

vignettes/04-synthesis-data.Rmd

@@ -13,7 +13,7 @@ The first contains all the canopy cover values for 2018, which was created in th
 The second is the cumulative growing degree days for all of 2018, which were calculated from the daily minimum and maximum temperatures in the weather vignette. 
 They are combined by their common column, the date. 
 
-```{r synth_setup}
+```{r synth_setup, message=FALSE}
 library(dplyr)
 library(ggplot2)
 library(jsonlite)
@@ -28,14 +28,14 @@ options(betydb_url = "https://terraref.ncsa.illinois.edu/bety/",
 ```{r get_trait_data, message = FALSE}
 trait_canopy_cover <- betydb_query(table     = "search", 
                       trait     = "canopy_cover", 
-                      date      = "~2017",
+                      date      = "~2018",
                       limit     =  "none")
 trait_canopy_cover_day = trait_canopy_cover %>% 
   mutate(day = as.Date(raw_date))
 ```
 
 ```{r get_weather_data}
-weather <- fromJSON('https://terraref.ncsa.illinois.edu/clowder/api/geostreams/datapoints?stream_id=46431&since=2017-01-01&until=2017-12-31', flatten = FALSE)
+weather <- fromJSON('https://terraref.ncsa.illinois.edu/clowder/api/geostreams/datapoints?stream_id=46431&since=2018-01-01&until=2018-12-31', flatten = FALSE)
 weather <- weather$properties %>% 
   mutate(time = ymd_hms(weather$end_time))
 daily_values = weather %>% 
@@ -74,24 +74,26 @@ We provide estimated values for the asymptote $c$ and initial canopy cover value
 The below provides better estimates for the $c$, $a$, and $b$ parameters, which are used to plot the model as an orange line on top of the black points which are actual values. 
 
 ```{r model_get_parameters}
-# TODO: replace with actual data
-# single_cultivar <- trait_weather_df %>% 
-#   filter(cultivar == "PI656026")
-single_cultivar <- data.frame(gdd_cum = c(5, 6.2, 7, 18, 25, 30, 36, 39, 40), mean = c(10, 10.11, 10.5, 20, 32.63432, 40, 50, 50.5, 50.5))
+single_cultivar <- trait_weather_df %>% 
+  filter(cultivar == "PI656026")
 
-c <- 51
-a <- 10
+c <- 90
+a <- 0.1
 y <- single_cultivar$mean[3]
 g <- single_cultivar$gdd_cum[3]
 b <- ((log((c/y) - 1)) - a)/g
 model_single_cultivar <- nls(mean ~ c / (1 + exp(a + b * gdd_cum)), 
                              start = list(c = c, a = a, b = b),
-                             data = single_cultivar, trace = TRUE)
+                             data = single_cultivar)
 summary(model_single_cultivar)
 coef(model_single_cultivar)
 
+single_c <- coef(model_single_cultivar)[1]
+single_a <- coef(model_single_cultivar)[2]
+single_b <- coef(model_single_cultivar)[3]
+
 single_cultivar <- single_cultivar %>% 
-  mutate(mean_predict = coef(model_single_cultivar)[1] / (1 + exp(coef(model_single_cultivar)[2] + coef(model_single_cultivar)[3] * gdd_cum)))
+  mutate(mean_predict = single_c / (1 + exp(single_a + single_b * gdd_cum)))
 ggplot(single_cultivar) +
   geom_point(aes(x = gdd_cum, y = mean)) +
   geom_line(aes(x = gdd_cum, y = mean_predict), color = "orange") +
@@ -103,9 +105,8 @@ We then calculate the inflection point for this cultivar's model.
 The maximum growth rate is the change in canopy cover per day at the rate of maximum growth. The growing degree day at which maximum growth is obtained is called the _inflection point_. This occurs near the midpoint of the y-axis, or $\frac{c - a}{2}$.
 
 ```{r}
-inf_y <- (as.numeric(coef(model_single_cultivar)[1]) - as.numeric(coef(model_single_cultivar)[2])) / 2
-
-inf_x <- ((log((as.numeric(coef(model_single_cultivar)[1]) / inf_y) - 1)) - as.numeric(coef(model_single_cultivar)[2])) / as.numeric(coef(model_single_cultivar)[3])
+inf_y <- (as.numeric(single_c) - as.numeric(single_a)) / 2
+inf_x <- ((log((as.numeric(single_c) / inf_y) - 1)) - as.numeric(single_a)) / as.numeric(single_b)
 
 ggplot(single_cultivar) +
   geom_point(aes(x = gdd_cum, y = mean)) +
@@ -119,26 +120,33 @@ We then use the parameters from a single cultivar to run a model for each of the
 These results are used to plot the model predictions, which are shown as an orange line. 
 We also calculated the inflection point from each cultivar's model, which will be used in the following section. 
 
-```{r model_all_cultivars, eval=FALSE}
+```{r model_all_cultivars}
 all_cultivars <- c(day = as.double(), cultivar = as.character(), mean = as.numeric(), 
-                   gdd_cum = as.numeric(), mean_predict = as.numeric())
+                   gdd_cum = as.numeric(), mean_predict = as.numeric(), 
+                   inf_y = as.numeric(), inf_x = as.numeric())
+
 for(each_cultivar in unique(trait_weather_df$cultivar)){
   each_cultivar_df <- filter(trait_weather_df, cultivar == each_cultivar)
-  each_cultivar_model <- nls(y ~ c / (1 + exp(a + b * g)), 
-                            start = list(c = c, a = a, b = b),
-                            data = each_cultivar_df)
+  each_cultivar_model <- nls(mean ~ c / (1 + exp(a + b * gdd_cum)), 
+                             start = list(c = c, a = a, b = b), 
+                             data = each_cultivar_df)
   model_c <- coef(each_cultivar_model)[1]
   model_a <- coef(each_cultivar_model)[2]
   model_b <- coef(each_cultivar_model)[3]
   each_cultivar_df <- each_cultivar_df %>% 
-      mutate(mean_predict = model_c / (1 + exp(model_a + model_b * g)), 
-             inf_point = ((log((model_c / 100) - 1)) - model_a) / model_b)
+    mutate(mean_predict = model_c / (1 + exp(model_a + model_b * gdd_cum)), 
+           inf_y = (as.numeric(model_c) - as.numeric(model_a)) / 2, 
+           inf_x = ((log((as.numeric(model_c) / inf_y) - 1)) - 
+                      as.numeric(single_a)) / as.numeric(single_b))
   all_cultivars <- rbind(each_cultivar_df, all_cultivars)
 }
+
 ggplot(all_cultivars) +
   geom_point(aes(x = gdd_cum, y = mean)) +
   geom_line(aes(x = gdd_cum, y = mean_predict), color = "orange") +
   facet_wrap(~cultivar, scales = "free_y") +
+  geom_hline(yintercept = inf_y, linetype = "dashed") +
+  geom_vline(xintercept = inf_x) +
   labs(x = "Cumulative growing degree days", y = "Canopy Height")
 ```
 
@@ -148,15 +156,9 @@ The last thing that we are going to do is assess the difference in this relation
 We are going to use the inflection point from the logistic growth model, which indicates when canopy cover stops increasing as quickly with increasingly more warm days. 
 The resulting inflection points for each cultivar are plotted as a histogram. 
 
-```{r plot_inflections, eval=FALSE}
-ggplot(all_cultivars) +
-  geom_point(aes(x = gdd_cum, y = mean)) +
-  geom_line(aes(x = gdd_cum, y = mean_predict), color = "orange") +
-  geom_vline(aes(xintercept = inf_point)) +
-  facet_wrap(~cultivar, scales = "free_y") +
-  labs(x = "Cumulative growing degree days", y = "Canopy Height")
-ggplot(data.frame(inf_points = unique(all_cultivars$inf_point))) +
-  geom_histogram(aes(x = inf_points)) +
+```{r plot_inflections, warning=FALSE}
+ggplot(data.frame(inf_points = unique(all_cultivars$inf_x))) +
+  geom_histogram(aes(x = inf_points), bins = 300) +
   xlim(min(all_cultivars$gdd_cum), max(all_cultivars$gdd_cum)) +
   labs(x = "Inflection points", y = "Number")
 ```