# p-hacking, or cheating on a p-value

Yesterday evening, I discovered some interesting slides on False-Positives, p-Hacking, Statistical Power, and Evidential Value, via  ‘s post on Twitter. More precisely, there was this slide on how cheating (because that’s basically what it is) to get a ‘good’ model (by targeting the p-value)

As mentioned by @david_colquhoun  one should be careful when reading the slides : some statistician might have a heart attack when they read

But still, there are interesting points in that slide.

# Bias of Hill Estimators

In the MAT8595 course, we’ve seen yesterday Hill estimator of the tail index. To be more specific, we did see see that if $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, then Hill estimators for $\alpha$ are given by

$\widehat{\alpha}_k = \left[\frac{1}{k}\sum_{i=0}^{k-1} \log X_{n-i,n} -\log X_{n-k,n}\right]^{-1}$
for $k\in\{1,2,\cdots,n\}$. Then we did say that $\widehat{\alpha}_k$ satisfies some consistency in the sense that $\widehat{\alpha}_k \overset{\mathbb{P}}{\rightarrow} \alpha$ if $k\rightarrow\infty$, but not too fast, i.e. $k/n\rightarrow0$ (under additional assumptions on the rate of convergence, it is possible to prove that $\widehat{\alpha}_k \overset{a.s.}{\rightarrow} \alpha$). Further, under additional technical conditions

$\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}(0,\alpha^2)$

In order to illustrate this point, consider the following code. First, let us consider a Pareto survival function, and the associated quantile function

> alpha=1.5
> S=function(x){ifelse(x>1,x^(-alpha),1)}
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root} The code here is obviously too complicated, since this power function can easily be inverted. But later on, we will consider a more complex survival function. Here are the survival function, and the quantile function, > u=seq(0,5,by=.01) > plot(u,Vectorize(S)(u),type="l",col="red") > u=seq(0,99/100,by=.01) > plot(u,Vectorize(Q)(u),type="l",col="blue",ylim=c(0,20)) Here, we need the quantile function to generate a random sample from this distribution, > n=500 > set.seed(1) > X=Vectorize(Q)(runif(n)) Hill plot is here > library(evir) > hill(X) > abline(h=alpha,col="blue") We can now generate thousands of random samples, and see how those estimators behave (for some specific $k$‘s). > ns=10000 > HillK=matrix(NA,ns,10) > for(s in 1:ns){ + X=Vectorize(Q)(runif(n)) + H=hill(X,plot=FALSE) + hillk=function(k) H$y[H$x==k] + HillK[s,]=Vectorize(hillk)(15*(1:10)) + } and if we compute the average, > plot(15*(1:10),apply(HillK,2,mean) we do get a series of estimators that can be considered as unbiased. So far, so good. Now, recall that being in the max-domain of attraction of the Fréchet distribution does not mean that $\overline{F}(x)=C x^{-\alpha}$, with $\alpha>0$, but is means that $\overline{F}(x)= x^{-\alpha} \mathcal{L}(x)$ for some slowly varying function $\mathcal{L}$, not necessarily constant! In order to understand what could happen, we have to be slightly more specific. And this can be done only by looking at second order regular variation property of the survival function. Assume, here that there is some auxilary function $a$ such that $\lim_{t\rightarrow\infty}\frac{\overline{F}(xt)/\overline{F}(t)-x^{-\alpha}}{a(t)}=x^{-\alpha}\frac{1-x^{-\beta}}{\beta}{}$ This (positive) constant $\beta$ is – somehow – related to the speed of convergence of the ratio of the survival functions to the power function (see e.g. Geluk et al. (2000) for some examples). To be more specific, assume that $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ then, the second order regular variation property is obtained using $a(t)=\beta t^{-\beta}$, and then, if $k$ goes to infinity too fast, then the estimator will be biased. More precisely (see Chapter 6 in Embrechts et al. (1997)), if $k=O(n^{2\beta/(\alpha+2\beta)})$, then, for some $\lambda>0$, $\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal L}{\rightarrow}\mathcal{N}\left(\frac{\alpha^3}{\beta-\alpha}\lambda,\alpha^2\right)$ The intuitive interpretation of this result is that if $k$ is too large, and if the underlying distribution is not exactly a Pareto distribution (and we do have this second order property), then Hill’s estimator is biased. This is what we mean when we say • if $k$ is too large, $\widehat{\alpha}_k$ is a biased estimator • if $k$ is too small, $\widehat{\alpha}_k$ is a volatile estimator (the later comes from properties of a sample mean: the more observations, the less the volatility of the mean). Let us run some simulations to get a better understanding of what’s going on. Using the previous code, it is actually extremly simple to generate a random sample with survival function $\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot x^{-\alpha}$ > beta=.5 > S=function(x){+ ifelse(x>1,.5*x^(-alpha)*(1+x^(-beta)),1) } > Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}

