Tag Archives: test

INF7100, statistiques

La seconde partie de mon intervention sur la science des données, dans le cadre du cours INF7100 portera sur les statistiques, univariées et multivarirées. Le plan sera le suivant

  • 201: De la Statistique aux Sciences de Données pdf video (14:24)
  • 211: Fonctions Usuelles en Statistique (fonction de réparition, densité, histogramme) pdf video (28:37)
  • 221: Indicateurs Statistique: Valeur Centrale (moyenne) pdf video (32:56)
  • 222: Indicateurs Statistique: Dispersion (variance, inégalités) pdf video (22:21)
  • 223: Indicateurs Statistique: Approximations (approximation normale) pdf video (18:42)
  • 224: Indicateurs Statistique: Quantiles pdf video (24:54)
  • 231: Inférence (statistique bayésienne) pdf video (39:33)
  • 241: Tests Statistiques (1) (tests, significativité, p-value) pdf video (43:41)
  • 242: Tests Statistiques (2) (erreurs) pdf video (16:51)
  • 261: Statistiques Bivariées pdf video (25:16)
  • 271: Statistiques Multivariées: Projections pdf video (29:06)
  • 272: Statistiques Multivariées: Clusters pdf video (32:21)
  • 281: Réseaux et Graphs pdf video (32:40)
  • 291: Données Chronologiques pdf video (29:01)

 

Probabilistic Foundations of Econometrics, part 4

This post is the fourth one of our series on the history and foundations of econometric and machine learning models. Part 3 is online here.

Goodness of Fit, and Model

In the Gaussian linear model, the determination coefficient – noted R^2 – is often used as a measure of fit quality. It is based on the variance decomposition formula \underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total variance}}=\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\text{residual variance}}+\underbrace{\frac{1}{n}\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2}_{\text{explained variance}} The R^2 is defined as the ratio of explained variance and total variance, another interpretation of the coefficient that we had introduced from the geometry of the least squares R^2= \frac{\sum_{i=1}^n (y_i-\bar{y})^2-\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}The sums of the error squares in this writing can be rewritten as a log-likelihood. However, it should be remembered that, up to one additive constant (obtained with a saturated model) in generalized linear models, deviance is defined by {Deviance}(\widehat{\beta}) = -2\log[\mathcal{L}] which can also be noted Deviance(\widehat{\mathbf{y}}). A null deviance can be defined as the one obtained without using the explanatory variables \mathbf{x}, so that \widehat{y}_i=\overline{y}. It is then possible to define, in a more general context (with a non-Gaussian distribution for y)R^2=\frac{{Deviance}(\overline{y})-{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}=1-\frac{{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}However, this measure cannot be used to choose a model, if one wishes to have a relatively simple model in the end, because it increases artificially with the addition of explanatory variables without significant effect. We will then tend to prefer the adjusted R^2,\bar R^2 = {1-(1-R^{2})\cdot{n-1 \over n-p}} = R^{2}-\underbrace{(1-R^{2})\cdot{p-1 \over n-p}}_{\text{penalty}}where p is the number of parameters of the model. Measuring the quality of fit will penalize overly complex models.

This idea will be found in the Akaike criterion, where AIC=Deviance+2\cdot p or in the Schwarz criterion, BIC=Deviance+log(n)\cdot p. In large dimensions (typically p>\sqrt{n}), we will tend to use a corrected AIC, defined by AIC_c=Deviance+2⋅p⋅n/(n-p-1) .

These criterias are used in so-called “stepwise” methods, introducing the set methods. In the “forward” method, we start by regressing to the constant, then we add one variable at a time, retaining the one that lowers the AIC criterion the most, until adding a variable increases the AIC criterion of the model. In the “backward” method, we start by regressing on all variables, then we remove one variable at a time, removing the one that lowers the AIC criterion the most, until removing a variable increases the AIC criterion from the model.

Another justification for this notion of penalty (we will come back to this idea in machine learning) can be the following. Let us consider an estimator in the class of linear predictors, \mathcal{M}=\big\lbrace m:~m(\mathbf{x})=s_h(\mathbf{x})^T\mathbf{y} \text{ where }S=(s(\mathbf{x}_1),\cdots,s(\mathbf{x}_n))^T\text{ is some smoothing matrix}\big\rbrace and assume that y=m_0 (x)+\varepsilon, with \mathbb{E}[\varepsilon]=0 and Var[\varepsilon]=\sigma^2\mathbb{I}, so that m_0 (x)=\mathbb{E}[Y|X=x] . From a theoretical point of view, the quadratic risk, associated with an estimated model \widehat{m}, \mathbb{E}\big[(Y-\widehat{m}(\mathbf{X}))^2\big], is written\mathcal{R}(\widehat{m})=\underbrace{\mathbb{E}\big[(Y-m_0(\mathbf{X}))^2\big]}_{\text{error}}+\underbrace{\mathbb{E}\big[(m_0(\mathbf {X})-\mathbb{E}[\widehat{m}(\mathbf{X})])^2\big]}_{\text{bias}^2}+\underbrace{\mathbb{E}\big[(\mathbb{E}[\widehat{m}(\mathbf{X})]-\widehat{m}(\mathbf{X}))^2\big]}_{\text{variance}} if m_0 is the true model. The first term is sometimes called “Bayes error”, and does not depend on the estimator selected, \widehat{m}.

The empirical quadratic risk, associated with a model m, is here: \widehat{\mathcal{R}}_n(m)=\frac{1}{n}\sum_{i=1}^n (y_i-m(\mathbf{x}_i))^2 (by convention). We recognize here the mean square error, “mse”, which will more generally give the “risk” of the model m when using another loss function (as we will discuss later on). It should be noted that:\displaystyle{\mathbb{E}[\widehat{\mathcal{R}}_n(m)]=\frac{1}{n}\|m_0(\mathbf{x})-m(\mathbf{x})\|^2+\frac{1}{n}\mathbb{E}\big(\|{Y}-m_0(\mathbf{X})\|^2\big)} We can show that:n\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]=\mathbb{E}\big(\|Y-\widehat{m}(\mathbf{x})\|^2\big)=\|(\mathbb{I}-\mathbf{S})m_0\|^2+\sigma^2\|\mathbb{I}-\mathbf{S}\|^2so that the (real) risk of \widehat{m} is: {\mathcal{R}}_n(\widehat{m})=\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]+2\frac{\sigma^2}{n}\text{trace}(\boldsymbol{S})So, if \text{trace}(\boldsymbol{S})\geq0 (which is not a too strong assumption), the empirical risk underestimates the true risk of the estimator. Actually, we recognize here the number of degrees of freedom of the model, the right-hand term corresponding to Mallow’s C_p, introduced in Mallows (1973) using not deviance but R^2.

