# Exotic link functions for GLMs

In my previous post on GLMs, I discussed power link functions. But there are much more links that can be used :

• The square root link (for the Poisson model)

Consider some random variable $Y$ with mean $\mu$ and variance $\sigma^2$. Using Taylor’s expansion,$$g(Y)\sim g(\mu)+(Y-\mu)g'(\mu)+\frac{1}{2}(Y-\mu)^2g''(\mu)$$we can write$$\mathbb{E}[g(Y)]\sim g(\mu)+\frac{\sigma^2}{2}g''(\mu)$$ $$\text{Var}[g(Y)]\sim [g'(\mu)]^2\sigma^2$$

Assume that $Y\sim\mathcal{P}(\lambda)$, a consider a square root transformation, $g(y)=\sqrt{y}$, then the second equality becomes $$\text{Var}[\sqrt{Y}]\sim \left[\frac{1}{2\sqrt{\mathbb{E}[Y]}}\right]^2\text{Var}[Y]=\frac{1}{4}$$

So, somehow, with a square-root transformation, we have variance stability, which might be interpreted as some homoscedasticity.

• The complementary log-log function for the Bernoulli model

Assume that the true variable of interest is a Poisson one, $N|\mathbf{X}=\mathbf{x}\sim\mathcal{P}(\lambda_{\mathbf{x}})$ where $\lambda_{\mathbf{x}}=\exp[\mathbf{x}^T\mathbf{\beta}]$Thus,$$\mathbb{P}[N=0|\mathbf{X}=\mathbf{x}]=\exp[-\lambda_{\mathbf{x}}]=\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]$$while$$\mathbb{P}[N>0|\mathbf{X}=\mathbf{x}]=1-\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]=H(\mathbf{x}^T\mathbf{\beta})$$where $H(s)=1-\exp[-\exp(s)]$. Let $Y=\mathbf{1}(N>0)$. The previous model seems like a Bernoulli regression with $H$ as link function,$$\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=H(\mathbf{x}^T\mathbf{\beta})$$

So, assume now that instead of observing $N$ we observe $Y=\boldsymbol{1}(N>0)$. In that case, running a Bernoulli regression with a complementary log-log link function would be the same (?) as running first a Poisson regression on the original data, and then use it on our binary variable, zero vs. non-zero. Let us generate some data, and see what’s going on. Let us compare $e^{\lambda_{\mathbf{x}}}$ and $p_{\mathbf{x}}$ obtained from a standard logistic regression

n=563 set.seed(1) base=data.frame(X1=rnorm(n),X2=rnorm(n)) lambda=base$X1+base$X2 base$Y=rpois(n,exp(lambda)) regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") What if $p_{\mathbf{x}}$ was obtained from a Bernoulli regression, with a cloglog link function ? regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") It looks like the fit is very good here ! Now, what if we have real data, like the dataset from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy (with 563 observations, and nine variables) base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) x=base$SEX base$SEX="M" base$SEX[x=="0"]="F" x=base$CHILDREN base$CHILDREN="YES" base$CHILDREN[x==0]="NO" regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") In that case the two models are very different. But actually, so is the second one regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") How can we interpret that ? Could it be because the Poisson model is not good ? Actually, if we run a zero-inflated model here, library(pscl) regZIP = zeroinfl(Y ~ . | ., data = base) summary(regZIP) Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.002274 0.048413 -0.047 0.963 X1 1.019814 0.026186 38.945 &lt;2e-16 *** X2 1.004814 0.024172 41.570 &lt;2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -4.90190 2.07846 -2.358 0.0184 * X1 -2.00227 0.86897 -2.304 0.0212 * X2 -0.01545 0.96121 -0.016 0.9872 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Hence, we reject here the Poisson distribution assumption, because of the inflation of zeros… It looks like the cloglog link can be used to check if the Poisson distribution is a good model, or not… # GLMs: link vs. distribution Usually, when I give a course on GLMs, I try to insist on the fact that the link function is probably more important than the distribution. In order to illustrate, consider the following dataset, with 5 observations x = c(1,2,3,4,5) y = c(1,2,4,2,6) base = data.frame(x,y) Then consider several model, with various distributions, and either an identity link (and in that case $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbf{x}^T\mathbf{\beta}$) or a log link function (so that $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=e^{\mathbf{x}^T\mathbf{\beta}}$) regNId = glm(y~x,family=gaussian(link="identity"),data=base) regNlog = glm(y~x,family=gaussian(link="log"),data=base) regPId = glm(y~x,family=poisson(link="identity"),data=base) regPlog = glm(y~x,family=poisson(link="log"),data=base) regGId = glm(y~x,family=Gamma(link="identity"),data=base) regGlog = glm(y~x,family=Gamma(link="log"),data=base) regIGId = glm(y~x,family=inverse.gaussian(link="identity"),data=base) regIGlog = glm(y~x,family=inverse.gaussian(link="log"),data=base One can also consider some Tweedie distribution, to be even more general library(statmod) regTwId = glm(y~x,family=tweedie(var.power=1.5,link.power=1),data=base) regTwlog = glm(y~x,family=tweedie(var.power=1.5,link.power=0),data=base) Consider the prediction obtained in the first case, with the linear link function library(RColorBrewer) darkcols = brewer.pal(8, "Dark2") plot(x,y,pch=19) abline(regNId,col=darkcols[1]) abline(regPId,col=darkcols[2]) abline(regGId,col=darkcols[3]) abline(regIGId,col=darkcols[4]) abline(regTwId,lty=2) The predictions are very very close, aren’t they ? In the case of the exponential prediction, we obtain plot(x,y,pch=19) u=seq(.8,5.2,by=.01) lines(u,predict(regNlog,newdata=data.frame(x=u),type="response"),col=darkcols[1]) lines(u,predict(regPlog,newdata=data.frame(x=u),type="response"),col=darkcols[2]) lines(u,predict(regGlog,newdata=data.frame(x=u),type="response"),col=darkcols[3]) lines(u,predict(regIGlog,newdata=data.frame(x=u),type="response"),col=darkcols[4]) lines(u,predict(regTwlog,newdata=data.frame(x=u),type="response"),lty=2) We can actually look closer. For instance, in the linear case, consider the slope obtained with a Tweedie model (that will include all the parametric familes mentioned here, actually) pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base))$coefficients[2,1:2] Vgamma = seq(-.5,3.5,by=.05) Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3,ylim=c(.965,1.03),xlab="power",ylab="slope")

