Monthly Archives: September 2013

2013, ACT2040, Actuarial science, Automne, Course, GLM, Ratemaking

Regression on variables, or on categories?

Leave a comment

I admit it, the title sounds weird. The problem I want to address this evening is related to the use of the stepwise procedure on a regression model, and to discuss the use of categorical variables (and possible misinterpreations). Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db2.txt",header=TRUE,sep=";")

First, let us change the reference in our categorical variable  (just to get an easier interpretation later on)

> db$X3=relevel(as.factor(db$X3),ref="E")

If we run a logistic regression on the three variables (two continuous, one categorical), we get

> reg=glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0758   0.1226   0.2805   0.4798   2.0345  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816    
X3B         -0.10558    0.32526  -0.325   0.7455    
X3C          0.63829    0.37838   1.687   0.0916 .  
X3D         -0.02776    0.33070  -0.084   0.9331    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Now, if we use a stepwise procedure, to select variables in the model, we get

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3

       Df Deviance    AIC
- X3    4   587.81 593.81
<none>      582.29 596.29
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=593.81
Y ~ X1 + X2

       Df Deviance    AIC
<none>      587.81 593.81
- X2    1   606.90 610.90
- X1    1   622.44 626.44

So clearly, we should remove the categorical variable if our starting point was the regression on the three variables.

Now, what if we consider the same model, but slightly different: on the five categories,

> X3complete = model.matrix(~0+X3,data=db)
> db2 = data.frame(db,X3complete)
> head(db2)
  Y       X1       X2 X3 X3A X3B X3C X3D X3E
1 1 3.297569 16.25411  B   0   1   0   0   0
2 1 6.418031 18.45130  D   0   0   0   1   0
3 1 5.279068 16.61806  B   0   1   0   0   0
4 1 5.539834 19.72158  C   0   0   1   0   0
5 1 4.123464 18.38634  C   0   0   1   0   0
6 1 7.778443 19.58338  C   0   0   1   0   0

From a technical point of view, it is exactly the same as before, if we look at the regression,

> reg = glm(Y~X1+X2+X3A+X3B+X3C+X3D+X3E,family=binomial,data=db2)
> summary(reg)

Call:
glm(formula = Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E, family = binomial, 
    data = db2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0758   0.1226   0.2805   0.4798   2.0345  

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.39528    0.86649  -6.227 4.77e-10 ***
X1           0.51618    0.09163   5.633 1.77e-08 ***
X2           0.24665    0.05911   4.173 3.01e-05 ***
X3A         -0.09142    0.32970  -0.277   0.7816    
X3B         -0.10558    0.32526  -0.325   0.7455    
X3C          0.63829    0.37838   1.687   0.0916 .  
X3D         -0.02776    0.33070  -0.084   0.9331    
X3E               NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 806.29  on 999  degrees of freedom
Residual deviance: 582.29  on 993  degrees of freedom
AIC: 596.29

Number of Fisher Scoring iterations: 6

Both regressions are equivalent. Now, what about a stepwise selection on this new model?

> step(reg)
Start:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D + X3E

Step:  AIC=596.29
Y ~ X1 + X2 + X3A + X3B + X3C + X3D

       Df Deviance    AIC
- X3D   1   582.30 594.30
- X3A   1   582.37 594.37
- X3B   1   582.40 594.40
<none>      582.29 596.29
- X3C   1   585.21 597.21
- X2    1   600.56 612.56
- X1    1   617.25 629.25

Step:  AIC=594.3
Y ~ X1 + X2 + X3A + X3B + X3C

       Df Deviance    AIC
- X3A   1   582.38 592.38
- X3B   1   582.41 592.41
<none>      582.30 594.30
- X3C   1   586.30 596.30
- X2    1   600.58 610.58
- X1    1   617.27 627.27

Step:  AIC=592.38
Y ~ X1 + X2 + X3B + X3C

       Df Deviance    AIC
- X3B   1   582.44 590.44
<none>      582.38 592.38
- X3C   1   587.20 595.20
- X2    1   600.59 608.59
- X1    1   617.64 625.64

Step:  AIC=590.44
Y ~ X1 + X2 + X3C

       Df Deviance    AIC
<none>      582.44 590.44
- X3C   1   587.81 593.81
- X2    1   600.73 606.73
- X1    1   617.66 623.66

What do we get now? This time, the stepwise procedure recommends that we keep one category (namely C). So my point is simple: when running a stepwise procedure with factors, either we keep the factor as it is, or we drop it. If it is necessary to change the design, by pooling together some categories, and we forgot to do it, then it will be suggested to remove that variable, because having 4 categories meaning the same thing will cost us too much if we use the Akaike criteria. Because this is exactly what happens here

> library(car)
> reg = glm(formula = Y ~ X1 + X2 + X3, family = binomial, data = db)
> linearHypothesis(reg,c("X3A=X3B","X3A=X3D","X3A=0"))
Linear hypothesis test

Hypothesis:
X3A - X3B = 0
X3A - X3D = 0
X3A = 0

Model 1: restricted model
Model 2: Y ~ X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.1446      0.986

So here, we should pool together categories A, B, D and E (which was here the reference). As mentioned in a previous post, it is necessary to pool together categories that should be pulled together as soon as possible. If not, the stepwise procedure might yield to some misinterpretations.

2013, ACT2040, Automne, Course, Datamining, GLM

ROC curves and classification

31 Comments

To get back to a question asked after the last course (still on non-life insurance), I will spend some time to discuss ROC curve construction, and interpretation. Consider the dataset we’ve been using last week,

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)

The first step is to get a model. For instance, a logistic regression, where some factors were merged together,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)