Statistical Tests

The most traditional test in econometrics is probably the significance test, corresponding to the nullity of a coefficient in a linear regression model. Formally, it is the test of H_0:\beta_k=0 against H_1:\beta_k\neq 0. The so-called Student test, based on the statistics t_k=\widehat{\beta}_k/se_{\widehat{β}_k}, allows to decide between the two alternatives, using the test p-value, defined by \mathbb{P}[|T|>|t_k|] avec T\overset{\mathcal{L}}{\sim} Std_\nu, where \nu is the number of degrees of freedom of the model (\nu=p+1 for the standard linear model). In large dimension, however, this statistic is of very limited interest, given a significant FDR (“False Discovery Ratio”). Classically, with a level of significance \alpha=0.05, 5% of the variables are falsely significant. Suppose that we have p=100 explanatory variables, but that 5 (only) are really significant. We can hope that these 5 variables will pass the Student test, but we can also expect that 5 additional variables (false positive test) will emerge. We will then have 10 variables perceived as significant, while only half are significant, i.e. an FDR ratio of 50%. In order to avoid this recurrent pitfall in multiple tests, it is natural to use the procedure of Benjamini & Hochberg (1995).

From a correlation to some causal effect

Econometric models are used to implement public policy evaluations. It is therefore essential to fully understand the underlying mechanisms in order to know which variables actually make it possible to act on a variable of interest. But then we move on to another important dimension of econometrics. Jerry Neyman was responsible for the first work on the identification of causal mechanisms, and then Rubin (1974) formalized the test, called the “Rubin causal model” in Holland (1986). The first approaches to the notion of causality in econometrics were based on the use of instrumental variables, models with discontinuity of regression, analysis of differences in differences, and natural or unnatural experiments. Causality is usually inferred by comparing the effect of a policy – or more generally of a treatment – with its counterfactual, ideally given by a random control group. The causal effect of the treatment is then defined as \Delta=y_1-y_0, i.e. the difference between what the situation would be with treatment (noted t=1) and without treatment (noted t=0). The concern is that only y=t\cdot y_1+(1-t)\cdot y_0 and t are observed. In other words, the causal effect of variable t  on t  is not observed (since only one of the two potential variables – y_0 or y_1  is observed for each individual), but it is also individual, and therefore a function of x-covariates. Generally, by making assumptions about the distribution of the triplet (Y_0,Y_1,T) , some parameters of the causal effect distribution become identifiable, based on the density of the observable variables (Y,T) . Classically, we will be interested in the moments of this distribution, in particular the average effect of treatment in the population, \mathbb{E}[\Delta] , or even just the average effect of treatment in the case of treatment \mathbb{E}[\Delta|T=1] . If the result (Y_0,Y_1) is independent of the processing access variable T, it can be shown that \mathbb{E}[\Delta]=\mathbb{E}[Y|T=1]- \mathbb{E} [Y|T=0]. But if this independence hypothesis is not verified, there is a selection bias, often associated with \mathbb{E}[Y_0|T=1]- \mathbb{E} [Y_0|T=0]. Rosenbaum & Rubin (1983) propose to use a propensity to be treated score, p(x)=\mathbb{P}[T=1|X=x] , noting that if variable Y_0\ is independent of access to treatment T conditionally to the explanatory variables X, then it is independent of T  conditionally to the score p(X) : it is sufficient to match them using their propensity score. Heckman et al (2003) thus proposes a kernel estimator on the propensity score, which simply provides an estimator of the effect of the treatment, provided that it is treated.