The slope here is always very very close to one ! Even more if we add a confidence interval

plot(Vgamma,Vpente[1,]) lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2)

Heuristically, for the Gamma regression, or the Inverse Gaussian one, because the variance is a power of the prediction, if the prediction is small (here on the left), the variance should be small. So, on the left of the graph, the error should be small with a higher power for the variance function. And that’s indeed what we observe here

erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(-.1,.04),xlab="power",ylab="error") abline(h=0,lty=2)

Of course, we can do the same with the exponential models

pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base))$coefficients[2,1:2] Vpente = Vectorize(pente)(Vgamma) plot(Vgamma,Vpente[1,],type="l",lwd=3) or, if we add the confidence bands, we obtain plot(Vgamma,Vpente[1,],ylim=c(0,.8),type="l",lwd=3,xlab="power",ylab="slope") lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2) lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2) So here also, the “slope” is rather similar… And if we look at the error we make on the left part of the graph, we obtain erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base),newdata=data.frame(x=1),type="response")-y[x==1] Verreur = Vectorize(erreur)(Vgamma) plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(.001,.32),xlab="power",ylab="error") So my point is that the distribution is usually not the most important point on GLMs, even if chapters of books on GLMs are distribution based… But as mentioned in an another post, if you consider a nonlinear transformation, like we have with GAMs, the story is more complicated… # Bailey (1963) and Poisson regression on two factors Consider the following dataset, from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy, with 563 observations, and nine variables : eight covariates, and the variable of interest, the number of extramarital affairs, over a year, base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) 'data.frame': 563 obs. of 9 variables:$ SEX : int 1 0 0 1 1 0 0 1 0 1 ... $AGE : num 37 27 32 57 22 32 22 57 32 22 ...$ YEARMARRIAGE: num 10 4 15 15 0.75 1.5 0.75 15 15 1.5 ... $CHILDREN : int 0 0 1 1 0 0 0 1 1 0 ...$ RELIGIOUS : int 3 4 1 5 2 2 2 2 4 4 ... $EDUCATION : int 18 14 12 18 17 17 12 14 16 14 ...$ OCCUPATION : int 7 6 1 6 6 5 1 4 1 4 ... $SATISFACTION: int 4 4 4 5 3 5 3 4 2 5 ...$ Y : int 0 0 0 0 0 0 0 0 0 0 ...

Let us focus on two categorical covariates, related to the importance of religion, and the occupation

