# Comparing quantiles for two samples

Recently, for a research paper, I got some samples, and I wanted to compare them. Not to compare their means (by construction, all of them were centered) but there dispersion. And not their variance, but more their quantiles. Consider the following boxplot type function, where everything here is quantile related (which is not the case for standard boxplot, see http://freakonometrics.hypotheses.org/4138, in French)

> boxplotqbased=function(x){
+ q=quantile(x[is.na(x)==FALSE],c(.05,.25,.5,.75,.95))
+ plot(1,1,col="white",axes=FALSE,xlab="",ylab="",
+ xlim=range(X),ylim=c(1-.6,1+.6))
+ polygon(c(q[2],q[2],q[4],q[4]),1+c(-.4,.4,.4,-.4))
+ segments(q[1],1-.4,q[1],1+.4)
+ segments(q[5],1,q[4],1)
+ segments(q[5],1-.4,q[5],1+.4)
+ segments(q[1],1,q[2],1)
+ segments(q[3],1-.4,q[3],1+.4,lwd=2)
+ xt=x[(x<q[1])|(x>q[5])]
+ points(xt,rep(1,length(xt)))
+ axis(1)
+ }

(one can easily adapt the code for lists, e.g.). Consider for instance temperature, when the (linear) trend is removed (see http://freakonometrics.hypotheses.org/1016 for a discussion on that series, in Paris),

from January 1st till December 31st. Let us remove now the seasonal cycle, i.e. we do have here the difference with the average seasonal temperature (with here upper and lower quantiles),

Seasonal boxplots are here (with Autumn on top, then Summer, Spring and Winter, below),

If we zoom in, we do have (where upper and lower segments are 95% and 5% quantiles, while classically, boxes are related to the 75% and the 25% quantiles)

Is there a (standard) test to compare quantiles – some of them perhaps ? Can we compare easily quantiles when have two (or more) samples ?

Note that this example on temperature could be related to other old posts (see e.g. http://freakonometrics.hypotheses.org/2190), but the research paper was on a very different topic.

Consider two (i.i.d.) samples $\{x_1,\cdots,x_m\}$ and $\{y_1,\cdots,y_n\}$, considered as realizations of random variables $X$ and $Y$. In all statistical courses, tests on the average are always considered, i.e.

$H_0:\mathbb{E}(X)=\mathbb{E}(Y)$

against

$H_1:\mathbb{E}(X)\neq\mathbb{E}(Y)$

Usually, the idea in courses is to start with a one sample test, and to test something like

$H_0:\mathbb{E}(X)=\mu_\star$

against

$H_1:\mathbb{E}(X)\neq\mu_\star$

The idea is to assume that samples are from Gaussian variables,

$T = \frac{\overline{x} - \mu_\star}{\widehat{\sigma}/\sqrt{n}}$
Under $H_0$$T$ has a Student t distribution. All that can be found in any Statistics 101 course. We can derive $p$value, computing probabilities that $T$ exceeds the observed values (for two sided tests, the probability that the absolute value of $T$ exceed the absolute value of the observed statistics). This test is closely related to the construction of confidence intervals for $\mu$. If $\mu_\star$ belongs to the confidence interval, then it might be a suitable value. The graphical representation of this test is related to the following graph

Here the observed value was 1,96, i.e. the $p$value (the area in red above) is exactly 5%.

To compare means, the standard test is based on

$T = {\overline{x} - \overline{y} \over \displaystyle\sqrt{{s_x^2 \over m} + {s_y^2 \over n}} }$

which has – under $H_0$ – a Student-t distribution, with $\nu$ degrees of freedom, where

$\nu = \frac{(s_x^2/m + s_y^2/n)^2}{(s_x^2/m)^2/(m-1) + (s_y^2/n)^2/(n-1)}.$

Here, the graphical representation is the following,

But tests on quantiles are rarely considered in statistical courses. In a general setting,define quantiles as

$Q_X(p)=\inf\left\{ x\in \mathbb R : p \le \mathbb P(X\leq x) \right\}$

one might be interested to test
$H_0:Q_X(p)=Q_Y(p)$
against
$H_1:Q_X(p)\neq Q_Y(p)$
for some $p\in(0,1)$. Note that we might be interested also to test if

$H_0:Q_X(p_k)= Q_Y(p_k)$
for all $k$, for some vector of probabilities $\boldsymbol{p}=(p_1,\cdots,p_d)\in(0,1)^d$.
One can imagine that this multiple test will be more complex. But more interesting, e.g. a test on boxplots (are the four quantiles equal ?).  Let us start with something a bit more simple: a test on quantiles for one sameple, and the derivation of a confidence interval for quantiles.

• Quantiles for one sample

The important idea here is that it should be extremely simple to get $p$values. Consider the following sample, and let us run a test to assess if the median can be zero.

> set.seed(1)
> X=rnorm(20)
> sort(X)
[1] -2.21469989 -0.83562861 -0.82046838 -0.62645381 -0.62124058 -0.30538839
[7] -0.04493361 -0.01619026  0.18364332  0.32950777  0.38984324  0.48742905
[13]  0.57578135  0.59390132  0.73832471  0.82122120  0.94383621  1.12493092
[19]  1.51178117  1.59528080
> sum(X<=0)
[1] 8

Here, 8 observations (out of 20, i.e. 40%) were below zero. But we do know the distribution of $N$ the number of observation below the target

$N=\sum_{i=1}^n \boldsymbol{1}(X_i\leq x_\star)$

It is a binomial distribution. Under $H_0$, it is a binomial distribution $\mathcal{B}(n,p_\star)$ where $p_\star$ is the probability target (here 50% since the test is on the median). Thus, one can easily compute the $p$-value,

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="white")
> abline(v=sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n,dbinom(n,size=20,prob=0.50),type="s")

which yields

Here, the $p$-value is

> 2*pbinom(sum(X<=0),20,.5)
[1] 0.5034447

Here the probability is easy to compute. But one can observe that there is some kind of disymmetry here. Actually, if the observed value was not 8, but 12, some minor changes should be done (to keep some symmetry),

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="grey")
> abline(v=20-sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n-1,dbinom(n,size=20,prob=0.50),type="s")

Based on those observations, one can easily write a code to test if the $p_\star$-quantile of a sample is $x_\star$. Or not. For a two sided test, consider

> quantile.test=function(x,xstar=0,pstar=.5){
+ n=length(x)
+ T1=sum(x<=xstar)
+ T2=sum(x< xstar)
+ p.value=2*min(1-pbinom(T2-1,n,pstar),pbinom(T1,n,pstar))
+ return(p.value)}

Here, we have

> quantile.test(X)
[1] 0.5034447

Now, based on that idea, due to the duality between confidence intervals and tests, one can easily write a function that computes confidence interval for quantiles,

> quantile.interval=function(x,pstar=.5,conf.level=.95){
+ n=length(x)
+ alpha=1-conf.level
+ r=qbinom(alpha/2,n,pstar)
+ alpha1=pbinom(r-1,n,pstar)
+ s=qbinom(1-alpha/2,n,pstar)+1
+ alpha2=1-pbinom(s-1,n,pstar)
+ c.lower=sort(x)[r]
+ c.upper=sort(x)[s]
+ conf.level=1-alpha1-alpha2
+ return(list(interval=c(c.lower,c.upper),confidence=conf.level))}
> quantile.interval(X,.50,.95)
$interval [1] -0.3053884 0.7383247$confidence
[1] 0.9586105

Because of the use of non-asymptotic distributions, we can not get exactly a 95% confidence interval. But it is not that bad, here.

• Comparing quantiles for two samples

Now, to compare quantiles for two samples… it is more complicated. Exact tests are discussed in Kosorok (1999) (see http://bios.unc.edu/~kosorok/…) or in Li, Tiwari and Wells (1996) (see http://jstor.org/…). For the computational aspects, as mentioned in a post published almost one year ago on http://nicebread.de/… there is a function to compare quantiles for two samples.

> install.packages("WRS")
> library("WRS")

Some multiple tests on quantiles can be performed here. For instance, on the temperature, if we compare quantiles for Winter and Summer (on only 1,000 observations since it can be long to run that function), i.e. 5%, 25%, 75% and 95%,

> qcomhd(Z1[1:1000],Z2[1:1000],q=c(.05,.25,.75,.95))
q   n1   n2      est.1      est.2 est.1_minus_est.2     ci.low     ci.up     p_crit p.value signif
1 0.05 1000 1000 -6.9414084 -6.3312131       -0.61019530 -1.6061097 0.3599339 0.01250000   0.220     NO
2 0.25 1000 1000 -3.3893867 -3.1629541       -0.22643261 -0.6123292 0.2085305 0.01666667   0.322     NO
3 0.75 1000 1000  0.5832394  0.7324498       -0.14921041 -0.4606231 0.1689775 0.02500000   0.338     NO
4 0.95 1000 1000  3.7026388  3.6669997        0.03563914 -0.5078507 0.6067754 0.05000000   0.881     NO

or if we compare quantiles for Winter and Summer

> qcomhd(Z1[1:1000],Z3[1:1000],q=c(.05,.25,.75,.95))
q   n1  n2      est.1     est.2 est.1_minus_est.2     ci.low       ci.up     p_crit p.value signif
1 0.05 1000 984 -6.9414084 -6.438318        -0.5030906 -1.3748624  0.39391035 0.02500000   0.278     NO
2 0.25 1000 984 -3.3893867 -3.073818        -0.3155683 -0.7359727  0.06766466 0.01666667   0.103     NO
3 0.75 1000 984  0.5832394  1.010454        -0.4272150 -0.7222362 -0.11997409 0.01250000   0.012    YES
4 0.95 1000 984  3.7026388  3.873347        -0.1707078 -0.7726564  0.37160846 0.05000000   0.539     NO

(the following graphs are then plotted)

Those tests are based on the procedure proposed in Wilcox, Erceg-Hurn,  Clark and Carlson (2013), online on http://tandfonline.com/…. They rely on the use of bootstrap samples. The idea is quite simple actually (even if, in the paper, they use Harrell–Davis estimator to estimate quantiles, i.e. a weighted sum of ordered statistics – as described in http://freakonometrics.hypotheses.org/1755 – but the idea can be understood with any estimator): we generate several bootstrap samples, and compute the median for all of them (since our interest was initially on the median)

>  Q=rep(NA,10000)
>  for(b in 1:10000){
+  Q[b]=quantile(sample(X,size=20,replace=TRUE),.50)
+  }

Then, to derive a confidence interval (with, say, 95% confidence), we compute quantiles of those median estimates,

> quantile(Q,c(.025,.975))
2.5%     97.5%
-0.175161  0.666113

We can actually visualize the distribution of that bootstrap median,

> hist(Q)

Now, if we want to compare medians from two independent samples, the strategy is rather similar: we bootstrap the two samples – independently – then compute the median, and keep in mind the difference. Then, we will look if the difference is significantly different from 0. E.g.

> set.seed(2)
> Y=rnorm(50,.6)
> QX=QY=D=rep(NA,10000)
> for(b in 1:10000){
+ QX[b]=quantile(sample(X,size=length(X),replace=TRUE),.50)
+ QY[b]=quantile(sample(Y,size=length(Y),replace=TRUE),.50)
+ D[b]=QY[b]-QX[b]
+ }

The 95% confidence interval obtained from the bootstrap difference is

> quantile(D,c(.025,.975))
2.5%      97.5%
-0.2248471  0.9204888

which is rather close to was can be obtained with the R function

> qcomhd(X,Y,q=.5)
q n1 n2    est.1     est.2 est.1_minus_est.2    ci.low     ci.up p_crit p.value signif
1 0.5 20 50 0.318022 0.5958735        -0.2778515 -0.923871 0.1843839   0.05    0.27     NO

(where the difference is here the oppositive of mine). And when testing for 2 (or more) quantiles, Bonferroni method can be used to take into account that those tests cannot be considered as independent.

# Chercher, et bloguer

Je vais tenter de poursuivre la discussion amorcée dans un précédant billet, intitulé enseigner, et bloguer. Car je suis à la fois enseignant et chercheur. Et comme tous les universitaires, je suis tiraillé entre ces deux activités

(comme je m’en suis expliqué à plusieurs reprises, j’en suis ravi, cela crée une espèce de balance, cf. /1240 par exemple, il y a quelques mois). Mais à ces deux activités s’ajoute l’activité de blogger. Ou peut-être pas, en fait. Car c’est peut-être là l’erreur de raisonnement (qui fait qu’à plusieurs reprises – depuis plusieurs années – on me demande comment j’ai le temps de tout faire): autant enseigner et chercher sont – je pense – deux activités relativement disjointes (j’enseigne finalement assez peu mes thèmes de recherche), autant – comme j’essayais de le montrer dans mon précédant billet – bloguer s’insère dans ces deux activités.

Par exemple, l’autre jour, j’avais besoin d’un test de comparaison de quantiles (sur plusieurs échantillons), alors que je travaillais sur un papier de recherche. Malheureusement, ce n’est pas un problème standard (la plupart des textbook de statistique n’évoquent pas le sujet). J’ai ouvert mon cahier, et j’ai pris des notes, de choses trouvées dans divers papiers (le principal ayant été publié cette année…). Puis quand j’ai commencé à coder (c’était un peu le but), j’ai fermé mon cahier, et j’ai commencé à écrire un billet sur le blog (qui sera bientôt en ligne, sur /4199). Une raison un peu personnelle est que cela me permet de garder une trace de ce code, et de clarifier un peu ce que j’ai pu glaner et lire (j’ai des milliers de fichiers de codes sur mon ordinateur, et j’avoue m’y perdre un peu). Une autre raison est que, quitte à passer une journée, autant rentabiliser ma “journée recherche” (oui, je réserve certaines de mes journées à la recherche, en me coupant un peu des activités administratives, et de la préparation de mes cours) un produisant du contenu publié (sous une forme moins conventionnelle et reconnue qu’une publication, mais en l’occurrence, mon gribouillage ne méritait pas un article académique). Rentabiliser n’est peut-être pas le mot correct toutefois. Mais il présente l’avantage de traduire le fait que mon employeur (et divers organismes subventionnaires) me payent pour produire de la recherche Ce qui signifie que je dois certes chercher, mais aussi transmettre, à un moment ou à un autre. Je n’apporte peut-être pas quelque chose de très novateur (au sens où je me suis contenté d’aller chercher des articles théoriques qui se penchaient sur un problème semblable, je l’ai codé, et j’ai tenté de l”illustrer). Mais… c’est un peu le processus de l’activité de recherche ça, non ?

En fait, à mes yeux, la grande différence entre ce que je publie sur mon blog et ce que je publie (ou que je tente – souvent désespérément – de faire publier) dans des revues académiques, c’est que j’estime que ce que je publie ici n’a pas la prétention d’une publication académique. Au sens où je ne perdrais pas de temps à essayer de faire croire que je vais révolutionner le monde (comme on doit tous essayer de le faire dans un article académique). Pour exprimer le fond de ma pensée, je ne me battrais pas pour qu’un billet de blog soit cité dans une revue académique. Si quelqu’un utilise mes codes dans un article de recherche, je trouve légitime d’avoir un citation (sinon je pense qu’on peut parler de plagiat) mais je ne me battrais pas sur ce terrain car je pense le combat vain. Soyons lucide: si j’ai vraiment l’impression qu’un billet mérite plus, j’en ferais une version que je mettrais en ligne sur arxiv.org.  Cela dit,  dit, arxiv est à peine mieux qu’un blog: on ne peut – malheureusement – pas de contenter de publier en ligne sur arxiv . Souvenez-vous (ça vous est arrivé forcément un jour ou l’autre): comme tout le monde, vous lisez un article sur arxiv  (ou sur n’importe quel site de dépot de preprints, comme hal ou ssrn) et vous le citez, dans un article que vous soumettez. Si l’article est accepté, lors de la relecture, on vous demande “et cet article il n’a pas été publié depuis ?“. Et comme tout le monde, vous mettez à jour la référence, en citant la version publiée. Sans vraiment le lire (à la rigueur, on vérifie que le théorème 2.4 porte toujours le même nom… si on est plus courageux, que les hypothèses mentionnées n’ont pas changé depuis). Oui, sur ce point, le monde des l’édition est d’une hypocrisie infinie. Dans un monde idéal, les blog et les preprints seraient cités, mais on ne vit pas dans le monde des Bisounours (ou Calinours comme on dit à l’ouest de l’océan Atlantique, c’était la minute culturelle du billet)…

Mais je m’égare. La différence fondamentale est qu’un billet n’est qu’une ébauche d’article. Je m’autorise à mettre en ligne des morceaux de réflexions non abouties, sans véritable conclusion. Bref, mon blog est un carnet. Un carnet de recherche, et d’enseignement (c’est probablement pour cela que j’ai été ravi de le publier sur la plateforme hypotheses.org qui propose de publier des carnets de recherche, avec en plus un conseil scientifique). Je renvois aux (nombreux) billets de Mélodie Faury qui a beaucoup écrit – sur son blog infusoir.hypotheses.org sur l’utilité des carnets de recherche, à plusieurs billets de  parus sur reflexivites.hypotheses.org ou de  sur penseedudiscours.hypotheses.org, ou à l’incroyable blog blog.homo-numericus.net tenu par Marin Dacos et Pierre Mounier.

Tous ceux qui m’ont côtoyé ces dernières années le savent: je remplis des carnets. Des carnets à la fois de recherche, et d’enseignement.

(au risque de tuer un mythe, un chercheur 2.0 ça peut être ça: pas de téléphone cellulaire et un cahier pour faire des dessins)

Mais ces cahiers ne permettent de garder des traces que de réflexions, de lectures, de théorèmes, de notes prises pendant des exposés, mais pas de codes informatiques, peu d’algorithmes. Bref, le blog m’apporte un complément important à mes cahiers papiers. L’autre principale différence est que le blog est public, et mes cahiers (relativement) privés. Et le monde de la recherche ayant sa part de secrets et de mystères, j’avoue ne pas publier exactement c de billets exactement sur mes thèmes de recherche. Dans le billet sur les quantiles mentionné auparavant, je prends une illustration sur des séries de température, mais le papier de recherche (sur lequel je travaillais lorsque les questions sont apparus) porte sur un tout autre sujet. Oui, compte tenu de l’esprit de concurrence qui existe dans le milieu de la recherche, il convient de garder une part de secret. Mais on peut rester confidentiel sur ses activités de recherche tout en communiquant…

L’autre intérêt d’un blog est qu’il va plus loin qu’un carnet: si je recopie un article dans mon cahier, c’est une activité de recherche. Si je le recopie sur mon blog, c’est du plagiat. En fait, le blog oblige à restituer avec ses propres mots des choses lues ici ou là. Car il serait temps de l’avouer: si on blogue en tant qu’enseignant, on blogue pour ses élèves; alors que si on blogue en tant que chercheur, on blogue avant tout pour soi. Quand Djalil publie un billet intitulé Euclidean kernel random matrices in high dimension sur son (superbe) blog djalil.chafai.net, il ne blogue pour 100,000 lecteur. Enfin, je pense. Et lorsque je blogue sur des sujets pointus, je blogue pour moi. Pour me forcer à assimiler des notions, et à les réécrire. J’ai presque envie de dire que les lecteurs, je m’en fous: s’il y a des gens intéressés par mes billets, c’est tant mieux, mais si personne ne lit mes billets “de recherche”, ça ne changera pas grand chose…

Maintenant, pour revenir sur la confidentialité, si elle est de mise tant que le papier n’est pas fini, il convient de noter qu’une fois le papier fini, c’est le contraire: il faut le faire savoir, haut et fort. Et là aussi, je dois avouer que le blog a de grandes vertus. Il donne une exposition incroyable,et permet de faire connaître ses travaux. Par exemple, en 6 semaines, 710 ont lu l’article sur la température à Montréal sur  mon blog /2909. Sans compter ceux qui ont lu la version reprise par les agrégateurs, comme r-bloggers.com/… ou architects.dzone.com/… (plus de 1300 personnes l’ont vu sur ce dernier) ou via un flux RSS. Mais en plus des billets sur mon blog, j’ai l’impression que mes articles académiques sont aussi beaucoup consultés (dans un moindre mesure, bien sûr !). Je suis toujours surpris par mon h-index, qui me semble plus élevés que beaucoup de chercheurs biens meilleurs que moi. Et je pense que le blog aide. Je pense que des papiers sont cités par beaucoup de monde parce qu’ils apparaissent tôt dans les recherches sur Google, mes blogs ayant tous été très bien référencés. Une étude (http://arxiv.org/0212043 mentionnée sur http://washingtonmonthly.com/…) disait que les articles cités étaient souvent peu lus, et effectivement, pour être cité, il suffit d’être bien référencé. Aussi, alors que j’ai longtemps cru que tenir un blog actif était contre-productif pour mes activités de recherche (et effectivement, le temps passé à écrire dans le blog se fait parfois au détriment d’autres taches) j’en suis de moins en moins convaincu. Je pense que des papiers ont été cités parce qu’ils sont mentionnés sur le blog. Soit parce j’y mets les textes en ligne (même si je n’en avais pas vraiment le droit), soit parce qu’ils étaient – juste – bien référencés, soit – en étant moins cynique – parce que pour certains papiers, je mets de codes en ligne pour mettre en œuvre les techniques proposées (je pense au papier sur le lissage spatial, /2416, ou tous les billets qui ont entouré celui sur la modélisation de la canicule d’août 2003, tag/heatwave).

Je me suis aussi rendu compte que j’ai gagné du temps dans des papiers car j’aboutissais sur des concepts connexes que j’avais pu aborder sur mon blog (souvent pour m’amuser). Par exemple, dans nos travaux sur les tremblements de terres avec Mathieu, on manipule des shapefiles (/tag/binar). Et je suis très content d’avoir joué un temps sur mon blog sur des problèmes de pistolets à eaux, qui m’ont permis de manipuler certaines fonctions sur les polygones, cf /tag/water. Enfin, plusieurs invitations à des séminaires de rechercher ou des conférences académiques ont été faites parce que le blog donne une relative notoriété, et visibilité.

Histoire de conclure en faisant une boucle, mon billet sur l’enseignement (et l’activité de blogger) venait d’une demande spécifique de Robin. Or Robin avait fait parti d’un projet génial, à savoir le blog statisfaction.wordpress.com lancé il y a quelques années par des doctorants (et post-doctorants) en statistique, sur Paris. L’idée était géniale car elle leur donnait une visibilité incroyable dans la communauté scientifique, avec des informations sur les séminaires, les cours, et leurs activités de recherche. Il n’est pas étonnant de voir autant de blogs scientifiques tenus pas des doctorants – ou lancés alors qu’ils étaient doctorants – je pense par exemple à Corey à Montréal, bayesianbiologist.com ou a plusieurs blogs d’économie, comme le défunt (ou disons resté très longtemps inactif) ecopublix.eu ou mafeco.fr, par Emmeline et  Jean-Edouard, ou encore rationalitelimitee.wordpress.fr lancé par Cyril quand il était encore doctorant. Les doctorants étant les principaux producteurs de recherche, il est logique de les voir apparaître aussi clairement dans la communauté des bloggers tenant des blogs de recherche (pour ne pas dire des carnets de recherche).

Mais bon, il serait temps de clôturer la séries des billets bloguer, c’est génial (surtout que je me rend compte que mon billet part un peu dans tous les sens, et n’a pas la profondeur d’analyse que l’ont peut retrouver dans les blogs d’ hypotheses.org mentionnés auparavant) La prochaine fois, promis, un billet sur bloguer, c’est nul !