version 0.0.2

DESCRIPTION

@@ -0,0 +1,32 @@
+Package: phenolocrop
+Type: Package
+Title: Time-Series Models to the Crop Phenology
+Version: 0.0.2
+Authors@R: c(
+  person("Shoji", "Taniguchi", email = "taniguchis532@affrc.go.jp", role = c("aut", "cre")),
+  person("The National Agriculture and Food Research Organization (NARO)", role = "cph")
+  )
+Maintainer: Shoji Taniguchi <taniguchis532@affrc.go.jp>
+Description: Fit a time-series model 
+    to a crop phenology data, such as time-series rice canopy height.
+    This package returns the model parameters as the summary statistics of crop phenology,
+    and these parameters will be useful to characterize the growth pattern of each cultivar and 
+    predict manually-measured traits, such as days to heading and biomass.
+    Please see Taniguchi et al. (2022) <doi:10.3389/fpls.2022.998803> for detail.
+    This package has been designed for scientific use.  
+    Use for commercial purposes shall not be allowed.
+Imports: purrr
+Suggests: dplyr, knitr, testthat (>= 3.0.0)
+License: CC BY-NC 4.0
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 7.2.1
+Depends: R (>= 4.2.0)
+Config/testthat/edition: 3
+VignetteBuilder: knitr
+NeedsCompilation: no
+Packaged: 2023-07-12 01:56:54 UTC; taniguchis532
+Author: Shoji Taniguchi [aut, cre],
+  The National Agriculture and Food Research Organization (NARO) [cph]
+Repository: CRAN
+Date/Publication: 2023-07-13 14:00:02 UTC

MD5

@@ -0,0 +1,16 @@
+cb1e80c1173719ac95fac479bd3274b4 *DESCRIPTION
+bfe7482b5ce7dfca2eb59d040860820a *NAMESPACE
+fdeac7a66702ac6572c0bc8c5043cfea *R/ShowStartupMessage.R
+d0bfe3d11b1baae616a20acd583840de *R/logisLateDicrease.R
+aeb0e2bd98dbf1ee8eda49b12c5273ae *R/phenololine.R
+4eaa75f2e49eab9f42002286b24fdbf9 *R/riceCH_eg.R
+708c6bf38e4d0c0084a8803256008299 *README.md
+bdad5015cf642a758478d88c2049c5aa *build/vignette.rds
+bf6f23db6d5c1214e505c8efc0c67e6f *data/riceCH_eg.rda
+5f56bf7be68b6bf1e22808a9636ba964 *inst/doc/logisLateDicr.R
+a74951c65b1527f17b3eda93cb495977 *inst/doc/logisLateDicr.Rmd
+2c8568957b4cd3984edcec20a6f1ccf3 *inst/doc/logisLateDicr.html
+83d0a972b8aecb0eb54b90df90e2d45e *man/logisLateDicr.Rd
+4bda39ce1eecdd64624e2f35c2f6b0c1 *man/phenololine.Rd
+cb5a0e65b00d4021e948180fa6551356 *man/riceCH_eg.Rd
+a74951c65b1527f17b3eda93cb495977 *vignettes/logisLateDicr.Rmd

NAMESPACE

@@ -0,0 +1,7 @@
+# Generated by roxygen2: do not edit by hand
+
+export(logisLateDicr)
+export(phenololine)
+importFrom(stats,lm)
+importFrom(stats,nls)
+importFrom(stats,predict)

R/ShowStartupMessage.R

@@ -0,0 +1,3 @@
+.onAttach <- function(libname, pkgname) {
+  packageStartupMessage("Package phenolocrop is designed for scientific use.\nUse for commercial purposes shall not be allowed.\nCopyright (C) 2023 The National Agriculture and Food Research Organization. All rights reserved.")
+}

R/logisLateDicrease.R