df=data.frame(y=base$Y, religion=as.factor(base$RELIGIOUS), occupation=as.factor(base$OCCUPATION), expo = 1) (E=xtabs(expo~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 8 4 16 9 0 2 23 3 11 17 56 36 6 3 29 1 10 12 39 25 2 4 38 7 12 21 59 44 2 5 13 1 3 10 19 19 3 (N=xtabs(y~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 13 3 13 7 0 2 1 1 13 10 25 43 10 3 15 0 12 11 34 35 1 4 24 1 3 15 11 9 10 5 6 0 0 6 11 7 0 The two tables above are the exposure (number of observations) and the number of extramarital affairs, here as contingency tables. Without any covariate, one can assume that $N\sim\mathcal{P}(\lambda\cdot E)$, where $\lambda$ would be sum(N)/sum(E) [1] 0.6305506 The idea with the margin method is to assume that $N_{i,j}=E_{i,j}\cdot\lambda_{i,j}$ where $\lambda_{i,j}=A_i\cdot B_j$. Bailey (1963) added two series of constraints : per row, $$\sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j$$ for any $i$ and similarly, for any $j$ $$\sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_j$$From the first series of constraints, write $$A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j}$$ and use the second series to write $$B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}$$Because we need $A_i$‘s to compute $B_j$‘s, and conversely, it is natural to consider some iterative procedure to solve it. Observe that we do not have unicity… Consider here some starting values for $A_i$‘s and $B_j$‘s A=rep(1,length(levels(df$religion))) B=rep(1,length(levels(df$occupation)))*sum(N)/sum(E) A [1] 1 1 1 1 1 B [1] 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 The predicted number of extramarital affairs would be $\hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j$ E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.5222025 0.6305506 5.0444050 2.5222025 10.0888099 5.6749556 0.0000000 2 14.5026643 1.8916519 6.9360568 10.7193606 35.3108348 22.6998224 3.7833037 3 18.2859680 0.6305506 6.3055062 7.5666075 24.5914742 15.7637655 1.2611012 4 23.9609236 4.4138544 7.5666075 13.2415631 37.2024867 27.7442274 1.2611012 5 8.1971581 0.6305506 1.8916519 6.3055062 11.9804618 11.9804618 1.8916519 sum(B*E[1,]) [1] 26.48313 sum(B*E[2,]) [1] 95.84369 apply(t(B*t(E)),1,sum) 1 2 3 4 5 26.48313 95.84369 74.40497 115.39076 42.87744 sum(A*E[,1]) [1] 107 sum(A*E[,2]) [1] 13 apply(A*E,2,sum) 1 2 3 4 5 6 7 107 13 44 64 189 133 13 From expressions above, observe that one can very easily write expressions of $A_i$‘s and $B_j$‘s as functions of $B_j$‘s and $A_i$‘s respectively A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) Let it iterate one thousand times for(i in 1:1000){ A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) } We obtain here A 1 2 3 4 5 1.5404346 1.0447195 1.4825650 0.6553159 0.6634763 B 1 2 3 4 5 6 7 0.4685515 0.2629769 0.8454435 0.7245310 0.4889697 0.7770553 1.6753750 E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044 4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997 5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148 That is our prediction, per category, of the number of affairs. Observe that here, sums per row are equal to observed numbers, apply(N,1,sum) 1 2 3 4 5 41 103 108 73 30 apply(E * A%*%t(B),1,sum) 1 2 3 4 5 41 103 108 73 30 as well as sums per colums apply(N,2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 apply(E * A%*%t(B),2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 Now, why should I mention that here, in the section on the Poisson regression in our course ? Because actually, this is exactly what we get if we run a Poisson regression on those two covariates reg=glm(y~religion+occupation,data=df,family=poisson) summary(reg) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.32604 0.21325 -1.529 0.126285 religion2 -0.38832 0.18791 -2.066 0.038783 * religion3 -0.03829 0.18585 -0.206 0.836771 religion4 -0.85470 0.19757 -4.326 1.52e-05 *** religion5 -0.84233 0.24416 -3.450 0.000561 *** occupation2 -0.57758 0.59549 -0.970 0.332083 occupation3 0.59022 0.21349 2.765 0.005699 ** occupation4 0.43588 0.20603 2.116 0.034381 * occupation5 0.04265 0.17590 0.242 0.808399 occupation6 0.50587 0.17360 2.914 0.003569 ** occupation7 1.27415 0.26298 4.845 1.27e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 First of all, observe that the total sum of predictions equals the total sum of observations yp = predict(reg,type="response") sum(yp) [1] 355 sum(df$y) [1] 355

But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey’s technique

xtabs(yp~df$religion+df$occupation) df$occupation df$religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586112 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450813 0.3898804 12.5342484 12.8899708 28.2722424 28.8008726 4.9677044 4 11.6678703 1.2063307 6.6483904 9.9707300 18.9053460 22.4055332 2.1957997 5 4.0413464 0.1744790 1.6827951 4.8070914 6.1639761 9.7955975 3.3347148 E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044 4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997 5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148

To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for $A_i$‘s

a=exp(coefficients(reg)[1]+c(0,coefficients(reg)[2:5])) a/a[1] religion2 religion3 religion4 religion5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 A/A[1] 1 2 3 4 5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072

but also for $B_j$‘s

b=exp(coefficients(reg)[1]+c(0,coefficients(reg)[6:11])) b/b[1] occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 1.0000000 0.5612551 1.8043769 1.5463210 1.0435773 1.6584203 3.5756477 B/B[1] 1 2 3 4 5 6 7 1.0000000 0.5612551 1.8043770 1.5463210 1.0435773 1.6584203 3.5756478

This will have major implications in non-life insurance models (for claims reserving).

# Variable explicative dans un intervalle

Suite a une question posée ce matin en cours, je vais faire un rapide billet pour expliquer comment extraire les valeurs inférieures et supérieures quand on a des intervalles, sous R. Commençons par générer des données,

n=200 set.seed(123) X=rnorm(n) Y=2+X+rnorm(n,sd = .3)

Supposons maintenant que l’on n’observe plus la vraie variable $x$ mais juste une classe (on va créer huit classes, avec chacune un huitième des observations)

Q=quantile(x = X,(0:8)/8) Q[1]=Q[1]-.00001 Xcut=cut(X,breaks = Q) B=data.frame(Y=Y,X=Xcut)

Par exemple, pour la premiere valeur, on a

as.character(Xcut[1]) [1] "(-0.626,-0.348]"

Pour extraire des informations sur ces bornes, on peut utiliser le petit code suivant, qui renvoie la borne inférieure, la borne supérieur, et le milieu de l’intervalle

extraire = function(x){ ax=as.character(x) lower1 = as.numeric( sub("\$$(.+),.*", "\\1", ax) ) lower2 = as.numeric( sub("\$(.+),.*", "\\1", ax) ) upper1 = as.numeric( sub("[^,]*,([^]]*)\$", "\\1", ax) ) upper2 = as.numeric( sub("[^,]*,([^]]*)\$$", "\\1", ax) ) lower = c(lower1,lower2) lower=lower[!is.na(lower)] upper = c(upper1,upper2) upper=upper[!is.na(upper)] mid = (lower+upper)/2 return(c(lower=lower,mid=mid,upper=upper)) }

On peut vérifier sur notre première observation

extraire(Xcut[1]) lower mid upper -0.626 -0.487 -0.348

Juste pour voir, on peut créer trois variables supplémentaires dans notre base (avec ces trois informations)

B2=Vectorize(function(i) extraire(Xcut[i]))(1:n) B$lower=B2[1,] B$mid =B2[2,] B$upper=B2[3,] et on peut comparer 4 régressions (i) on régresse sur nos 8 classes, i.e. nos 8 facteurs (ii) on régresse sur la borne inférieure de l’intervalle, (iii) sur la valeur “moyenne” de l’intervalle (iv) sur la borne supérieure regF=lm(Y~X,data=B) regL=lm(Y~lower,data=B) regM=lm(Y~mid,data=B) regU=lm(Y~upper,data=B) On peut comparer les prévisions avec nos quatre modeles plot(B$Y,predict(regF),ylim=c(0,4)) points(B$Y,predict(regM),col="red") points(B$Y,predict(regU),col="blue") points(B$Y,predict(regL),col="purple") abline(a=0,b=1,lty=2) Pour aller plus loin, on peut aussi comparer les AIC de nos modèles, AIC(regF) [1] 204.5653 AIC(regM) [1] 201.1201 AIC(regL) [1] 266.5246 AIC(regU) [1] 255.0687 Si l’utilisation des bornes inférieures et supérieures n’est pas concluante, ici, on notera qu’utiliser la valeur moyenne de l’intervalle donne des résultats un peu meilleurs que d’utiliser 8 facteurs. # Monte Carlo techniques to create counterfactuals In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample $\{y_1,\cdots,y_n\}$, and more generally, in econometrics $\{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}$. But let’s get back to statistics (without covariates) to illustrate. We assume that observations $y_i$ are realizations of an underlying random variable $Y_i$. We assume that $Y_i$ are i.id. random variables, with (unkown) distribution $F_{\theta}$. Consider here some estimator $\widehat{\theta}$ – which is just a function of our sample $\widehat{\theta}=h(y_1,\cdots,y_n)$. So $\widehat{\theta}$ is a real-valued number like . Then, in mathematical statistics, in order to derive properties of the estimator $\widehat{\theta}$, like a confidence interval, we must define $\widehat{\theta}=h(Y_1,\cdots,Y_n)$, so that now, $\widehat{\theta}$ is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, $\widehat{\theta}$ is . There are two strategies here. In classical statistics, we use probability theorem, to derive properties of $\widehat{\theta}$ (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, $\{y_1,\cdots,y_n\}$, and that’s a pity. But maybe we can create another one $\{y_1^{(1)},\cdots,y_n^{(1)}\}$, as realizations of $F_{\theta}$, and another one $\{y_1^{(2)},\cdots,y_n^{(2)}\}$, anoter one $\{y_1^{(3)},\cdots,y_n^{(3)}\}$, etc. From those counterfactuals, we can now get a collection of estimators, $\widehat{\theta}^{(1)}$,$\widehat{\theta}^{(2)}$, $\widehat{\theta}^{(3)}$, etc. Instead of using mathematical tricks to calculate $\mathbb{E}(\widehat{\theta})$, compute $$\frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}$$That’s what we’ve seen last friday. I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students, 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)] X=Davis$height[Davis$sex=="M"] We can visualize its distribution (density and cumulative distribution) u=seq(155,205,by=.5) par(mfrow=c(1,2)) hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) lines(u,dnorm(u,178,6.5),col="black") Xs=sort(X) n=length(X) p=(1:n)/(n+1) plot(Xs,p,type="s",col="blue") lines(u,pnorm(u,178,6.5),col="black") Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a $\mathcal{N}(178,6.5^2)$. In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) Y=rnorm(n,178,6.5) hist(Y,col=rgb(1,0,0,.3)) lines(density(Y),col="red",lwd=2) Ys=sort(Y) plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205)) polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA) lines(Xs,p,type="s",col="blue",lwd=2) lines(Ys,p,type="s",col="red",lwd=2) We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution $\widehat{F}$ obtained from $\{y_1,\cdots,y_n\}$ (the blue curve), and some conterfactual, $\widehat{F}^{(s)}$ obtained from $\{y_1^{(s)},\cdots,y_n^{(s)}\}$ generated from $F_{\theta_0}$ – where $\theta_0$ is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance. d=rep(NA,1e5) for(s in 1:1e5){ d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic } ds=density(d) plot(ds,xlab="",ylab="") dks=ks.test(X,"pnorm",178,6.5)$statistic id=which(ds$x&gt;dks) polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA) abline(v=dks,col="red") If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic $\widehat{d}^{(1)}$, $\widehat{d}^{(2)}$, etc, and compare it with the one observe on our sample $\widehat{d}$. The proportion of samples where the test-statistics exceeded the one observed mean(d&gt;dks) [1] 0.78248 is the computational version of the $p$-value ks.test(X,"pnorm",178,6.5) One-sample Kolmogorov-Smirnov test data: X D = 0.068182, p-value = 0.8079 alternative hypothesis: two-sided I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations). # Combining automatically factor levels in R Each time we face real applications in an applied econometrics course, we have to deal with categorial variables. And the same question arise, from students : how can we combine automatically factor levels ? Is there a simple R function ? I did upload a few blog posts, over the pas years. But so far, nothing satistfying. Let me write down a few lines about what could be done. And if some wants to write a nice R function, that would be awesome. To illustrate the idea, consider the following (simulated dataset) n=200 set.seed(1) x1=runif(n) x2=runif(n) y=1+2*x1-x2+rnorm(n,0,.2) LB=sample(LETTERS[1:10]) b=data.frame(y=y,x1=x1, x2=cut(x2,breaks= c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2), labels=LB)) str(b) 'data.frame': 200 obs. of 3 variables:$ y : num 1.345 1.863 1.946 2.481 0.765 ... $x1: num 0.266 0.372 0.573 0.908 0.202 ...$ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ... table(b$x2)[LETTERS[1:10]] A B C D E F G H I J 11 12 23 34 23 36 12 32 3 14 There is one (continuous) dependent variable $y$, one continuous covariable $x_1$ and one categorical variable $x_2$, with here ten levels. We can plot the data using plot(b$x1,y,col="white",xlim=c(0,1.1)) text(b$x1,y,as.character(b$x2),cex=.5)

The output of a linear regression yield the following predictions

for(i in 1:10){ p=function(x) predict(lm(y~x1+x2,data=b),newdata=data.frame(x1=x,x2=LETTERS[i])) u=seq(-1,1.065,by=.01) v=Vectorize(p)(u) lines(u,v)}

the slope for $x_1$ is the same, we simply add a different constant for each level. As we can see, some levels are very very close, so it seems legitimate to combine them into one single category. Here is the output of the linear regression,

summary(lm(y~x1+x2,data=b)) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 0.843802 0.119655 7.052 3.23e-11 *** x1 1.992878 0.053838 37.016 &lt; 2e-16 *** x2A 0.055500 0.131173 0.423 0.6727 x2H 0.009293 0.121626 0.076 0.9392 x2F -0.177002 0.121020 -1.463 0.1452 x2B -0.218152 0.130192 -1.676 0.0955 . x2D -0.206970 0.121294 -1.706 0.0896 . x2G -0.407417 0.129999 -3.134 0.0020 ** x2C -0.526708 0.123690 -4.258 3.24e-05 *** x2J -0.664281 0.128126 -5.185 5.54e-07 *** x2E -0.816454 0.123625 -6.604 3.94e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 0.2014 on 189 degrees of freedom Multiple R-squared: 0.8995, Adjusted R-squared: 0.8942 F-statistic: 169.1 on 10 and 189 DF, p-value: &lt; 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -60.74443 BIC(lm(y~x1+x2,data=b)) [1] -21.16463

Here the reference category is “I”. And it looks like we could actually combine that category with several others. One strategy here would be to select all categories that seem to be not significantly different, and to run a (multiple) test

library(car) linearHypothesis(lm(y~x1+x2,data=b), c("x2A = 0", "x2H = 0", "x2F = 0"))   Hypothesis: x2A = 0 x2H = 0 x2F = 0   Model 1: restricted model Model 2: y ~ x1 + x2   Res.Df RSS Df Sum of Sq F Pr(&gt;F) 1 192 8.4651 2 189 7.6654 3 0.79971 6.5726 3e-04 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

It seems that we can combine those four categories together.

Here, we can see what’s going on when we change the reference category (actually, loop on all categories)

P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERS[1:10] plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) text(1:10,0,LETTERS[1:10]) text(0,1:10,LETTERS[1:10]) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERS[i],names(p)]=p p=P[LETTERS[i],] idx=which(p&gt;.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p&gt;.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)} We are glad to see that it is symmetric : if “H” should be combined with “I”, “I” should also be combined with “H”. Here black points are related with the 10% $p$-value, and white points the 5% $p$-value. This graph is actually hard to read… And actually, this reminds us of Bertin (1967). Here, we can predefine manually some ordering (we will see below how it might be automatised) LETTERSord=c("I","A","H","F","B","D","G","C","J","E") P=matrix(NA,nlevels(b$x2),nlevels(b$x2)) colnames(P)=rownames(P)=LETTERSord plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5), ylim=c(0,10.5)) ct=c(3,3,2,1,1) abline(v=.5+c(0,cumsum(ct)),lty=2) abline(h=.5+c(0,cumsum(ct)),lty=2) text(1:10,0,LETTERSord) text(0,1:10,LETTERSord) for(i in 1:nlevels(b$x2)){ #levels(b$x2)=LETTERS[1:10] b$x2=relevel(b$x2,LETTERSord[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,3) P[LETTERSord[i],names(p)]=p p=P[LETTERSord[i],] idx=which(p&gt;.05) points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2) idx=which(p&gt;.1) points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2) }

Here we get the following

It looks like we have our combined categories…

Actually, it is possible to use another strategy. We start from some level, say “A”. Then, we merge it with all non-significantly different levels. If “B” is not one of them, we use it as the new reference. Etc.

for(i in 1:nlevels(b$x2)){ if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p&gt;.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }} The final categories are table(b$x2)   A+I+H B+D+F C+G E J 46 82 35 23 14

