Tag Archives: significance

Regression

Estimates on training vs. validation samples

Leave a comment

Before moving to cross-validation, it was natural to say “I will burn 50% (say) of my data to train a model, and then use the remaining to fit the model”. For instance, we can use training data for variable selection (e.g. using some stepwise procedure in a logistic regression), and then, once variable have been selected, fit the model on the remaining set of observations. A natural question is usually “does it really matter ?”.

In order to visualize this problem, consider my (simple) dataset

MYOCARDE=read.table(
  "http://freakonometrics.free.fr/saporta.csv",
  head=TRUE,sep=";")

Let us generate 100 training samples (where we keep about 50% of the observations). On each of them, we use a stepwise procedure, and we keep the estimates of the remaining variables (and their standard deviation actually)

n=nrow(MYOCARDE)
M=matrix(NA,100,ncol(MYOCARDE))
colnames(M)=c("(Intercept)",names(MYOCARDE)[1:7])
S1=S2=M1=M2=M
for(i in 1:100){
idx = which(sample(0:1,size=n, replace=TRUE)==1)
reg=step(glm(PRONO=="DECES"~.,data=MYOCARDE[idx,]))
nm=names(reg$coefficients)
M1[i,nm]=reg$coefficients
S1[i,nm]=summary(reg)$coefficients[,2]
f=paste("PRONO=='DECES'~",paste(nm[-1],collapse="+"),sep="")
reg=glm(f,data=MYOCARDE[-idx,])
M2[i,nm]=reg$coefficients
S2[i,nm]=summary(reg)$coefficients[,2]
}

Then, for the 7 covariates (and the constant) we can look at the value of the coefficient in the model fitted on the training sample, and the value on the model fitted on the validation sample (of course, only when they were remaining)

for(j in 1:8){
idx=which(!is.na(M1[,j]))
plot(M1[idx,j],M2[idx,j])
abline(a=0,b=1,lty=2,col="gray")
segments(M1[idx,j]-2*S1[idx,j],M2[idx,j],M1[idx,j]+2*S1[idx,j],M2[idx,j])  
segments(M1[idx,j],M2[idx,j]-2*S2[idx,j],M1[idx,j],M2[idx,j]+2*S2[idx,j])  
}

For instance, with the intercept, we have the following

 

where horizontal segments are confidence intervals of the parameter on the model fitted on the training sample, the vertical on the validation sample. The green part means some sort of consistency, while the red one means that actually, the coefficient was negative with one model, positive with the other one. Which is odd (but in that case, observe that coefficients are rarely significant).

We can also visualize the joint distribution of the two estimators,

for(j in 1:8){
library(ks)
idx = which(!is.na(M1[,j]))
Z = cbind(M1[idx,j],M2[idx,j])
H = Hpi(x=Z)
fhat = kde(x=Z, H=H)
image(fhat$eval.points[[1]],
fhat$eval.points[[2]],fhat$estimate)
abline(a=0,b=1,lty=2,col="gray")
abline(v=0,lty=2)
abline(h=0,lty=2)
}

which are here, almost on the diagonal,

meaning that the intercept on the two samples is (more or less) the same. We can then look at other parameters (which is actually more interesting).

On that variable, it seems that it is significant on the training dataset (somehow, it is consistent with the fact that it is remaining in the model after the stepwise procedure) but not on the validation sample (or hardly significant).

Others are much more consistent (with some possible outliers)

 

 

On the next one, we have again significance on the training sample, but not on the validation sample,

 

 

and probably more interesting

where the two are very consistent.

M1-statistique

Tests, Power and Significance

3 Comments

In the mathematical statistics course today, we started talking about tests, and decision rules. To illustrate all the concepts introduced today, we considered the case where we have a sample  with . And we want to test

  against 

In the course, we’ve seen that we could use a test based on the order statistics .  The test would be

i.e. if  we choose , and if , we choose .

From the definition of the first order risk,

we can easily get that

Thus, the power is then

To visualize it, use the following parameters

n=5
alpha=.1
theta0=1

Then

C1=theta0*(1-alpha)^(1/n)
theta=seq(0,2,by=.01)
P1=(1-(theta0/theta)^n*(1-alpha))*(theta>C1)
plot(theta,P1,type="l",lwd=2,col="blue",xlab="",ylab="Power")

Note that, so far, we did never consider the maximum of our sample. Assume that the maximum is , then we can compute the -value,

Here it is

PV=(1-theta^n)*(theta<=1)
plot(theta,PV,type="l",lwd=2,col="blue",xlab="",ylab="p-value")

Now, why not consider another test, based on the minimum (since we have the distribution of the minimum of a sample from a uniform distribution). The test is the same as before

but here, the threshold is

The power of the test is here

This test has the same significance level (by construction), but the power of the test is clearly lower than the one we got using the maximum of our sample, when 

C2=theta0*(1-alpha^(1/n))
P2=(1-(theta0/theta)*(1-alpha^(1/n)))^n*(theta>C2)
lines(theta,P2,type="l",lwd=2,col="red")