To be continued next time, we’ll introduce “machine learning techniques” (references mentioned above are online here)

Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
Davis=data.frame(Y=Davis$weight * 2.204622,
                 X1=Davis$sex)

On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') 
 Y=Davis$Y

car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)

 solve(t(X)%*%X)
            [,1]       [,2]
[1,] 0.008928571 0.00000000
[2,] 0.000000000 0.01136364

et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

 solve(t(X)%*%X) %*% (t(X)%*%Y)
         [,1]
[1,] 125.4272
[2,] 167.3258

ce qui correspond effectivement à la sortie de R,

 reg=lm(Y~0+X1,data=Davis)
 summary(reg)
 
Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1F  125.427      1.960   64.00   <2e-16 ***
X1M  167.326      2.211   75.68   <2e-16 ***

Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes

x=Y[X[,1]==1]
y=Y[X[,2]==1]
nx=length(x)
ny=length(y)

On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M   est le poids moyen des hommes. Et en effet

 mean(y)
[1] 167.3258

C’est finalement très naturel, ou intuitif.

On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici

 sqrt(var(y)/ny)
[1] 2.794391
 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)
[1] 2.794391

car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2

 s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2))
 sqrt(s2/ny)
[1] 2.210863

et en effet, on retombe sur la valeur donnée dans le tableau de régresion

 sqrt(s2/nx)
[1] 1.959721

(pareil pour l’autre coefficient).

On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).

 X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y)
          [,1]
[1,] 125.42724
[2,]  41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)
[1] 41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.427      1.960   64.00   <2e-16 ***
X1M           41.899      2.954   14.18   <2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -6.4475, df = 286, p-value = 4.826e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.62603 -16.30035
sample estimates:
mean of x mean of y 
 143.8626  167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Call:
lm(formula = Y ~ X1, data = Davis)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-57.458 -22.184  -5.512  17.809 118.912 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 143.4382     2.7820   51.56   <2e-16 ***
X1M           0.9645     4.1940    0.23    0.818    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 29.44 on 198 degrees of freedom
Multiple R-squared:  0.000267,	Adjusted R-squared:  -0.004782 
F-statistic: 0.05289 on 1 and 198 DF,  p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y
 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -0.22998, df = 198, p-value = 0.8183
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.235209  7.306165
sample estimates:
mean of x mean of y 
 143.4382  144.4027

et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.

Dominance stochastique, ordre 1 ou ordre 2 ?

Mercredi dernier, on commençait le cours avec une présentation rapide du test de Wilcoxon-Mann-Whitney (parfois appelé test de la somme des rangs, de Wilcoxon). Ce test a été proposé par Frank Wilcoxon en 1945, et (presque en même temps) par Henry Mann et Donald Whitney (en 1947), dans un article intitulé On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other. Le titre est clair. Ce test (non paramétrique) a vocation à tester l’hypothèse selon laquelle la distribution dans deux échantillons est la même. Ou, pour reprendre les notations du derniers billets, la distribution est la même dans deux groupes.