with the following regression output

summary(lm(y~x1+x2,data=b))   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 0.86407 0.03950 21.877 &lt; 2e-16 *** x1 1.99180 0.05323 37.417 &lt; 2e-16 *** x2B+D+F -0.21517 0.03699 -5.817 2.44e-08 *** x2C+G -0.50545 0.04528 -11.164 &lt; 2e-16 *** x2E -0.83617 0.05128 -16.305 &lt; 2e-16 *** x2J -0.68398 0.06131 -11.156 &lt; 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 0.2008 on 194 degrees of freedom Multiple R-squared: 0.8975, Adjusted R-squared: 0.8948 F-statistic: 339.6 on 5 and 194 DF, p-value: &lt; 2.2e-16 AIC(lm(y~x1+x2,data=b)) [1] -66.76939 BIC(lm(y~x1+x2,data=b)) [1] -43.68117

Which is consistent with the group we got before. But actually, if we change the order, we can get different combinations. For instance, if we go from “J” to “A”, instead of “A” to “J”, we obtain

for(i in nlevels(b$x2):1){ #levels(b$x2)=LETTERS[1:10] if(LETTERS[i]%in%levels(b$x2)){ b$x2=relevel(b$x2,LETTERS[i]) p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4] names(p)=substr(names(p),3,nchar(p)) idx=which(p&gt;.05) mix=c(LETTERS[i],names(p)[idx]) b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep="")) }} table(b$x2) E G+C I+A+B+D+F+H J 23 35 128 14 with different information criteria here AIC(lm(y~x1+x2,data=b)) [1] -36.61665 BIC(lm(y~x1+x2,data=b)) [1] -16.82675 I guess it would be necessary to run randomly the order we go through the levels. Last, but not least, one can use regression trees (even if it not per se in the syllabus of the course). The problem is that there is another explanatory variable that might interphere. So I would suggest (1) to fit a linear model $y=\beta_0+\beta_1x_1+u_i$, to calculate the residuals, $\widehat{u}_i$ (2) to run a regression tree, to explain $\widehat{u}_i$ with categorical variable $x_2$ (I did explain how trees are build when the explanatory variable is a categorical one in a previous post) library(rpart) library(rpart.plot) b$e=residuals(lm(y~x1,data=b)) arbre=rpart(e~x2,data=b) prp(arbre,type=2,extra=1)

Observe that the leaves have the same groups as the one we got.

arbre n= 200   node), split, n, deviance, yval * denotes terminal node   1) root 200 22.563500 7.771561e-18 2) x2=G,C,J,E 72 4.441495 -3.232525e-01 4) x2=J,E 37 1.553520 -4.578492e-01 * 5) x2=G,C 35 1.509068 -1.809646e-01 * 3) x2=I,A,H,F,B,D 128 6.366628 1.818295e-01 6) x2=F,B,D 82 2.983381 1.048246e-01 * 7) x2=I,A,H 46 2.030229 3.190993e-01 *