If we use the code above. Here, with

> n=500
> set.seed(1)
> X=Vectorize(Q)(runif(n))

the Hill plot becomes

> library(evir)
> hill(X)
> abline(h=alpha,col="blue")

But it’s based on one sample, only. Again, consider thousands of samples, and let us see how Hill’s estimator is behaving,

so that the (empirical) mean of those estimator is

# Likelihood Based Methods, for Extremes

This week, in the MAT8595 course, we will start the section on inference for extreme values. To start with something simple, we will use maximum likelihood techniques on a Generalized Pareto Distribution (we’ve seen Monday Pickands-Balkema-de Hann theorem).

• Maximum Likelihood Estimation

In the context of parametric models, the standard technique is to consider the maximum of the likelihood (or the log-likelihod).Let $\boldsymbol{\theta}$ denote the parameter (with $\boldsymbol{\theta}\in\Theta\subset\mathbb{R}^k$). Given some – stnardard – technical assumptions, such as $\boldsymbol{\theta}_0\notin\partial\Theta$, or $f(x;\boldsymbol{\theta})­>0$ on some neighbourhood of $\boldsymbol{\theta}_0$, then

$\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0\big)\ \xrightarrow{\mathcal{L}}\ \mathcal{N}(0,\,I^{-1})$

where $I$ denotes Fisher information matrix (see any textbook for mathematical statistics courses). Consider here some i.i.d. sample, from a Generalized Pareto Distribution, with parameter $\boldsymbol{\theta}=(\xi,\sigma)$, so that

$F_{(\xi,\sigma)}(x) = \begin{cases} 1 - \left(1+ \frac{\xi x}{\sigma}\right)^{-1/\xi} &, \xi \neq 0 \\ 1 - \exp \left(-\frac{x}{\sigma}\right) &, \xi = 0 \end{cases}$

If we solve (numerically) the first order condition of the maximum likelihood, we get an estimator  $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ which satisfies

$\sqrt{n}\left(\left[\begin{array}{c}\widehat{\xi}_n\\\widehat{\sigma }_n\end{array}\right]-\left[\begin{array}{c}\xi_0\\\sigma_0 \end{array}\right]\right)\rightarrow \mathcal{N}\left(\left[\begin{array}{c}0\\0\end{array}\right],\left[\begin{array}{cc}(1+\xi_0)^2 & \sigma_0[1+\xi_0]\\ \sigma_0 [1+\xi_0] & 2\sigma^2_0(1+\xi_0) \end{array}\right]\right)$

The idea of this asymptotic normality is the following : if the true distribution of the sample is a GPD with parameter $\boldsymbol{\theta}_0$, then, if $n$ is large enough, then $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ will have a joint normal distribution. So if we generate a lot of sample (sufficently large, say 200 observations), then the scatterplot of the estimator should the same as the scatterplot of a Gaussian distribution,

