In this vignette we show how to use the do.fit = FALSE flag to debug a WHAM model.

# devtools::install_github("timjmiller/wham", dependencies=TRUE)
library(wham)

Load the ASAP3 data file.

wham.dir <- find.package("wham")
asap3 <- read_asap3_dat(file.path(wham.dir,"extdata","ex7_SNEMAYT.dat"))

Prepare the WHAM model (m3 from example 1).

input <- prepare_wham_input(asap3, recruit_model=2, model_name="Ex 1: SNEMA Yellowtail Flounder",
                              NAA_re = list(sigma="rec+1", cor="iid"))

Try to fit the model… uh oh.

mod <- fit_wham(input, do.osa = F, do.retro = F)
#> Order of parameters:
#>  [1] "log_catch_sig_scale"  "log_index_sig_scale"  "log_N1_pars"         
#>  [4] "log_NAA_sigma"        "trans_NAA_rho"        "log_NAA"             
#>  [7] "logR_proj"            "mean_rec_pars"        "logit_q"             
#> [10] "q_re"                 "q_repars"             "q_prior_re"          
#> [13] "logit_selpars"        "selpars_re"           "sel_repars"          
#> [16] "catch_paa_pars"       "index_paa_pars"       "F_devs"              
#> [19] "log_F1"               "M_a"                  "M_re"                
#> [22] "M_repars"             "log_b"                "Ecov_re"             
#> [25] "Ecov_beta"            "Ecov_process_pars"    "Ecov_obs_logsigma"   
#> [28] "Ecov_obs_sigma_par"   "Ecov_obs_logsigma_re"
#> Not matching template order:
#>  [1] "mean_rec_pars"        "logit_q"              "q_prior_re"          
#>  [4] "q_re"                 "q_repars"             "log_F1"              
#>  [7] "F_devs"               "log_N1_pars"          "log_NAA_sigma"       
#> [10] "trans_NAA_rho"        "log_NAA"              "logR_proj"           
#> [13] "logit_selpars"        "selpars_re"           "sel_repars"          
#> [16] "catch_paa_pars"       "index_paa_pars"       "M_a"                 
#> [19] "M_re"                 "M_repars"             "log_b"               
#> [22] "log_catch_sig_scale"  "log_index_sig_scale"  "Ecov_re"             
#> [25] "Ecov_beta"            "Ecov_process_pars"    "Ecov_obs_logsigma"   
#> [28] "Ecov_obs_logsigma_re" "Ecov_obs_sigma_par"  
#> Your parameter list has been re-ordered.
#> (Disable this warning with checkParameterOrder=FALSE)
#> iter: 1  Error in iterate(par) : 
#>   Newton dropout because inner gradient had non-finite components.
#> Warning in stats::nlminb(model$par, model$fn, model$gr, control = list(iter.max
#> = 1000, : NA/NaN function evaluation
#> iter: 1  Error in iterate(par) : 
#>   Newton dropout because inner gradient had non-finite components.
#> Error in ff(x, order = 1) : 
#>   inner newton optimization failed during gradient calculation
#> outer mgc:  NaN
#> Error in stats::nlminb(model$par, model$fn, model$gr, control = list(iter.max = 1000, : NA/NaN gradient evaluation

What’s wrong? It looks like the likelihood function is returning NaN. It is often easier to diagnose problems like this using the unoptimized model, i.e. look at everything using the initial parameter values. WHAM includes a do.fit = F flag in fit_wham to return the unoptimized model object returned by TMB::MakeADFun. Let’s see how it works.