From this model, we can predict a probability, not a  variable,

> S=predict(reg,type="response")

Let https://latex.codecogs.com/gif.latex?\widehat{S} denote this variable (actually, we can use the score, or the predicted probability, it will not change the construction of our ROC curve). What if we really want to predict a  variable. As we usually do in decision theory. The idea is to consider a threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png, so that

  • if https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-05.png, then  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png will be https://latex.codecogs.com/gif.latex?1, or “positive” (using a standard terminology)
  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-06.png, then  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png will be https://latex.codecogs.com/gif.latex?0, or “negative“

Then we derive a contingency table, or a confusion matrix

     observed value https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png
predicted
value
https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png
 “positive“ “négative“
“positive“ TP FP
“négative“ FN TN

where TP are the so-called true positive, TN  the true negative, FP are the false positive (or type I error) and FN are the false negative (type II errors). We can get that contingency table for a given threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png

> roc.curve=function(s,print=FALSE){
+ Ps=(S>s)*1
+ FP=sum((Ps==1)*(Y==0))/sum(Y==0)
+ TP=sum((Ps==1)*(Y==1))/sum(Y==1)
+ if(print==TRUE){
+ print(table(Observed=Y,Predicted=Ps))
+ }
+ vect=c(FP,TP)
+ names(vect)=c("FPR","TPR")
+ return(vect)
+ }
> threshold = 0.5
> roc.curve(threshold,print=TRUE)
        Predicted
Observed   0   1
       0   5 231
       1  19 745
      FPR       TPR 
0.9788136 0.9751309

Here, we also compute the false positive rates, and the true positive rates,

  • TPR = TP / P = TP / (TP + FN) also called sentivity, defined as the rate of true positive: probability to be predicted positve, given that someone is positive (true positive rate)
  • FPR = FP / N = FP / (FP + TN) is the rate of false positive: probability to be predicted positve, given that someone is negative (false positive rate)

The ROC curve is then obtained using severall values for the threshold. For convenience, define

> ROC.curve=Vectorize(roc.curve)

First, we can plot https://latex.codecogs.com/gif.latex?(\widehat{S}_i,Y_i) (a standard predicted versus observed graph), and visualize true and false positive and negative, using simple colors

> I=(((S>threshold)&(Y==0))|((S<=threshold)&(Y==1)))
> plot(S,Y,col=c("red","blue")[I+1],pch=19,cex=.7,,xlab="",ylab="")
> abline(v=threshold,col="gray")

And for the ROC curve, simply use

> M.ROC=ROC.curve(seq(0,1,by=.01))
> plot(M.ROC[1,],M.ROC[2,],col="grey",lwd=2,type="l")

This is the ROC curve. Now, to see why it can be interesting, we need a second model. Consider for instance a classification tree

> library(tree)
> ctr <- tree(Y~X1+X2+X3bis,data=db)
> plot(ctr)
> text(ctr)

To plot the ROC curve, we just need to use the prediction obtained using this second model,

> S=predict(ctr)

All the code described above can be used. Again, we can plot https://latex.codecogs.com/gif.latex?(\widehat{S}_i,Y_i) (observe that we have 5 possible values for https://latex.codecogs.com/gif.latex?\widehat{S}_i, which makes sense since we do have 5 leaves on our tree). Then, we can plot the ROC curve,

An interesting idea can be to plot the two ROC curves on the same graph, in order to compare the two models

> plot(M.ROC[1,],M.ROC[2,],type="l")
> lines(M.ROC.tree[1,],M.ROC.tree[2,],type="l",col="grey",lwd=2)

The most difficult part is to get a proper interpretation. The tree is not predicting well in the lower part of the curve. This concerns people with a very high predicted probability. If our interest is more on those with a probability lower than 90%, then, we have to admit that the tree is doing a good job, since the ROC curve is always higher, comparer with the logistic regression.

2013, ACT2040, Actuarial science, Automne, Computer, Course

Nice tutorials to discover R

3 Comments

A series of tutorials, in R, by Anthony Damico. As claimed on http://twotorials.com/, “how to do stuff in r. two minutes or less, for those of us who prefer to learn by watching and listening“. So far,

(I guess there is no need now to make my own….)

2013, ACT2040, Actuarial science, Automne, Course, GLM, Tests

Logistic regression and categorical covariates

8 Comments

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
     Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 ** 
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543    
X3D          0.08004    0.34916   0.229 0.818672    
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality . Which means that for someone with characteristics , we predict the following probability

where  denotes the cumulative distribution function of the logistic distribution

For someone with characteristics , we predict the following probability

For someone with characteristics , we predict the following probability

(etc.) Here, if we accept  (against ), it means that modality  cannot be considerd as different from .

A natural idea can be to change the reference modality, and to look at the -values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k])
+ reg = glm(Y~X1+X2+X3,family=binomial,data=db)
+ M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4]
+ } 
> M
             A            B            C            D            E
A           NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05
B 4.771853e-04           NA 4.841204e-04 0.0009474491 4.743636e-01
C 9.225428e-01 0.0004841204           NA 0.7506242347 9.194193e-05
D 8.186724e-01 0.0009474491 7.506242e-01           NA 1.730589e-04
E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589           NA

and if we simply want to know if the -value exceeds – or not – 5%, we get the following,

> M.TF = M>.05
> M.TF
      A     B     C     D     E
A    NA FALSE  TRUE  TRUE FALSE
B FALSE    NA FALSE FALSE  TRUE
C  TRUE FALSE    NA  TRUE FALSE
D  TRUE FALSE  TRUE    NA FALSE
E FALSE  TRUE FALSE FALSE    NA

