Estuarine Bayesian Single-station Estimation method for ecosystem metabolism

Usage

ebase(
  dat,
  Z,
  interval,
  ndays = 1,
  aprior = c(4, 2),
  rprior = c(300, 150),
  bprior = c(0.251, 0.125),
  bmax = 0.502,
  doave = TRUE,
  maxinterp = 43200/interval,
  n.iter = 10000,
  update.chains = TRUE,
  n.burnin = n.iter * 0.5,
  n.chains = 3,
  n.thin = 10,
  progress = NULL,
  model_file = NULL
)

Arguments

dat: input data frame
Z: numeric as single value for water column depth (m) or vector equal in length to number of rows in dat that will be averaged within ndays
interval: timestep interval in seconds
ndays: numeric for number of days in dat for optimizing the metabolic equation, see details
aprior: numeric vector of length two indicating the mean and standard deviation for the prior distribution of the a parameter, see details
rprior: numeric vector of length two indicating the mean and standard deviation for the prior distribution of the R parameter, see details
bprior: numeric vector of length two indicating the mean and standard deviation for the prior distribution of the b parameter, see details
bmax: numeric value for the upper limit on the prior distribution for bprior, set as twice the default value of the mean
doave: logical indicating if the average dissolved oxygen concentration is used as the starting value for the estimation (default), otherwise the first observation will be used if FALSE, see details
maxinterp: numeric value for minimum number of continuous observations that must not be interpolated within a group defined by ndays to assign as NA in output, see details
n.iter: number of MCMC iterations, passed to jags
update.chains: logical to run metab_update if chains do not converge
n.burnin: number of MCMC chains to delete, passed to jags
n.chains: number of MCMC chains to run, passed to jags
n.thin: number of nth iterations to save for each chain, passed to jags
progress: character string of path where progress is saved to as 'log.txt', use NULL to suppress (default)
model_file: NULL to use the model file included with the package or a path to a model text file can be used

Value

A data frame with metabolic estimates for areal gross production (P, O2 mmol m-2 d-1), respiration (R, O2 mmol m-2 d-1), and gas exchange (D, O2 mmol m-2 d-1, positive values as ingassing, negative values as outgassing). Additional parameters estimated by the model that are returned include a and b. The a parameter is a constant that represents the primary production per quantum of light with units of (mmol m-2 d-1)/(W m-2) and is used to estimate gross production (Grace et al., 2015). The b parameter is a constant used to estimate gas exchange in (cm hr-1)/(m2 s-2) (provided as 0.251 in eqn. 4 in Wanninkhof 2014). Observed dissolved oxygen (DO_obs, mmol m-3), modeled dissolved oxygen (DO_mod, mmol m-3), and delta dissolved oxygen of the modeled results (dDO, mmol m-3 d-1) are also returned. Note that delta dissolved oxygen is a daily rate.

95% credible intervals for a, b, and R are also returned in the corresponding columns alo, ahi, blo, bhi, Rlo, and Rhi, for the 2.5th and 97.5th percentile estimates for each parameter, respectively. These values indicate the interval within which there is a 95% probability that the true parameter is in this range. Note that all values for these parameters are repeated across rows, although only one estimate for each is returned based on the number of days defined by ndays.

Model fit can also be assessed using the converge and rsq columns. The values in these columns apply to each group in the grp column as specified with the ndays argument. The converge column indicates "Check convergence" or "Fine" if the JAGS estimate converged at that iteration (repeated across rows for the group). The n.chains argument can be increased if convergence is not achieved. Similarly, the rsq column shows the r-squared values of the linear fit between the modeled and observed dissolved oxygen (repeated across rows for the group). These values can also be viewed with fit_plot.

Details

Required input data are time series for dissolved oxygen (mg L-1), water temperature (C), salinity (psu), total PAR (W m-2), and wind speed (m s-1). See the exdat example data file for a representation of the required data. Data are typically from continuously monitored water quality and weather parameters are hourly of sub-hourly time steps. Oxygen concentrations are converted to mmol/m3 prior to metabolic estimation. Water column depth is also required. This can be supplied as a single value or a vector of length equal to the number of rows in dat.

The metabolic estimates are based on a mass balance equation in Grace et al. 2015 with the gas exchange estimate from Wanninkhof 2004. It is similar to that provided by the BASEmetab R package at https://github.com/dgiling/BASEmetab, with modifications to estimate different parameters. The equation optimized in the JAGS model is:

$$ Z\frac{dC_d}{dt} = aPAR - R + bU_{10}^2\left(\frac{Sc}{600} \right)^{-0.5} \left(C_{Sat} - C_d \right )$$

More simply:

$$ Z\frac{dC_d}{dt} = P - R + D$$

