Update logistic model equation and code

Author: KristinaRiemer

vignettes/04-synthesis-data.Rmd

@@ -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      = "~2018",
+                      date      = "~2017",
                       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=2018-01-01&until=2018-12-31', flatten = FALSE)
+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 <- weather$properties %>% 
   mutate(time = ymd_hms(weather$end_time))
 daily_values = weather %>% 
@@ -60,61 +60,79 @@ trait_weather_df <- full_join(trait_canopy_cover_day, daily_values, by = "day")
 
 We are interested in how growing degree days affects canopy cover. 
 To investigate this, we are going to model and plot their relationship. 
-We want to know the relationship for each cultivar, so we'll start of by determining the parameters of the model for one of the cultivars in our dataset. 
-We are using a logistic growth model here because it is appropriate for the shape of the GDD-cover relationship. (From Bolker? ... need to cite)
-
-This needs to be flushed out ... and the variable names should be consistent b/w equations and code
+We are using a logistic growth model here because it is appropriate for the shape of the GDD-cover relationship.
 
 The logistic growth model is specified as 
 
-$Y=\frac{z}{1+e^{-k(t-t_0)}}$
-
-$Y = \frac{c}{1+e^{a + b * \textrm{gdd}}}$
+$$y = \frac{c}{1+e^{a + b * \textrm{x}}}$$
 
-Where $Y$ is canopy cover, $c$ is the maximum canopy cover, $a$ is ... $b$ is ...
+where $y$ is the response variable canopy cover, $x$ is the predictor growing degree days, $c$ is the asymptote or maximum canopy cover, $a$ is the initial value for canopy cover, and $b$ is the steepness of the curve. (reference)
 
-From this we can find the asymptote and maximum growth rate. 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_ $\textrm{gdd}_i=\frac{\log{a}}{b}$.
-<!--
-or maybe start from the bolker form of the model (I added )
-if $Y=\frac{e^{a+b(t-t_0)}}{1+e^{a+b(t-t_0)}}$
+We want to know the relationship for each cultivar, so we'll start of by determining the parameters of the model for one of the cultivars in our dataset. 
+We provide estimated values for the asymptote $c$ and initial canopy cover value $a$, and provide canopy cover $y$ with corresponding growing degree days $x$ for one measurement of the chosen cultivar. 
 
-then $t-t_0$ is the number of days after planting.
+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. 
 
-then the inflection point is at $x=\frac{-a}{b}$ and the maximum rate of growth, $\max{\frac{dY}{dt}}=\frac{b}{4}$
--->
 ```{r model_get_parameters}
-single_cultivar <- trait_weather_df %>% 
-  filter(cultivar == "PI656026")
-cap     <- 150
-initial <- 25
-mean    <- single_cultivar$mean[15]
-gdd_cum <- single_cultivar$gdd_cum[15]
-rate    <- ((log((cap/mean) - 1)) - initial)/gdd_cum
-model_single_cultivar <- nls(mean ~ cap / (1 + exp(initial + rate * gdd_cum)), 
-                             start = list(cap = cap, initial = initial, rate = rate),
+# 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))
+
+c <- 51
+a <- 10
+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)
+summary(model_single_cultivar)
+coef(model_single_cultivar)
+
 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)))
+ggplot(single_cultivar) +
+  geom_point(aes(x = gdd_cum, y = mean)) +
+  geom_line(aes(x = gdd_cum, y = mean_predict), color = "orange") +
+  labs(x = "Cumulative growing degree days", y = "Canopy Height")
+```
+
+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])
+
+ggplot(single_cultivar) +
+  geom_point(aes(x = gdd_cum, y = mean)) +
+  geom_line(aes(x = gdd_cum, y = mean_predict), color = "orange") +
+  geom_hline(yintercept = inf_y, linetype = "dashed") +
+  geom_vline(xintercept = inf_x) +
+  labs(x = "Cumulative growing degree days", y = "Canopy Height")
 ```
 
 We then use the parameters from a single cultivar to run a model for each of the rest of the cultivars. 
 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}
+```{r model_all_cultivars, eval=FALSE}
 all_cultivars <- c(day = as.double(), cultivar = as.character(), mean = as.numeric(), 
                    gdd_cum = as.numeric(), mean_predict = 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(mean ~ cap / (1 + exp(initial + rate * gdd_cum)), 
-                            start = list(cap = cap, initial = initial, rate = rate),
+  each_cultivar_model <- nls(y ~ c / (1 + exp(a + b * g)), 
+                            start = list(c = c, a = a, b = b),
                             data = each_cultivar_df)
-  model_cap     <- coef(each_cultivar_model)[1]
-  model_initial <- coef(each_cultivar_model)[2]
-  model_rate    <- coef(each_cultivar_model)[3]
+  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_cap / (1 + exp(model_initial + model_rate * gdd_cum)), 
-             inf_point = ((log((model_cap / 100) - 1)) - model_initial) / model_rate)
+      mutate(mean_predict = model_c / (1 + exp(model_a + model_b * g)), 
+             inf_point = ((log((model_c / 100) - 1)) - model_a) / model_b)
   all_cultivars <- rbind(each_cultivar_df, all_cultivars)
 }
 ggplot(all_cultivars) +
@@ -130,7 +148,7 @@ 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}
+```{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") +