Unable to retrieve parameters from a dynamic model

Hello,

TL;DR

I apologize if these are too many questions for a single post – I’ll be grateful if you can help me to solve question 1, and I could keep the rest of the questions for future posts if you think that would be more appropriate. My questions, in order of decreasing importance, are:

  1. Why am I unable to retrieve the parameters of this model? Model (soy_dev.stan) plus fake data and sampling statement (soy_dev_fake_data.R) are attached
  2. Is my approach to modifying mu in the model block of soy_dev.stan correct? Or would it be better instead to create a custom function in the functions block of soy_dev.stan to later determine mu?
  3. In this example, we are simulating fake data for a single soybean genotype, but do you think that this model can be scaled into a hierarchical model, where parameters dOptFiFl, dOptSdPm, dOptSdPm, Pcrit and Psen can vary between several soybean genotypes?
  4. How would you call this model? From my ignorance this is a “dynamic model” but my definition may be too broad!

Longer story

I am interested in modeling how a soybean crop develops. Development in this species depends on temperature and day length (i.e. photoperiod) – it is faster at an optimum temperature (around 30 °C) and when day length shortens. Inspired by Ben Bolker’s Ecological Models and Data in R, Chapter 11, Dynamic models, I tried to write a dynamic model with a daily time step. For this, I coded developmental stages, which are discrete events, into a continuous scale (called R) as follows:

  • R = 0: plant emergence
  • R = 1: end of juvenile phase
  • R = 2: flowering initiation
  • R = 3: flowering
  • R = 4: start of seed growth
  • R = 5: end of seed growth (physiological maturity)

Temperature and photoperiod are observed daily (weather.csv).

Under ideal temperature and photoperiod conditions, development progreses at a rate (with units 1/day) equal to:

  • 1/dOptEmJuv: between emergence (R = 0) and the end of the juvenile stage (R = 1)
  • 1/dOptJuvFi: between the end of the juvenile phase (R = 1) and flowering induction (R = 2)
  • 1/dOptFiFl: between flowering induction (R = 2) and flowering (R = 3)
  • 1/dOptFlSd: between flowering (R = 3) and the start of seed growth (R = 4)
  • 1/dOptSdPm: between the start (R = 4) and the end of seed growth (R = 5)

Under non-optimal temperature and photoperiod conditions, these developmental rates are penalized by temperature (a, in the following figure) and photoperiod multipliers (b) that look like this:

The parameters of the equation describing the temperature multiplier (i.e. Tmin, Topt and Tmax) don’t change much among soybean genotypes, so I fixed them as constants (see lines 24:32 in soy_dev_fake_data.R, and lines 3:19 in soy_dev.stan).

In turn, the parameters describing the photoperiod multiplier can vary between soybeans, and I am therefore interested in their estimation. There’s a threshold photoperiod (Pcrit), below which developmental rate is maximum (multiplier = 1), and a photoperiod sensitivity (Psen). For photoperiods above Pcrit, the developmental rate is reduced at a rate equal to Psen.

Parameters dOptEmJuv and dOptJuvFi are assumed constant and are thus fixed as 3.7 days (transformed data block in soy_dev.stan).

To test the model I generated some fake data, starting from the following parameters

pars <- c(Topt      = 30,
          Pcrit     = 13,
          Psen      = 0.29,
          dOptEmJuv = 3.7,
          dOptJuvFi = 3.7,
          dOptFiFl  = 16,
          dOptFlSd  = 17,
          dOptSdPm  = 39,
          sigma     = 0.1)

and simulating for three emergence dates. Simulations look like the following figure:

In summary, I am interested in estimating parameters dOptFiFl, dOptSdPm, dOptSdPm, Pcrit and Psen.

I am not proficient coding Stan scripts so I used brms::make_stancode() to make an initial Stan script (soy_dev_to_be_edited.stan) which I later modified (soy_dev.stan). Comparing both files, you’ll notice that

  • I modified how mu is estimated in the model block (lines 87:104)
  • I added a one-dimensional array called time, so that I can define when mu = 0 at time = 1 (line 50)

There are no warnings after fitting this model, but I am not able to retrieve the original parameter values (e.g. b_dOptSdPm), whereas the posterior distributions of other parameters resemble that of their priors (e.g. b_dOptFiFl, b_dOptFlSd):

Inference for Stan model: soy_dev.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                 mean se_mean    sd    2.5%     25%     50%     75%   97.5%
b_dOptFiFl[1]   14.95    0.06  5.90    3.51   11.07   14.85   18.87   26.53
b_dOptFlSd[1]   15.97    0.07  7.00    2.52   11.09   15.96   20.86   29.47
b_dOptSdPm[1]   13.13    0.01  0.45   12.19   12.88   13.18   13.42   13.91
b_Pcrit[1]      11.09    0.01  0.38   10.28   10.89   11.12   11.34   11.71
b_Psen[1]        0.16    0.00  0.02    0.13    0.15    0.16    0.17    0.20
sigma            0.94    0.00  0.03    0.89    0.92    0.94    0.97    1.01
lp__          -921.31    0.46 15.61 -952.03 -931.78 -921.46 -910.43 -891.56
              n_eff Rhat