mod <- fit_wham(input, do.fit = F)
#> Order of parameters:
#>  [1] "log_catch_sig_scale"  "log_index_sig_scale"  "log_N1_pars"         
#>  [4] "log_NAA_sigma"        "trans_NAA_rho"        "log_NAA"             
#>  [7] "logR_proj"            "mean_rec_pars"        "logit_q"             
#> [10] "q_re"                 "q_repars"             "q_prior_re"          
#> [13] "logit_selpars"        "selpars_re"           "sel_repars"          
#> [16] "catch_paa_pars"       "index_paa_pars"       "F_devs"              
#> [19] "log_F1"               "M_a"                  "M_re"                
#> [22] "M_repars"             "log_b"                "Ecov_re"             
#> [25] "Ecov_beta"            "Ecov_process_pars"    "Ecov_obs_logsigma"   
#> [28] "Ecov_obs_sigma_par"   "Ecov_obs_logsigma_re"
#> Not matching template order:
#>  [1] "mean_rec_pars"        "logit_q"              "q_prior_re"          
#>  [4] "q_re"                 "q_repars"             "log_F1"              
#>  [7] "F_devs"               "log_N1_pars"          "log_NAA_sigma"       
#> [10] "trans_NAA_rho"        "log_NAA"              "logR_proj"           
#> [13] "logit_selpars"        "selpars_re"           "sel_repars"          
#> [16] "catch_paa_pars"       "index_paa_pars"       "M_a"                 
#> [19] "M_re"                 "M_repars"             "log_b"               
#> [22] "log_catch_sig_scale"  "log_index_sig_scale"  "Ecov_re"             
#> [25] "Ecov_beta"            "Ecov_process_pars"    "Ecov_obs_logsigma"   
#> [28] "Ecov_obs_logsigma_re" "Ecov_obs_sigma_par"  
#> Your parameter list has been re-ordered.
#> (Disable this warning with checkParameterOrder=FALSE)

This runs without an error. The optimization failed because the likelihood was NaN, and now we can see which of the likelihood components was responsible. To do so, we need to look at the REPORTed objects with "nll" in their name. We can use the $report() function from TMB.

therep = mod$report()

See all of the REPORTed objects.

names(therep)
#>  [1] "pred_log_catch"        "q_prior_re"            "selAA"                
#>  [4] "sigma_q"               "FAA_tot"               "log_SPR0"             
#>  [7] "Ecov_beta"             "log_FXSPR"             "Ecov_obs_sigma_par"   
#> [10] "NAA_devs"              "mean_rec_pars"         "QAA"                  
#> [13] "nll_catch_acomp"       "M_a"                   "pred_catch_paa"       
#> [16] "MAA"                   "log_Y_FXSPR_static"    "SSB"                  
#> [19] "F"                     "nll_agg_indices"       "nll"                  
#> [22] "log_YPR_FXSPR"         "selpars"               "nll_Ecov_obs"         
#> [25] "nll_Ecov_obs_sig"      "Ecov_re"               "log_SPR_FXSPR"        
#> [28] "log_SPR0_static"       "pred_indices"          "nll_NAA"              
#> [31] "log_FXSPR_static"      "pred_IAA"              "logit_q_mat"          
#> [34] "pred_CAA"              "nll_agg_catch"         "nll_sel"              
#> [37] "log_FXSPR_iter"        "nll_q"                 "logit_selpars"        
#> [40] "Fbar"                  "nll_Ecov"              "log_FXSPR_iter_static"
#> [43] "Ecov_obs_sigma"        "rho_q"                 "q_re"                 
#> [46] "Ecov_process_pars"     "sel_repars"            "nll_index_acomp"      
#> [49] "ZAA"                   "pred_log_indices"      "M_re"                 
#> [52] "selpars_re"            "Ecov_x"                "NAA"                  
#> [55] "pred_catch"            "log_SSB_FXSPR"         "M_repars"             
#> [58] "log_SSB_FXSPR_static"  "pred_NAA"              "FAA"                  
#> [61] "log_SPR_FXSPR_static"  "nll_q_prior"           "log_Y_FXSPR"          
#> [64] "pred_index_paa"        "Ecov_out"              "nll_M"                
#> [67] "log_YPR_FXSPR_static"  "q"

Now just get the objects with "nll" in their name, and sum over all individual values.

sapply(grep("nll",names(therep),value=T), function(x) sum(therep[[x]]))
#>  nll_catch_acomp  nll_agg_indices              nll     nll_Ecov_obs 
#>        2844.2380        1562.8441              NaN           0.0000 
#> nll_Ecov_obs_sig          nll_NAA    nll_agg_catch          nll_sel 
#>           0.0000         346.2991              NaN           0.0000 
#>            nll_q         nll_Ecov  nll_index_acomp      nll_q_prior 
#>           0.0000           0.0000        3921.8316           0.0000 
#>            nll_M 
#>           0.0000