I guess that it should be possible to put all that in an R function, to suggest combinations of level that might improve the regression.

# Traitement des valeurs manquantes : remplacer les NA par une constante ?

Un rapide billet pour répondre à une question posée à la fin du cours de ce matin (ST5100), par Jean-Pierre Liégeois, jeune lecteur du Var (pour préserver un peu d’anonymat)

Dans mon stage, quand on avait des valeurs manquantes, on me disait de remplacer par -1, puis de rajouter une indicatrice comme quoi la variable vaut -1. Ça permet de ne supprimer ni variable, ni observations. On peut faire ça ?

Si je formalise un peu, on va simuler ici des données, disons $x_1$ et $x_2$, on génère ensuite des données suivant un modèle, de la forme $y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon$. Une proportion $\alpha$ de $x_1$ sera transformée en NA.  Ce que suggérais Jean-Pierre, c’est de remplacer les valeurs manquantes par -1, puis d’ajuster un modèle $$y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon$$. Côté code, c’est assez simple. Par défaut, la stratégie de R est de supprimer les valeurs manquantes. Si 50% des données de $x_1$ sont manquantes, la moitié des lignes sont supprimées

n=1000 x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.05 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1) base$x1[indice]=NA reg=lm(y~x1+x2,data=base) Au lieu de générer un unique échantillon, on va en simuler 10,000, et regarder la distribution de $\widehat{\beta}_1$, m=10000 B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red")

Bien sur, avec un taux de valeurs manquantes plus faibles – disons $\alpha=5\%$ – on perd moins d’observations, et donc l’estimateur a une variance plus faible.

Tentons maintenant la stratégie consistant à remplacer les valeurs manquantes par des valeurs numériques fixes, et de rajouter une indicatrice,

