dataesr/scientific-paragraphs-categorization huggingface dataset:
Theophylline pharmacokinetics The data: this is again a well-known dataset, distributed in NONMEM, Monolix and R, in nonlinear mixed effect models, and is a classic example of pharmacokinetic data. The data was collected in a study by Upton of the kinetics of the anti-asthmatic drug theophylline and analysed by  Boeckmann et al. (1994). Twelve subjects were given oral doses of the anti-asthmatic drug theophylline, then serum concentrations (in mg/L) were measured at 11 time points over the next 25 hours. In package saemix, we removed the data at time 0 to avoid some unexplained non-zero values in a supposedly single-dose study; in addition we transformed the doses to body doses in mg instead of doses by weight mg/kg as in the original dataset, and we added a dummy variable for gender with values 0 and 1 to illustrate fitting categorical covariate model.
The resulting dataset is available under the name theo. saemix in the saemix package. The data is shown in Figure  4. In the following code, we use the dataset from the library to create the saemixData object: R> data("theo. saemix") R> theo. data <-saemixData(name. data = theo. saemix, header = TRUE, + sep = " ", na = NA, name. group = c("Id"), + name. predictors = c("Dose", "Time"), name. response = c("Concentration"), + name. covariates = c("Weight", "Sex"), + units = list(x = "hr", y = "mg/L", covariates = c("kg", "-")), + name. X = "Time") R> plot(theo. data,type = "b", col = "DarkRed", main = "Theophylline data") The data includes 2 predictors, Dose and Time; we use the name. X = "Time" argument to specify which of the predictor will be used as the independent variable in the plots. We read the two covariates in the data object, and we specify units for the independent and dependent variables, as well as for the covariates.
model <-saemixModel(model = model1cpt, + description = "One-compartment model with first-order absorption", + psi0 = matrix(c(1, 20, 0. 5, 0. 1, 0, -0. 01), ncol = 3, byrow = TRUE, + dimnames = list(NULL, c("ka", "V", "CL"))), transform. par = c(1, 1, 1), + covariate. model = matrix(c(0, 0, 1, 0, 0, 0), ncol = 3, byrow = TRUE)) The covariate model is specified as a matrix, and is reported in the model object as: Covariate model: ka V CL [1,] 0 0 1 [2,] 0 0 0 This indicates that the first covariate in the dataset (Weight in our example) will be assumed to have a relationship with parameter CL. In the code, psi0 now includes a second line, which is used to specify initial values for the parameter-covariate relationships. The same initial values would be used for all covariates, but we could specify different starting values for each parameter-covariate relationship by adding more lines to the matrix psi0.
The penalty using n tot is implemented in the R package nlme and SPSS procedure MIXED, while the number of subject, N , is used in Monolix and SAS proc NLMIXED. An optimal BIC penalty based on two terms proportional to log N and log n tot and that is consistent with the random effect structure in a mixed effects model was recently proposed by  Delattre, Poursat, and Lavielle (2014). Other criteria have been proposed and studied for nonlinear mixed effect models (see review in  Bertrand, Comets, Chenel, and Mentré (2012) ); they can be computed from the estimate of the log-likelihood, which can be accessed through elements of the results object. For instance, the log-likelihood computed by importance sampling can be obtained as: R> logLik(theo. fit) LL by is "-172. 41 (df=8)" We should also note that a conditional AIC (cAIC) has been proposed by  Vaida and Blanchard (2005)  for linear mixed effects models.
fit, plot. type = "observations. vs. predictions") R> plot(saemix. fit, plot. type = "npde") R> plot(saemix. fit, plot. type = "vpc") Overall, these plots show good model adequacy, and individual profiles indicates that the model is able to reproduce observed data reasonably well. Some points to note may be the departure of the npde from the theoretical normal distribution seen in Figure  7 , as well as a slight tendency in the plots of npde versus time. The VPC plot in Figure  8  shows that the model predicts the median evolution very well; the variability may be slightly overestimated but with only 12 subjects in the dataset this may be expected. Corresponding code for nlme and nlmer: again, we can cast the same model for use in nlme and nlmer. For nlme, the relationship between parameters and covariates is set in the fixed argument within the function call: R> startvec.
Evaluating the performances of saemix Data and model: In this section, we use simulated data to evaluate the performance of saemix. The simulated datasets were generated by  Plan et al. (2012)  in a paper comparing different software to estimate parameters in nonlinear mixed effect models, and were obtained directly from the authors of that paper. The link is provided in Acknowledgments. The model used for this simulation is a sigmoid E max model, a standard model in doseresponse studies where the effect of a drug in response to a dose d, E(d), is a sigmoid function corresponding to the following equation: E(d) = E 0 + E max d γ i d γ i + ED γ 50 ( 12 ) The saemix package This model involves 4 parameters, the initial effect E 0 , the maximum effect E max , the concentration at which half the maximum effect is achieved EC 50 and the sigmoidicity factor γ which controls the nonlinearity of the model through the curvature.
1. We also used either the default settings of the algorithm, or tuned the settings by increasing the number of chains to 5 and increasing the number of iterations to 300 and 150 respectively in the burn-in and convergence phases of the algorithm. Combining starting values and tuning parameters, we therefore performed 4 successive parameter estimations for each dataset. In all settings, we used the same random seed for all runs. Evaluation: the performance of both algorithms was evaluated by computing the relative bias and root mean square error (RMSE) over the 100 simulations, using the simulation parameters as true values. Denoting Ψ k 0 the true value of parameter Ψ k , and Ψk l the value estimated in the lth simulated dataset, the relative bias for this parameter was given by: RB(Ψ k ) = 100 l=0 Ψk l -Ψ k 0 Ψ k 0 ( 13 ) and the relative RMSE by: RRM SE(Ψ k ) = 100 l=0 ( Ψk l -Ψ k 0 ) 2 Ψ k 2 0 (14) Results: tables 2 and 3 shows the relative bias and root mean square error obtained for both methods in the two scenarios, in the four settings.
In the rich sampling scenario, all parameters are estimated without bias (less than 5% relative bias) irrespective of whether initial parameters were set to the true value or not, or whether the algorithm was tuned. The relative root mean square errors are less than 20% for the fixed effects, and between 20 and 50% for the variability parameters. Figure  9  shows a boxplot of the biases observed in the 100 datasets, when running the algorithm with the default settings and false starting values: for most of the datasets and all parameters except the covariance the bias is less than 20% in the majority of the datasets.
Comparing saemix to nlme/nlmer Simulation study: the datasets from section 4. 3 have been fit using nlme by  Plan et al. (2012) , who report convergence for only 5% (rich design) and 16% (sparse design) of the datasets. We had the same issue when trying to fit the same data (results not shown). We attempted to use nlmer from lme4 and had the same convergence problems; with nlmer we also had a technical issue, in that the derivative with respect to γ is not defined for a dose of 0, which required some tinkering with the data to bypass. Because we could fit so few datasets, we could not compare meaningfully the performance of saemix to those of the other two modelling packages. In  (Plan et al. (2012) ), the rate of convergence decreased as γ increased, suggesting it was related to model nonlinearity.
In this section, we therefore decided to use the datasets simulated with γ = 1 and we fit them using an E max model corresponding to Equation 12 with γ = 1 (fixed, not estimated). The datasets for the sparse design were derived from the datasets with a rich design (scenario R1A) by using for each dataset the same sampling times as in scenario S3A analysed in the previous section. Initial conditions for nlmer were set to the true values used in the simulation, given in Table  1 , as for saemix, except for cov(E 0 , E max ) and σ which could not be initialised in the other packages. As we encountered convergence issues again, we changed the initial parameter estimates within 20% of the original values; this was repeated a maximum of 10 times per dataset. Results: tables 4 and 5 shows the bias and RMSE for the parameters estimated by each of three packages for the rich and sparse scenarios respectively.
Computation times We performed two series of simulations to evaluate the scalability of the algorithm. In the first, we investigated the effect of the size of the datasets by simulating datasets using the E max model from section 4. 3, and varying the number of samples and the number of subjects. In the second, we assessed the effect of increasing model complexity on the runtimes. For this second simulation, we simulated three models with one, two and three exponential terms (respectively 2, 4 and 6 parameters). Each dataset was simulated assuming 100 subjects and 6 sample times, which had been optimised by the PFIM software  (Retout, Comets, Samson, and Mentré (2007) ) to ensure reasonable identifiability of the most complex model. The model with two exponentials is similar to the E max model in terms of model complexity.
In each simulation, 100 datasets were generated, and the runtimes for the fit of the 100 datasets were collected through the Rprof() function in R, and summarised. Table  6  shows the total run times in seconds reported by Rprof(); the same information is displayed in Figure  11 , with the graph on the left showing the evolution of run times with the number of subjects, stratified by the number of samples, and the graph on the right showing the influence of model complexity. Number of subjects Runtimes for 100 simulations (s) q q q q q q q q q q q q q q q 3 samples 6 samples 12 samples  the dataset. In Figure  11 , there is a slight bend in the curve between a number of subjects equal to 30 and a number of subjects higher than 100. This is because when there are less than 50 subjects in the datasets, more chains are generated by the algorithm in order to improve model stability, increasing the number of 'individuals' considered in the algorithm.
As a result, these models have been extensively studied over the past two decades, with many developments in methods for parameter estimation, experimental design, or model evaluation. Parameter estimation can be performed through maximum likelihood (ML), restricted maximum likelihood (REML), or Bayesian approaches. The saemix package provides an implementation of one of these recent estimation algorithms, SAEM, in R, to obtain maximum likelihood estimates. It offers an interesting alternative to linearisation-based estimation methods as implemented for instance in nlme and lme4. The different packages all have their own specificities and limitations. The oldest, nlme, is included in the base packages distributed with R, and as lme4, can handle multiple levels of variability. They both implement maximum and restricted maximum likelihood estimation.
In addition, for the moment saemix can only handle one response. Finally, the current version of the package is better suited to models with closedform solutions, although there is no limitations from an algorithmic point of view on the models that can be fit. Using a system of differential equations (ODE) can be implemented quite easily within R using deSolve for instance (code available on request), and this approach could be used to implement a model with Michaelis-Menten elimination, for which no closed form solution of the model exists. We however recommend against integrating ODE models with ode(). Indeed, runtimes are exceedingly long even with the theophylline example, which is a very simple model first-order elimination model with only 12 subjects. They would in addition increase linearly with the number of subjects, quickly becoming intractable.
It is also possible with deSolve to implement the structural model in Fortran to improve the computational time but this is beyond the scope of this paper. The SAEM algorithm has been shown to have good statistical properties in various theoretical  (Kuhn and Lavielle (2005) ) and practical applications  (Girard and Mentré (2005) ;  Panhard and Samson (2009) ;  Comets, Lavenu, and Lavielle (2011) ). We show here that the R implementation saemix performs well on reference datasets used to compare different softwares by  Plan et al. (2012). Most of the scenarios run with nlme showed various degrees of convergence issues in that paper, whether started from the true or from altered parameters. Estimates were obtained for only 5% of the datasets simulated under a rich design, and 16% of those simulated under a sparse design in  (Plan et al.
(2012) ). Convergence issues were related to the nonlinearity between dose and response, with a γ of 3. Of note, FOCE in NONMEM or SAS did not show the same convergence problems, suggesting the R implementation of the FOCE algorithm could probably be improved. We attempted here to fit the same datasets with nlmer (data not shown), but encountered even more convergence problems. Note also we could not fit the original dataset directly with nlmer, as the gradient of the model is undefined for a dose of 0. We had to change the dataset by assuming a dose of 0. 01 instead of 0 to avoid this numerical problem. nlme for the same datasets also showed major convergence issues. In the present paper, we therefore compared saemix to nlme and nlmer on a less challenging setting using an E max model instead of the Hill model.
However, nlmer still failed in over half of the datasets we attempted to fit, despite trying different initial parameter estimates. nlme fared better in terms of convergence, but still exhibited more bias than saemix. Both packages were run with the default settings, which might explain in part the poor results we found here. The SAEM algorithm involves repeated evaluations of the model function, and runtimes increase with the number of subjects in the dataset. One possibility to improve the speed of saemix will be to implement a stopping criterion to move into the second phase of the algorithm when the estimates of the parameters stabilise, as has been done in Monolix  (Lavielle (2014) ). A sizeable portion of the time spent in the main algorithm is devoted to evaluations of the model predictions.
This is a much more flexible structure to include for instance complex dose regimens or other inputs to the system. The structure of the saemixModel object is also more easy to interface with other outputs (C programs, specialised libraries) which may not comply with R standards. The future of saemix will indeed be tied in with larger developments that are ongoing in the pharmacometrics community. There is currently a European consortium called DDMoRe working on developing a general language for pharmacometrics, encompassing structural models from pharmacokinetics to system biology and statistical models ranging from frequentist to Bayesian approaches  (Harnisch, Matthews, Chard, and  Karlsson on behalf of the DDMoRe consortium Partners and Contributors (2013)). Two authors of this paper  (Comets and Lavielle)  are part of this endeavour and saemix will be extended to support this language.
In conclusion, we have developed the saemix algorithm in an R package to estimate the parameters of nonlinear mixed effect models. We do not intend to suggest that it is better than the other modelling packages available in R, merely that it can be used to complement the existing array of estimation methods. It is very easy to use, as we show in the current paper how to cast models in nlme and nlmer in terms of saemix, so that it can be used to compare the results for different estimation methods and evaluate the robustness of the analysis. However, the underlying algorithm, SAEM, has been shown to have superior statistical properties and good operational performances when compared to gradient-based algorithms. We hope to extend saemix to more complicated models by using the MLX-TRAN language and solver, which are being externalised to C/C++ so that the models can be called and model predictions obtained from different software.
Other extensions will include more complex data, involving simultaneous modelling of several responses, left-censored data, as well as allow the user to model the interindividual variability through other parameter distributions and write his or her own residual error model. These more advanced features are already available in the Monolix software, which is implemented in MATLAB and has many specific features adapted to PK/PD applications as well as an extensive library of PK/PD models and a modelling language where ODE models can be easily coded. However, the objective in developing the saemix was not to provide a complete rewrite of Monolix, but to offer the SAEM estimation algorithm to the R community, to be used in conjunction with or as an alternative to other packages for nonlinear mixed effect models.
Relative bias and RRMSE for the three packages in the scenario with rich sampling. A first result concerns convergence rates: while saemix provided estimates for all datasets, we could not obtain parameter estimates for more than half of the datasets with nlmer despite retrying a number of times. With nlme, all datasets could be fit in the rich scenario but 7% failed with a sparse design. In terms of estimation errors, saemix provided unbiased estimates in the rich design, while both nlme and nlmer had difficulties to estimate ED 50 and its variability. In the sparse scenario, both saemix and nlme tended to slightly overestimate σ and underestimate variabilities. saemix tended to underestimate the variability on E 0 and ED 50 , while nlme strongly underestimated the variabilities on E rmmax and ED 50.
Acknowledgments We thank  Elodie Plan  and her co-authors Alan Maloney,  France Mentré ,  Mats Karlsson ,  Julie Bertrand , for kindly providing us the datasets used in Plan et al. (2012) to evaluate saemix in Section 4. 3. The datasets can be downloaded from: https://www. dropbox. com/sh/haag3i20lfukyhb/C6kam13Ngs.
Using the object called Orange in this R session as the data. The following SaemixData object was successfully created: Object of class SaemixData longitudinal data for use with the SAEM algorithm Dataset Orange Structured data: circumference ~age | Tree Predictor: age (Days) Tree age circumference 1 1 118 30 2 1 484 58 3 1 664 87 4 1 1004 115 5 1 1231 120 6 1 1372 142 R> orange <-saemixData(name. data = Orange,name. group = "Tree", + name. predictors = "age", name. response = "circumference", + units = list(x = "Days", y = "mm")).
The results of nlmer were not meaningful considering estimates could be obtained in only 16 datasets out of 100, but they also indicated similar issues. saemix nlme nlmer E 0 -2. 51 (5. 84) 0. 29 (4. 69) 2. 12 (5. 46) E max 1. 09 (9. 43) -1. 92 (7. 98) -2. 56 (8. 57) ED 50 -2. 54 (14. 41) -8. 05 (13. 48) -6. 34 (12. 36) σ 1. 12 (6. 86) 1. 94 (6. 26) -0. 95 (6. 11) ω 2 E 0 Emax ω 2 ω 2 ED 50 1. 09 (42. 38) 7. 72 (31. 12) -1. 84 (19. 21) -8. 5 (18. 22) -1. 66 (18. 32) 0. 76 (32) -2. 13 (37. 46) -27. 6 (34. 83) -16. 6 (30. 95) Failed runs 0 0 56.
Materials and Methods In this study, the climatological data cover the period from 1931 to 2014, recorded by weather observatory of N'zerekore (7. 7° N, 8. 8° W), and extracted from the well-preserved archives of Guinea National Weather Service. Among the Guinea observatories, N'zerekore meteorological data are the subject of national, regional and global uses, so they were used for an international study  [12]. As for the variability of agroclimatic parameters, the daily rainfall data from 1961 to 2013 of the synoptic observatory of N'zerekore were used. The Informatics Centre of the National Weather Service has developed robust PC uses CLICOM software for managing its climate data. Other software, such as Climsoft, Excel, Word and INSTAT+, were used for data processing and presentation of results.
This ranking is done according to several criteria described below. Finally, the first observation in the sorted list is selected as the next observation. The filling process is repeated until the gap is filled or there are not target observations left. Each target selected by the Short-term scheduler is sent to the Instrument Control System, and the information on the success of the observation is stored in the database for use in subsequent scheduler runs. Besides, the sorted list of objects can be provided to the operator for override in case of need. Fig. 3  shows an example of the most common situation where the Short-term scheduler needs to repair the mid-term plan to select the next target to be observed. The ranking of the targets is key in the process of repairing the mid-term plan by filling a gap between the last executed observation and the next observations recommended by the mid-term plan.
CAST in real operation CAST is in real operation since September 2016 for the CARMENES survey. The used version includes some additional functionality that is specific to the observing program and that is intended to adapt its performance to the actual working conditions. We list here some of these additional spe-cific features that extend beyond the general design described above: • During the survey it is possible to change the priority of an object to increase the chances of it being selected by CAST. • It is possible to disable an object to make it unavailable for CAST. • CAST can include the observations of telluric standard stars during the nautical twilight in the planning. • Bright stars can also be planned during nautical twilight updating the Night constraint of the affected objects with the desired Sun altitude.
Exoplanet yield simulation In this section we analyze the impact of our scheduler algorithms in the scientific results of the CARMENES survey. We estimate how many planets would be detected in a radialvelocity survey following our planning scheme. As described in the introduction, the CARMENES survey has been designed to discover a statistically significant sample of exoplanets orbiting M dwarfs  [Quirrenbach et al. , 2014]. We have shown that CAST can optimize the survey to obtain over 67 observations per target distributed during the three years of operations. To evaluate the exoplanet yield, we have performed a number of simulations that are described below. We simulated 100 hypothetical planet scenarios using the stellar properties available for our 309 potential targets in the CARMENES database: right ascension, declination, magnitude (which are input parameters of the scheduling algorithms) and spectral type.
Results We tested 100 different random exoplanet scenarios and obtained the average number of planet detections for 50 scheduling trials, thus adding to a total of 5000 simulations. In particular, we computed the average number of detected planets, the average number of detected planetary systems (considering that a system is detected if all its planets are found), and the average number of false positive detections. In order to assess the scientific return of the simulated survey, we analysed the number of planets that we expected to detect and their properties. Each of the scenarios contained a different number of planets orbiting the stars selected from the Carmencita database according to the adopted planet distribution functions. On average, a total of 505 +16 -17 planets were generated for each scenario.