The likelihood components that are equal to 0 are not used in the model (no random effects on M, Ecov, selectivity, etc.). nll is the total likelihood and is NaN. The one troublesome component is nll_agg_catch. This is the total (aggregate) catch from the fishery in each year. It could be an issue with the catch data or the model-predicted catch, since both are in the likelihood calculation.

Search the WHAM .cpp file for “nll_agg_catch”. We find that on line 888, this depends on agg_catch (the catch data) and pred_catch (model-predicted catch).

The catch data looks ok.

input$data$agg_catch
#>             [,1]
#>  [1,] 14549.0000
#>  [2,] 17088.0000
#>  [3,]  5732.0000
#>  [4,]  3436.0000
#>  [5,]  5223.0000
#>  [6,]  8085.0000
#>  [7,]  9883.0000
#>  [8,]  8021.0000
#>  [9,]  6607.0000
#> [10,] 15764.0000
#> [11,] 22211.0000
#> [12,] 11225.0000
#> [13,]  4817.0000
#> [14,]  4620.0000
#> [15,]  2652.0000
#> [16,]  2782.0000
#> [17,]  8349.0000
#> [18,] 17916.0000
#> [19,]  6430.0000
#> [20,]  2695.0000
#> [21,]   771.0008
#> [22,]   735.0004
#> [23,]   343.0004
#> [24,]   759.0001
#> [25,]  1222.0000
#> [26,]  1087.0000
#> [27,]  1403.0000
#> [28,]  1397.0000
#> [29,]  1449.0000
#> [30,]   945.0005
#> [31,]   666.0005
#> [32,]   619.0002
#> [33,]   346.0002
#> [34,]   396.0000
#> [35,]   502.0004
#> [36,]   583.0006
#> [37,]   453.0006
#> [38,]   290.9995
#> [39,]   390.0002
#> [40,]   563.0002
#> [41,]   645.9999
#> [42,]   625.0000
#> [43,]   337.0005
#> [44,]   151.9999

The predicted catch is calculated on lines 879-880. pred_catch depends on NAA (numbers-at-age), FAA (F at age), ZAA (Z at age), and waa (weight-at-age data).

NAA looks ok.

therep$NAA
#>            [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
#>  [1,] 134809.00 11710.97 25945.08 17508.26 11476.00 13803.70
#>  [2,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [3,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [4,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [5,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [6,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [7,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [8,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#>  [9,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [10,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [11,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [12,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [13,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [14,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [15,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [16,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [17,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [18,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [19,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [20,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [21,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [22,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [23,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [24,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [25,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [26,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [27,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [28,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [29,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [30,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [31,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [32,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [33,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [34,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [35,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [36,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [37,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [38,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [39,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [40,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [41,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [42,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [43,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47
#> [44,]  22026.47 22026.47 22026.47 22026.47 22026.47 22026.47

FAA looks ok - no issues with F or selectivity.

therep$FAA[,1,] # middle dim is n.fleets = 1
#>            [,1]       [,2]       [,3]       [,4]       [,5] [,6]
#>  [1,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [2,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [3,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [4,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [5,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [6,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [7,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [8,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#>  [9,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [10,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [11,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [12,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [13,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [14,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [15,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [16,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [17,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [18,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [19,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [20,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [21,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [22,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [23,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [24,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [25,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [26,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [27,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [28,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [29,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [30,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [31,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [32,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [33,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [34,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [35,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [36,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [37,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [38,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [39,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [40,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [41,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [42,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [43,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1
#> [44,] 0.0501933 0.07552447 0.09057153 0.09692876 0.09922137  0.1

ZAA looks ok - no issues with M either.