Autrement dit, on va chercher à tester, ici, https://latex.codecogs.com/gif.latex?H_0:F=G, au sens où https://latex.codecogs.com/gif.latex?H_0:\forall%20x,F(x)=G(x). Dans le test de Kolmogorov Smirnov, l’hypothèse alternative est https://latex.codecogs.com/gif.latex?H_1:\exists%20x,F(x)\neq%20G(x).

Ici, comme le notent d’ailleur Mann et Whitney, “the test is shown to be consistent with respect to the class of alternatives https://latex.codecogs.com/gif.latex?F(x)%3EG(x) for every https://latex.codecogs.com/gif.latex?x“. Cette propriété est un peu particulière, et correspond à la notion de dominance stochastique à l’ordre 1. https://latex.codecogs.com/gif.latex?F\preceq_1%20Gcorrespond à https://latex.codecogs.com/gif.latex?F(x)\geq%20G(x) https://latex.codecogs.com/gif.latex?\forall%20x, ou encore, pour tout fonction https://latex.codecogs.com/gif.latex?u(\cdot) croissante, si https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)] avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G.

Mais regardons un peu le test. On considère deux populations https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_{n_x}\} et https://latex.codecogs.com/gif.latex?\boldsymbol{y}=\{y_1,\cdots,y_{n_y}\}. L’idée est de réordonner les https://latex.codecogs.com/gif.latex?(n_x%20+%20n_y) éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{z}=\boldsymbol{x}\cup\boldsymbol{y} par ordre croissant. On note https://latex.codecogs.com/gif.latex?\boldsymbol{r} le vecteur des rangs associés. Soit https://latex.codecogs.com/gif.latex?S_x la somme des https://latex.codecogs.com/gif.latex?n_x (premiers) rangs des éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{x}, i.e. https://latex.codecogs.com/gif.latex?S_x=r_1+\cdots+r_{n_x}.

Comme l’ont montré Wilcoxon, Mann et Whitney, sous https://latex.codecogs.com/gif.latex?H_0, https://latex.codecogs.com/gif.latex?S_x suit une “loi connue”, où classiquement, “connue” veut dire tabulée pour de petits échantillons, et pour des échantillons plus grands.

Histoire d’illustrer, considerons deux échantillons, tirés suivant deux lois https://latex.codecogs.com/gif.latex?F et https://latex.codecogs.com/gif.latex?G avec https://latex.codecogs.com/gif.latex?F\preceq_1%20G, par exemple deux lois normales de moyennes différentes (et de même variance)

> u=seq(-3,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,.5,1),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 1. Les densité sont ici

>plot(u,dnorm(u,0,1),col="red",type="l")
>lines(u,dnorm(u,.5,1),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,.5,1))
+  P[i,1:2]=c(sum(rank(x)[1:100]),sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Sous https://latex.codecogs.com/gif.latex?H_0, la somme aurait eu la distribution avec les pointillés. On en est manifestement assez loin, autrement dit, on va rejeter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G.

Intuitivement, c’est normal, car effectivement, les rangs dans un cas sont – en moyenne – plus faibles que pour l’autre,

Au lieu de simuler des lois dominées, à l’ordre 1, on peut tenter des lois dominées à l’ordre 2. Rappelons que https://latex.codecogs.com/gif.latex?F\preceq_2%20G signifie

https://latex.codecogs.com/gif.latex?\int_{-\infty}^x%20F(t)dt%20\geq\int_{-\infty}^x%20G(t)dt,\forall%20x

ou encorehttps://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]

avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est unne fonction croissante et concave.

En gestion des risques (ou sur les inégalités), on va supposer en plus que https://latex.codecogs.com/gif.latex?\mathbb{E}[X]=\mathbb{E}[Y], ce qui correspond à une notion de “mean preserving spread” (on peut relire le document sur la notion de
single crossing, par Alain Chateauneuf). Une caractérisation est

https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est concave. Une caractérisation alternative est https://latex.codecogs.com/gif.latex?Y=X+\varepsilon avec https://latex.codecogs.com/gif.latex?\mathbb{E}[\varepsilon\vert%20X]=0 (on pourra relire Rothschild et Stiglitz, Increasing Risk: A Definition)