@@ -0,0 +1,109 @@
+#' Three-parameter logistic with the decrease in the late growth phase
+#'
+#' @description Apply a time-series model, three-parameter logistic with the decrease in the
+#' late growth phase, to the time-series trait data.
+#' This function was originally developed for the time-series data of rice canopy height.
+#' Fitting the time-series model is done by the "two-step procedure".
+#' For more information, see Taniguchi et al. (2022).
+#'
+#' @param dat data.frame including date and trait (e.g. canopy height).
+#' @param x Column name (character) for the date after sowing or plantingl.
+#' @param y Column name (character) for the trait.
+#' @param returnModels Logical value whether to return the time-series model object. Default is F.
+#' @param start Start values to estimate 'd0', 'r' and 'a'. Default is 'c(d0 = 50, r = 0.05, a = 0.0001)'.
+#' @param upper Upper bounds to estimate 'd0', 'r' and 'a'. Default is 'c(d0 = 200, r = Inf, a = 1)'.
+#' @param lower Lower bounds to estimate 'd0', 'r' and 'a'. Default is 'c(d0 = 0, r = 0, a = 0)'.
+#'
+#' @return
+#' logisLateDicr function returns the vector of estimated parameter values.
+#' If returnModels = TRUE, this function also returns the cubic polynomial
+#' regression object and logistic with the decrease in the
+#' late growth phase regression object.
+#'
+#' @examples
+#' library(phenolocrop)
+#' riceCH_eg |>
+#'    logisLateDicr("x", "height")
+#'
+#' @references
+#' S. Taniguchi et al., (2022) Prediction of heading date, culm length, and
+#' biomass from canopy-height-related parameters derived from time-series UAV observations of rice.
+#' Frontiers in Plant Science, 13:998803.
+#'
+#' @importFrom stats lm
+#' @importFrom stats nls
+#' @importFrom stats predict
+#'
+#' @export
+
+
+logisLateDicr <- function(dat, x, y,
+                          returnModels = FALSE,
+                          start = c(d0 = 50, r = 0.05, a = 0.0001),
+                          upper = c(d0 = 200, r = Inf, a = 1),
+                          lower = c(d0 = 0, r = 0, a = 0)){
+
+  x <- dat |> purrr::pluck(x)
+  y <- dat |> purrr::pluck(y)
+
+  ### 1st step: cubic polynomial regression
+
+  dat_poly <- data.frame(y = as.numeric(y),
+                         x1 = x,
+                         x2 = x^2,
+                         x3 = x^3)
+
+  lm_poly <- lm(y ~ x1 + x2 + x3, data = dat_poly)
+  x_max <- max(x)
+  x_new_vec <- 1:x_max
+  y_pred <- predict(lm_poly,
+                    newdata = list(x1 = x_new_vec,
+                                   x2 = x_new_vec^2,
+                                   x3 = x_new_vec^3))
+
+  d1 <- findMaxima(y_pred)
+  if(is.null(d1)){
+    stop(message =
+           "The cubic polynomial did not detect maximum time point 'd1'")
+  }
+
+  ### 2nd step: 3-parameter logistic with decrease in the late growth phase
+
+  K <- max(y)
+  x0 <- x
+  x0 <- x0 - d1
+  id_zero <- which(x < d1)
+  x0[id_zero] <- 0
+  dat_logi <- data.frame(y = y,
+                         x = x,
+                         x0 = x0)
+
+  nls_logi <- tryCatch(nls(y ~ K/(1 + exp(r * (d0 - x))) - a * x0^2,
+                           data = dat_logi,
+                           start = start,
+                           algorithm = "port",
+                           upper = upper,
+                           lower = lower),
+                           error = function(e){
+                             stop(message =
+                                    "Failed to apply nls. Modifying the start, upper, and lower values may solve the problem.")
+                           })
+
+  nls_logi_summ <- summary(nls_logi)
+  param <- nls_logi_summ$coefficients[, "Estimate"]
+  param <- c(K = K, param, d1 = d1)
+  if(returnModels == FALSE){
+    return(param)
+  }else{
+    return(list(param, lm_poly, nls_logi))
+  }
+}
+
+findMaxima <- function(x_vec){
+  n <- length(x_vec) - 2
+  for(i in 2:(n + 1)){
+    if(x_vec[i - 1] < x_vec[i] && x_vec[i + 1] < x_vec[i]){
+      return(i)
+    }
+  }
+}

R/phenololine.R

@@ -0,0 +1,43 @@
+#' Return time-series curve given the time-series model and parameter values
+#'
+#' @description phenololine function returns the predicted values given the
+#' model name and model parameters.
+#' @param param Vector of model parameter values.
+#' @param x Vector of time (e.g. vector of dates).
+#' @param method Character name of the time-series model. At present, only "logisLateDicr" is accepted.
+#' @return phenololine function returns the trait values given x.
+#' @details
+#' If method = "logisLateDicr", param should be the vector of c(K, d0, r, a, d1).
+#' @examples
+#' library(phenolocrop)
+#' y <- riceCH_eg |>
+#'    logisLateDicr("x", "height") |>
+#'    phenololine(x = 1:160, method = "logisLateDicr")
+#' plot(1:160, y, type = "l")
+#'
+#' @export
+
+phenololine <- function(param, x, method){
+  if(method == "logisLateDicr"){
+    if(length(param) != 5){
+      stop(message =
+             "\"logisLateDicr\" method need 5 parameters (K, d0, r, a, d1).")
+    }
+    K <- param[1]
+    d0 <- param[2]
+    r <- param[3]
+    a <- param[4]
+    d1 <- param[5]
+
+    x0 <- x
+    x0 <- x0 - d1
+    id_zero <- which(x < d1)
+    x0[id_zero] <- 0
+
+    y <- K/(1 + exp(r * (d0 - x))) - a * x0^2
+
+    return(y)
+  }else{
+    stop(message = "Method name is invalid.")
+  }
+}

