# 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}\|^2$$so 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_i$$On 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(&gt;|t|) X1F 125.427 1.960 64.00 &lt;2e-16 *** X1M 167.326 2.211 75.68 &lt;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_i$$autrement dit $\widehat{{\beta}}_M$ est le poids moyen des hommes. Et en effet  mean(y)  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)  2.794391 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)  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)  2.210863 et en effet, on retombe sur la valeur donnée dans le tableau de régresion  sqrt(s2/nx)  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)  41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis) summary(reg2)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 125.427 1.960 64.00 &lt;2e-16 *** X1M 41.899 2.954 14.18 &lt;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(&gt;|t|) (Intercept) 143.4382 2.7820 51.56 &lt;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, , au sens où . Dans le test de Kolmogorov Smirnov, l’hypothèse alternative est .

Ici, comme le notent d’ailleur Mann et Whitney, “the test is shown to be consistent with respect to the class of alternatives for every “. Cette propriété est un peu particulière, et correspond à la notion de dominance stochastique à l’ordre 1. correspond à , ou encore, pour tout fonction  croissante, si avec  et .

Mais regardons un peu le test. On considère deux populations  et . L’idée est de réordonner les  éléments de par ordre croissant. On note  le vecteur des rangs associés. Soit  la somme des  (premiers) rangs des éléments de , i.e. .

Comme l’ont montré Wilcoxon, Mann et Whitney, sous , 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 et avec , 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 (et celle de ).

> 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 , 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, . 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 signifie

ou encore

avec et . où est unne fonction croissante et concave.

En gestion des risques (ou sur les inégalités), on va supposer en plus que , 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

et . où est concave. Une caractérisation alternative est avec (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 (et celle de ).

> 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, . 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  "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
 > X=(echantillon=="pile")*1 > X  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") # 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 , 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, # 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…

# 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é ), on va chercher à tester . Dans un test statistique, on espère rejeter l’hypothèse nulle .

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

# 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. # 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.

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

# Introduction aux tests nonparamétriques

Aujourd’hui, on va commencer la partie du cours sur les tests, avec une introduction par les tests non-paramétriques. On parlera (dans la continuité de ce qu’on a fait jusqu’à présent) de tests d’ajustement de lois, mais aussi de tests plus divers, comme des tests d’indépendance (entre observations successives).

Le test le plus classique est le test des signes. L’idée est de regarder puis de compter le nombre de fois où (différence positive) et le nombre de fois où  (différence négative).

> (X=rnorm(20))
  0.81352495  0.78076709 -0.58654359  0.52805014 -0.31101282 -0.65022034 -0.93619729 -0.16483729 -1.97396753 -0.22506985  1.14959851 -0.38895448 -0.09298421
> (S=c("N","P")[1+1*(diff(X)>0)])
 "N" "N" "P" "N" "N" "N" "P" "N" "P" "P" "N" "P" "P" "N" "P" "N" "P" "P" "N"

On verra ce test en cours. On peut par exemple compter le nombre de fois où la différence était positive,

> sim=function(a=0){
+ n=100
+ epsilon=X=rnorm(n)
+ for(t in 2:n) X[t]=a*X[t-1]+epsilon[t]
+ S=c("N","P")[1+1*(diff(X)>0)]
+ return(sum(S=="P"))}

(je fais ici une fonction un peu plus compliqué pour voir ce qui se passe si les données ne sont pas indépendantes, justement). Si on génère plein d’échantillons, on peut voir que le nombre de variations positives suit une loi normale, si les observations sont effectivement indépendantes,

> library(ggplot2)
> set.seed(1)
> vectS=unlist(lapply(rep(0,1000),sim))
> ggplot(data=data.frame(vectS = vectS),
+ aes(x = vectS)) +
+ geom_histogram(aes(y = ..density..),
+ binwidth = 1, fill = "light blue",
+ col = "black") +
+ stat_function(fun=function(x) dnorm(x,
+ (n-1)/2+.5,sqrt((n+1)/12)),
+ col = "red", size = 1) + ylab("Densite") +
+ ggtitle("Histogramme") 