B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") Ce qui ne change pas grand chose, on en conviendra… y compris si le taux de valeurs manquantes passe à 5%, On peut se demander ce qui se passe si le shift n’est plus de 1 mais de 10 (a priori, c’est arbitraire, la variable $x_1$ pouvant être plus ou moins dispersée… -1 pour une variable entre 0 et 1, ou entre 0 et 1000, ça n’est pas tout à fait pareil). Mais non, par exemple avec toujours 5% de valeurs manquantes, on a Si on regarde notre échantillon, en particulier le nuage de points $(x_1,y)$, on observe ici, les valeurs manquantes sont choisies au hasard, de manière totalement indépendante, x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha)) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr) (ici avec 1/3 de valeurs manquantes, en rouge). Mais on pourrait supposer que les valeurs manquantes sont les plus grandes valeurs de $x_1$, par exemple, x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha),prob = x1^3) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr) On peut se demander ce que ça donnerait sur l’estimateur $\widehat{\beta}_1$ Ça ne change pas grand chose, mais on a plus de variance, si on regarde bien. Dernier essai : que se passe-t-il si les variables $x_1$ et $x_2$ sont maintenant corrélées, B=rep(NA,m) library(mnormt) r=.8 S = matrix(c(1,r,r,1),2,2) for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red")

Cette fois, on a un estimateur biaisé (de l’ordre de 10% sur cet exemple numérique). Manifestement, cette technique n’est pas très concluante…

Je pourrais ajouter que cette méthode ne revient pas à la première, même si la distribution des estimateurs est proche

set.seed(1) x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg1=lm(y~x1+x2+I(x1==(-1)),data=base) coefficients(reg1) (Intercept) x1 x2 I(x1 == (-1))TRUE 1.0988005 1.7454385 -0.5149477 3.1000668 base$x1[indice]=NA reg2=lm(y~x1+x2,data=base) coefficients(reg2) (Intercept) x1 x2 1.1123953 1.8612882 -0.6548206

