dataesr/scientific-paragraphs-categorization huggingface dataset:
Wrapper mleHetGP feeds lower-level objective hetGP:::logLikHet and derivative hetGP:::dlogLikeHet into optim after finding replicates with find_reps. Binois et al. (2018)  show that l is maximized when ∆ n = Λ n and g = 0. In other words, smoothing the latent noise variables (  8 ) is unnecessary at the MLE. Optimizing the objective naturally smooths the latent ∆ n values, leading to GP predictions that interpolate those smoothed values. However, smoothing is useful as a device in three ways: (1) theoretically: connecting SK to Goldberg's latent representation; (2) numerically: annealing to avoid local minima; and (3) practically: yielding a smooth solution when the numerical optimizer is stopped prematurely, which may be essential in big data (large unique-n) contexts.
Package structure The main functions in hetGP are either related to GP fitting or their use in sequential design tasks, which are the subject of Section 4. Inferential interface functions mleHomGP, mleHetGP, mleHomTP and mleHetTP return homoskedastic/heteroskedastic GP/TP models in the form of S3 objects. They are completed with S3 methods: • plot, print, predict; • update to add one or several new observations, processing them to identify replicates. Warm-started re-estimation of the hyperparameters can be invoked with argument maxit > 0. • strip to save a bare representation of the model by removing all pre-computed quantities that enable fast prediction like inverse covariance matrices. Conversely, the rebuild method recomputes all missing quantities. Note that Cholesky decomposition is used in the hyperparameter optimization for speed, while the generalized inverse via MASS's ginv can optionally be used when rebuilding.
Require: previous hetGP model, Xnew, Znew 1: if Xnew is a replicate then 2: (Optional) Use leave-one-out prediction to re-estimate δ i and empirical estimate σ i 3: else 4: Predict λ n+1 based on previous (n) fit 5: (Optional) Combine λ n+1 prediction with empirical estimate 6: end if 7: if maxit > 0 then 8: Fit heteroskedastic GP on D N +1 with initial hyperparameters θ Hom , θ (g) , g and ∆ n+1. 9: else 10: Update internal quantities (e. g. , inverse covariance matrices) 11: end if Algorithm 2: Pseudo code for the update of hyperparameters procedure in mleHetGP.
R> Y <-apply(X, 1, sirEval) The code below fits a hetGP model. By default, lengthscales for the variance GP are linked to those from the mean GP, requiring that the former be a scalar multiple k > 1 of the latter. That can be undone, however, as we do below primarily for illustrative purposes. Bayes factor surfaces. Bayes factor (BF) calculation for model selection is known to be sensitive to hyperparameter settings, which is further complicated (and obscured) by Monte Carlo evaluation that interjects a substantial source of noise. To study BF surfaces in such settings,  Franck and Gramacy (2018)  propose treating expensive BF calculations, via MCMC say, as a (stochastic) computer simulation experiment. The idea is that BF calculations at a space-filling design in the (input) hyperparameter space can be used to map out the space and better understand the sensitivity of model selection to those settings.
If a random request comes in for a product that is missing a key item, a replenishment order is executed, which is filled after a random delay. Holding items in inventory is expensive, so there is a balance between inventory costs and revenue. Hong and Nelson built a MATLAB simulator for this setup, which was reimplemented by  Xie, Frazier, and Chick (2012). Binois et al. (2018)  describe an out-of-sample experiment based on this latter implementation in its default setting (originally from Hong and Nelson), specifying item cost structure, product makeup (their items) and revenue, and distribution of demand and replenishment time, under target stock vector inputs b ∈ {0, 1,. , 20} 8 for eight items. Here we provide code that can be used to replicate results from that experiment, which involved a uniform design of size n tot = 2000 in 8D space, with ten replicates for a total design size of N tot = 20000.
R> data("ato") In order to create an out-of-sample test setting, random training-test partitions were constructed by randomly selecting n = 1000 unique locations and a uniform number of replicates a i ∈ {1,. , 10} thereupon. The ato data object also contains one such random partition, which is subsequently coded into the unit cube [0, 1] 8. Further detail is provided in the package documentation for the ato object. Actually the object provides two test sets. One is a genuine out-of-sample test set, where the test sites do not intersect with any of the training sets. The other is replicate based, involving replicates 10a i not selected for training. The training set is large, which makes MLE calculations slow, so the ato object provides a fitted model for comparison. In the examples section of the ato documentation, code is provided to reproduce that fit or to create a new one based on a new random training-test partition.
R> c(n = nrow(Xtrain), N = length(unlist(Ztrain)), time = out$time) n N time. elapsed 1000. 000 5594. 000 8583. 767 The main reason for storing these objects is to enable a fast illustration of prediction and outof-sample comparison, in particular to have something to compare with a more thoughtful sequential design scheme outlined in Section 4. The code below performs predictions at the held-out test locations and then calculates a proper score (Gneiting and Raftery 2007, Equation (  27 )) against the ten replicates observed at each of those locations. Higher scores are better. A function scores in hetGP computes scores and root mean squared error (RMSE). R> sc <-scores  (out, Xtest, matrix(unlist(Ztest) , byrow = TRUE, ncol = 10)) A similar calculation for the held-out training replicates is shown below.
IMSPE I n+1 ≡ IMSPE(x 1 ,. , xn , x n+1 ) = x∈D σ 2 n+1 (x) dx. ( 11 ) This formula is written in terms of unique-n inputs for reasons that we shall clarify shortly. It could, however, instead be phrased in terms of full-N equations. The idea would be to choose the next design location, x N +1 , which could be a new unique location (x n+1 ) or a repeat of an existing one (i. e. , one of x1 ,. , xn ), by minimizing I n+1. The presentation in Binois et al. is more careful, but also more cumbersome, with the notation in this regard. In general, the integral in Equation (  11 ) requires numerical evaluation  (Seo et al. 2000; Gramacy and Lee 2009; Gauthier and Pronzato 2014; Gorodetsky and Marzouk 2016; Pratola, Harari, Bingham, and Flowers 2017b). For example, the tgp package (option Ds2x=TRUE) sums over a reference set, as well as averaging over the posterior distribution of (treed) Gaussian process parameterization.
However, conditional on GP hyperparameters, and when the study region D is an easily integrable domain such as a hyperrectangle, the requisite integration may be calculated in closed form. Although examples exist in the literature (e. g. ,  Ankenman et al. 2010; Anagnostopoulos and Gramacy 2013; Burnaev and Panov 2015; Leatherman, Santner, and Dean 2017)  for particular cases (and approximations), we are not aware of examples providing the level of specificity and generality in derivation, or implementation in code as provided in hetGP.
R> XX <-matrix(seq(0, 1, by = 0. 005)) Below we shall refer to these surfaces as "green" and "blue," respectively, referencing the colors from Figure  6. The "knots" are shown as red open circles. Code below implements the closed form IMSPE of Equation (  12 ) for a generic variance function r, like one of our splines from above. The implementation uses some hetGP functions such as Wij and cov_gen. (We shall illustrate the intended hooks momentarily; these lowlevel functions are of value here in this toy illustration, and as a window into the spirit of implementation in the package. ) R> IMSPE. r <-function(x, X0, theta, r) { + x <-matrix(x, nrow = 1) + Wijs <-Wij(mu1 = rbind(X0, x), theta = theta, type = "Gaussian") + K <-cov_gen(X1 = rbind(X0, x), theta = theta) + K <-K + diag(apply(rbind(X0, x), 1, r)) + return(1 -sum(solve(K) * Wijs)) + } The next step is to apply this function on a grid for each of the two choices for r(x), green and blue.
, xn rather than a new xn+1 when r(x N +1 ) ≥ k n (x N +1 ) ⊤ K -1 n W n K -1 n k n (x N +1 ) -2w ⊤ n+1 K -1 n k n (x N +1 ) + w n+1,n+1 tr(B k * W n ) -σ 2 n (x N +1 ), where k * = argmin k∈{1,. ,n} IMSPE(x k ) and B k = (K -1 n ). ,k (K -1 n ) k,. νλ k a k (a k +1) -(Kn) -1 k,k. Basically, this relationship says that IMSPE will prefer replication when the variance is "large enough. " Rather than read tea leaves more deeply than that, let us see it in action in our toy example. The code below utilizes some of hetGP's internals to enable evaluation of the right-hand side of the inequality above, in other words, treating it as an equality. R> rx  X0, rn, theta, Ki, kstar, Wijs) { + x <-matrix(x, nrow = 1) + kn1 <-cov_gen(x, X0, theta = theta) + wn <-Wij(mu1 = x, mu2 = X0, theta = theta, type = "Gaussian") + a <-kn1 %*% Ki %*% Wijs %*% Ki %*% t(kn1) -2 * wn %*% Ki %*% t(kn1) + a <-a + Wij(mu1 = x, theta = theta, type = "Gaussian") + Bk <-tcrossprod(Ki[, kstar], Ki[kstar,]) / + (2 / rn[kstar] -Ki[kstar, kstar]) + b <-sum(Bk * Wijs) + sn <-1 -kn1 %*% Ki %*% t(kn1) + return(a / b -sn) + } Calling that function with the XX grid defined above, with matching covariance structure details, commences as follows.
The second attempts to adapt to minimize IMSPE regardless of computational cost: Adapt: h N +1 ∼ Unif{a ′ 1 ,. , a ′ n } with a ′ i := max(0, a * i -a i ), and a * i ∝ r(x i )(K -1 n W n K -1 n ) i,i comes from a criterion in the SK literature  (Ankenman et al. 2010). The code below duplicates the example above with the adapt heuristic. Alternatively, a horizon call can be provided target and previous_ratio arguments to implement the target heuristic instead. First, reinitialize the design. R> X <-seq(0, 1, length = 10) R> Y <-fr(X) R> mod. a <-mleHetGP(X = X, Z = Y, lower = 0. 0001, upper = 1) R> h <-rep  (NA, 500)  Next, loop to obtain N = 500 observations under the adaptive horizon scheme. The left panel of Figure  8  shows the adaptively selected horizon over the iterations of sequential design.
Contour finding. A related problem is contour finding, also known as level-set estimation. The objective is to identify a region of inputs of interest: Γ f = x ∈ R d : Y (x) > T with T a given threshold, often zero without loss of generality. We describe briefly here the criteria defined by  Lyu et al. (2018)  for both GPs and TPs that are implemented in hetGP. Several others can be found in the literature  (Chevalier, Ginsbourger, Bect, and Molchanov 2013; Bogunovic et al. 2016; Azzimonti, Ginsbourger, Chevalier, Bect, and Richet 2016) , with a selection of implementations provided by the KrigInv package  (Chevalier, Picheny, and Ginsbourger 2014b; Chevalier, Picheny, Ginsbourger, and Azzimonti 2018). The simplest criterion is maximum contour uncertainty (MCU), implemented by crit_MCU in the hetGP package.
For illustration, let us consider finding the zero contour in our simple 1D example. The code below illustrates the crit_cSUR criteria, the hetGP package's generic crit_optim interface. Otherwise the setup and for loop are identical to those for our IMSPE and EI-based examples. R> X <-seq(0, 1, length = 10) R> X <-c(X, X, X) R> Y <-fr(X) R> mod <-mleHetGP(X = X, Z = Y) R> R> for(i in 1:470) { + opt <-crit_optim(mod, crit = "crit_cSUR", h = 5, ncores = ncores) + X <-c(X, opt$par) + Ynew <-fr(opt$par) + Y <-c(Y, Ynew) + mod <-update(mod, Xnew = opt$par, Znew = Ynew, ginit = mod$g * 1. 01) + if(i %% 25 == 0) { + mod2 <-mleHetGP(X = list(X0 = mod$X0, Z0 = mod$Z0, mult = mod$mult), + Z = mod$Z, lower = 0. 0001, upper = 1) + if(mod2$ll > mod$ll) mod <-mod2 + } + } R> R> p <-predict  (mod, matrix(xgrid, ncol = 1 )) R> pvar <-p$sd2 + p$nugs Again, we see that using the biasing lookahead procedure reduces the number of unique designs, which is beneficial in terms of speed as well as accuracy in general.
Summary and discussion We have introduced the R package hetGP for heteroskedastic Gaussian process regression, sequential experimental design, and Bayesian optimization. Although the package is designed for dealing with noise and changing signal-to-noise ratios in the setting of Gaussian process regression, it offers a full-featured GP regression approach. Ordinary homoskedastic (and noise-free) GP modeling is supported. When the data is observed with noise, the package implements a Woodbury trick to make inference more efficient, decomposing matrices sized new sequential design acquisition strategies toward novel surrogate modeling and prediction applications.
Acknowledgments The work of MB is supported by the  U. S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research , under Contract No. DE-AC02-06CH11357. RBG gratefully acknowledges funding from a  DOE   LAB 17-1697  via subaward from  Argonne National Laboratory for SciDAC/DOE Office of Science ASCR  and  High Energy Physics , and partial support from  National Science Foundation  grants  DMS-1849794 ,  DMS-1821258  and  DMS-1621746. Many thanks to  D. Austin Cole  for comments on early drafts and to  Gail Pieper  for her useful language editing.
2% per °C) in the range of temperature considered (figure 1, pink line), which is consistent with some experimental data  (Cochard et al 2007)  PEER COMMUNITY IN FOREST AND WOOD SCIENCES 5 The osmotic potential PT of leaf cells increases with leaf temperature according to van't Hoff relation: Π W = -𝑐𝑅𝑇 ( 10 ) where c is the concentration of solutes in the cells and R the gas constant. Therefore, if P20 is the osmotic potential at 20°C and T is expressed in °C I have: Π W = Π [a Wl[:m. 79 [nm. 79 (11) As the total leaf water potential Yleaf remains constant when temperature is changing, the leaf turgor pressure Pleaf must also increase with temperature as: 𝜓 &'() = Π &'() + 𝑃 &'() (12) This may impact in return leaf stomatal conductance and transpiration under the hypothesis that these variables are determined by Pleaf.
Acknowledgements I was supported in part by  IDEX-ISITE  initiative  16-IDEX-0001  ( CAP2025 ) and the  ANR-18-CE20-0005  Hydrauleaks. Version 2 of this preprint has been peer-reviewed and recommended by  Peer Community In Forest and Wood Sciences  (https://doi. org/10. 24072/ pci. forestwoodsci. 100002 ).
Here, we focused on all aspects downstream of these steps, as these upstream steps are largely dataset-and technology-dependent, and have been previously reviewed and benchmarked. Instead, we set out to optimise statistical analysis workflows that were originally developed for bulk RNA-seq studies in the sc-eQTL mapping setting. In particular, we take advantage of our unique, population-scale, dataset of matched (i. e. from the same individuals) single-cell and bulk expression data from a homogenous cell population (induced pluripotent cells; iPSCs) to evaluate the effect of different statistical choices. Because we have matched samples, we are able to measure performance both in terms of sc-eQTL yield (i. e. the number of significant sc-eQTL) and in terms of concordance with bulk-eQTL, which can be considered as the gold standard for eQTL mapping.
We also validate these findings on simulated population-scale scRNA-seq datasets, where the eQTL effects are known, and performance can be measured in terms of power and false discovery rates. Using both empirical and simulated datasets, we identify best practices in terms of aggregation and normalisation, the type and number of expression covariates to include as fixed effects, methods to account for single-cell sampling variation in the model, and methods to guide multiple testing correction using information from bulk RNA sequencing data. Finally, we demonstrate how applying these best practices increases power to detect sc-eQTL and reduce the number of false discoveries.
We chose to explore the dr aggregation level because in the datasets we considered some donors are present in multiple runs  [12, 13]. In all cases (i. e. using any of the aggregation methods), aggregated expression values were only calculated for samples (i. e. donors or donor-run combinations) with at least 5 cells. RNA pre-processing steps to obtain count matrices to perform eQTL mapping, including gene expression quantification, cell and gene-level quality control (QC) and cell type annotation. These steps are not optimised/tested in this work (shown in grey). (b) Different approaches tested to perform eQTL mapping using scRNA-seq profiles. Starting from one gene x cell count matrix obtained as in a, counts were aggregated per sample (i. e. donor, or donor-run combination), either by summing the data first at samplelevel and then normalising using methods designed for bulk RNA-seq (i.
We found that the normalised counts are highly correlated (Pearson's R 2 mean: 0. 88, median 0. 93, Methods). On the other hand, summed count values (both dr-sum and d-sum) were obtained directly from the raw count data (i. e. non-normalised). Normalisation was then applied on the resulting pseudo-bulk counts, using methods typically used for bulk RNA-seq data. In particular, we perform TMM normalisation on the aggregated counts, as implemented in edgeR  [28] , one of the best-established methods for bulk RNA-seq normalisation (Fig. 1b , Table  1 ), followed by log transformation. We test these approaches on two empirical datasets and two simulated datasets. The main analyses are performed on expression data gathered from iPSCs, for which we have bulk RNAsequencing  [29]  and Smart-Seq2 single-cell RNA-sequencing  [12]  data.
This unique data allowed us to compare approaches to discovery of eQTL in single-cell expression data, using bulk eQTL results from the same samples as the gold standard for evaluating approaches for single-cell data. Specifically, we assessed the numbers of eGenes discovered in single-cell data and the extent of replication of bulk eQTL effects in single-cell analyses. Additionally, we consider one cell type (midbrain floor plate progenitor; FPP) from a large 10X single-cell RNA-sequencing differentiation study, differentiating iPSCs towards dopaminergic neurons  [13]. Lastly, to support our findings, we test our approaches on simulated single-cell eQTL datasets generated using SplatPop (Methods), based on either the Smart-Seq2 iPSC expression dataset or the 10X differentiating neuron dataset.
Smart-Seq2 datasets We first focused on the single-cell data from Cuomo et al. [12] , from which we selected the iPSC data (day0) to compare to bulk data from the HipSci consortium  [2, 29]. We selected the iPSC data because the homogeneous expression profiles of these cells make it an ideal cell type for performing this kind of study. From the HipSci resource we selected 87 healthy donors of European descent for which we had quality controlled: 1) genotype data, 2) Smart-Seq2 sc-RNAseq, and 3) bulk RNA-seq (i. e. matched bulk, from hereon "m-bulk"). Additionally we used a superset of 526 samples (including the previous 87) for which we had quality controlled: 1) genotype data, and 2) bulk RNAseq data, for reference (i. e. all bulk, "a-bulk"). We (re)processed the raw RNA-seq data from the single-cell study to match the bulk processing as much as possible (Methods).
Replication fractions remain consistent when considering all HVGs for eQTL (Table  S1 ). Moreover, we see higher replication rates considering a-bulk as compared to m-bulk, indicating that some of the effects found in the single-cell data can only be picked up from bulk datasets with more samples. When specifically looking at the effects that get replicated, we observe that they are highly overlapping (86%) across aggregation strategies. This result indicates that the same effects that get replicated in d-sum, are also replicated in dr-mean, but since there are more effects found in dr-mean the fraction is lower. While the percentage of bulk-eQTL replicated by sc-eQTL mapping in a matched dataset is a powerful performance metric, a limitation of empirical benchmarks like this is that the ground truth is not known.
To validate our results on a dataset with completely known true eQTL effects, we simulated single-cell expression profiles for genes on chromosome 2 for a population with known eQTL effects applied to 35% of genes (total expressed genes on chr2: n=1,255, eQTL Genes=439). To ensure the simulations reflected real data, expression statistics from the iPSC Smart-Seq2 dataset were used to estimate key parameters used in the simulations (Fig. S2 ; see Methods for details). We then performed aggregation, normalisation, and eQTL mapping as described for the empirical study on 10 replicate simulated datasets and quantified performance in terms of power (fraction of true eQTL detected), empirical FDR (i. e. fraction of false eQTL detected at an FDR < 5%), and effect size correlation (see Methods for details).
In both cases, we observed a high number of eGenes (1,668 for the first approach, corresponding to ~60% genes tested and 840 using the second approach, i. e. 31% of the genes tested). However, the replication rates were very low, ranging from 10 to 15% when considering the two approaches and the two sets of bulk results (Methods), likely suggesting inflation and a high rate of false positives, confirmed by poor correlation of effect sizes between these results and those obtained using bulk RNA-seq (Fig. S4 ). To further confirm this result, we tested the first approach on one of the simulated datasets described above. This analysis resulted in a 75. 6% false discovery rate, a power of just 18%, and a correlation between the absolute value of simulated and estimated effect sizes of 0. 832 (Fig.
10X datasets To check the generality of the results obtained from Smart-Seq2 data, we proceeded to assess the same aggregation/normalisation strategies on data generated with the 10X Chromium technology  [32]. Recent large-scale single-cell studies have predominantly used 10X because of the lower cost and higher cell throughput of the droplet-based scRNA-seq technology, as compared to plate-based technologies such as Smart-Seq2  [33]. Though cell throughput is higher, most often the number of reads quantified is lower for 10X studies, and reads are from the '3 or '5 end of transcripts (or both), but not from full-length transcripts as in Smart-Seq2. Specifically, we selected a midbrain floor plate progenitor (FPP) cell population from a recent differentiation study of iPSCs toward dopaminergic neurons  [13]  (n=174, Methods).
Unfortunately, we did not have matching bulk RNA-seq data to assess replication, and thus could only assess differences in terms of cis-eQTL discovery power. In general, we observed very similar effects as found in the Smart-Seq2 dataset (Table  3 ). Again the donor-run methods (dr-mean, dr-sum and dr-median) outperformed donor methods (d-mean, d-sum and d-median), and as for Smart-Seq2 we observe that the single-cell normalisation method with mean aggregation outperformed the sum-based aggregation followed by bulk-like normalisation, supporting the previous observation (both from the simulated and the empirical Smart-Seq2 iPSC data) that single-cell-specific normalisation and treating different runs as separate entities is beneficial. Moreover, for 10X data, single-cell normalisation with median aggregation performed especially poorly, likely due to the greater sparsity in 10X data compared to Smart-Seq2 data.