Dans le cas des lois normales, on a cette notion de ‘mean preserving spread’ lorsque les moyennes sont identiques, mais pas les variances. On a alors plus de dispersion,

> u=seq(-4,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,0,1.5),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 2, avec du single crossing. Les densités sont ici

> plot(u,dnorm(u,0,1),col="red",type="l")
> lines(u,dnorm(u,0,1.5),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,0,1.5))
+  P[i,1:2]=c(sum(rank(x)[1:100]), sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Autrement dit, ici les sommes des rangs sont identiques, avec les deux échantillons,  et donc on aura tendance à accepter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G. L’idée naturelle est que avec les dans le premiers cas, effectivement, on avait plus de petits rangs avec le petit échantillon. Mais plus ici: le fait d’avoir plus de dispersion fait qu’on a en même temps plus de petits, et plus de grands. Et en moyenne, ça ne change rien

Interpréter un test de Wald

Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l’interprétation d’un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données

> p=.5

En l’occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face)

> n=20
> set.seed(1)
> echantillon=sample(c("pile","face"),
+    size=n,replace=TRUE,
+    prob=c(p,1-p))
> echantillon
[1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile"

Classiquement, on va recoder en https://latex.codecogs.com/gif.latex?\{0,1\}

> X=(echantillon=="pile")*1
> X
[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1

La vraisemblance est ici donnée par

> L=function(p) prod(dbinom(X,1,p))
> logL=function(p) sum(log(dbinom(X,1,p)))
> xp=seq(0,1,by=.01)
> yL=Vectorize(L)(xp)
> plot(xp,yL,type="l")

mais on va plutôt travailler sur la log-vraisemblance

> ylogL=Vectorize(logL)(xp)
> plot(xp,ylogL,type="l")

Continue reading Interpréter un test de Wald

Niveau et puissance de test (paramétrique)

Lors du dernier cours de statistique, nous avons vu comme se construisent des tests, et comment on compare leur puissance. Considérons pour cela des observations \{x_1,\cdots,x_n\} tirées suivant une loi uniforme \mathcal{U}([0,\theta]). Et on se pose la question de savoir si l’hypothèse H_0:\theta\leq\theta_0 est rejetée, ou acceptée, avec \theta_0=1 par exemple. L’hypothèse alternative sera ici le complémentaire, autrement dit H_1:\theta>\theta_0.

Comme on l’a vu, la construction d’un test se fait en utilisant une statistique de test (ou plus simplement un estimateur du paramètre inconnu), avec \psi(\mathbf{x})=\mathbf{1}(\widehat{\theta}(\mathbf{x}))\leq c. On prend la décision D_1 (consistant à rejeter H_0) si \psi(\mathbf{x})=1, et la décision D_0 (consistant à ‘accepter’ H_0) si \psi(\mathbf{x})=0.

Le risque d’erreur de première espèce est basé sur \mathbb{E}[\psi(\mathbf{X})]. Si l’hypothèse H_0 est une hypothèse simple H_0:\theta=\theta_0, le risque d’erreur de première espèce est https://latex.codecogs.com/gif.latex?\alpha=\mathbb{E}[\psi(\mathbf{X})\vert%20X_1,\cdots,X_n\sim%20F_{\theta_0}], soit \alpha=\mathbb{P}[\widehat{\theta}(\mathbf{X})>c \vert X_1,\cdots,X_n\sim F_{\theta_0}] pour un seuil c. On va déterminer le seuil c de manière à ce que le risque de première espèce soit 5%.

Pour cela, il nous fait un estimateur \widehat{\theta} dont on connaîtrait la loi. Commençons par l’estimateur du maximum de vraisemblance. Rappelons que l’estimateur du maximum de vraisemblance est \widehat{\theta}=\max\{x_i\}. Nous avons vu – dans la section sur les statistiques d’ordre – que \mathbb{P}[\max\{X_i\}]=\left(\frac{x}{\theta}\right)^n

pour x\leq\theta. Aussi, pour un seuil c, \mathbb{P}[\widehat{\theta}>c]=1-\left(\frac{c}{\theta}\right)^n

pour c\leq\theta. Le risque de première espèce \alpha=\underset{\theta\leq\theta_0}{\max}\left\lbrace\mathbb{P}\big[\widehat{\theta}\geq c\vert X\sim F_{\theta}\big]\right\rbrace

soit \alpha=\underset{\theta\leq\theta_0}{\max}\left\lbrace1-\left(\frac{c}{\theta}\right)^n\right\rbrace

qui atteindra son maximum en \theta_0. Aussi \alpha=1-\left(\frac{c}{\theta_0}\right)^n

soit c=\theta_0\cdot[1-\alpha]^{1/n}. Les régions d’acceptation et de rejet sont les suivantes,

Continue reading Niveau et puissance de test (paramétrique)

Statistical Tests: Asymptotic, Exact, ou based on Simulations?

This morning, in our mathematical statistics course, we’ve been discussing the ‘proportion test‘, i.e. given a sample of Bernoulli trials, with , we want to test

against 

A natural test (which can be related to the maximum likelihood ratio test) is  based on the statistic

The test function is here

To get the bounds of the acceptance region, we need the distribution of , under . Consider here a numerical application

n=20
p=.5
set.seed(1)
echantillon=sample(0:1,size=n,
            prob=c(1-p,p),
            replace=TRUE)
  • the asymptotic distribution

The first (and standard idea) is to use the central limit theorem, since

So, under ,

Then  while . The acceptance region is then between the two red lines, below,

T=sqrt(n)*(mean(echantillon)-.5)/
  sqrt(mean(echantillon)*
  (1-mean(echantillon)))
u=seq(-3,3,by=.01)
v=dnorm(u)
plot(u,v,type="l",lwd=2)
abline(v=qnorm(.025),col="red")
abline(v=qnorm(.975),col="red")
abline(v=T,col="blue")

  • the exact distribution

Here we use the fact that

Using transformation of the ‘density’, we can (at least numerically) compute the (exact) distribution of

 

u=seq(-3,3,by=.01)
v=sqrt(.5*(1-.5))*n*dbinom(round(
  (sqrt(.5*(1-.5))*u/sqrt(n)+.5)*n),
  size=n,prob=.5)/sqrt(n)

Here I used a round value, it guess it would be better with a floor function, but here the graph looks symmetric (which is something I like)

abline(v=sqrt(n)*(qbinom(.025,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
abline(v=sqrt(n)*(qbinom(.975,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
lines(u,v,type="s")

  • distribution based on Monte Carlo simulations

Probably more interesting, here we do not use the fact that we might know the distribution of the mean. We just generate random samples, under , and then compute ,

T=rep(NA,1000)
for(i in 1:1000){
x=sample(0:1,size=n,
         prob=c(1-.5,.5),
         replace=TRUE)
m=mean(x)
T[i]=(m-.5)/sqrt(m*(1-m))*sqrt(n)}
lines(density(T),lwd=2)
abline(v=quantile(T,.025),col="red")
abline(v=quantile(T,.975),col="red")

Where does that 2 come from in the likelihood ratio test?

This afternoon, in class, we’ve seen Wald test, the likelihood-ratio test, and finally the score test. All of them rely on the same idea

and then, use that if   with , we can write

Or – slightly more interesting – if  , then

Then one can get that

Based on that property, we can derive Wald statistics,

that can be visualized below

The score test is a test on the square of the slope

The idea for the likelihood ratio test is to consider

Observe that  can be written, using Taylor’s expansion

for some . The first term is null, since the maximum likelihood estimator is precisely at the maximum of the (log) likelihood. So

That’s more or less where the 2 comes from. Then observe that

and therefore

This test will be discussed further next week (since it is related to Neyman-Pearson’s theorem), but also, that result can be used to derive confidence intervals. With a log-likelihood as follows

it is possible to get a confidence interval for the parameter by looking for‘s such that

We will discuss that idea later on, in the context of profile likelihood.

Tests, Power and Significance

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…).

Dear reader, who are you?

While working on an updated (not to say revised) version of Blogging in Academia, a Personal Experience, I noticed that I now have a better understanding of me as a blogger (this academic alter-ego, to use Vanessa Paz Dennen’s words), the role of blogs within academia, and the community of academic-bloggers. But tonight, in the train, I was asking myself about the readers. I mean You. I know from the stats of my blog that a few thousand people will read this post, that you are 32% American, 21% French, 6% Brazilian and 3% Swedish. That it is quite unlikely that you use a Linux computer. But I don’t know much more…

Since I have always claimed that my blog was a place to experiment, I wanted to try one! Dear reader, if you have 10 seconds, can you please post a comment, with your age, your gender, and your background (I’d like to know if you’re a student, a colleague, an economist, a mathematician, a biologist, a journalist, an astronaut, a lion tamer, an actuary, a modeler, a model, a robot… ). If you want to spend more time explaining how you ended up on my blog, I’d be glad to know… But at least, your age, your gender, and your ‘academic background’ (since I have some sort of academic blog), that would be great! If you come frequently on my blog, you should start to know me, and I thought it would be great if I could know more about you…

Puissance et tests statistiques (simples)

Pour comprendre la notion de puissance dans un test statistique, reprenons Statistical power analysis for the behavioral sciences de Jacob Cohen. Comme il le rappelle Fisher (dans The Design of Experiments) ne souhaite pas “prouver” qu’une hypothèse (appelée H_0) est valide, mais il espère, en réalité, rejeter l’hypothèse en question (c’est ce que rappellent d’ailleurs Hubbard et Bayarri en 2003). Autrement dit, supposons que l’on cherche à faire un test médical, et à détecter si une personne est malade (ou pas). Si une personne est généralement malade dès lors qu’un taux de quelque chose est faible (noté https://latex.codecogs.com/gif.latex?\theta), on va chercher à tester https://latex.codecogs.com/gif.latex?H_0:\theta_\star%3E\theta_0. Dans un test statistique, on espère rejeter l’hypothèse nulle https://latex.codecogs.com/gif.latex?H_0.

The power of a statistical test of a null hypothesis is the probability that it will lead to the rejection of the null hypothesis.

Continue reading Puissance et tests statistiques (simples)

Tests basés sur la vraisemblance – score

Une autre grandeur intéressant est le score, qui est la dérivée de la vraisemblance. Intuitivement (c’est l’idée de la condition du premier ordre),   et  seront proches si les dérivées en ces points sont proches. En  la dérivée est nulle, donc on va se demander ici, tout simplement, si la dérivée en  est proche de 0. Ou pas.

Continue reading Tests basés sur la vraisemblance – score

Tests basés sur la vraisemblance – Rapport de Vraisemblance

Chose promise, chose due. J’avais dit qu’on parlerait du test de rapport de vraisemblance. L’idée – visuelle – est d’avoir une lecture dans l’autre sens : au lieu de se demander si  et  sont proches, on va se demander si   et  sont proches.

Si la fonction de vraisemblance est suffisamment régulière, on se pose la même question.

Continue reading Tests basés sur la vraisemblance – Rapport de Vraisemblance

Tests basés sur la vraisemblance – Wald

Petit rappel préliminaire. Si on dispose d’un échantillon , i.id de loi , où le paramètre  est inconnu, alors on peut calculer le maximum de la vraisemblance, ce qui nous donnera un estimateur intéressant (cf. cours de stat). Illustrons avec un jeu de pile ou face. Reprenons l’échantillon de mon précédent billet,

> X=c(0, 0, 1, 1, 0, 1, 1, 1, 1, 0,
      0, 0, 1, 0, 1, 0, 1, 1, 0, 1)

On peut ici tracer la fonction de vraisemblance

ou (un peu plus malin) la fonction de log-vraisemblance

> p=seq(0,1,by=.01)
> logL=function(p)
+ {sum(log(dbinom(X,size=1,prob=p)))}
> plot(p,Vectorize(logL)(p),
+ type="l",col="red",lwd=2)
> p0=.5
> abline(v=p0,col="blue")

La valeur correspondant au trait bleu correspondant au cas que l’on va chercher à tester, de pièce non pipée, soit  (de manière très générale)

Continue reading Tests basés sur la vraisemblance – Wald

Les tests et la logique, modus tollens

Quand on apprend la logique, on apprend la notion de modus tollens, correspondant à une notion de contraposition. Si on a une proposition du genre , alors la proposition contraposée est .  Et on apprend que les deux propositions sont équivalentes (je fais de la logique classique). Par exemple si  correspond à “feu” et  à “fumée“,  signifie que tout feu fait de la fumée. Si cette affirmation est vrai, alors il n’y a pas de fumée sans feu, et donc , autrement dit, s’il n’y a pas de fumée, c’est qu’il n’y a pas de feu.

Continue reading Les tests et la logique, modus tollens