Comme je le disais (lors de la discussion qui a suivi le cours) une méthode plus prometteuse est l’imputation. L’idée est de prédire une valeur pour les $x_1$ manquants. On pourrait être tenté de prendre  $\overline{x}_1$, la moyenne des $x_1$ observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de $x_1$, ici, donc on doit pouvoir faire mieux ! On sait aussi que $x_1$  et $x_2$ sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si $x_2$ est grand, on sait que le $x_1$ (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, $x_1=\alpha_0+\alpha_2x_2+\eta_i$, calibré sur les valeurs non-manquantes. Puis on utilise $\widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2$ pour les valeurs manquantes. C’est simpliste, mais pourquoi pas ? On estime alors le modèle sur cette nouvelle base.

for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red")

sur l’exemple numérique, on obtient

base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg3=lm(y~x1+x2,data=base) coefficients(reg3) (Intercept) x1 x2 1.1593298 1.8612882 -0.6320339

Cette méthode a au moins pu corriger du biais…

Après, si on regarde attentivement, on a exactement la même valeur qu’avec la première méthode qui consiste à supprimer les lignes avec des valeurs manquantes !

summary(reg3)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.15933 0.06649 17.435 &lt; 2e-16 *** x1 1.86129 0.21967 8.473 &lt; 2e-16 *** x2 -0.63203 0.20148 -3.137 0.00176 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 1.051 on 997 degrees of freedom Multiple R-squared: 0.1094, Adjusted R-squared: 0.1076 F-statistic: 61.23 on 2 and 997 DF, p-value: &lt; 2.2e-16   summary(reg2)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.11240 0.06878 16.173 &lt; 2e-16 *** x1 1.86129 0.21666 8.591 &lt; 2e-16 *** x2 -0.65482 0.20820 -3.145 0.00172 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 1.037 on 797 degrees of freedom (200 observations deleted due to missingness) Multiple R-squared: 0.1223, Adjusted R-squared: 0.12 F-statistic: 55.5 on 2 and 797 DF, p-value: &lt; 2.2e-16

Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d’imputation. Par exemple prendre la moyenne des $k$ valeurs de $x_1$ (observées) pour les $k$ individus ayant des $x_2$ les plus proches du $x_2$ de l’individu ayant $x_1$ manquant :

kpp=function(i,basena,k=5){ x2=basena$x2[i] sb=basena[!is.na(basena$x1),] idx=rank(abs(sb$x2-x2)) mean(sb[which(idx&lt;=k),"x1"]) } Sur notre base simulée on obtient base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg4=lm(y~x1+x2,data=base) coefficients(reg4) (Intercept) x1 x2 1.197944 1.804220 -0.806766 Si on regarde ce que ça donne sur nos 10,000 simulations, on a (c’est un peu long, car j’ai codé ça très rapidement, et pas du tout de manière optimale) for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg=lm(y~x1+x2,data=base0) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") Le biais semble ici plus faible que sans imputation… autrement dit, les méthodes d’imputation me semblent plus robustes que cette stratégie visant à remplacer des NA par une valeur arbitraire, et de rajouter une indicatrice dans la régression. # 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) [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(&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.

# Modélisation “appliquée”

Depuis la rentrée, je donne un cours de modélisation appliquée, ou pour être plus précis, intitulé “modèles linéaires appliqués“. Ça a été l’occasion de me rendre compte que le mot “appliqué” pouvait prêter à confusion. Dans ma tête, on allait appliquer une théorie (la théorie des modèles linéaires) à cas concrets. Généralement, je présente les grandes lignes de la théorie : les hypothèses, les principaux résultats, et je fais des démonstrations “avec les mains” pour expliquer de manière intuitive certains passages (genre la perte des degrés de liberté dans le test de Student, le signe de la covariance des estimateurs dans le modèle simple, etc) en renvoyant à des livres de référence pour les preuves complètes. Et j’essaye, autant que faire se peut, d’illustrer sur des données, avec du code R (par exemple vendredi dernier, c’était sur la notion de “modalité de référence” quand on régresse sur une variable catégorielle, montrant que tous les modèles sont équivalent et que ce n’est vraiment qu’une histoire de conventions). J’ai aussi donné un gros devoir, avec des vraies données, 2400 observations, 34 variables (quantitatives mais aussi qualitatives). Bref, un vrai cas d’étude. Et j’ai même proposé de revenir sur leurs modélisations en fin de cours…

Et j’ai été surpris quand un groupe d’élèves est venu me voir à la pause, pour me demander s’il y aurait des “chiffres” : si au lieu d’écrire $${\displaystyle {\hat {\beta }}_{1}={\frac {\sum x_{i}\sum y_{i}-n\sum x_{i}y_{i}}{\left(\sum x_{i}\right)^{2}-n\sum x_{i}^{2}}}={\frac {\sum (x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum (x_{i}-{\bar {x}})^{2}}}}$$ou encore $${\displaystyle {\hat {\beta }}_{0}={\frac {\sum y_{i}-{\hat {\beta }}_{1}\sum x_{i}}{n}}={\bar {y}}-{\hat {\beta }}_{1}{\bar {x}}}$$il y aurait des calculs, à la calculatrice, avec des chiffres au lieu des $x$ et des $y$. J’ai été très supris, parce qu’on est en 2018, que tout le monde parle de “big data” et revenir à la calculatrice pour calculer des sommes de 10 nombres, c’est possible. En tous cas, certains profs le font encore. J’ai dit que si ça pouvait les rassurer, je pouvais chercher des bases dans des vieux livres, mais que c’était pas ma manière de comprendre le terme “appliqué” dans le nom du cours. Je me suis souvenu du livre de Robert Pindyck et Daniel Rubinfeld, Econometric Models and Economic Forecasts, qui était le livre de référence de la SOA il y a 20 ans, qui avait des tableaux de ce genre. Mais je ne l’ai pas trouvé en ligne. Par contre, je suis tombé sur la Table 4.3. dans le livre de John Neter, Applied Linear Regression Models, paru en 1974 (repris dans la Table 3.4. de Applied Linear Statistical Models)

On a ici juste 11 observations, un modèle linéaire simple. Il est possible de faire tous les calculs à la main. Mais je me suis dit que ce n’était pas un cadeau : on ne peut pas faire grand chose sans faire un minimum de maths… Faisons juste une régression linéaire sur ces données, pour “voir”

base = data.frame(x=c(125,100,200,75,150,175,75,175,125,200,100), y=c(160,112,124,28,152,156,42,124,150,104,136)) plot(base,xlab="Size of Minumum Deposit", ylab="Number of New Accounts",pch=19,ylim=c(30,180)) abline(lm(y~x,data=base),col="red")

La droite de régression (qui minimise la somme des carrés des erreurs) est la courbe (ou ici “droite”, on est sur un modèle linéaire) rouge

summary(lm(y~x,data=base))   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 50.7225 39.3979 1.287 0.23 x 0.4867 0.2747 1.772 0.11

On ne voit rien, littéralement ! Surtout que sur cet exemple, la question qui nous brule les lèvres, c’est “où est le maximum ?” car on voit bien qu’on a une courbe qui semble concave, croissante au début, décroissante sur la fin… On peut faire du lissage non-paramétrique pour s’en convaincre

with(base, scatter.smooth(x, y, lpars = list(col = "red"), xlab="Size of Minumum Deposit", ylab="Number of New Accounts", ylim=c(30,180),pch=19))

Peut-aller plus loin, pour répondre à la question “où est le maximum ?”. En particulier, peut-on suggérer une valeur, voir trouver un intervalle de confiance ?

On peut continuer avec un modèle paramétrique – l’essentiel de ce qu’on voit dans les cours de statistique tourne autour des modèles paramétrique. Mais si on est plus pragmatique, on peut considérer un modèle quadratique, autrement dit, notre prévision sera une parabole. Ça pourrait être cohérent avec ce qu’on voit quand on trace les points, et puis souvent, les paraboles ont un maximum

with(base, scatter.smooth(x, y, lpars = list(col = "red",lty=2), xlab="Size of Minumum Deposit", ylab="Number of New Accounts", ylim=c(30,180),pch=19)) reg=lm(y~x+I(x^2),data=base) vx=seq(min(base$x),max(base$x),length=251) vy=predict(reg,newdata=data.frame(x=vx)) lines(vx,vy,col="red")

On peut voir que le modèle pourrait sembler réaliste, non seulement visuellement, mais en plus si on regarde la sortie de la régression (même si on pourrait discuter ici l’interprétation de la significativité…)

summary(reg)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) -3.255e+02 5.589e+01 -5.824 0.000394 *** x 6.569e+00 8.744e-01 7.513 6.84e-05 *** I(x^2) -2.203e-02 3.143e-03 -7.011 0.000111 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Or on peut montrer que pour une parabole de la forme $y=\beta_2x^2+\beta_1x+\beta_0$ l’optimum (qui sera un maximum si $\beta_2<0$ – la courbe part en $-\infty$ à gauche et à droite) sera obtenu en $x^\star=\theta=-\beta_1/2\beta_2$. Un estimateur naturel de $\theta$ est alors $\widehat{\theta}=-\widehat{\beta}_1/2\widehat{\beta}_2$ où les estimateurs sont obtenus ici par minimisation de la somme des carrés des erreurs. Mais comment avoir la variance de cet estimateur ? On peut naturellement tenter la $\Delta$-méthode (décrite sur wikipedia), en considérant $$\displaystyle{\theta=g(\beta_1,\beta_2)=\frac{-\beta_1}{2\beta_2}}$$Comme $$\displaystyle{\frac{\partial \theta}{\partial \beta_1}=\frac{-1}{2\beta_2}}$$et $$\displaystyle{\frac{\partial \theta}{\partial \beta_2}=\frac{\beta_1}{2\beta_2^2}}$$la variance (asymptotique) vaut ici $$\left[ \frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right] \left[ \begin{array}{cc} \sigma_{1}^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{array} \right] \left[\frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right]^{\text{ T}}$$soit $$\frac{\sigma_1^2\beta_2^2-2\beta_1\beta_2\sigma_{12}+\beta_1^2\sigma_2^2}{4\beta_2^2}$$On obtient alors un estimateur de cette variance asymptotique en remplaçant les valeurs inconnues par leurs estimateurs. Ici, on aurait, numériquement

beta=coefficients(reg) n=nrow(base) sigma=n*vcov(reg)[2:3,2:3] dg = c(-.5/beta[3],.5*beta[2]/beta[3]^2) theta=-beta[2]/(2*beta[3]) theta [1] 149.0676 s2=t(dg) %*% sigma %*% dg s2 [,1] [1,] 94.59911 sqrt(s2) [,1] [1,] 9.726207

Autrement dit, si on suppose la normalité de notre estimateur (ce que nous dit d’ailleurs la $\Delta$-méthode, tout du moins asymptotiquement), on a l’intervalle de confiance suivant

theta+qt(c(.025,.975),n-3)*sqrt(s2) [1] 126.6389 171.4963 vx=theta vy=predict(reg,newdata=data.frame(x=vx)) points(vx,vy,pch=19,cex=1.5,col="red") arrows(vx-qt(.975,n-3)*sqrt(s2),vy, vx+qt(.975,n-3)*sqrt(s2),vy,code=3,angle=90,length=.1, col="red",lwd=1.2)

On peut aussi tenter une autre stratégie, basée sur l’utilisation de la vraisemblance profilée. Je renvoie à mon précédent billet pour expliquer l’idée de la méthode, malheureusement pas assez enseignée dans les cours de statistique… Pour notre modèle quadratique, on maximisait la log-vraisemblance, classiquement, sous une hypothèse de résidus Gaussiens et homoscédastiques

logL = function(pm){ b0=pm[1] b1=pm[2] b2=pm[3] b3=pm[4] -sum(log(dnorm( (base$y-(b0+b1*base$x+b2*base$x^2)) ,0,b3))) } Ici, la première astuce est d’introduire $\theta$ (là où se trouve le maximum de la parabole) comme un des paramètres du modèle – par exemple à la place de $\beta_2$ logL = function(pm){ b0=pm[1] b1=pm[2] theta=pm[3] b3=pm[4] -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3))) }