R/riceCH_eg.R

@@ -0,0 +1,12 @@
+#' Rice canopy height example data
+#'
+#' This data was a simulated rice CH data to show the usage of this package.
+#'
+#' @format
+#' \describe{
+#'   \item{date}{Date of observation}
+#'   \item{height}{CH}
+#'   \item{id}{Cultivar ID}
+#'   \item{x}{Dayes after sowing}
+#' }
+"riceCH_eg"

README.md

@@ -0,0 +1,17 @@
+## phenolocrop
+
+### Overview
+
+`phenolocrop` has been developed for the time-series analysis of crop traits. In time-series analysis, one data is composed of serial observation and interpreted by means of "time-series models." Therefore, `phenolocrop` is aimed at applying appropriate time-series model to the time-series observation of crop, which is called phenology.
+
+### Installation
+
+    install.packages("phenolocrop")
+
+### Target traits and models
+
+At the latest version (0.0.1), `phenolocrop` package offers a model for rice canopy height data.
+
+-   **Trait**: Rice canopy height.
+
+    -   **Model**: Three-parameter logistic with the decrease in the late growth phase. Use `logisLateDicr` function.

inst/doc/logisLateDicr.R

@@ -0,0 +1,24 @@
+## ---- include = FALSE---------------------------------------------------------
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+
+## ----setup--------------------------------------------------------------------
+library(phenolocrop)
+
+## -----------------------------------------------------------------------------
+riceCH_eg
+
+## -----------------------------------------------------------------------------
+riceCH_eg |>
+  logisLateDicr("x", "height")
+
+## -----------------------------------------------------------------------------
+x_vec <- 1:max(riceCH_eg$x)
+y <- riceCH_eg |>
+  logisLateDicr("x", "height") |>
+  phenololine(x = x_vec, method = "logisLateDicr")
+plot(x_vec, y, type = "l", ylab = "rice CH", xlab = "Dayes after sowing")
+points(riceCH_eg$x, riceCH_eg$height)
+

inst/doc/logisLateDicr.Rmd

@@ -0,0 +1,68 @@
+---
+title: "logisLateDicr"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{logisLateDicr}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+## 1. Package loading
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(phenolocrop)
+```
+
+## 2. Data
+
+The data frame should include the columns of the objective trait and time point.
+In the example rice CH data, date column is for the measuring date, height is for CH (m),
+id is for the cultivar id, and x is for the days after sowing.
+The following CH time-series data was generated by the computer simulation.
+
+```{r}
+riceCH_eg
+```
+
+## 3. Estimate model parameters
+
+```{r}
+riceCH_eg |>
+  logisLateDicr("x", "height")
+```
+
+## 4. Plot time-series model and measured CH data
+
+```{r}
+x_vec <- 1:max(riceCH_eg$x)
+y <- riceCH_eg |>
+  logisLateDicr("x", "height") |>
+  phenololine(x = x_vec, method = "logisLateDicr")
+plot(x_vec, y, type = "l", ylab = "rice CH", xlab = "Dayes after sowing")
+points(riceCH_eg$x, riceCH_eg$height)
+```
+
+## 5. Time-series model
+
+Time-series model of `logisLateDicr` for the rice CH is as the following.
+
+$$
+\mathrm{CH}=
+\begin{cases}
+\frac{K}{1+\mathrm{exp}(r(d_0-x))} & (x\le d_1)\\
+\frac{K}{1+\mathrm{exp}(r(d_0-x))}-a(x-d_1)^2 & (x>d_1)
+\end{cases}
+$$
+
+## 6. Reference
+
+S. Taniguchi et al., (2022) Prediction of heading date, culm length, and
+biomass from canopy-height-related parameters derived from time-series UAV observations of rice.
+Frontiers in Plant Science, 13:998803.