Gross production is provided by aPAR, respiration is provided by R, and gas exchange is provided by the remainder. The likelihood of the parameters a, R, and b given the observed data are estimated from the JAGS model using prior distributions shown in the model file. At each time step, the change in oxygen concentration between time steps is calculated from the equation using model inputs and parameter guesses, and then a finite difference approximation is used to estimate modeled oxygen concentration. The first modeled value starts at the mean oxygen concentration for all measurements in the optimization period. The estimated concentration at each time step is also returned for comparison to observed as one measure of model performance.

The prior distributions for the a, R, and b parameters are defined in the model file included with the package as normal distributions with mean and standard deviations provided by the aprior, rprior, and bprior arguments. The default values were chosen based on approximate values from national syntheses of metabolic estimates. An additional constraint is that all the prior distributions are truncated to be positive values as required by the core metabolism equation above. The upper limit for b is set as two times 0.251, as given in eqn. 4 in Wanninkhof 2014. Units for each parameter are (mmol m-2 d-1)/(W m-2) for a, mmol m-2 d-1 for R, and (cm hr-1)/(m2 s-2) for b.

The ndays argument defines the model optimization period as the number of days that are used for optimizing the above mass balance equation. By default, this is done each day, i.e., ndays= 1 such that a loop is used that applies the model equation to observations within each day, evaluated iteratively from the first observation in a day to the last. Individual parameter estimates for a, R, and b are then returned for each day. However, more days can be used to estimate the unknown parameters, such that the loop can be evaluated for every ndays specified by the argument. The ndays argument will separate the input data into groups of consecutive days, where each group has a total number of days equal to ndays. The final block may not include the complete number of days specified by ndays if the number of unique dates in the input data includes a remainder when divided by ndays, e.g., if seven days are in the input data and ndays = 5, there will be two groups where the first has five days and the second has two days. The output data from ebase includes a column that specifies the grouping that was used based on ndays.

Missing values in the input data are also interpolated prior to estimating metabolism. It is the responsibility of the user to verify that these interpolated values are not wildly inaccurate. Missing values are linearly interpolated between non-missing values at the time step specified by the value in interval. This works well for small gaps, but can easily create inaccurate values at gaps larger than a few hours. The interp_plot function can be used to visually assess the interpolated values. Records at the start or end of the input time series that do not include a full day are also removed. A warning is returned to the console if gaps are found or dangling records are found.

The maxinterp argument specifies a minimum number of observations that must not be interpolated within groups defined by ndays that are assigned NA in the output (except Date and DateTimeStamp). Groups with continuous rows of interpolated values with length longer than this argument are assigned NA. The default value is half a day, i.e., 43200 seconds divided by the value in interval.

The doave argument can be used to define which dissolved oxygen value is used as the starting point in the Bayesian estimation for the optimization period. The default setting (doave = TRUE) will use the average of all the dissolved oxygen values in the optimization period defined by ndays. For example, the average of all dissolved oxygen values in each 24 hour period will be used if doave = TRUE and ndays = 1. The first dissolved oxygen observation of the time series in the optimization period will be used as the starting point if doave = F. In this case, the simulated dissolved oxygen time series will always start at the first observed dissolved oxygen value for each optimization period.

References

Grace, M.R., Giling, D.P., Hladyz, S., Caron, V., Thompson, R.M., Nally, R.M., 2015. Fast processing of diel oxygen curves: Estimating stream metabolism with BASE (BAyesian Single-station Estimation). Limnology and Oceanography: Methods 13, e10011. https://doi.org/10.1002/lom3.10011

Wanninkhof, R., 2014. Relationship between wind speed and gas exchange over the ocean revisited. Limnology and Oceanography: Methods 12, 351–362. https://doi.org/10.4319/lom.2014.12.351

Examples

# get one day of data
dat <- exdat[as.Date(exdat$DateTimeStamp, tz = 'America/Jamaica') == as.Date('2012-06-01'), ]

# run ebase, use more chains and iterations for a better fit, update.chains as T
ebase(dat, interval = 900, Z = 1.85, n.chains = 2, n.iter = 50, 
 update.chains = FALSE)