The first column is obtained when  is the reference, and then, we see which parameter should be considered as null. The interpretation is the following:

  •  and  are not different from 
  •  is not different from 
  •  and  are not different from 
  •  and  are not different from 
  •  is not different from 

Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand)

> library(car)
> db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 ** 
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *  
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

  Res.Df Df  Chisq Pr(>Chisq)
1    996                     
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories,  and , let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis)
> reg=glm(Y~X1+X2+X3bis,family=binomial,data=db)
> summary(reg)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.39439    1.02791  -4.275 1.91e-05 ***
X1           0.51378    0.11138   4.613 3.97e-06 ***
X2           0.20807    0.07234   2.876  0.00402 ** 
X3bisBE      1.94905    0.36852   5.289 1.23e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 398.31  on 996  degrees of freedom
AIC: 406.31

Number of Fisher Scoring iterations: 7

Here, all the categories are significant. So we do have a proper model.

2013, ACT2040, Actuarial science, Automne, Course, GLM, Premium principle, Ratemaking

Régression de Poisson

Leave a comment

Mercredi, on finira les arbres de classification, et on commencera la modélisation de la fréquence de sinistre. Les transparents sont en ligne.

Comme annoncé lors du premier cours, je suggère de commencer la lecture du Practicionner’s Guide to Generalized Linear Models. Le document correspond au minimum attendu dans ce cours.

Academe, Research, Statistics

De la significativité (statistique)

11 Comments

On va encore dire que je fais mon numéro d’anti-journaliste primaire, mais je voulais revenir sur une histoire, relayée par certains journaux (français) qui me laisse songeur sur les liens entre les journalistes et le milieu scientifique. Je ne reviendrais pas sur l’information que j’avais relayée alors (que @Tournyol a reprise sur Twitter ce week-end, et qui a beaucoup circulé ces dernières heures sur Twitter) où plusieurs journalistes traduisaient “augmenté de 50%” par “doublé” (j’étais revenu sur ces raccourcis en juillet dernier). Je voulais revenir sur cette information, évoquée en une du Parisien il y a 15 jours, en bas à gauche

 

où on apprend

Je n’invente rien: le titre est “la cigarette électronique aide mieux que le patch à arrêter de fumer“. Bon, les plus fidèles lecteurs de mon blog me feront remarquer que je pourrais citer d’autres sources quand j’essaye de relayer une étude scientifique. Soit. Le même jour, Le Monde en parle, sur son site

Ici, on parle de “réduire le tabac” (ça doit être la consommation), mais l’idée est la même: la cigarette électronique est plus efficace que le patch. Mais cette fois, si on prend le temps de lire le détail de l’article (ce qui peut m’arriver à l’occasion, par excès de zèle), on trouve ce paragraphe surprenant,

Sur un test, mené sur plusieurs personnes, il y a une plus grande efficacité, mais statistiquement non significative ! Ben, comment dire… c’est qu’il n’y a pas d’efficacité alors. Avant de m’énerver, on peut alors voir ce que dit l’article en question ?

On retrouve clairement écrit qu’il n’y a pas d’effet statistiquement significatif, et si on regarde un peu l’analyse, on voit qu’effectivement, il y a un léger mieux avec la cigarette électronique – en témoigne la distribution du temps avant la rechute –

mais comme le répètent à plusieurs reprises les chercheurs, l’effet n’est pas statistiquement significatif. On peut aussi aller voir un peu plus loin dans la revue, car après l’article, il y a un commentaire.

Le commentaire est très éclairant, mais je retiendrais un paragraphe,

Je me souviens avoir évoqué il y a quelques mois des études qui disaient que les études qui concluaient qu’il n’y avaient pas de différence significatives entre deux produits avaient plus de mal à être publiées. Et je trouve bien que ce genre d’étude paraissent. Mais je trouve inquiétant la lecture qui en est faite ensuite par les quotidiens. J’ai beaucoup de mal à comprendre ce qui s’est passé ici. Par exemple, pour l’article du Monde, il n’y a pas d’auteur à l’article. Je trouve cela gênant. Car je me demande s’il s’agit d’un article écrit par quelqu’un qui n’a pas compris ce qu’il a lu ? ou est-ce un publi-reportage d’un fabriquant de ces cigarettes électroniques ? J’ai beaucoup de mal à comprendre comment on est passé de cet article scientifique en anglais, à ces phrases en français dans des journaux non-spécialisés.

Cela dit, si on peut commencer à faire des articles d’une page dans un quotidien sur des effets non-significatifs, je peux en suggérer des pistes d’articles ! Dans mes modèles de régressions, j’ai toujours plein de variables non significatives. Si je pouvais les garder pour raconter des histoires, on pourrait rigoler.

2013, ACT2121, Course, Exams, Probability

Loi de Poisson

Leave a comment

Suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur le thème 6 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

Pour rappels, l’intra du 27 septembre portera sur les thèmes 1-6, c’est à dire sur ce qui a été abordé dans les feuilles mises en ligne.

Geometry, Maths, Optimization

Newton-Raphson avec un dessin

1 Comment

Dans un précédant billet, publié il a quelques mois maintenant, je trouvais dommage que l’on oublie de citer des personnes qui ont introduit des concepts fondamentaux. En l’occurence, je pensais à Bruno de Finetti. En traînant sur Google Books l’autre jour, j’ai eu l’impression inverse, en trouvant dommage que l’on ne se souvienne que de la personne qui a introduit un concept, et pas celui qui en a compris la portée générale. Mais je devrais peut-être revenir un peu sur le contexte.