> library(evir)
> n=200
> param=matrix(NA,1000,2)
> for(s in 1:1000){
+ x=rgpd(n,xi=1/1.5,beta=1)
+ param[s,]=gpd(x,0)$par.ests + } > m=apply(param,2,mean) > S=var(param) > library(mnormt) > x=seq(min(param[,1])-.05,max(param[,1])+.05,length=101) > y=seq(min(param[,2])-.05,max(param[,2])+.05,length=101) > vx=rep(x,each=length(y)) > vy=rep(y,length(x)) > vz=dmnorm(cbind(vx,vy),m,S) > z=matrix(vz,length(y),length(x)) > COL=rev(heat.colors(100)) > image(x,y,z,col=COL) > points(param) and to get a 3d representation > x=seq(min(param[,1])-.05,max(param[,1])+.05,length=31) > y=seq(min(param[,2])-.05,max(param[,2])+.05,length=31) > vx=rep(x,each=length(y)) > vy=rep(y,length(x)) > vz=dmnorm(cbind(vx,vy),m,S) > z=matrix(vz,length(y),length(x)) > persp(x,y,t(z),shade=TRUE,col="green",theta=-30,phi=20,ticktype="detailed", + xlab="xi",ylab="sigma") With 200 observations, if the true underlying distribution is a GPD, then, indeed, the joint distribution of $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ seems to be normal. That would be interesting to generate some confidence intervals for instance, or define some tests. To go further, see any standard textbook on statistical mathematics, e.g. Casella & Berger (2002). • Delta Method Another important property is the so called delta-method (we’ve seen Monday in class that it was obtained easily using a first order Taylor expansion). The idea is that if $\widehat{\boldsymbol{\theta}}_n$ is asymptotically normal, and if $h:\Theta\subset\mathbb R^k \rightarrow \mathbb R^m$ is sufficently smooth, then $h(\widehat{\boldsymbol{\theta}}_n)$ will also be asymptotically Gaussian. More precicely (see also the header of this blog) $\sqrt{n}\big(h(\widehat{\boldsymbol{\theta}}_n)-h(\boldsymbol{\theta}_0)\big)\ \xrightarrow{\mathcal{L}}\ \mathcal{N}(0,\,\nabla h((\boldsymbol{\theta}_0)'\cdot I^{-1}\cdot \nabla h((\boldsymbol{\theta}_0))$ From this property, we can get the normality of $\widehat{\alpha}_n=\widehat{\xi}_n^{-1}$ (which is another parametrization used in extreme value models), or on any quantile, $\widehat{Q}_u=F^{-1}_{\widehat{\boldsymbol{\theta}}_n}(u)=h_u(\widehat{\xi}_n,\widehat{\sigma}_n)$. Let us run some simulation, one more time to check that we actually have a joint normality. > library(evir) > n=200 > param=riskm=matrix(NA,1000,2) > for(s in 1:1000){ + x=rgpd(n,xi=1/1.5,beta=1) + param[s,]=gpd(x,0)$par.ests
+ xihat=param[s,1]
+ shat=param[s,2]
+ q=shat * (.01^(-xihat) - 1)/xihat
+ tvar=q+(shat + xihat * q)/(1 - xihat)
+ riskm[s,]=c(1/xihat,q)
+ }
> m=apply(riskm,2,mean)
> S=var(riskm)
> library(mnormt)
> x=seq(min(riskm[,1])-.05,max(riskm[,1])+.05,length=101)
> y=seq(min(riskm[,2])-.05,max(riskm[,2])+.05,length=101)
> vx=rep(x,each=length(y))
> vy=rep(y,length(x))
> vz=dmnorm(cbind(vx,vy),m,S)
> z=matrix(vz,length(y),length(x))
> image(x,y,t(z),col=COL)
> points(riskm)

As we can see bellow, with samples of size 200, we cannot use this asymptotical result: it looks like we do not have enought data. But if we run the same code with

> n=5000

We get the joint normality of $\widehat{\alpha}_n$ and $\widehat{Q}_n(u)$. This is what we can get from this result, called delta-method in statistical textbooks. See again Casella & Berger (2002) for more details.

• Profile Likelihood

Another interesting tool is the concept of profile likelihood. This would be interesting here since the main interest is the tail index $\xi$, $\sigma$ being here some kind of auxilary parameter. See Venzon & Moolgavkar (1988) for more details. Here, we will plot

But more generally, it is possible to consider

where  is the set of interesting parameters. Then (under standard suitable conditions) we can prove that

which can be used to derive confidence intervals. In the GPD case, for each $\xi$, we have to find an optimal $\sigma^\star(\xi)$. We compute the (profile) likelihood i.e. $\mathcal{L}(\xi,\sigma^\star(\xi))$. And we can compute the maximum of this profile likelihood. This two-stage optimization is, in general, not equivalent with the (global) maximization of the likelihood, as computed below

>  n=500
>  set.seed(1)
>  x=rgpd(n,xi=1/1.5,beta=1)
>  loglikelihood=function(xi,beta){
+  sum(log(dgpd(x,xi,mu=0,beta))) }
>  XIV=(1:300)/100;L=rep(NA,300)
>  for(i in 1:300){
+  XI=XIV[i]
+  profilelikelihood=function(beta){
+  -loglikelihood(XI,beta) }
+  L[i]=-optim(par=1,fn=profilelikelihood)$value } > plot(XIV,L,type="l") > XIV[which.max(L)] [1] 0.67 > gpd(x,0)$par.ests
xi      beta
0.6730145 0.9725483

We are not far away. Actually, if we want to compute the maximum of the profile likelihood (and not only compute the values of the profile likelihood on a grid, as before), we use

>  PL=function(XI){
+  profilelikelihood=function(beta){
+  -loglikelihood(XI,beta) }
+  return(optim(par=1,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(0,3)))$minimum
[1] 0.6731025

$objective [1] 822.5574 Observe that, indeed, we are not far away from the maximum likelihood estimator of $\xi$ (I believe that it’s mainly a computational issue here, and theat the two are similar, here… actually, I’d be glad to hear about cases where maximum of the profile likelihood is not the same as the maximum of the likelihood). The interesting point is that we can use this technique to compute a confidence interval, and even visualize it on a graph > up=OPT$objective
>  abline(h=-up)
>  abline(h=-up-qchisq(p=.95,df=1),col="red")
>  I=which(L>=-up-qchisq(p=.95,df=1))
>  lines(XIV[I],rep(-up-qchisq(p=.95,df=1),length(I)),
+  lwd=5,col="red")
>  abline(v=range(XIV[I]),lty=2,col="red")

The vertical lines are the lower and the upper bound of a 95% confidence interval for parameter $\xi$.

To go further, see Murphy, S.A & van der Vaart, A.W. (2000). On Profile Likelihood.

# Pariwise and global dependence

I was just contacted by some researchers willing to test if a multivariate copula is – or not – Gaussian. They use a test proposed in an unpublished paper by Malevergne and Sornette, stating that one should simply test for pairwise normality. In this paper, the following result is mentioned

In other words, it is stated that (or at least something closed to this)

Unfortunately, this result is (probably) not correct (and if it is valid, it is nontrivial). It should be possible to construct a counterexample (thanks to Roger Nelsen (web) for the idea) by letting the correlation in each pair be close to -1 (i.e. global pairwise contercomonotonicity). Then the correlation matrix of the triplet would fail to be positive definite (which is a requierment for Gaussian vectors). This idea can be related to the “compatibility” condition in Joe (1997).

My colleague Mihai Gradinaru (web) mentioned the book by Wlodzimierz Bryc, Normal distribution: characterizations with applications, but I did not find any result in it (I looked very quickly). Anyway, I am still working on this… any comment is welcome !