Why not consider a test based on ? The problem is that we need the distribution (more specifically the survival distribution) of . We can compute it, numerically. But that might be painful. An alternative is to consider some approximation, based on the central limit theorem, i.e.

Our test is based on , and to get the same significance as before, use

The power of the test is then

Here it is

mu=2*(theta0/2)
s2=2^2*(theta0^2/12)/n
C3=qnorm(1-alpha,mu,sqrt(s2))
(P=1-pnorm(C3,theta,sqrt(s2)))*(theta>C3)
lines(theta,P)

Observe here that the test based on the maximum is not more powerful than the one based on the average (I just wonder if it could be due to the Gaussian approximation…).

Statistics

Un peu plus près des étoiles (***)

1 Comment

Il y a eu un gros buzz, récement autour du papier de Valen Johnson paru dans PNAS. L’article a été repris un peu partout (http://nature.com/news/, http://blogs.scientificamerican.com/absolutely-maybe/, http://arstechnica.com/science/ ou encore http://passeurdesciences.blog.lemonde.fr/ qui a repris l’information, en français). Et plusieurs personnes m’ont fait suivre des liens, en me demandant mon avis, par courriel ou via twitter. Je ne vais pas revenir sur l’étude (pour l’instant) ni sur les mauvaises lectures de l’étude, mais plutôt sur le buzz qui a eu lieu autour. Car manifestement, on ne vit pas dans la même temporalité avec les journalistes. Par ‘on‘ j’entends les scientifiques au sens large. Je veux dire par là que quand j’ai écris mon billet, il y a quelques semaines sur l’origine de la -value, ça m’a pris du temps. Et depuis, je travaille sur les suites du billet (en particulier la lecture bayésienne des tests, mais aussi les tests multiples). Faire une revue de littérature, se replonger sur l’origine d’un concept, ça prend du temps. Et je le prendrais. Je l’avais fait il y a quelques mois sur les martingales et la notion d’efficience de marché, en revenant sur les travaux d’Eugène Fama (dont certains semblaient alors ignorer l’existence). Heureusement, c’était la fin de ma session, et j’avais envie de me changer les idées.

Pour revenir rapidement sur l’article paru dans PNAS (et c’est pour cela que je fais un parallèle avec mon billet sur les martingales), il n’y a rien de bien nouveau pour ceux qui ont suivi la statistique depuis une soixantaine d’années (depuis les travaux de Fisher, en fait). Je veux dire par là qu’il y a eu des dizaines, des centaines d’articles pour critiquer l’usage systématique de la -value. On pourra penser au papier de Ziliak et McCloskey, The Cult of Statistical Significance par exemple, ou surtout le livre The Significance Test Controversy par Denton Morrison et Ramon Henkel, publié en 1970, et toutes les références mentionnées. Comme le populaire http://psychclassics.yorku.ca/Rozeboom/ de William Rozeboom, publié en 1960, avec le paragraphe que tout le monde connait,

ce paragraphe est souvent cité par Bruno Lecoutre, à qui on doit le joli Uses, Abuses and Misuses of Significance Tests in the Scientific Community: Won’t the Bayesian Choice be Unavoidable? (tiens, ce n’est pas ce que dit le papier dans PNAS ?). Si je ne devais mentionner que deux références à relire en ces temps de buzz, je penserais à Harold Jefreys qui publiait en 1961 son Theory of Probabilities, dont toute une section parle justement de la significativité

Et en 1982, Sir David Cox publiait Statistical Significance Tests, avec là aussi plusieurs paragraphes qu’il serait bon de relire,

Quant à l’alternative Bayésienne1, je renvoie au très joli papier de Bradley Efron, publié en 1998, sur Fisher in the 21st century. Bref, critiquer l’utilisation abusive de la -value n’est pas nouveau, loin de là. Et quand je dis abusive… @mc_hankins a posté sur twitter plusieurs exemples d’utilisations de -value, comme

suivi quelques lignes plus loin par un

Manifestement, pour certains, 5%, c’est 5% ! C’est un peu pour comprendre ce fanatisme que j’avais fait mon précédant billet…

Maintenant, pour revenir à mon premier paragraphe, et au buzz autour de cet article, je suis convaincu que les journalistes et les scientifiques ne vivent pas dans la même temporalité. Pas seulement dans le temps que l’on met pour écrire un article, mais aussi sur le fait que dans mon étude sur les martingales, j’ai relu Fama, Samuelson ou encore LeRoy, des papiers vieux de 20 ou 30 ans. Pour revenir sur la notion de significativité, et comprendre le mythe des 5%, je me suis replongé dans les travaux de Fisher, qui datent de 60 ou 70 ans. Sans émettre de jugement de valeurs, on relit les scientifiques 20, 30 ou 50 ans après, mais je relis rarement un journaliste 5 ans après2 . Et c’est normal, je pense. Je veux dire par là que les journalistes sont là pour informer, pas éduquer. C’est aussi pour ça que je crois que les blogs académiques ont une raison d’être, car la nature a horreur du vide… Il manque de lieux pour apprendre des choses (et pas seulement se tenir informer) et les blogs académiques sont parfait pour remplir cette mission.

Mais je m’égare car je voulais parler du buzz récent, autour de la significativité statistique, et faire deux commentaires.

Le premier est que je trouve que la statistique telle que Fisher la pratiquait a pris trop de place dans l’enseignement de la statistique (là dessus, je suis assez d’accord avec les critiques). En particulier les résultats asymptotiques. Et il serait temps d’enseigner davantage la statistique bayésienne, le plus tôt possible. Je n’ai pas suivi de cours de statistique bayésienne, et je le regrette. Cela dit, dans de nombreux domaines (je pense en particulier à des conférences en hydrologie auxquelles j’ai participé) tout le monde fait du bayésien, et de manière naturelle. Mais c’est un long débat, et je ne vais pas le commencer maintenant.

Mon second est une forme de réponse à cette critique sur l’importance de la significativité. Regardons un peu l’évolution des logiciels de statistique pour mieux comprendre. Un des plus anciens logiciels pour faire des régression, TSP, a été créé au milieu des années 60. Par exemple, dans la version TSP 4.1 datant de 1987, une régression, ça ressemblait à ça (à en croire Renfro (2009))

   METHOD OF ESTIMATION = ORDINARY LEAST SQUARES

   DEPENDENT VARIABLE : DIST

            SUM OF SQUARED RESIDUALS = 11353.52
    STANDARD ERROR OF THE REGRESSION = 15.37959
          MEAN OF DEPENDENT VARIABLE = 42.98
                  STANDARD DEVIATION = 25.76938
                           R-SQUARED = 0.6510794
                  ADJUSTED R-SQUARED = 0.6438102
             DURBIN-WATSON STATISTIC = 1.6762
                   F-STATISTIC(1,48) = 89.56711
          LOG OF LIKELIHOOD FUNCTION = -206.5784
              NUMBER OF OBSERVATIONS = 50

  VARIABLE      ESTIMATED  STANDARD 
              COEFFICIENT     ERROR   T-STATISTIC
     C          -17.5791     6.7584      -2.601  
     SPEED        3.9324    0.4155        9.464

car en 1987, je faisais certes de l’informatique, mais les régressions linéaires m’intéressaient peu. Bref, je n’ai jamais pu tester ce logiciel, et j’en crois ce que j’en ai lu. Mais toutes les sorties de cette époque semblaient ressembler à celle-ci. Aujourd’hui, avec R, la même régression se lit

    lm(formula = dist ~ speed, data = cars)

    Residuals:
       Min      1Q  Median      3Q     Max 
   -29.069  -9.525  -2.272   9.215  43.201 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791  6.7584  -2.601   0.0123 *  
speed         3.9324  0.4155   9.464 1.49e-12 ***
---
      Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Residual standard error: 15.38 on 48 degrees of freedom
    Multiple R-squared:  0.6511,	Adjusted R-squared:  0.6438 
    F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Vous avez-vu la différence ? Regardez tranquillement, on a le temps…

Oui… les étoiles ! Les fameuses étoiles ! Quand j’enseigne les méthodes de régression, la seule chose qui intéresse et fascine les étudiants, ce sont les étoiles. C’est expliqué sous le tableau, les étoiles, c’est la significativité. En dessous de 5% de significativité, on a le droit a une étoile, et en dessous de 0.1%, on a carrément le droit à trois étoiles ! Dois-je reprocher à mes étudiants de ne retenir que cette histoire stellaire de mon cours de prévision (predictive modeling comme on dit) ? De plus en plus, on veut des formations professionnelles, orientées vers la pratique, vers les vraies applications. Alors dans les cours, on enseigne des recettes de cuisine… et une d’entre elles est basée sur le fait que pour qu’un modèle soit valide, avec des variables significatives, il faut garder les variables avec des étoiles. Vivement un cours où on aura le temps d’apprendre des trucs inutiles quand on va ensuite travailler en entreprise, des trucs qui prennent du temps à être compris (oui, les concepts, c’est souvent subtil, c’est plus difficile à expliquer qu’une recette simple: pas d’étoile, tu recommences), un cours où on prendrait le temps de s’interroger sur le sens des recettes que l’on apprend par cœur. Parce que mine de rien, toute une génération de personnes qui ont suivi des cours de stats n’ont retenu qu’un truc: regarder les étoiles. Et les gens qui regardent les étoiles, c’est dur de leur demander du jour au lendemain de regarder ailleurs…

1.pour les plus courageux, qui liront rapidement le billet, je rappelle que vendredi matin, Xiao-Li Meng donnera un exposé depuis le Fields Institute à Toronto, intitulé sobrement “Who is crazier: Bayes or Fisher?” et l’exposé sera retransmis en direct sur http://fields.utoronto.ca.
2.bon, je devrais peut-être nuancer un peu, car j’ai eu un grand plaisir à découvrir il y a quelques années les écrits d’Albert Londres (surtout au bagne ou en Afrique, avec terre d’ébène) ou de Georges Orwell (mais je ne sais pas si ses témoignage sont des essais ou du journalisme…)