En travaillant avec un étudiant cet été, on devait résoudre des problèmes d’optimisation (en l’occurence l’estimation de paramètres dans un mélange de lois par maximum de vraisemblance) et je me souviens encore lui avoir dit, “c’est facile, tu programmes une descente de gradient, la méthode de Newton-Raphson, et c’est bon, on réglera l’histoire des contraintes plus tard” (oui, on avait un ensemble – simple – de contraintes qu’on pouvait intégrer en reparamétrant le programme d’optimisation). La méthode de Newton repose sur l’idée que pour trouver le zéro d’une fonction, on utilise une suite définie par récurence,

La fameuse illustration de cette formule étant un dessin de la forme

De là, si on se souvient que la condition du premier ordre pour cherche un extremum d’une fonction revient à annuler la dérivée, on peut utiliser cette méthode pour chercher un extremum, avec une suite définie par

Cette méthode se généralise en dimension quelconque, avec la méthode dite de la descente de gradient

Bref, ces méthodes sont très géométriques. Et j’avais toujours cru que les vieux mathématiciens étaient tous des géomètres. Aussi, j’espérais trouver plein de dessins dans les textes de Newton et Raphson. Mais non….

La méthode de Newton a été introduite par Isaac Newton dans De analysi per aequationes numero terminorum infinitas, écrit en 1669 et publié en 1711. Mais c’est surtout dans De metodis fluxionum et serierum infinitarum, écrit en 1671, traduit (du latin) et publié sous le titre Methods of Fluxions en 1736 que l’on va trouver l’idée de la méthode dite de Newton pour trouver les zéros d’une fonction à valeurs réelles. C’est ce texte que l’on peut trouver en ligne, numérisé. Mais Newton n’a aucunement proposé une méthode générale, car il ne l`appliquait qu’aux seuls polynômes. A sa décharge, la notion de dérivée n’était pas (clairement) définie à cette époque… donc imaginer une descente en suivant la tangente, en 1710, ne serait peut-être pas raisonnable.

C’est Thomas Simpson qui proposa de généraliser cette méthode pour mettre en place une méthode de calcul itératif des solutions d’une équation non linéaire, en utilisant les dérivées (qu’il appelait fluxions, comme Newton).  En 1690,  Raphson publie une description de sa méthode dans Analysis aequationum universalis. Comme Newton, il va proposer une méthode de calcul récursif des approximations successives d’un zéro d’un polynôme (avec une petite nuance, car Newton essayait de construire une suite de polynômes, comme l’expliquent Ypma (1995) ou Deuflhard (2005))

Maintenant, si on regarde les publications, le dessin avec la tangente, on ne le voit pas ! Loin de là ! On voit des pages de calculs algébriques,

Il m’a fallu un certain temps pour découvrir un traité abordant ce problème avec des dessins. Pour les amateurs de mathématiques marseillaises, on peut découvrir un traité rafraîchissant, datant de 1768, publié d’ailleurs par un homme qui allait devenir par la suite maire de Marseille (à l’époque où les hommes politiques faisaient des maths, comme j’en parlais dans un vieux billet, avec une preuve du théorème de Pythagore proposé par un mathématicien amateur, qui allait devenir ensuite Président des États-Unis d’Amérique), monsieur Jean-Raymond Pierre  Mouraille.

Ce traité est passionnant, avec un style d’une autre époque

Mourailles se positionne clairement dans la lignée de Newton, avec quelques bémols toutefois,

Le plus intéressant, ce sont les pages d’annexes, remplies de dessins, où enfin, on peut découvrir le dessin avec la suite définie par récurrence, et la descente le long de la tangente,

Bref, tout le monde pense à ces dessins quand il entend parler de la méthode de Newton-Raphson, même si aucune représentation géométrique n’est présentée.

Je ne sais pas si Jean-Raymond Pierre  Mouraille a le premier a avoir représenté la recherche de zéros d’une fonction de cette manière, mais il est dommage qu’on n’en parle pas plus souvent…

Probability

Monty Hall (oh no, not again)

3 Comments

Quite frequently, someone on the internet discovers the Monty Hall paradox, and become so enthusiastic that it becomes urgent to publish an article – or a post – about it. The latest example can be http://www.bbc.co.uk/news/magazine-24045598. I won’t blame them, I did the same a few years ago (see http://freakonometrics.hypotheses.org/776, or http://freakonometrics.hypotheses.org/775, in French).

My point today is that the Monty Hall paradox raise an important question, about information. How comes that something to sounds like non-informative can actually be extremely informative. I will not get back on the blue eyes paradox (see http://freakonometrics.hypotheses.org/1963, in French) or the exam paradox (see http://freakonometrics.hypotheses.org/2328, in French one more time), which are related to information, but not with a probabilistic approach. I will stay close to Monty Hall’s paradox today.

This morning, in my probability class, we were looking at a simple exercise (I say simple because it is only the second course of the session). The problem was the following

Consider an urn , with 15 blue balls, and 10 red balls, and an urn , with 10 blue balls, and 15 red balls. We select randomly one urn (with probability 50% for each urn).
We draw a ball, which turns out to be blue, and we put it back in the urn, Now, we draw a (second) ball. What is the probability that this (second) ball is blue?

Please, take your time to read that carefully…

Ready? Your first thought should be that since we put back the ball, after the first draw, it does not change the probabilities, right? So, why did we say that? It is necessary? (about the last question, yes, when something is mentioned in an exercise, we should use it).

Let’s forget about this second ball story, as an introduction to this problem. What was, actually, the probability for the first ball to be blue? Trivially, it was

i.e.

Let us run a code to get that, using simulations:

> n=1000000
> set.seed(1)

First, let us draw the urn, randomly

> urn=sample(1:2,size=n,replace=TRUE)

Then, let us draw the first, and the second ball,

> urns=matrix(c(15,10,10,15),2,2)
> colnames(urns)=c("blue","red")
> sample.urn=(urns[urn,])
> prob.urn=sample.urn/apply(sample.urn,1,sum)
> u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
> u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]

The probability that the first ball was blue is here

> sum(u1=="blue")/n
[1] 0.499953

and for the second one

> sum(u2=="blue")/n
[1] 0.499221

So, indeed, the probability to have a blue ball is 50%. Now, what was the question? Given that the first ball was blue, what it the probability that the second one is blue? Here, on our simulations, it is

> sum(u2[u1=="blue"]=="blue")/sum(u1=="blue")
[1] 0.5194088

Which is close to 52%.And if you run more simulations, you get

> f=function(seed){
+ set.seed(seed)
+ urns=matrix(c(15,10,10,15),2,2)
+ colnames(urns)=c("blue","red")
+ sample.urn=(urns[urn,])
+ prob.urn=sample.urn/apply(sample.urn,1,sum)
+ u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
+ u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]
+ return(sum(u2[u1=="blue"]=="blue")/
+ sum(u1=="blue"))
+ }
> Vectorize(f)(1:20)
 [1] 0.5194088 0.5200931 0.5203338 0.5192104 0.5196960 0.5206121 0.5195453
 [8] 0.5184580 0.5203755 0.5200154 0.5196557 0.5179276 0.5188652 0.5204724
[15] 0.5197437 0.5209244 0.5205770 0.5208725 0.5206228 0.5190711

The probability is always close to 52%, and is (significantly) different from 50%.

Still not convinced that we have some information here that should be used? Imagine that in the first urn, we add 1 blue ball, and 24 red balls; and the opposite in the second one. In that case, if we say that the first ball was blue, it means that it is very likely that the urn chosen was the second one. Let’s look at by it running some simulations

> set.seed(1)
> urns=matrix(c(1,24,24,1),2,2)
> colnames(urns)=c("blue","red")
> sample.urn=(urns[urn,])
> prob.urn=sample.urn/apply(sample.urn,1,sum)
> u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
> u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]

As before, the probability that the second ball is blue is 50% (because of the symmetry actually)

> sum(u2=="blue")/n
[1] 0.500362

But if I tell you that the first one was blue, the probability that the second one is blue becomes

> sum(u2[u1=="blue"]=="blue")/sum(u1=="blue")
[1] 0.9236433

So even if – somehow – we do not change much by replacing the ball in its urn, we do have here some information, since it was mentioned that the ball was blue. And we should use it. Again, the important point is that the sentence was not “we draw a ball and we put it back”, but “we draw a blue ball, and we put it back”. Now, it we do the maths, everything become simple, and clear (as usual).

The question is here to compute

and according to Bayes formula, it is

Now, to compute those two probabilities, we have to condition on the urn,

Given the urn, since we replace the ball,

i.e.

So if we substitute numerical probabilities to get a blue ball in the previous formula, we get

which not the same as

Here, we get

> {(15/25)^2+(10/25)^2}/((15/25)+(10/25))
[1] 0.52

which confirms our empirical 52%, and note that in the second case (where there was only 1 blue ball in one urn, and 24 in the second one)

> {(24/25)^2+(1/25)^2}/((24/25)+(1/25))
[1] 0.9232

which again is close to the empirical 92.3% we got.

I strongly believe that the mis-intuition we might have is close to the one we can observe in Monty Hall paradox. And unless you write things properly, it is difficult to conclude anything….

PS [48  hours later] thanks @mikeandallie for the animated version of my post

2013, ACT2121, Probability

Variables aléatoires discrètes

Leave a comment

Suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveau jeu d’exercices, sur les thèmes 4-5 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

  • Formule de la probabilité totale, et formule de Bayes, #4, et lois discrètes #5 ACT2121-A2013-45.pdf

On fera des exercices sur la loi de Poisson la semaine prochaine, et l’intra du 27 septembre portera sur les thèmes 1-6. On commencera début octobre les variables continues. A suivre…

 

2013, ACT2040, Automne, Computer, econometrics, Insurance, Ratemaking, Regression

Non-observable vs. observable heterogeneity factor

5 Comments

This morning, in the ACT2040 class (on non-life insurance), we’ve discussed the difference between observable and non-observable heterogeneity in ratemaking (from an economic perspective). To illustrate that point (we will spend more time, later on, discussing observable and non-observable risk factors), we looked at the following simple example. Let  denote the height of a person. Consider the following dataset

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")

There is a small typo in the dataset, so let us make manual changes here

> Davis[12,c(2,3)]=Davis[12,c(3,2)] 

Here, the variable of interest is the height of a given person,

> X=Davis$height 

If we look at the histogram, we have

> hist(X,col="light green", border="white",proba=TRUE,xlab="",main="")

Can we assume that we have a Gaussian distribution ?

Maybe not… Here, if we fit a Gaussian distribution, plot it, and add a kernel based estimator, we get

> (param <- fitdistr(X,"normal")$estimate) 
> f1 <- function(x) dnorm(x,param[1],param[2]) 
> x=seq(100,210,by=.2) 
> lines(x,f1(x),lty=2,col="red") 
> lines(density(X))

 

If you look at that black line, you might think of a mixture, i.e. something like

(using standard mixture notations). Mixture are obtained when we have a non-observable heterogeneity factor: with probability , we have a random variable  (call it type [1]), and with probability , a random variable  (call it type [2]). So far, nothing new. And we can fit such a mixture distribution, using e.g.



> library(mixtools) 
> mix <- normalmixEM(X)
 number of iterations= 335 
> (param12 <- c(mix$lambda[1],mix$mu,mix$sigma)) 
[1] 0.4002202 178.4997298 165.2703616 6.3561363 5.9460023  


If we plot that mixture of two Gaussian distributions, we get


> f2 <- function(x){ param12[1]*dnorm(x,param12[2],param12[4])
+ (1-param12[1])*dnorm(x,param12[3],param12[5]) }
> lines(x,f2(x),lwd=2, col="red") lines(density(X))


Not bad. Actually, we can try to maximize the likelihood with our own codes,


> logdf <- function(x,parameter){
+ p <- parameter[1]
+ m1 <- parameter[2]
+ s1 <- parameter[4]
+ m2 <- parameter[3]
+ s2 <- parameter[5]
+ return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2)))
+ }
> logL <- function(parameter) -sum(logdf(X,parameter))
> Amat <- matrix(c(1,-1,0,0,0,0,
+ 0,0,0,0,1,0,0,0,0,0,0,0,0,1), 4, 5)
> bvec <- c(0,-1,0,0)
> constrOptim(c(.5,160,180,10,10), logL, NULL, ui = Amat, ci = bvec)$par

[1]   0.5996263 165.2690084 178.4991624   5.9447675   6.3564746


Here, we include some constraints, to insurance that the probability belongs to the unit interval, and that the variance parameters remain positive. Note that we have something close to the previous output.


Let us try something a little bit more complex now. What if we assume that the underlying distributions have the same variance, namely




In that case, we have to use the previous code, and make small changes,


> logdf <- function(x,parameter){
+ p <- parameter[1]
+ m1 <- parameter[2]
+ s1 <- parameter[4]
+ m2 <- parameter[3]
+ s2 <- parameter[4]
+ return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2)))
+ }
> logL <- function(parameter) -sum(logdf(X,parameter))
> Amat <- matrix(c(1,-1,0,0,0,0,0,0,0,0,0,1), 3, 4)
> bvec <- c(0,-1,0)
> (param12c= constrOptim(c(.5,160,180,10), logL, NULL, ui = Amat, ci = bvec)$par)

[1]   0.6319105 165.6142824 179.0623954   6.1072614


This is what we can do if we cannot observe the heterogeneity factor. But wait… we actually have some information in the dataset. For instance, we have the sex of the person. Now, if we look at histograms of height per sex, and kernel based density estimator of the height, per sex, we have




So, it looks like the height for male, and the height for female are different. Maybe we can use that variable, that was actually observed, to explain the heterogeneity in our sample. Formally, here, the idea is to consider a mixture, with an observable heterogeneity factor: the sex,




We now have interpretation of what we used to call class [1] and [2] previously: male and female. And here, estimating parameters is quite simple,


>  (pM <- mean(sex=="M"))
[1] 0.44
>  (paramF <- fitdistr(X[sex=="F"],"normal")$estimate)
      mean         sd 
164.714286   5.633808 
>  (paramM <- fitdistr(X[sex=="M"],"normal")$estimate)
      mean         sd 
178.011364   6.404001


And if we plot the density, we have


> f4 <- function(x) pM*dnorm(x,paramM[1],paramM[2])+(1-pM)*dnorm(x,paramF[1],paramF[2])
> lines(x,f4(x),lwd=3,col="blue")




What if, once again, we assume identical variance? Namely, the model becomes


Then a natural idea to derive an estimator for the variance, based on previous computations, is to use




The code is here


> s=sqrt((sum((height[sex=="M"]-paramM[1])^2)+sum((height[sex=="F"]-paramF[1])^2))/(nrow(Davis)-2))
> s
[1] 6.015068


and again, it is possible to plot the associated density,


> f5 <- function(x) pM*dnorm(x,paramM[1],s)+(1-pM)*dnorm(x,paramF[1],s)
> lines(x,f5(x),lwd=3,col="blue")


Now, if we think a little about what we’ve just done, it is simply a linear regression on a factor, the sex of the person,




where .  And indeed, if we run the code to estimate this linear model,


> summary(lm(height~sex,data=Davis))

Call:
lm(formula = height ~ sex, data = Davis)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.7143  -3.7143  -0.0114   4.2857  18.9886 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 164.7143     0.5684  289.80   <2e-16 ***
sexM         13.2971     0.8569   15.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.015 on 198 degrees of freedom
Multiple R-squared:  0.5488,	Adjusted R-squared:  0.5465 
F-statistic: 240.8 on 1 and 198 DF,  p-value: < 2.2e-16


we get the same estimators for the means and the variance as the ones obtained previously. So, as mentioned this morning in class, if you have a non-observable heterogeneity factor, we can use a mixture model to fit a distribution, but if you can get a proxy of that factor, that is observable, then you can run a regression. But most of the time, that observable variable is just a proxy of a non-observable one…



	

	
	




	
	

				

			Academe, Internet, Miscellaneous
		

			
Blogging in Academia

		

			 			Leave a comment
						

	


		

		
In a few weeks, I will attend the World Social Science Forum, and participate to a panel committee, on “Minor forms of academic communication: revamping the relationship between science and society?” The Forum will take place at the Palais des Congrès in Montréal.





First developed by physicists, the open access movement has significantly widened in scope and has been taken up by the European Commission and the G8. If the motivations behind this are essentially to do with innovation, competiveness and the economic efficiency of state investments, open access also introduces a major change in the relationship between science and society. Once citizens have access to the results of social science research, in real time and in their entirety, the whole nature of the relationship between science and society is renewed. However, the debate is focused on the major forms of academic communication: journal articles and books. And yet, in ways almost invisible to the academy, so-called minor forms of academic communication are developing in the interstices, creating a kind of “permanent virtual seminar” and intermeshing with promising heuristic, pedagogic and societal possibilities. Over the last decade, research blogs have embodied a new experience of academic communication, allowing for an experimentation of formats, schedules and interactions that differ from those of the traditional academic process, at both its early and later stages. These blogs dovetail with societal questions, open up new frontiers and step outside the academic ivory tower. From “just-in-time sociology” to probing interpretations of the “Arab Spring”, from analysis of contemporary visual culture to knowledge of contemporary movements such as rap, tattoos and religious conversions, academic blogs place the social sciences at the heart of the society they study. Academic blogs have already found a following. Through this readership a democratisation of access to science is underway. This momentum brings opportunities and reveals new ethical, epistemological and scientific questions.


See also  http://oep.hypotheses.org/319 and http://oep.hypotheses.org/322



For an academic, imparting knowledge revolves around two main activities. As a researcher, an academic must produce articles destined for a very limited readership and adhering to a very rigid process, one that involves severe constraints of time and form. As a teacher, an academic must think of his or her students; he or she must make knowledge accessible and arouse curiosity while also seeking to instil the notion of scientific rigour: outlining hypotheses; identifying the results that may be obtained if these hypotheses are borne out; questioning the validity of the said hypotheses. I will come back to the experience of the Freakonometrics blog, which originated as a response to a pedagogical difficulty and is based on two recent ventures. The first was Freakonomics, Stephen Dubner and Steven Levitt’s blog (then books), which aimed to explain economic behaviours using surprising and sometimes provocative examples. The second was the recent explosion of “data visualisation”, an offshoot of open data and big data that showed statistics could be elegant as well as informative. The Freakonometrics blog emerged from a desire to explain in very practical terms how econometric modelling works, by providing (or explaining how to find) data and sharing codes with which to produce graphics. This experience provided the opportunity to publish research studies in an unconventional form. If the rigour expected is the same as that of an article submitted to a peer-reviewed journal, the blog format allows authors to include anecdotes and exploit the rich potential of online documents (animations, links, etc.). Blog posts also extend the notion of reproducible research, the idea being not to impress the reader by making him or her think (s)he is reading something groundbreaking (which we all seek to do in an academic article), but instead to instil the idea of do-it-yourself by allowing all readers to reproduce the analysis.



My paper will discuss the creation of a blog on the Hypotheses platform in relation to my position as a “young academic” and drawing on my experience of academic blogging. I would like to put forward the hypothesis that for a young academic, blogging is a means to liberate oneself from the rules of the academic world. It offers unrivalled editorial freedom and potential academic recognition. I would like to show that, in turn, this academic liberation emancipates knowledge itself, allowing it to reach sectors and readerships beyond those originally intended. I also wish to point out the practical aspects of this interplay between the liberation of the young academic and the liberation of knowledge. First, I will show why for a young academic, blogging is a way to liberate (or at least distance) oneself from the rules of the academic establishment. Secondly, I will attempt to show that by freeing themselves from the academic sphere, young academics impart knowledge, skills and qualities that are useful to everyone.



In a context of financial difficulties and waning influence, the social sciences now place more importance on producing authority than producing knowledge. The new tools of digital micro-publication, which compete with traditional publications yet lack institutional legitimacy, have met with strong resistance. The characteristics that make them such formidable tools for research, communication and academic discussion, and for collective and collaborative work in particular, remain largely unrecognised. In the absence of suitable incentives and training programs, use of these tools is developing virally and falls far short of full potential. Will such tools continue to develop as best they can at the sidelines of the academy? While the humanities shun their social responsibilities, academic blogging will remain an art rather than a science.


See also http://www.arhv.lhivic.org/index.php/2008/09/15/807-why-blog


    

  • A blind spot? Digital infrastructures for digital publishing, and for academic blogging in particular, Marin Dacos
    • 



After several centuries of development, knowledge technologies today form a highly organised ecosystem, structured around books and journals and with its own clearly identified professions, infrastructures and actors. From publishers to librarians, authors to booksellers, a book industry has emerged and encourages the circulation of ideas. With the rise of the network, these roles are slowly being redefined and new actors are rapidly emerging. The 2006 ACLS report (“Our Cultural Commonwealth: The Report of the American Council of Learned Societies Commission on Cyberinfrastructure for the Humanities and Social Sciences”) is one of the first signs of recognition of the need for digital infrastructures. These infrastructures are not simply confined to “noble” publications i.e. books and journals. They also concern the so-called minor forms of academic communication. Yet developing such infrastructures requires much more than simply installing a server under a desk. On the contrary, digital infrastructures necessitate the creation of platforms, which in turn entail the emergence of new teams and new professions – those of digital publishing. These platforms are often developed or bought up by predatory multinationals (for example, Mendeley absorbed by Elsevier). Academic-led alternatives do exist (Zotero for bibliographies, Hypotheses for blogs), yet the academic community has failed to fully recognise the associated opportunities and risks. The academy has every interest in making sure it does not become marginalised within its own infrastructures. The alternative is to reproduce the vagaries of the extraordinarily concentrated global publishing system, which has stripped the research sector of some of its intellectual and budgetary initiative-taking capacities.


See also :  ”Scholarly blogs, a space on the side for academic dialogue” by Marin Dacos and Pierre Mounier. Part 1 and Part 2.


(to be continued…)

	

	
	




	
	

				

			2013, ACT2040, Actuarial science, Automne
		

			
Introduction à la régression logistique et aux arbres

		

			 			Leave a comment
						

	


		

		
Pour le second cours ACT2040, on va finir l’introduction (et les rappels de statistique inférentiel) puis attaque la première grosse section, sur la régression logistique et aux arbres de classification. base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/…






> baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",
+ header=TRUE,sep=",")
> tail(baseavocat)
     CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE  LOSS
1335   34204        2      2       2        2        1     26 0.161
1336   34210        2      1       2        2        1     NA 0.576
1337   34220        1      2       1        2        1     46 3.705
1338   34223        2      2       1        2        1     39 0.099
1339   34245        1      2       2        1        1     18 3.277
1340   34253        2      2       2        2        1     30 0.688






On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire.


Les transparents sont en ligne sur le blog,




Des compléments théoriques sur les arbres peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/….

	

	
	




	
	

				

			Computer, Regression
		

			
Linear regression from a contingency table

		

			 			21 Comments
						

	


		

		
This morning, Benoit sent me an email, about an exercise he found in an econometric textbook, about linear regression. Consider the following dataset,




Here, variable X denotes the income, and Y the expenses. The goal was to fit a linear regression (actually, in the email, it was mentioned that we should try to fit an heteroscedastic model, but let us skip this part). So Benoit’s question was more or less: how do you fit a linear regression from a contingency table?


Usually, when I got an email on Saturday morning, I try to postpone. But the kids had their circus class, so I had some time to answer. And this did not look like a complex puzzle… Let us import this dataset in R, so that we can start playing


> df=read.table("http://freakonometrics.free.fr/baseexo.csv",sep=";",header=TRUE)
> M=as.matrix(df[,2:ncol(df)])
> M[is.na(M)]<-0
> M
      X14 X19 X21 X23 X25 X27 X29 X31 X33 X35
 [1,]  74  13   7   1   0   0   0   0   0   0
 [2,]   6   4   2   7   4   0   0   0   0   0
 [3,]   2   3   2   2   4   0   0   0   0   0
 [4,]   1   1   2   3   3   2   0   0   0   0
 [5,]   2   0   1   3   2   0   6   0   0   0
 [6,]   2   0   2   1   0   0   1   2   1   0
 [7,]   0   0   0   2   0   0   1   1   3   0
 [8,]   0   1   0   1   0   0   0   0   2   0
 [9,]   0   0   0   0   1   1   0   1   0   1


The first idea I had was to use those counts as weights. Weighted least squares should be perfect. The dataset is built from this matrix,


> W=as.vector(M)
> x=df[,1]
> X=rep(x,ncol(M))
> y=as.numeric(substr(names(df)[-1],2,3))
> Y=rep(y,each=nrow(M))
> base1=data.frame(X1=X,Y1=Y,W1=W)


Here we have


> head(base1,10)
   X1 Y1 W1
1  16 14 74
2  23 14  6
3  25 14  2
4  27 14  1
5  29 14  2
6  31 14  2
7  33 14  0
8  35 14  0
9  37 14  0
10 16 19 13


The regression is the following,


> summary(reg1)

Call:
lm(formula = Y1 ~ X1, data = base1, weights = W1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.35569    2.03022   2.145    0.038 *  
X1           0.68263    0.09016   7.572 3.04e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.892 on 40 degrees of freedom
Multiple R-squared:  0.589,	Adjusted R-squared:  0.5787 
F-statistic: 57.33 on 1 and 40 DF,  p-value: 3.038e-09


It looks like the output is the same as what Benoit found, so we should be happy. Now, I had a second thought. Why not create the implied dataset. Using replicates, we should be able to create the dataset that was used to get this contingency table,


> vX=vY=rep(NA,sum(W))
> sumW=c(0,cumsum(W))
> for(i in 1:length(W)){
+ if(W[i]>0){
+ vX[(1+sumW[i]):sumW[i+1]]=X[i]
+ vY[(1+sumW[i]):sumW[i+1]]=Y[i]
+ }}
> base2=data.frame(X2=vX,Y2=vY)


Here, the dataset is much larger, and there is no weight,


> tail(base2,10)
    X2 Y2
172 31 31
173 33 31
174 37 31
175 31 33
176 33 33
177 33 33
178 33 33
179 35 33
180 35 33
181 37 35


If we run a linear regression on this dataset, we obtain


> summary(reg2)

Call:
lm(formula = Y2 ~ X2, data = base2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.35569    0.95972   4.538 1.04e-05 ***
X2           0.68263    0.04262  16.017  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.731 on 179 degrees of freedom
Multiple R-squared:  0.589,	Adjusted R-squared:  0.5867 
F-statistic: 256.5 on 1 and 179 DF,  p-value: < 2.2e-16


If we compare the two regressions, we have


> rbind(coefficients(summary(reg1)),
+ coefficients(summary(reg2)))
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 4.3556857 2.03021637  2.145429 3.804237e-02
X1          0.6826296 0.09015771  7.571506 3.038443e-09

(Intercept) 4.3556857 0.95972279  4.538483 1.036711e-05
X2          0.6826296 0.04261930 16.016913 2.115373e-36


The estimators are exactly the same (which does not surprise me), but standard deviation (ans significance levels) are quite different. And to be honest, I find that surprising. Which approach here is the most legitimate (since they are finally not equivalent)?

	

	
	




	
	

				

			2013, ACT2121, Course
		

			
Probabilités

		

			 			Leave a comment
						

	


		

		
Ce vendredi aura lieu – au local SH-2560 – le premier cours ACT2121, de préparation pour l’examen P de la SOA (probability). Je renvoie à http://beanactuary.org/exams/… pour davantage d’information sur l’examen.


Le plan de cours est en ligne, ainsi que deux jeux d’exercices, sur les thèmes 1-3 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)