Si on cherche juste le maximum de la vraisemblance, on obtient

 optim(par=c(-323,6.5,140,5),logL) $par [1] -325.521244 6.569328 149.070597 13.693488$value [1] 44.39629

ce qui est cohérent avec notre calcul précédant (et c’est normal, c’est la propriété du cours de statistique qui dit qu’on peut reparamétrer le modèle, ça ne change rien). La seconde astuce est celle de la vraisemblance profilée : on dit que dans notre paramètre multivarié, un est plus important que les autres. Et les autres sont juste maximisés, à $\theta$ donné. Techniquement, la log-vraisemblance profilée du paramètre $\theta$ est

logL=function(theta){ f=function(pm){ b0=pm[1] b1=pm[2] b3=pm[3] -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3))) } optim(par=c(-323,6.5,5),f)$value }

On peut d’ailleurs la tracer

v1 = seq(130,190) v2 = Vectorize(logL)(v1) plot(v1,-v2,type="l",xlim=range(base$x),xlab="",ylab="log vraisemblance") Le maximum est atteint, ici en opt$minimum [1] 149.068

Ce qui est cohérent avec nos calculs. On peut ensuite un résultat sur le test du rapport de vraisemblance, qui nous dit que la différence entre la log-vraisemblance au maximum et en la vraie valeur suit la moitié du quantile d’une loi du chi-deux à 1 degré de liberté

points(opt$minimum,-opt$objective,pch=19,cex=1.5,col="red") ref=-opt$objective-.5*qchisq(.95,1) abline(h=ref,lty=2,col="red") On peut alors alors tracer l’intervalle de confiance sur le graphique initial with(base, scatter.smooth(x, y, lpars = list(col = "red",lty=2), xlab="Size of Minumum Deposit", ylab="Number of New Accounts", ylim=c(30,180),pch=19)) reg=lm(y~x+I(x^2),data=base) vx=seq(min(base$x),max(base$x),length=251) vy=predict(reg,newdata=data.frame(x=vx)) lines(vx,vy,col="red") vx=opt$minimum vy=predict(reg,newdata=data.frame(x=vx)) points(vx,vy,pch=19,cex=1.5,col="red") arrows(min(v1[id]),vy,max(v1[id]),vy,code=3,angle=90,length=.1, col="red",lwd=1.2)

Voilà au moins deux méthodes, basées sur deux résultats du cours de statistique. Après, comme la majorité des techniques de statistique, ce sont des résultats asymptotiques, et rien ne garantie leur validité avec seulement 11 observations (“rien” est un peu fort… on pourrait faire des développements sur ce point).

Une autre solution est d’utiliser des simulations. On suppose que les observations, c’est un modèle, et du bruit. On peut prendre comme modèle le modèle non-paramétrique (lissage local) et supposer le bruit Gaussien. Pour générer d’autres échantillons, on va garder nos observations en $x$, par contre, pour $y$, on va utiliser $\widehat{y}+\varepsilon$$\varepsilon$ sera tiré suivant une loi normale

reg=loess.smooth(x = base$x, y= base$y, evaluation = 501) pred=function(x) reg$y[which.min(abs(x-reg$x))] ypred=Vectorize(pred)(base$x) epsilon = base$y-ypred simu = function(graph=FALSE){ newbase = data.frame(x= base$x, y= ypred + rnorm(n,mean(epsilon),sd(epsilon))) reg=loess.smooth(x = newbase$x, y= newbase$y, evaluation = 501) lines(reg$x,reg$y,col=rgb(1,0,0,.2)) reg$x[which.max(reg$y)]} for(i in 1:20) simu(TRUE) lines(loess.smooth(x = base$x, y= base\$y, evaluation = 501),lwd=1.2,col="red")

On récupère ainsi plein de modèles, estimés sur ces nouvelles données. Compte tenu de l’asymétrie des données, on va oublier le modèle quadratique, et utiliser un modèle non-paramétrique, là encore. Et on calcule numériquement le maximum. Et on répète 10,000 fois.

V=Vectorize(simu)(1:1e5) hist(V,probability = TRUE, col="light blue", border="white",ylim=c(0,.07)) lines(density(V),col="blue") quantile(V,c(.025,.975))

On a ici la distribution empirique du maximum, observée sur nos 10,000 échantillons simulés. On peut même obtenir un intervalle de confiance, en prenant les quantiles empiriques

vy=pred(mean(V)) arrows(quantile(V,.025),vy,quantile(V,.975),vy,code=3,angle=90,length=.1, col="red",lwd=1.2)

Bref, j’ai l’impression que pour faire des applications, il vaut mieux bien connaître la théorie… non ?

# STT5100, Modèles linéaires appliqués

La semaine prochaine commencera le premier cours sur les modèles linéaires appliqués, STT5100. Pas de notes de cours ou de slides/powerpoint, a priori je ferais le cours au tableau, en interagissant autant que possible avec R. Je mettrais sur github les codes et les bases de données qui seront utilises pour illustrer. Il y a aussi des examens des années passées (avec les corrections). Le cours aura lieu les vendredi matin, salle SH-3140. Le blog (via le tag STT5100) permettra d’apporter certains compléments, ici ou la, a priori sur des points techniques.