therep$ZAA
#>            [,1]      [,2]      [,3]      [,4]      [,5]     [,6]
#>  [1,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [2,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [3,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [4,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [5,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [6,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [7,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [8,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#>  [9,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [10,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [11,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [12,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [13,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [14,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [15,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [16,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [17,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [18,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [19,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [20,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [21,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [22,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [23,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [24,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [25,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [26,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [27,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [28,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [29,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [30,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [31,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [32,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [33,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [34,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [35,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [36,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [37,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [38,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [39,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [40,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [41,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [42,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [43,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099
#> [44,] 0.4551933 0.4115245 0.3865715 0.3719288 0.3552214 0.331099

Ah, here is the problem. The weight at age data has an entry of -99. This means that pred_catch is negative on line 880, and we take the log of a negative number on line 883.

input$data$waa[1,,]
#>          [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
#>  [1,]   0.210 0.296 0.348 0.374 0.382 0.428
#>  [2,]   0.203 0.308 0.352 0.396 0.439 0.457
#>  [3,]   0.218 0.289 0.376 0.432 0.435 0.481
#>  [4,]   0.228 0.303 0.408 0.498 0.499 0.557
#>  [5,]   0.215 0.283 0.381 0.504 0.513 0.542
#>  [6,]   0.234 0.292 0.383 0.536 0.662 0.656
#>  [7,]   0.189 0.301 0.364 0.475 0.590 0.662
#>  [8,]   0.206 0.281 0.384 0.500 0.682 0.925
#>  [9,]   0.140 0.262 0.342 0.474 0.596 0.650
#> [10,]   0.226 0.263 0.353 0.499 0.660 0.833
#> [11,]   0.175 0.261 0.339 0.496 0.668 0.819
#> [12,]   0.182 0.237 0.295 0.388 0.487 0.656
#> [13,]   0.183 0.260 0.365 0.408 0.504 0.608
#> [14,]   0.186 0.284 0.331 0.463 0.587 0.642
#> [15,]   0.247 0.268 0.353 0.404 0.520 0.631
#> [16,]   0.270 0.293 0.396 0.493 0.611 0.821
#> [17,]   0.061 0.216 0.275 0.489 0.735 0.957
#> [18,]   0.204 0.255 0.290 0.366 0.613 0.884
#> [19,]   0.090 0.214 0.294 0.378 0.664 0.798
#> [20,]   0.110 0.303 0.375 0.447 0.631 0.918
#> [21,]   0.122 0.314 0.420 0.439 0.640 1.040
#> [22,]   0.078 0.247 0.321 0.387 0.480 0.622
#> [23,]   0.076 0.216 0.325 0.401 0.579 0.758
#> [24,]   0.102 0.335 0.368 0.457 0.604 0.740
#> [25,]   0.139 0.251 0.396 0.466 0.584 0.768
#> [26,]   0.160 0.287 0.367 0.494 0.567 0.726
#> [27,]   0.131 0.309 0.400 0.557 0.629 1.695
#> [28,]   0.185 0.321 0.444 0.561 0.667 0.752
#> [29,]   0.145 0.360 0.419 0.567 0.684 0.824
#> [30,]   0.164 0.330 0.438 0.574 0.764 0.751
#> [31,]   0.095 0.313 0.413 0.572 0.722 0.945
#> [32,]   0.136 0.295 0.436 0.540 0.581 0.799
#> [33,]   0.102 0.295 0.415 0.511 0.634 0.795
#> [34,]   0.110 0.251 0.373 0.475 0.607 0.783
#> [35,]   0.111 0.268 0.363 0.472 0.628 0.834
#> [36,]   0.151 0.266 0.388 0.461 0.574 1.077
#> [37,]   0.105 0.281 0.367 0.502 0.601 0.753
#> [38,]   0.099 0.281 0.414 0.470 0.573 0.702
#> [39,]   0.130 0.280 0.412 0.491 0.572 0.717
#> [40,]   0.120 0.251 0.365 0.475 0.554 0.709
#> [41,]   0.154 0.254 0.351 0.453 0.565 0.683
#> [42,] -99.000 0.345 0.396 0.459 0.565 0.689
#> [43,]   0.046 0.281 0.388 0.476 0.548 0.685
#> [44,]   0.000 0.270 0.349 0.434 0.521 0.629

If we fix this issue with the data file, the model runs fine!