b_dOptFiFl[1]  9135    1
b_dOptFlSd[1]  8718    1
b_dOptSdPm[1]  2254    1
b_Pcrit[1]     2017    1
b_Psen[1]      2872    1
sigma          3317    1
lp__           1134    1

Samples were drawn using NUTS(diag_e) at Tue Dec  3 15:56:10 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

My questions are:

  1. Why am I unable to retrieve the parameters of this model? Model (soy_dev.stan) plus fake data and sampling statement (soy_dev_fake_data.R) are attached
  2. Is my approach to modifying mu in the model block of soy_dev.stan correct? Or would it be better instead to create a custom function in the functions block of soy_dev.stan to later determine mu?
  3. In this example, we are simulating fake data for a single soybean genotype, but do you think that this model can be scaled into a hierarchical model, where parameters dOptFiFl, dOptSdPm, dOptSdPm, Pcrit and Psen can vary between several soybean genotypes?
  4. How would you call this model? From my ignorance this is a “dynamic model” but my definition may be too broad!

Thanks in advance.

soy_dev.stan (4.0 KB) soy_dev_fake_data.R (5.4 KB) soy_dev_to_be_edited.stan (3.3 KB) weather.csv (100.3 KB)

sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 10 (buster)

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  grid      methods  
[8] base     

other attached packages:
 [1] brms_2.10.0        Rcpp_1.0.3         rstan_2.19.2       StanHeaders_2.19.0
 [5] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.3        purrr_0.3.3       
 [9] readr_1.3.1        tidyr_1.0.0        tibble_2.1.3       tidyverse_1.3.0   
[13] ggplot2_3.2.1     

loaded via a namespace (and not attached):
 [1] nlme_3.1-140         matrixStats_0.55.0   fs_1.3.1            
 [4] xts_0.11-2           lubridate_1.7.4      threejs_0.3.1       
 [7] httr_1.4.1           tools_3.6.1          backports_1.1.5     
[10] R6_2.4.1             DT_0.10              DBI_1.0.0           
[13] lazyeval_0.2.2       colorspace_1.4-1     withr_2.1.2         
[16] tidyselect_0.2.5     gridExtra_2.3        prettyunits_1.0.2   
[19] processx_3.4.1       Brobdingnag_1.2-6    compiler_3.6.1      
[22] cli_1.1.0            rvest_0.3.5          shinyjs_1.0         
[25] xml2_1.2.2           labeling_0.3         colourpicker_1.0    
[28] scales_1.1.0         dygraphs_1.1.1.6     ggridges_0.5.1      
[31] callr_3.3.2          digest_0.6.23        base64enc_0.1-3     
[34] pkgconfig_2.0.3      htmltools_0.4.0      dbplyr_1.4.2        
[37] fastmap_1.0.1        htmlwidgets_1.5.1    rlang_0.4.2         
[40] readxl_1.3.1         rstudioapi_0.10      shiny_1.4.0         
[43] farver_2.0.1         generics_0.0.2       zoo_1.8-6           
[46] jsonlite_1.6         gtools_3.8.1         crosstalk_1.0.0     
[49] inline_0.3.15        magrittr_1.5         loo_2.1.0           
[52] bayesplot_1.7.1      Matrix_1.2-17        munsell_0.5.0       
[55] abind_1.4-5          lifecycle_0.1.0      stringi_1.4.3       
[58] pkgbuild_1.0.6       plyr_1.8.4           parallel_3.6.1      
[61] promises_1.1.0       crayon_1.3.4         miniUI_0.1.1.1      
[64] lattice_0.20-38      haven_2.2.0          hms_0.5.2           
[67] zeallot_0.1.0        ps_1.3.0             pillar_1.4.2        
[70] igraph_1.2.4.2       markdown_1.1         shinystan_2.5.0     
[73] codetools_0.2-16     reshape2_1.4.3       stats4_3.6.1        
[76] rstantools_2.0.0     reprex_0.3.0         glue_1.3.1          
[79] modelr_0.1.5         vctrs_0.2.0          httpuv_1.5.2        
[82] cellranger_1.1.0     gtable_0.3.0         assertthat_0.2.1    
[85] mime_0.7             xtable_1.8-4         broom_0.5.2         
[88] coda_0.19-3          later_1.0.0          rsconnect_0.8.15    
[91] shinythemes_1.1.2    bridgesampling_0.7-2
1 Like

Potential answers to this question will widen our ability to tackle agronomic questions through modelling physiological processes using Stan, so I will be following this thread up close!

1 Like