#> Warning: executing %dopar% sequentially: no parallel backend registered
#> Loading required package: rjags
#> Loading required package: coda
#> Linked to JAGS 4.3.0
#> Loaded modules: basemod,bugs
#> 
#> Attaching package: ‘R2jags’
#> The following object is masked from ‘package:coda’:
#> 
#>     traceplot
#> 
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:stats’:
#> 
#>     filter, lag
#> The following objects are masked from ‘package:base’:
#> 
#>     intersect, setdiff, setequal, union
#>          DateTimeStamp       Date grp    Z  DO_obs   DO_mod DO_modlo DO_modhi
#> 1  2012-06-01 00:00:00 2012-06-01   1 1.85 140.625 166.1132 166.1129 166.1136
#> 2  2012-06-01 00:15:00 2012-06-01   1 1.85 159.375 165.3936 165.3218 165.5435
#> 3  2012-06-01 00:30:00 2012-06-01   1 1.85 150.000 164.7099 164.5632 165.0197
#> 4  2012-06-01 00:45:00 2012-06-01   1 1.85 150.000 164.0038 163.7843 164.4666
#> 5  2012-06-01 01:00:00 2012-06-01   1 1.85 150.000 163.2883 162.9970 163.9011
#> 6  2012-06-01 01:15:00 2012-06-01   1 1.85 156.250 162.5434 162.1827 163.2977
#> 7  2012-06-01 01:30:00 2012-06-01   1 1.85 168.750 161.8078 161.3771 162.7057
#> 8  2012-06-01 01:45:00 2012-06-01   1 1.85 168.750 161.0710 160.5706 162.1118
#> 9  2012-06-01 02:00:00 2012-06-01   1 1.85 165.625 160.3364 159.7662 161.5203
#> 10 2012-06-01 02:15:00 2012-06-01   1 1.85 165.625 159.6248 158.9833 160.9575
#> 11 2012-06-01 02:30:00 2012-06-01   1 1.85 150.000 158.8925 158.1814 160.3683
#> 12 2012-06-01 02:45:00 2012-06-01   1 1.85 159.375 158.2529 157.4661 159.8949
#> 13 2012-06-01 03:00:00 2012-06-01   1 1.85 131.250 157.5133 156.6576 159.2960
#> 14 2012-06-01 03:15:00 2012-06-01   1 1.85 118.750 156.7307 155.8090 158.6431
#> 15 2012-06-01 03:30:00 2012-06-01   1 1.85 121.875 155.9822 154.9924 158.0327
#> 16 2012-06-01 03:45:00 2012-06-01   1 1.85 115.625 155.2356 154.1777 157.4242
#> 17 2012-06-01 04:00:00 2012-06-01   1 1.85 118.750 154.5369 153.4080 156.8748
#> 18 2012-06-01 04:15:00 2012-06-01   1 1.85 131.250 153.7836 152.5871 156.2577
#> 19 2012-06-01 04:30:00 2012-06-01   1 1.85 156.250 153.1179 151.8488 155.7485
#> 20 2012-06-01 04:45:00 2012-06-01   1 1.85 162.500 152.4392 151.0983 155.2228
#> 21 2012-06-01 05:00:00 2012-06-01   1 1.85 134.375 151.7478 150.3361 154.6812
#> 22 2012-06-01 05:15:00 2012-06-01   1 1.85 150.000 151.0579 149.5753 154.1409
#> 23 2012-06-01 05:30:00 2012-06-01   1 1.85 150.000 150.6449 149.0761 153.9395
#> 24 2012-06-01 05:45:00 2012-06-01   1 1.85 146.875 150.1249 148.4763 153.6063
#> 25 2012-06-01 06:00:00 2012-06-01   1 1.85 150.000 149.4519 147.7335 153.0814
#> 26 2012-06-01 06:15:00 2012-06-01   1 1.85 150.000 148.7726 146.9885 152.5384
#> 27 2012-06-01 06:30:00 2012-06-01   1 1.85 137.500 148.1076 146.2624 151.9981
#> 28 2012-06-01 06:45:00 2012-06-01   1 1.85 140.625 147.4427 145.5408 151.4450
#> 29 2012-06-01 07:00:00 2012-06-01   1 1.85 150.000 146.8061 144.8514 150.9111
#> 30 2012-06-01 07:15:00 2012-06-01   1 1.85 153.125 146.2361 144.2331 150.4362
#> 31 2012-06-01 07:30:00 2012-06-01   1 1.85 159.375 145.9501 143.9422 150.1419
#> 32 2012-06-01 07:45:00 2012-06-01   1 1.85 168.750 146.1791 144.1881 150.3371
#> 33 2012-06-01 08:00:00 2012-06-01   1 1.85 168.750 146.6068 144.6472 150.7043
#> 34 2012-06-01 08:15:00 2012-06-01   1 1.85 168.750 146.5963 144.5981 150.8008
#> 35 2012-06-01 08:30:00 2012-06-01   1 1.85 162.500 147.0069 145.0675 151.0652
#> 36 2012-06-01 08:45:00 2012-06-01   1 1.85 165.625 147.5823 145.7171 151.4481
#> 37 2012-06-01 09:00:00 2012-06-01   1 1.85 168.750 148.4424 146.6538 152.1146
#> 38 2012-06-01 09:15:00 2012-06-01   1 1.85 171.875 149.4843 147.7819 152.9427
#> 39 2012-06-01 09:30:00 2012-06-01   1 1.85 171.875 151.3225 149.5145 154.7566
#> 40 2012-06-01 09:45:00 2012-06-01   1 1.85 168.750 153.3085 151.0880 156.7725
#> 41 2012-06-01 10:00:00 2012-06-01   1 1.85 168.750 153.5644 151.0514 157.3022
#> 42 2012-06-01 10:15:00 2012-06-01   1 1.85 175.000 153.7757 151.1354 157.5760
#> 43 2012-06-01 10:30:00 2012-06-01   1 1.85 178.125 154.1029 151.4212 157.8242
#> 44 2012-06-01 10:45:00 2012-06-01   1 1.85 178.125 155.4090 152.7609 158.7264
#> 45 2012-06-01 11:00:00 2012-06-01   1 1.85 178.125 156.6951 154.1425 159.5303
#> 46 2012-06-01 11:15:00 2012-06-01   1 1.85 175.000 158.0054 155.5798 160.3096
#> 47 2012-06-01 11:30:00 2012-06-01   1 1.85 162.500 159.3459 157.0527 161.1048
#> 48 2012-06-01 11:45:00 2012-06-01   1 1.85 165.625 160.7278 158.5634 162.4647
#> 49 2012-06-01 12:00:00 2012-06-01   1 1.85 159.375 162.2002 160.1511 164.1497
#> 50 2012-06-01 12:15:00 2012-06-01   1 1.85 159.375 164.0740 162.0282 166.2938
#> 51 2012-06-01 12:30:00 2012-06-01   1 1.85 162.500 165.6557 163.5179 168.1274
#> 52 2012-06-01 12:45:00 2012-06-01   1 1.85 165.625 167.2160 164.9953 169.9521
#> 53 2012-06-01 13:00:00 2012-06-01   1 1.85 153.125 168.7745 166.5460 171.7842
#> 54 2012-06-01 13:15:00 2012-06-01   1 1.85 140.625 169.6355 167.3747 172.8035
#> 55 2012-06-01 13:30:00 2012-06-01   1 1.85 153.125 171.0881 168.7735 174.4969
#> 56 2012-06-01 13:45:00 2012-06-01   1 1.85 153.125 172.4370 170.0801 176.0689
#> 57 2012-06-01 14:00:00 2012-06-01   1 1.85 146.875 174.0277 171.4462 177.9162
#> 58 2012-06-01 14:15:00 2012-06-01   1 1.85 143.750 175.6070 172.6233 179.7471
#> 59 2012-06-01 14:30:00 2012-06-01   1 1.85 131.250 177.2261 173.7999 181.6251
#> 60 2012-06-01 14:45:00 2012-06-01   1 1.85 168.750 178.4056 174.7342 182.9782
#> 61 2012-06-01 15:00:00 2012-06-01   1 1.85 175.000 179.2437 175.4161 183.9356
#> 62 2012-06-01 15:15:00 2012-06-01   1 1.85 150.000 180.4956 176.4630 185.3485
#> 63 2012-06-01 15:30:00 2012-06-01   1 1.85 190.625 182.2572 177.9473 187.3224
#> 64 2012-06-01 15:45:00 2012-06-01   1 1.85 181.250 183.6345 179.0203 188.8803
#> 65 2012-06-01 16:00:00 2012-06-01   1 1.85 200.000 184.6039 179.8324 189.9584
#> 66 2012-06-01 16:15:00 2012-06-01   1 1.85 200.000 185.6548 180.7627 191.1090
#> 67 2012-06-01 16:30:00 2012-06-01   1 1.85 196.875 186.5095 181.4785 192.0541
#> 68 2012-06-01 16:45:00 2012-06-01   1 1.85 203.125 186.9395 181.8541 192.5207
#> 69 2012-06-01 17:00:00 2012-06-01   1 1.85 206.250 187.4098 182.2662 193.0298
#> 70 2012-06-01 17:15:00 2012-06-01   1 1.85 206.250 188.1033 182.8536 193.7891
#> 71 2012-06-01 17:30:00 2012-06-01   1 1.85 206.250 188.6004 183.2891 194.3237
#> 72 2012-06-01 17:45:00 2012-06-01   1 1.85 209.375 188.6085 183.3209 194.3162
#> 73 2012-06-01 18:00:00 2012-06-01   1 1.85 209.375 188.5256 183.2893 194.2035
#> 74 2012-06-01 18:15:00 2012-06-01   1 1.85 206.250 188.6583 183.4488 194.3237
#> 75 2012-06-01 18:30:00 2012-06-01   1 1.85 206.250 188.5367 183.3447 194.1819
#> 76 2012-06-01 18:45:00 2012-06-01   1 1.85 203.125 188.1723 183.0367 193.7743
#> 77 2012-06-01 19:00:00 2012-06-01   1 1.85 203.125 187.7099 182.6758 193.2488
#> 78 2012-06-01 19:15:00 2012-06-01   1 1.85 203.125 187.1669 182.2310 192.6415
#> 79 2012-06-01 19:30:00 2012-06-01   1 1.85 203.125 186.5392 181.7460 191.9323
#> 80 2012-06-01 19:45:00 2012-06-01   1 1.85 196.875 185.8227 181.1773 191.1287
#> 81 2012-06-01 20:00:00 2012-06-01   1 1.85 190.625 185.0886 180.6043 190.3032
#> 82 2012-06-01 20:15:00 2012-06-01   1 1.85 193.750 184.3134 179.9767 189.4386
#> 83 2012-06-01 20:30:00 2012-06-01   1 1.85 196.875 183.5295 179.3348 188.5668
#> 84 2012-06-01 20:45:00 2012-06-01   1 1.85 184.375 182.7656 178.7223 187.7144
#> 85 2012-06-01 21:00:00 2012-06-01   1 1.85 175.000 182.0113 178.1231 186.8720
#> 86 2012-06-01 21:15:00 2012-06-01   1 1.85 175.000 181.2682 177.5392 186.0415
#> 87 2012-06-01 21:30:00 2012-06-01   1 1.85 171.875 180.5434 176.9812 185.2295
#> 88 2012-06-01 21:45:00 2012-06-01   1 1.85 103.125 179.8355 176.4470 184.4340
#> 89 2012-06-01 22:00:00 2012-06-01   1 1.85 150.000 179.1410 175.9309 183.6531
#> 90 2012-06-01 22:15:00 2012-06-01   1 1.85 162.500 178.4413 175.4041 182.8700
#> 91 2012-06-01 22:30:00 2012-06-01   1 1.85 131.250 177.7537 174.8948 182.0980
#> 92 2012-06-01 22:45:00 2012-06-01   1 1.85 156.250 177.0818 174.2971 181.3448
#> 93 2012-06-01 23:00:00 2012-06-01   1 1.85 181.250 176.3702 173.5784 180.5506
#> 94 2012-06-01 23:15:00 2012-06-01   1 1.85 184.375 175.6662 172.8639 179.7671
#> 95 2012-06-01 23:30:00 2012-06-01   1 1.85 175.000 174.9048 172.1102 178.9239
#> 96 2012-06-01 23:45:00 2012-06-01   1 1.85 162.500 174.1385 171.3528 178.0762
#>            dDO          converge       rsq         a       alo     ahi
#> 1           NA Check convergence 0.3337845        NA        NA      NA
#> 2  -69.0843195 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 3  -65.6339597 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 4  -67.7906180 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 5  -68.6816694 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 6  -71.5085041 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 7  -70.6179486 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 8  -70.7346970 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 9  -70.5248981 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 10 -68.3166482 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 11 -70.2986658 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 12 -61.3973806 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 13 -70.9988852 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 14 -75.1372681 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 15 -71.8481405 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 16 -71.6805952 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 17 -67.0751906 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 18 -72.3177654 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 19 -63.8998122 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 20 -65.1627638 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 21 -66.3661385 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 22 -66.2376453 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 23 -39.6492267 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 24 -49.9160980 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 25 -64.6110604 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 26 -65.2111881 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 27 -63.8354209 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 28 -63.8297456 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 29 -61.1170448 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 30 -54.7185449 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 31 -27.4601724 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 32  21.9861213 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 33  41.0565129 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 34  -1.0084103 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 35  39.4180148 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 36  55.2452710 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 37  82.5640792 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 38 100.0276378 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 39 176.4594406 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 40 190.6617040 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 41  24.5619357 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 42  20.2839792 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 43  31.4162646 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 44 125.3832655 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 45 123.4710070 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 46 125.7855831 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 47 128.6853792 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 48 132.6611765 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 49 141.3545197 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 50 179.8836518 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 51 151.8472491 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 52 149.7849452 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 53 149.6136261 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 54  82.6613664 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 55 139.4500937 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 56 129.4867233 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 57 152.7067891 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 58 151.6199004 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 59 155.4263160 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 60 113.2370317 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 61  80.4525927 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 62 120.1847577 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 63 169.1141198 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 64 132.2180272 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 65  93.0620388 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 66 100.8938990 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 67  82.0504339 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 68  41.2761192 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 69  45.1544122 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 70  66.5738715 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 71  47.7219609 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 72   0.7781609 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 73  -7.9638485 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 74  12.7371630 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 75 -11.6712338 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 76 -34.9805048 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 77 -44.3935632 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 78 -52.1281917 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 79 -60.2603594 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 80 -68.7797421 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 81 -70.4720725 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 82 -74.4232741 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 83 -75.2544767 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 84 -73.3332098 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 85 -72.4146708 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 86 -71.3311257 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 87 -69.5793019 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 88 -67.9670364 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 89 -66.6641178 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 90 -67.1715898 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 91 -66.0129340 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 92 -64.5036942 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 93 -68.3103327 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 94 -67.5862398 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 95 -73.0925296 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#> 96 -73.5703052 Check convergence 0.3337845 0.9122054 0.6136367 1.05631
#>            b       blo       bhi          P         Plo        Phi       R
#> 1         NA        NA        NA         NA          NA         NA      NA
#> 2  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 3  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 4  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 5  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 6  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 7  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 8  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 9  0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 10 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 11 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 12 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 13 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 14 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 15 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 16 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 17 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 18 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 19 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 20 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 21 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 22 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 23 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 24 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 25 0.3482113 0.2079293 0.4644247   1.366788   0.9194324   1.582704 146.991
#> 26 0.3482113 0.2079293 0.4644247   4.761712   3.2031837   5.513937 146.991
#> 27 0.3482113 0.2079293 0.4644247   9.589559   6.4508562  11.104456 146.991
#> 28 0.3482113 0.2079293 0.4644247  13.778103   9.2684715  15.954679 146.991
#> 29 0.3482113 0.2079293 0.4644247  18.760264  12.6199508  21.723891 146.991
#> 30 0.3482113 0.2079293 0.4644247  25.968968  17.4692151  30.071378 146.991
#> 31 0.3482113 0.2079293 0.4644247  80.596389  54.2168508  93.328489 146.991
#> 32 0.3482113 0.2079293 0.4644247 125.259486  84.2615280 145.047175 146.991
#> 33 0.3482113 0.2079293 0.4644247 148.164204  99.6694350 171.570233 146.991
#> 34 0.3482113 0.2079293 0.4644247  61.329090  41.2558203  71.017466 146.991
#> 35 0.3482113 0.2079293 0.4644247 177.175377 119.1851286 205.164405 146.991
#> 36 0.3482113 0.2079293 0.4644247 207.200618 139.3829817 239.932841 146.991
#> 37 0.3482113 0.2079293 0.4644247 230.854864 155.2950936 267.323833 146.991
#> 38 0.3482113 0.2079293 0.4644247 252.106210 169.5907840 291.932330 146.991
#> 39 0.3482113 0.2079293 0.4644247 259.336958 174.4548778 300.305345 146.991
#> 40 0.3482113 0.2079293 0.4644247 254.685471 171.3258418 294.919045 146.991
#> 41 0.3482113 0.2079293 0.4644247  37.873248  25.4771744  43.856221 146.991
#> 42 0.3482113 0.2079293 0.4644247  98.629170  66.3474263 114.209972 146.991
#> 43 0.3482113 0.2079293 0.4644247 151.845713 102.1459706 175.833323 146.991
#> 44 0.3482113 0.2079293 0.4644247 331.842845 223.2292822 384.265246 146.991
#> 45 0.3482113 0.2079293 0.4644247 354.659383 238.5778710 410.686194 146.991
#> 46 0.3482113 0.2079293 0.4644247 371.589915 249.9669688 430.291303 146.991
#> 47 0.3482113 0.2079293 0.4644247 378.445899 254.5789602 438.230352 146.991
#> 48 0.3482113 0.2079293 0.4644247 383.317836 257.8562917 443.871926 146.991
#> 49 0.3482113 0.2079293 0.4644247 391.430383 263.3135677 453.266041 146.991
#> 50 0.3482113 0.2079293 0.4644247 406.310734 273.3235169 470.497094 146.991
#> 51 0.3482113 0.2079293 0.4644247 329.484034 221.6425199 381.533805 146.991
#> 52 0.3482113 0.2079293 0.4644247 338.147705 227.4705348 391.566107 146.991
#> 53 0.3482113 0.2079293 0.4644247 375.006885 252.2655497 434.248064 146.991
#> 54 0.3482113 0.2079293 0.4644247 247.432677 166.4469185 286.520502 146.991
#> 55 0.3482113 0.2079293 0.4644247 333.826892 224.5639421 386.562720 146.991
#> 56 0.3482113 0.2079293 0.4644247 322.143061 216.7042783 373.033152 146.991
#> 57 0.3482113 0.2079293 0.4644247 364.270987 245.0435567 421.816178 146.991
#> 58 0.3482113 0.2079293 0.4644247 365.086651 245.5922502 422.760695 146.991
#> 59 0.3482113 0.2079293 0.4644247 380.429946 255.9136201 440.527825 146.991
#> 60 0.3482113 0.2079293 0.4644247 290.949437 195.7204588 336.911760 146.991
#> 61 0.3482113 0.2079293 0.4644247 236.740870 159.2545846 274.139672 146.991
#> 62 0.3482113 0.2079293 0.4644247 292.558719 196.8030163 338.775266 146.991
#> 63 0.3482113 0.2079293 0.4644247 365.086651 245.5922502 422.760695 146.991
#> 64 0.3482113 0.2079293 0.4644247 332.724644 223.8224644 385.286346 146.991
#> 65 0.3482113 0.2079293 0.4644247 258.411070 173.8320365 299.233191 146.991
#> 66 0.3482113 0.2079293 0.4644247 261.453275 175.8785150 302.755984 146.991
#> 67 0.3482113 0.2079293 0.4644247 243.817303 164.0148715 282.333995 146.991
#> 68 0.3482113 0.2079293 0.4644247 179.732593 120.9053569 208.125593 146.991
#> 69 0.3482113 0.2079293 0.4644247 185.596553 124.8500184 214.915904 146.991
#> 70 0.3482113 0.2079293 0.4644247 220.846451 148.5624759 255.734354 146.991
#> 71 0.3482113 0.2079293 0.4644247 190.755075 128.3201341 220.889336 146.991
#> 72 0.3482113 0.2079293 0.4644247 117.455569  79.0118658 136.010445 146.991
#> 73 0.3482113 0.2079293 0.4644247 101.693420  68.4087344 117.758293 146.991
#> 74 0.3482113 0.2079293 0.4644247 131.983200  88.7845421 152.833058 146.991
#> 75 0.3482113 0.2079293 0.4644247 103.434972  69.5802692 119.774964 146.991
#> 76 0.3482113 0.2079293 0.4644247  66.972601  45.0521862  77.552502 146.991
#> 77 0.3482113 0.2079293 0.4644247  46.933729  31.5721213  54.348018 146.991
#> 78 0.3482113 0.2079293 0.4644247  37.586664  25.2843902  43.524364 146.991
#> 79 0.3482113 0.2079293 0.4644247  18.583905  12.5013143  21.519671 146.991
#> 80 0.3482113 0.2079293 0.4644247   7.076433   4.7602870   8.194323 146.991
#> 81 0.3482113 0.2079293 0.4644247   2.072227   1.3939781   2.399584 146.991
#> 82 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 83 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 84 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 85 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 86 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 87 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 88 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 89 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 90 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 91 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 92 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 93 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 94 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 95 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#> 96 0.3482113 0.2079293 0.4644247   0.000000   0.0000000   0.000000 146.991
#>         Rlo      Rhi          D        Dlo        Dhi
#> 1        NA       NA         NA         NA         NA
#> 2  125.7794 161.9538  19.184960  11.456023  25.587969
#> 3  125.7794 161.9538  25.568126  15.283144  33.999484
#> 4  125.7794 161.9538  21.578308  12.916025  28.604481
#> 5  125.7794 161.9538  19.929863  11.942181  26.343469
#> 6  125.7794 161.9538  14.700219   8.816695  19.378630
#> 7  125.7794 161.9538  16.347746   9.811919  21.489761
#> 8  125.7794 161.9538  16.131762   9.690742  21.141346
#> 9  125.7794 161.9538  16.519890   9.930138  21.599063
#> 10 125.7794 161.9538  20.605152  12.394580  26.872113
#> 11 125.7794 161.9538  16.938420  10.197978  22.035284
#> 12 125.7794 161.9538  33.405797  20.123816  43.368932
#> 13 125.7794 161.9538  15.643014   9.437872  20.253658
#> 14 125.7794 161.9538   7.987005   4.820836  10.322266
#> 15 125.7794 161.9538  14.071891   8.496525  18.138899
#> 16 125.7794 161.9538  14.381850   8.686641  18.506382
#> 17 125.7794 161.9538  22.901849  13.837814  29.418221
#> 18 125.7794 161.9538  13.203085   7.983243  16.931709
#> 19 125.7794 161.9538  28.776299  17.404817  36.840966
#> 20 125.7794 161.9538  26.439838  16.006725  33.792989
#> 21 125.7794 161.9538  24.213595  14.673073  30.882855
#> 22 125.7794 161.9538  24.451308  14.829497  31.123878
#> 23 125.7794 161.9538  73.639882  44.676027  93.699632
#> 24 125.7794 161.9538  54.646170  33.246041  69.359798
#> 25 125.7794 161.9538  26.093702  15.908954  33.031602
#> 26 125.7794 161.9538  21.588541  13.169283  27.296858
#> 27 125.7794 161.9538  19.305863  11.784646  24.364402
#> 28 125.7794 161.9538  15.127819   9.238392  19.063037
#> 29 125.7794 161.9538  15.164154   9.263495  19.079918
#> 30 125.7794 161.9538  19.792676  12.090823  24.888746
#> 31 125.7794 161.9538  15.593244   9.526806  19.597655
#> 32 125.7794 161.9538  62.405789  38.120675  78.445309
#> 33 125.7794 161.9538  74.781296  45.751611  93.996967
#> 34 125.7794 161.9538  83.796302  51.353577 105.416451
#> 35 125.7794 161.9538  42.738902  26.263299  53.681022
#> 36 125.7794 161.9538  41.994085  25.823159  52.801991
#> 37 125.7794 161.9538  68.879633  42.348794  86.877958
#> 38 125.7794 161.9538  79.935871  49.225248 100.988521
#> 39 125.7794 161.9538 214.103958 132.095243 271.108231
#> 40 125.7794 161.9538 245.029633 152.360843 310.205200
#> 41 125.7794 161.9538 154.557284  97.067439 195.372907
#> 42 125.7794 161.9538  85.887143  54.273152 108.132300
#> 43 125.7794 161.9538  53.265328  33.751215  67.003071
#> 44 125.7794 161.9538  47.107147  29.851761  59.387787
#> 45 125.7794 161.9538  20.752931  13.168159  26.291562
#> 46 125.7794 161.9538   8.104365   5.135319  10.346373
#> 47 125.7794 161.9538   6.613003   4.182866   8.514473
#> 48 125.7794 161.9538   9.096291   5.738446  11.823169
#> 49 125.7794 161.9538  17.066430  10.740801  22.398956
#> 50 125.7794 161.9538  73.464973  46.161010  97.392591
#> 51 125.7794 161.9538  98.424328  61.942567 131.722394
#> 52 125.7794 161.9538  85.945395  54.299460 115.769957
#> 53 125.7794 161.9538  48.769275  30.899989  66.192861
#> 54 125.7794 161.9538  52.481802  33.268575  72.016358
#> 55 125.7794 161.9538  71.146733  45.109660  98.130340
#> 56 125.7794 161.9538  64.398329  40.929472  89.641764
#> 57 125.7794 161.9538  65.227524  41.528765  91.629027
#> 58 125.7794 161.9538  62.401116  39.790209  88.675748
#> 59 125.7794 161.9538  54.099690  34.545765  77.833447
#> 60 125.7794 161.9538  65.530023  41.886180  95.666380
#> 61 125.7794 161.9538  59.087378  37.934942  87.301941
#> 62 125.7794 161.9538  76.774034  49.468642 114.402700
#> 63 125.7794 161.9538  94.765422  61.363951 142.817965
#> 64 125.7794 161.9538  58.869658  38.574796  90.845819
#> 65 125.7794 161.9538  60.744653  39.785413  94.793223
#> 66 125.7794 161.9538  72.191390  47.380240 113.555564
#> 67 125.7794 161.9538  54.966951  36.241014  87.264128
#> 68 125.7794 161.9538  43.619179  28.908929  70.173104
#> 69 125.7794 161.9538  44.930061  29.825451  72.592116
#> 70 125.7794 161.9538  49.306163  32.716146  79.747001
#> 71 125.7794 161.9538  44.521504  29.706525  72.977157
#> 72 125.7794 161.9538  30.974980  20.735393  51.170400
#> 73 125.7794 161.9538  30.564412  20.493345  50.568181
#> 74 125.7794 161.9538  38.571503  25.897213  63.802679
#> 75 125.7794 161.9538  21.964197  14.815297  36.648980
#> 76 125.7794 161.9538  15.304416  10.295917  25.376820
#> 77 125.7794 161.9538  17.929131  12.037390  29.547916
#> 78 125.7794 161.9538  12.967133   8.670707  21.137188
#> 79 125.7794 161.9538  16.925382  11.192232  27.002297
#> 80 125.7794 161.9538  12.671995   8.361338  20.022435
#> 81 125.7794 161.9538  14.545390   9.559475  22.708354
#> 82 125.7794 161.9538   9.307894   6.112317  14.419680
#> 83 125.7794 161.9538   7.770169   5.086935  11.921037
#> 84 125.7794 161.9538  11.324513   7.383326  17.191631
#> 85 125.7794 161.9538  13.023810   8.468346  19.595823
#> 86 125.7794 161.9538  15.028369   9.739557  22.394811
#> 87 125.7794 161.9538  18.269243  11.831664  27.053704
#> 88 125.7794 161.9538  21.251934  13.788230  31.355759
#> 89 125.7794 161.9538  23.662333  15.335765  34.655631
#> 90 125.7794 161.9538  22.723510  14.619793  32.821818
#> 91 125.7794 161.9538  24.867023  16.099940  35.944534
#> 92 125.7794 161.9538  27.659117  17.794560  39.518071
#> 93 125.7794 161.9538  20.616836  13.345950  29.445360
#> 94 125.7794 161.9538  21.956408  14.069242  30.956977
#> 95 125.7794 161.9538  11.769772   7.562459  16.555044
#> 96 125.7794 161.9538  10.885887   6.975266  15.230043