dataesr/scientific-paragraphs-categorization huggingface dataset:
Outlier sample detection To detect a few sample outliers, we compare the standard rule of "6 SDs from the mean" (6SD) used in e. g. EIGENSOFT to the statistic we propose in section 2. 4. Our statistic identifies only isolated samples or isolated pairs that seems to be outliers driving structure of PC17-20 of 1000G (Figure  2 ). All but one outliers are distantly-related pairs that disappear if using a more stringent threshold on relatedness (i. e. using a KING-relatedness cutoff of ~0. 0442 instead of ~0. 0884, see tutorial in section "code availability"). In contrast, rule 6SD identifies a lot of outliers, of which some are part of a relatively large cluster (Figure  S1 ). We recall that, in theory, if all PCs are normally distributed, after correcting for multiple testing of 2500 individuals and 20 PCs, 6SD has a small probability of 0.
4). the "White British" subset of UKBB homogeneous samples, to the use of the robust Mahalanobis distance we propose. Here, we visually choose a threshold of 5 on the log-distance and show that this gives a similar subset of individuals than the "White British" subset reported by the UK Biobank (Figure  S2 ). Moreover, when using this threshold, only 3 out of 10,936 people of self-reported Asian ancestry (1 "Chinese" and 2 "Indian") are kept, and 1 "African" out of 7622 people with Black background is kept (Table  S1 ). In contrast, 416,492 out of 431,090 "British" (96. 6%) and 12,620 out of 12,759 "Irish" (98. 9%) are kept. 3. 3 Projecting onto the PCA space from a reference dataset We use 60% of individuals in the 1000G data (section 2. 6) to compute K=20 PCs. Then, we project the remaining 40% individuals using three methods: 1/ simply multiplying the genotypes of these individuals by the previously computed loadings; 2/ correcting the simple projections using asymptotic shrinkage factors as determined by R package hdpca v1.
1. 3 (Dey and Lee 2019) , with all eigenvalues derived from the Genetic Relationship Matrix (GRM) computed with bed_tcrossprodSelf, one of the new functions of R package bigsnpr; and 3/ the OADP projection (section 2. 5). When simply projecting using loadings, there is negligible shrinkage for PC1 and PC2, a small shrinkage for PC3 and PC4, and a large shrinkage for PC5 to PC8 (Figure  3 ). In contrast, there is no visible shrinkage when projecting new individuals with OADP (Figure  3 ). Simple projection is affected even more by this shrinkage for PC9 to PC20, while OADP still appears free of this bias (Figure  S3 ). We show the same results when projecting the full UK Biobank data onto PCA computed using 1000G data (Figure  S6 ). When correcting projected PC scores with asymptotic shrinkage factors, bias is smaller than with simple projection, yet, there is a visible bias for PC7-8 (Figure  S4 ).
Capturing subtle population structure in the UK Biobank We recomputed PCA in the UK Biobank after restricting the individuals included in the computations: we randomly subsampled UKBB data to use only 10,000 British individuals (out of 431,029) and 5000 Irish individuals (out of 12,755), while keeping all individuals with other or unknown self-reported ancestry. We further removed all pairs of related individuals reported by the UKBB (i. e. both individuals in each pair). This resulted in 48,942 individuals that we used to compute 50 PCs, which took less than 3 hours using function bed_autoSVD (that converged after 4 iterations of automatic LD removal). We show that we are able to capture more PCs (at least 40 instead of 16-18) that display visual population structure (Figures  4  and S12 ).
2018) , resulting PCs are always very similar (absolute correlation larger than 0. 99 for the 20 computed PCs, results not shown). This justifies performing PCA with mean imputation directly on PLINK bed files with a few missing values; this has the advantage to be much faster than having to impute genotyped data using snp_fastImpute, which took 4 days for 406,545 individuals and 240,444 variants, or based on external reference datasets.
Discussion In this work, we have compiled different pitfalls that can arise with Principal Component Analysis (PCA) of genetic data. Then, we have investigated possible solutions to these pitfalls and selected the ones that we found most advantageous, both with respect to properties such as accuracy and robustness, but also computational efficiency and ease of use. We then implemented these solutions in R packages bigsnpr and  bigutilsr. The new functions we provide in R package bigsnpr can be directly applied to genotypes stored as PLINK bed/bim/fam files with some missing values. This contrasts with previous releases of bigsnpr that could only use format "bigSNP". This data format can store both genotype calls and dosages, but requires conversion from other formats and imputation of missing values using functions provided in the package  (Privé et al.
After selecting a subset of individuals, we apply several steps of outlier detection, one for outlier variants that capture long-range LD variation (automatic), and one for detecting outlier samples (semi-automatic and visual). To make these steps more computationally efficient, we explored solutions for not recomputing PCA from scratch when removing a few samples or a few variants. Using educated guesses in R package PRIMME based on low-rank approximations of the updated PCA seemed to be a promising approach but did not reduce computation time by much, so we did not pursue this idea  (Brand 2003; Wu et al. 2017). Once PCA is done, one should check the PC scores (scores of each individual for each PCA dimension) and PC loadings (weights for each variant for each PCA dimension). We differentiate PCs in 3 broad types: the ones capturing LD structure, the ones capturing population structure, and noise.
g. PC1-16 in figure  S9 ) and excluding PCs that do not seem to capture any population structure (e. g. PC17-20 in figure  S9 ). When analysing a dataset that is composed mainly of one population (e. g. British people in the UK Biobank), we found that it is useful to subset these individuals to reduce the imbalance between the different population sizes. Likewise, previous works have shown that uneven population sizes can distort PCs In contrast, a common strategy in genetic analyses is to restrict the analysis to a homogeneous sample to reduce risk of confounding due to population stratification. For that purpose, we show that using the Mahalanobis distance on PC scores can efficiently achieve this goal, which we used in previous analyses  (Privé et al. 2019). When the homogeneous sample is not predominant in the dataset, one solution is to compute the center and covariance of the robust Mahalanobis distance using only the population of interest, and then computing the distances for all individuals using these robust estimates.
Software and code availability All code used for this paper is available at https://github. com/privefl/paper4-bedpca/ tree/master/code. Package bigsnpr can be installed from GitHub (see https://github. com/ privefl/bigsnpr). A tutorial on the steps to perform PCA on 1000G data is available at https: //privefl. github. io/bigsnpr/articles/bedpca. html.
Acknowledgements Authors thank  Rounak Dey  and  Seunggeun Lee  for helpful discussions about PCA projection. This research has been conducted using the UK Biobank Resource under Application Number  25589. F. P. , J. M. and B. V. are supported by the  Danish National Research Foundation  ( Niels Bohr Professorship  to J. M. ), and also acknowledge the  Lundbeck Foundation Initiative for Integrative Psychiatric Research, iPSYCH  ( R248-2017-2003 ).
Effect of the number of particles A dataset of T = 100 observations was generated from the model with parameters a = 0. 5 and σ X = 1. We study how the meeting times and the efficiency vary as a function of N , the number of particles. We set the initial distribution to U[0, 1] over a and U[0, 5] over σ X , and the proposal covariance of the Normal random walk proposals to 0. 2 2 I, corresponding to acceptance rate for the exact algorithm of approximately 36. 6%. In the following we consider a grid of values for the number of particles, varying N between 50 and 250. We estimate large quantiles of the distribution of the meeting time over 20,000 repetitions of coupled PMMH, with the results shown in Figure  3a. As expected, increasing N generally reduces the meeting time at the cost of more computation per iteration.
We also see that for case (3) that the larger value of k can ameliorate the efficiency of the estimators for small numbers of particles. We also examine the inefficiency weighted by the cost of obtaining each estimator, i. e. N IF[H k:m ], and compare this to the inefficiency of the serial algorithm using N V as , with the notation of Section 1. Here, V as was estimated using the spectrum0. ar function in R's CODA package  [Plummer et al. , 2006] , averaging over 10 estimators obtained through running the serial algorithm for 500,000 iterations and discarding the first 10% as burn-in. Figure  4  shows the results of this procedure, showing ±2 sample standard errors for the inefficiency estimates. Figure  4a  demonstrates that despite the lower cost of obtaining unbiased estimators for lower values of N , the initial decline in inefficiency is still significant.
Model, data and target distribution The model aims at capturing the activation of neurons of rats as their whiskers are being moved with a periodic stimulus. The experiment involves M = 50 repeated experiments, and T = 3000 measurements (one per millisecond) during each experiment. The activation of a neuron is recorded as a binary variable for each time and each experiment. These activation variables are then aggregated by summing over the M experiments at each time step, yielding a series of variables Y t taking values between 0 and M ; see  Zhang et al. [2018]  for an alternative analysis that avoids aggregating over experiments. Letting Bin(•; n, p) denote the binomial distribution for n trials with success probability p, the model for neuron activation is given by X 0 ∼ N (0, 1) and, for t ≥ 1, X t |{X t-1 = x} ∼ N (•; ax, σ 2 X ), Y t |{X t = x} ∼ Bin(•; M, s(x)) where s(x) := (1 + exp(-s)) -1.
We use a random sample of T = 2, 000 participants over three waves (one a year) in the run up to the United Kingdom's European Union membership referendum on 23 rd June 2016, regressing the binary outcome asking participants how they would vote in an EU referendum against how they perceive the general economic situation in the UK has changed over the previous 12 months (graded 1-5, with 1='Got a lot worse', 5='Got a lot better'). A detailed description of the data is provided in Appendix A. 6. We allow for correlations between waves through modelling the perturbations ( i1 , i2 , i3 ) ∼ N (0, Σ ρ ) with a generic correlation matrix Σ ρ. In total, we have five unknown parameters θ = (β 1 , β 2 , ρ 2,1 , ρ 3,1 , ρ 3,2 ), with β 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 100 200 300 400 500 5 10 20 600 700 800 N meeting time block pm (a) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e+03 1e+05 ∼ N (0, 10 2 ) and ρ s,t i.
Experimental setup of the acoustic acquisitions Acoustic acquisitions were carried out by B&K type 4189 microphones behind the plate as shown on Fig. 2. The microphone are pre-polarized with a bandwidth of 7 Hz ≤ f ≤ 20 KHz and its measurement range for acoustic levels is between 15. 2 dB and 146 dB. Acoustic measurements were performed with a sampling frequency of 15 KHz. Acoustic signals were acquired by LabVIEW software installed on a National Instruments N I P X I 1036 station equipped with a N I P X I 4496 dynamic card having 16 outputs sampled simultaneously with a sampling frequency of 204. 8 KHz and a resolution of 24 bits.
6) (Table  1 ). When stratifying the general population from the Society Islands by age ( 16 years and >16 years), there was a decline in adults in the two consecutive cross-sectional studies conducted in 2014 and 2015, from 35. 4% (22. 2-50. 5%) to 21. 3% (18. 2-24. 5%) (Figure  1 ). A decline in adults was still observed, albeit with larger uncertainty, when the two datasets were standardised according to the age distribution of the population, with age-adjusted seroprevalence decreasing from 32. 0% (16. 7-62. 1%) to 26. 0% (20. 1-33. 9%) (Table  2 ). In Fiji, in the subset of individuals who were aged over 16 years (n = 122), there was a decrease in seroprevalence by MIA from 24% (17-33%) in 2015 to 7. 3% (3. 4-13%) 2017 (Figure  1 ). There were two seroconversions in the collected samples over this period but 23 seroreversions (McNemar's test, p<0.
Results from the cross-sectional surveys in the Society Islands, French Polynesia, in 2014 and 2015 show a decline in seroprevalence by MIA against all 4 DENV serotypes and ZIKV. However, the 2014 sample included more individuals that tested positive for >1 DENV serotype and are assumed to be a higher risk group. We used a bootstrap method with 10,000 iterations which estimated seroprevalence from a sample of the 2014 dataset, taken with replacement, weighted by the exposure distribution to other DENV viruses in the 2015 survey. After adjusting for the sample bias, there was no evidence of a decline in seroprevalence for DENV-1, DENV-3, or DENV-4, which had circulated in the years preceding the 2014 sample collection (World Health Organisation, 2019), but there remained strong evidence that ZIKV seroprevalence declined between 2014-15.
Serological analysis French polynesia Serum samples collected from blood donors between July 2011 and October 2013 and samples collected from the general population and schoolchildren in 2014 were all tested for presence of IgG antibodies against ZIKV and each of the four DENV serotypes using a recombinant antigen-based indirect ELISA as reported previously  (Aubry et al. , 2017; Aubry et al. , 2015a). Samples collected from the general population in 2015 and from schoolchildren in 2018 were tested by microsphere immunoassay (MIA) using the same recombinant antigens as for the ELISA  (Cao-Lormeau et al. , 2016; Aubry et al. , 2017; Kama et al. , 2019). Recombinant antigens used in both assays comprised domain III of the envelope glycoprotein of ZIKV, DENV-1, DENV-2, DENV-3, or DENV-4 strains (respective GenBank accession no.
3%) for DENV (k = 0. 80). In this study, we tested what proportion of samples for the 45 participants in the full dataset (i. e. 135 samples in total) that were seropositive or seronegative by MIA had the same result by the neutralization assay. We found that 54/68 (79. 4%)   samples that were positive to ZIKV by MIA were also positive in the neutralization assay, and 42/67 (62. 7%) who were seronegative were also negative by neutralization assay (k = 0. 42). We also calculated the proportion positive by neutralization assay that had the same result by MIA. We found that 54/79 (68. 4%) samples that were positive to ZIKV in the neutralization assay were also positive by MIA, and 42/56 (75%) who were seronegative were also negative by MIA.
Statistical analysis For data from Fiji, where serial samples were collected from the same individual, changes in seroprevalence between studies were tested using McNemar's test. In French Polynesia, chi-squared tests were performed to test for evidence of a change in seroprevalence between two cross-sectional surveys. Changes in mean log titre between groups were analyzed using a t-test. To analyse the potentially non-linear relationship between DENV neutralization titres and seroprevalence by MIA and neutralization test (Figure  3 ), we used a generalized additive model via the mcgv package in R  (Wood, 2019). The model was of the form g(E(y))=b + f(x), where y was the binary outcome variable, x was the predictor (i. e. titre), g was the link function, b was the intercept, and f was a smooth function represented by a penalized regression spline.
Mean DENV titre was calculated as the mean of log titres against the four DENV serotypes for each participant. All data and code used in the analysis are available at: https://github. com/a-henderson91/zika-sero-pacific/ (Henderson and Kucharski, 2019 ; copy archived at https://github. com/elifesciences-publications/zika-sero-pacific/ settings). French The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements We are grateful to the  Minister for Education of French Polynesia  and to the directors, teachers, nurses and schoolchildren from the elementary and junior high schools selected for the serosurvey in 2018. We greatly thank all the participants and community leaders in Fiji who generously contributed to the study over the three visits. We would like to acknowledge the work of the field teams:  Jessica Paka ,  Amele Ratevono ,  Warren Fong ,  Manisha Prakash ,  Jonetani Bola ,  Mosese Ligani , and  Taina Naivalu  (2017);  Meredani Taufa ,  Adi Kuini Kadi ,  Jokaveti Vubaya ,  Colin Michel ,  Mereani Koroi ,  Atu Vesikula , and  Josateki Raibevu  (2015);  Dr. Kitione Rawalai ,  Jeremaia Coriakula ,  Ilai Koro ,  Sala Ratulevu ,  Ala Salesi ,  Meredani Taufa , and  Leone Vunileba  (2013).
To validate a new myelin estimation technique, histological studies are powerful and considered as the standard criterion, but its significance is somewhat limited as histological staining is evaluated based on optical density contrasts and therefore does not directly quantify myelin  [6]. Recently, myelin measurement using quantitative synthetic MRI (SyMRI) was proposed  [7]. This approach was originally developed for the simultaneous measurement of longitudinal T1 relaxation rate (R1), transverse T2 relaxation rate (R2), proton density (PD), and local radiofrequency field B1, using multislice, multiecho, and multidelay acquisition. The estimated B1 field is used for the correction of local variations in flip angle. Full head coverage is obtained in approximately 6 minutes  [7]. From these absolute parameters, it is possible to create any contrast-weighted image that is clinically useful, including T1-weighted (T1w) or T2-weighted (T2w) images, using SyMRI software  [8] [9] [10].
Study Participants Thirty-seven patients with relapsing-remitting MS were prospectively included in this study, from March 2017 to July 2018. These patients were diagnosed with MS according to the standard McDonald criteria  [40]. Of these patients, 14 did not have any plaque or had only small (<5 mm maximal length) plaques on brain MRI and 2 had diffuse extensive WM abnormalities. Therefore, these 16 patients were excluded and the remaining 21 patients were included in the analysis. Disability was assessed using the Expanded Disability Status Scale (EDSS) score. The demographic and clinical data of the patients with MS are provided in Table  1. The institutional review board of Juntendo University Hospital approved this study on 18 March 2016 and written informed consent was obtained from all participants.
Acquisition and Processing of SyMRI Data MR relaxometry was performed with the QRAPMASTER (quantification of relaxation times and proton density by multi-echo acquisition of a saturation recovery using turbo spin-echo readout) pulse sequence, which is a multi-slice, multi-echo, and multi-saturation delay acquisition sequence. Two sets of echo times (TE) and 4 sets of delay times were used to generate 8 complex images in each section in order to quantify R1, R2, and PD. The TE were 22 and 99 ms, the delay times were 170, 620, 1970, and 4220 ms and the repetition time (TR) was 4250 ms. An in-plane resolution of 0. 8 × 0. 8 mm with a slice thickness/gap of 4. 0/1. 0 mm was used for 30 slices. The field-of-view (FOV) was set to 230 × 186 mm with a matrix size of 320 × 208. The acquisition time was 5 min and 8 s.
With the assumption that all the R1, R2, and PD values of MVF, excess parenchymal water volume fraction, cellular volume fraction, and free water volume fraction contribute to the effective R1, R2, and PD in each acquisition voxel; a model was produced to estimate partial volumes of these four compartments, as described by Warntjes et al. [7]. This was done by running Bloch equations and optimizing model parameters in a spatially normalized and averaged brain from a group of healthy controls. Damage to myelin, even in the NAWM of patients with MS  [18] , is supposed to cause deviation from the R1, R2, and PD in the healthy state, leading to a decrease in the calculated SyMRI MVF. Using this model, MVF maps were created from R1, R2, and PD maps via SyMRI software (version 8. 04; SyntheticMR, Linköping, Sweden).
Processing of the T1w/T2w Ratio Synthetic T1w and T2w images were created using post-processing TR of 500 ms and TE of 10 ms and TR of 4500 ms and TE of 100 ms, respectively, on SyMRI software based on R1, R2, and PD maps. To reduce the effect of intensity scaling, we used a calibration algorithm based on 2 anatomical masks (eye and temporal muscle), as proposed by Ganzetti et al. [41]. We recorded the modes as reference values for the eyes as follows: 28. 2 for T1w images and 99. 9 for T2w images. For the temporal muscle, these values were 58. 6 for T1w images and 21. 1 for T2w images. After calibrating the T1w and T2w images, their ratio was calculated to produce T1w/T2w ratio images. The T1w/T2w ratio was further linearly calibrated, as reported in a previous study  [13] , to be comparable to SyMRI MVF.
Acquisition and Processing of Radial Diffusivity Whole-brain diffusion-weighted imaging was performed using spin-echo planar imaging employing a b-value of 1000 s/mm 2 along 64 uniformly distributed motion-probing directions in the anteroposterior phase-encoding direction. Standard and reverse phase-encoded blipped images with no diffusion weighting (b = 0 s/mm 2 ) were also acquired. Other imaging parameters were TR, 3300 ms; TE, 70 ms; flip angle, 90 • ; FOV, 229 × 229 mm; matrix size, 130 × 130; resolution, 1. 8 × 1. 8 mm; slice thickness, 1. 6 mm; acquisition time, 3 min and 55 s. All datasets were free from severe artifacts, such as gross geometric distortion, signal dropout, or bulk motion. Diffusion-weighted imaging data were then corrected for susceptibility-induced geometric distortions, eddy current distortions, and inter-volume subject motion using TOPUP and EDDY tools  [42].
Image Analysis Synthetic T2w images and maps of SyMRI MVF , MTsat, RD, and T1w/T2w (Figure  1 ) were analyzed by using ITK-SNAP software (3. 8. 0-beta version) by a radiologist (L. S. ) with 3 years of experience. Plaques were defined as a WM area of more than 5 mm in diameter in the supratentorial area, with abnormally high intensity on T2w images. Plaques were semi-automatically segmented on synthetic T2w images using an adaptive brush tool for automatic and objective adjustment to lesion boundaries, choosing an adaptive algorithm with high granularity and low smoothness to undersegment plaques slightly in order to mitigate possible misregistration and the partial volume effect. Periplaque regions were manually segmented with a small round-shaped brush (approximately 1/5th of plaque diameter).