Classifieur Naïf Bayésien

Introduction

Le classifieur naïf bayésien est l'une des méthodes les plus simples en apprentissage supervisé basée sur le théorème de Bayes.
il est peu utilisé par les praticiens du data mining au détriment des méthodes traditionnelles que sont les arbres de décision ou les régressions logistiques.
Un avantage de cette méthode est la simplicité de programmation, la facilité d'estimation des paramètres et sa rapidité (même sur de très grandes bases de données). Malgré ses avantages, son peu d'utilisation en pratique vient en partie du fait que ne disposant pas d'un modèle explicite simple (l'explication de probabilité conditionnelle à priori), l'intérêt pratique d'une telle technique et remise en question.

Le classifieur naïf bayésien : L'indépendance conditionnelle

On considère deux variables aléatoires X et Y. Et on essaye de prédire Y à partir de X. Imaginons par exemple que :
Y = “Acheter” ou “ Ne pas acheter” X = “Cher” ou “ Moins cher”

	Acheter	Ne pas acheter
Cher	50	60
Pas cher	120	90

La règle de Bayes indique que :
\(P( Y = Acheter | X =cher ) =P( Y = Acheter et X = cher ) / P( X = cher)\) ou encore \[P( Y = Ne pas acheter | X = cher ) = P( Y = Ne pas acheter et X = cher )/ P( X = cher )\]
La règle de décision sera de dire : Si la première proabilité est plus grande, alors l'achat sera effectué.
Le classifieur naïf bayésien découle directement de cette formule. Ainsi donc, \(P( Y = Acheter | X =cher ) =P(X= cher et Y = Acheter ) * P( Y= Acheter) P( X = cher)\)
Le terme naïf désigne ici le fait que les descripteurs sont conditionnellement indépendants. Ce qui en pratique n'est presque jamais vrai. Mais cela ne remet pas en cause les résultats comme on peut le lire ici : http://www.cs.iastate.edu/~honavar/bayes-lewis.pdf
En résumé :
Pour mettre en place un classifieur naïf de Bayes :

• On détermine un ensemble d'apprentissage
• On détermine des probabilités à priori de chaque classe (par exemple en observant les effectifs)
• On applique la règle de Bayes : \(P(y=c/X=x)= P(X=x/Y=c) * P(Y=c) /(PX=x)\) pour obtenir la probabilité à posteriori des classes au point x.
• On choisit la classe la plus probable.

Dans la formule de Bayes :
\(P(Y=c)\) : La proabilité à priori
\(P(X=x/Y=c)\) : La probailité conditionnelle d'appartenance à la classe c, qui peut être interprété comme la vraisemblance de cette classe.

Mise en oeuvre sous R

Le logiciel R fournit à travers plusieurs packages des fonctions pour mettre en oeuvre une classification à l'aide d'un bayésien naïf.
Nous allons utiliser le package e1071.
Pour une application, nous allons utiliser un jeu de données de G. Saporta qui contient les données relatives aux pronostics de survie suite à un infacrtus des patients selon des cractéristiques observées telles que : la fréquence cardiaque, le nombre de pulsations, etc…

rm(list = ls())
setwd("F:/NaivesBayes")
library(e1071)

## Loading required package: class

# Les données
survie <- read.table(file = "Data/survie.txt", header = T)
head(survie)

##   FRCAR INCAR INSYS PRDIA PAPUL PVENT REPUL  PRONO
## 1    90  1.71  19.0    16  19.5  16.0   912 SURVIE
## 2    90  1.68  18.7    24  31.0  14.0  1476  DECES
## 3   120  1.40  11.7    23  29.0   8.0  1657  DECES
## 4    82  1.79  21.8    14  17.5  10.0   782 SURVIE
## 5    80  1.58  19.7    21  28.0  18.5  1418  DECES
## 6    80  1.13  14.1    18  23.5   9.0  1664  DECES

# Echantillon d'apprentissage & de test
set.seed(1)
n = dim(survie)[1]
index = sample(n, 0.7 * n)
Appren = survie[index, ]
Test = survie[-index, ]
# Modélisation
nb.model <- naiveBayes(PRONO ~ ., data = Appren)

l'implémentation du modèle est simple en utilisant la fonction. Les résultats en sortie sont les suivants

## 
## Naive Bayes Classifier for Discrete Predictors
## 
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
## 
## A-priori probabilities:
## Y
##  DECES SURVIE 
## 0.3878 0.6122 
## 
## Conditional probabilities:
##         FRCAR
## Y         [,1]  [,2]
##   DECES  92.16 14.87
##   SURVIE 87.13 14.69
## 
##         INCAR
## Y         [,1]   [,2]
##   DECES  1.366 0.3519
##   SURVIE 2.353 0.6219
## 
##         INSYS
## Y         [,1]  [,2]
##   DECES  15.01 3.834
##   SURVIE 27.93 9.500
## 
##         PRDIA
## Y         [,1]  [,2]
##   DECES  21.42 5.409
##   SURVIE 15.03 4.560
## 
##         PAPUL
## Y         [,1]  [,2]
##   DECES  28.55 7.695
##   SURVIE 20.93 5.519
## 
##         PVENT
## Y         [,1]  [,2]
##   DECES  11.87 4.216
##   SURVIE  8.75 3.997
## 
##         REPUL
## Y          [,1]  [,2]
##   DECES  1777.0 630.9
##   SURVIE  766.3 292.9

La qualité du modèle dépend de sa capacité à bien classer dans le jeu de données test

PRONOTH <- predict(object = nb.model, newdata = Test)
Test.mod <- cbind(Test, PRONOTH)
head(Test.mod, 5)

##   FRCAR INCAR INSYS PRDIA PAPUL PVENT REPUL  PRONO PRONOTH
## 2    90  1.68  18.7    24  31.0  14.0  1476  DECES   DECES
## 3   120  1.40  11.7    23  29.0   8.0  1657  DECES   DECES
## 5    80  1.58  19.7    21  28.0  18.5  1418  DECES   DECES
## 6    80  1.13  14.1    18  23.5   9.0  1664  DECES   DECES
## 9    78  2.16  27.7    15  20.5  11.5   759 SURVIE  SURVIE

tail(Test.mod, 5)

##    FRCAR INCAR INSYS PRDIA PAPUL PVENT REPUL  PRONO PRONOTH
## 50    75  1.21  16.1    19    24     4  1587  DECES   DECES
## 51    80  2.41  30.9    19    24     7   797 SURVIE  SURVIE
## 64   100  1.76  17.6    23    33     2  1500 SURVIE   DECES
## 70    87  2.51  28.8    16    24    20   765  DECES  SURVIE
## 71   100  2.31  23.1     8    12     1   416 SURVIE  SURVIE

# Taux de bien classé
(Confusion = table(Test.mod$PRONO, Test.mod$PRONOTH))

##         
##          DECES SURVIE
##   DECES      9      1
##   SURVIE     2     10

# En pourcetages:
round(prop.table(Confusion), 2)

##         
##          DECES SURVIE
##   DECES   0.41   0.05
##   SURVIE  0.09   0.45

Soit un taux de bien classé de 86.3636 %
Comparons la classification naïve bayésienne à des méthodes traditionnelles comme une régession logistique

log.model <- glm(PRONO ~ ., data = Appren, family = "binomial")
log.model

## 
## Call:  glm(formula = PRONO ~ ., family = "binomial", data = Appren)
## 
## Coefficients:
## (Intercept)        FRCAR        INCAR        INSYS        PRDIA  
##   -36.27863      0.22251     -0.22697      1.25918     -0.41212  
##       PAPUL        PVENT        REPUL  
##    -0.49393      0.26358      0.00785  
## 
## Degrees of Freedom: 48 Total (i.e. Null);  41 Residual
## Null Deviance:       65.4 
## Residual Deviance: 20.2  AIC: 36.2

LOG.PRONOTH <- predict(object = log.model, newdata = Test, type = "response")
Log.test <- cbind(Test, LOG.PRONOTH = ifelse(LOG.PRONOTH >= 0.5, "SURVIE", "DECES"))
# Taux de bien classé
(LOG.Confusion = table(Log.test$PRONO, Log.test$LOG.PRONOTH))

##         
##          DECES SURVIE
##   DECES      9      1
##   SURVIE     2     10

# En pourcetages:
round(prop.table(LOG.Confusion), 2)

##         
##          DECES SURVIE
##   DECES   0.41   0.05
##   SURVIE  0.09   0.45

Soit un taux de bien classé de 86.3636 %
On peut par ailleurs vérifier si les erreurs sont commises au même endroit, ce qui indiquerait le même comportement des algorithmes : Ce qui est bien le cas

(ConfMethode = table(Log.test$LOG.PRONOTH, Test.mod$PRONOTH))

##         
##          DECES SURVIE
##   DECES     11      0
##   SURVIE     0     11

A faire :

• Tester si ce comportement est identique sur une volumétrie plus importante.
• D'autres classifieurs naïfs
• Implémenter manuellement un classifieur naïf
• Détailler un exemple numérique sur un exemple

What are my chances to talk to this girl? Fisher or Bayes

Robert Mathews said that : "Ronald Fisher gave scientists a mathematical machine for turning baloney into breakthroughs, and ukes into funding. It is time to pull the plug.". He's right.
In one previous life, I wrote a thesis in Philosophy. But, a specific area, Epistemology also called
theory of knowledge, because, It questions what knowledge is and how it can be acquired,
and the extent to which any given subject or entity can be known.

My thesis deal about : The tradition, since Cournot in applying mathematical modelling to social sphere and more specific, how the climate modelling interact with interdiscplinary source of knowledge (mathematics, physics,geography, philosophy).

After reading this : http://www.academia.edu/1075253/Climate_Change_Epistemic_Trust_and_Expert_Trustworthiness 
It seems that the use of bayesian statistics is misunderstood.

Assume two thesis :

• According to the Bayesian's statistic is the science who deal about the degree of of proofs in the observations. That means that, Bayesian statistics self-contained paradigm providing tools and techniques for all statistical problems.
• In the classical frequentist view point of statistical theory, a statistical  procedure is judged by averaging its performance over all possible data

However, the bayesian approach gives prime importance to how a given procedure performs for
the actual data observed in a given situation.

The core of this theory have been formalised by Popper (Karl) with two main principles :

1. Knowledge cannot start from nothing — from a tabula rasa – nor yet from observation. The advance of knowledge consists, mainly, in the modification of earlier knowledge. Although we may sometimes, for example in archaeology, advance through a chance observation,the significance of the discovery will usually depend upon its power to modify our earlier theories.
2. Any probability is a degree of belief about something; It's not a property.That means that any scientific model produce data absolutely or conditionnaly on probability. It's possible to measure how the data can modify the degrees of belief (baye's principle).

One of the result of frequentist theory, may be the most criticized, is the fisher s p-value.


A problem can be to evaluate how the diploma is important to get in the first job

Hypothese :

- The diploma has no effect on the first job. In frequentist, we compute the p-value that can be
interpreted as the probability to  observe a difference at least as  important observed in the data if our hypothesis is true.

Read more: http://www.answers.com/topic/bayesian-statistics#ixzz2UKZDkq3v

So, let us talk about the p-value  and this problem :
If you cross seven times a beautiful girl each day for  10 days at the same place, can you conclude that she is always there? And then, tomorrow, youalways have a chance to speak with her.
We're going to examine both approach : Fisher (based on p_value) or Jeffreys (Bayesian)

• With Fisher, we have :

H0 : p=0.5 and H1 :p=0.5
where p is "the probability to cross the girl". If we reject H0 (fisher) we can conclude that the girl is always there.
the p_value is :

> 2*pbinom(3,10,0.5)

0.34375
which is greater than 5% (the arbitrary threshold that everybody likes). indeed, we cannot conclude... Thank's Fisher; I can't trust you. 

With Bayes, we take our two hypothesis again and we suppose p(H0)=p(H1)=0.5. We have the same chance to cross her every day. With bayesian approach, wa need to compute the likelihood,  wich require, a priori distribution, that can be sum up (in our example) by "your a priori belief of cross again and to talk to her"


Let s compute this :
> (p <- c(0.5,0.6,0.7))
[1] 0.5 0.6 0.7
> (apriori <- c(0.5,0.3,0.2))
[1] 0.5 0.3 0.2
> (vraisemblance <- dbinom(7,10,p))
[1] 0.1171875 0.2149908 0.2668279
> (loi.jointe <- apriori*vraisemblance)
[1] 0.05859375 0.06449725 0.05336559
> (p.y <- sum(loi.jointe))
[1] 0.1764566
> (aposteriori <- loi.jointe/p.y)
[1] 0.3320576 0.3655134 0.3024290
 
So, the conditionnal probability p(H0/y)=is 33% and P(H1/y) = 66% .
There are two out of three chance that the girl is always there. In other words;
I have 66% of chances to cross her tomorrow.
Bayesian approach is hopeful. I like it; I can take my time with the girl next door.
 
 
 
 
 

Analyse discriminante linéaire ou Regression logistique

Supposons que l'on dispose d'iris de Paris (en population >100khabts) et qu'on veuille pouvoir les classer selon leurs caractéristiques sociodémos :

• Population
• taux de chômage
• Etudiants
• CSP
• etc...

Une fois, les iris classés, on se demande si l'on peut transporter cette typologie à une autre grande ville (Lyon) par exemple : Il faudrait alors pouvoir utiliser un modèle d'affectation des iris selon leurs caractéristiques respectives à des classes prédéfinies. Le choix de la méthode dépend en général dans ce cas de ... pas grand chose comme indiqué ici. L'analyse discriminante ou une régression logistique, ou un SVM, pour ne citer que ces trois méthodes, ferait très bien l'affaire; Mais les résultats sont-ils identiques ou mieux, quelle est la méthode qui commet le moins d'erreur et pourquoi?

On va tester deux méthodes ici :

• L'analyse discriminante 
• La régression logistique multinomiale

Le programme :
# Etape 1 : Réduction de la dimension par une ACP pour diminuer le nombre de variables
# Etape 2 : CAH sur les facteurs de  l'ACP --> Définition des classes
# Etape 3: Analyse discriminante pour trouver les règles d'affectation aux classes
# Etape 4: Reg. log. pour trouver les mêmes règles d'affectation
# Etape 5: Comparaison des résultats des méthodes sur un échantillon test et mesure du taux de mal classés




1. Réduction de la dimension

# Etape 1 : Réduction de la dimension
require(ade4)
Vm.acp <-dudi.pca(df=Vm,center=TRUE,scale=TRUE,scan=FALSE)
# Eboulis des valeurs propres
screeplot(Vm.acp,type ="l", main = 'choix du nombre d axes')
nbaxes =3;
Vm.acp <-dudi.pca(df=Vm,center=T,scale=T,nf=nbaxes,scannf=FALSE)
s.corcircle(Vm.acp$co,xax=1,yax=3,clabel=0.7,sub="C. de corr des variables",possub="topright")
summary(Vm.acp)
New.vm<-Vm.acp$li

2. CAH sur les facteurs de l'ACP

# CAH sur les facteurs de l'ACP
preclus<-dist(New.vm)**2 #  La CAH op?re sur les distances
cah.vm<-hclust(d=preclus,method="ward")
# Visualisation de la CAH
par(mfrow=c(1,1))
nbclasses.fi =3
groups<-cutree(cah.vm,k=nbclasses.fi)
plot(cah.vm)
rect.hclust(cah.vm,k=nbclasses.fi,border="blue")

Les résultats de la classif

# Rattachement des groupes pour chacun des iris
> discrim.vm<-cbind(Vm,cluster=as.factor(groups))
> table(discrim.vm$cluster)

  1   2   3 
180 483 250 

3. Analyse discriminante pour affectation dans les classes

> index <- sample(nrow(discrim.vm), nrow(discrim.vm)*.70)
> #Echantillon d'apprentissage
> appren <- discrim.vm[index, ]
> #Echantillon de test
> test <- discrim.vm[-index, ]
> vm.disc= f.lda(appren[,-29],groups=appren$cluster)
> # Matrice de confusion
> pred.vm <-predict(vm.disc,newdata=test)
> # Taux  d'erreur
> Tx_err <- function(y,ypred){
+ mc <- table(y,ypred)
+ error <- 100*(mc[1,2]+mc[2,1])/sum(mc)
+ print(mc)
+ print(paste(round(error,2),"%",sep =""))
+ }
> Tx_err(test$cluster,pred.vm$class)
   ypred
y     1   2   3
  1  47   4   1
  2   3 142   2
  3   0  11  64
[1] "2.55%"
 
Une erreur de prévision de 2% avec une analyse discriminante
 
4. Avec une régression logistique multinomiale
 
# Avec une régression logistique multinomiale:
> # Avec une régression logistique multinomiale:
> library(nnet)
> logi=multinom(cluster~., data = appren)
# weights:  90 (58 variable)
initial  value 702.013252 
iter  10 value 113.648468
iter  20 value 75.506364
iter  30 value 69.704088
iter  40 value 65.256664
iter  50 value 62.358589
iter  60 value 61.152375
iter  70 value 60.246843
iter  80 value 59.528363
iter  90 value 58.957304
iter 100 value 58.452446
final  value 58.452446 
stopped after 100 iterations
> pp = predict(logi,newdata=test)
> Tx_err(test$cluster,pp)
   ypred
y     1   2   3
  1  51   0   1
  2   1 140   6
  3   1  12  62
[1] "0.36%"
 
 
Verdict : La reg. log. fait beaucoup mieux que l'analyse discriminante.
 
Certains auteurs  avaient déjà commencé à traiter la question. Ils voyaient l'analyse discriminante linéaire comme un cas particulier de la régression logistique (ce avec quoi je ne suis pas totalement d'accord, car aucune de ces deux méthodes ne s'affranchit de l'hypothèse forte de normalité et utilise les moments d'ordre deux à chaque fois)
Quelques conclusions néanmoins  intéressantes :

1. Lorsque le nombre de paramètre à estimer est important, le temps que met une reg log peut être 1,5 fois plus important que l'ADL du fait de l'algo itératif
2. La rég. log. peut être plus précise sur de petits échantillons, car du fait des modalités de référence, le nombre de paramètres à estimer est plus faible que si l'on utilise une ADL
 

GT
 
 

ggmap : Interesting toolbox for spatial analysis

ggmap is a new tool which enables such visualization by combining the spatial information of static maps from Google Maps, OpenStreetMap, Stamen Maps or CloudMade Maps with the layered grammar of graphics implementation of ggplot2

The library is developped by David Kahle and Hadley Wickham and in the latest R/Journal (Volume 5/1, June 2013), there is a whitepaper, very interesting that should to be read.

Let's use it to see where we can always buy tobacco at night in Paris.

Data are from data.ratp.fr, and give all authorized dealers by ratp. There's long and lat in this file, but this geocoding use lambert. So to show how ggmap work, we can extract adress, run geocoding using api/google and then ouput thematic map.

> require(ggmap)
> setwd("C:\\Users\\guibertt\\Desktop\\Geocodage - A faire")
> fic<-"commerces.csv"
> commerces<-read.csv(fic,header=T,sep=";", dec=",")
> head(commerces[,c(1,3:6)])
  DEA_CODE    ADRESSE_LIVRAISON DEA_CODE_POSTAL_LIVRAISON DEA_COMMUNE_LIVRAISON INSEE_LIVRAISON
1   510001       11   R. Mozart                     92230         Gennevilliers           92036
2   510002     81   Bd Voltaire                     92600    Asnières-sur-Seine           92004
3   510003 60   Av. Jean Moulin                     92390 Villeneuve-la-Garenne           92078
4   510004     2   Av. Michelet                     93400            Saint-Ouen           93070
5   510006      31   R. d'Anjou                     92600    Asnières-sur-Seine           92004
6   510007 45   R. Jules Larose                     92230         Gennevilliers           92036

We use the api google (limited at 2500 requests by day/ip adress) to geocode our places
>ad<-as.vector(tabac$adresse2)
>system.time(gc <- geocode(ad,output='latlona',messaging=FALSE))
> head(gc)
       lon      lat                                        address
1 2.267043 48.85066         108 avenue mozart, 75116 paris, france
2 2.297897 48.84552 56 rue de la croix nivert, 75015 paris, france
3 2.276529 48.84276             20 rue cauchy, 75015 paris, france
4 2.293242 48.83890  157 rue de la convention, 75015 paris, france
5 2.285362 48.83406       37 boulevard victor, 75015 paris, franc
>tabacp<-cbind(tabac,gc)

But, don't forget that ggmap is just a ggplot and we can get this, for example using OpenStreetMap

So we can now plot our places
png("paris4.png", width=800,height=600)
map <- get_map(location = 'paris',zoom=13,maptype="roadmap",color="color",source="google")
mymap = ggmap(map, darken = c(.3,'white'))
mymap+
  stat_bin2d(
aes(x = lon, y = lat, colour = typecommerce, fill = typecommerce),
size = .5, bins = 30, alpha = 1/2,
data = tabacp)
dev.off()

And at the end, with wrap, to split
png("paris5.png", width=800,height=600)
map <- get_map(location = 'paris',zoom=13,maptype="roadmap",color="bw",source="google")
mymap = ggmap(map, darken = c(.8,'white'))
mymap +  stat_bin2d(
aes(x = lon, y = lat, colour = TCO_LIBELLE, fill = TCO_LIBELLE),
size = .5, bins = 30, alpha = 1/2,
data = tabacp)+facet_wrap(~ TCO_LIBELLE)
dev.off()

How to read quickly large dataset in R?

Here, or there, I read many techniques to import a large dataset in R.

The option read.table or read.csv doesn't work anyway because, as discusshere, R load in memory. And sometimes, when we try to load a big dataset, we got this message :

Warning messages: 

1: Reached total allocation of 8056Mb: see help(memory.size)

2: Reached total allocation of 8056Mb: see help(memory.size) 

Many techniques can be used to load a large dataset. I found some there, or there. But there is two techniques that I never think before. 

Suppose that we have a large dataset with 10 millions rows

Comparing the methods for loading in R. 

- Using read.table

read.csv() performs a lot of analysis of the data it is reading, to determine the data types. So we can help R, by reading the first rows, determine the data type of the columns, and then, read the big data and provide the type of each columns and/or squeeze some of them if it doesn't need for analysis anyway;

Example

First we try to read a big data file (10 millions rows)

> system.time(df <-read.table(file="bigdf.csv",sep =",",dec=".")) Timing stopped at: 160.85 0.75 161.97 

 I let this run for a long period but no answer.

With this new method, we load the first rows, determine the data type and then, run read.table with indications of datatype.

> system.time (ds <- read.table("bigdf.csv", nrows=100, dec=".",sep=",")) user system elapsed 0 0 0 > classes <-sapply(ds, class) > classes V1 V2 V3 V4 "integer" "factor" "factor" "factor"

system.time(ds<-read.table("bigdf.csv",dec=".",sep=","colClasses=classes))

user  system elapsed 

234     432    128

As we see, this technique is not very interesting. It's also longer.

- We can use the package sqldf.

> require(sqldf)

> f <- file("bigdf.csv")

> system.time(SQLf <- sqldf("select * from f", dbname = tempfile(),

+                           file.format = list(header = T, row.names = F)))

Le chargement a nécessité le package : tcltk

   user  system elapsed 

  53.64    4.17   58.20 

Less of 1 minute  to import 10 millions rows  of an object of

> print(object.size(SQLf), units="Mb")

267 Mb

-          We can aslo used package read.table

> require(data.table)

Le chargement a nécessité le package : data.table

data.table 1.8.8  For help type: help("data.table")

> system.time(DT <- fread("bigdf.csv"))

   user  system elapsed 

 133.11    0.56  133.93 

But DT is a data.table format and a bit of transformation is require for use the table as dataframe using ddply from plyr package.

So. The point is : the package Sqldf is very useful to read quickly a large dataset in R. 10 millions rows in Less of a minute.

How logistic regression work ?

Discussing with a non statistician colleague, it seems that the logistic regression is not intuitive; Some basics questions like :
 - Why don't use the linear model?
 - What's logistic function?
 - How can we compute by hand, step by step to listen what is dealing by the glm function?

This post aims to answer that questions and may be this helps.

Suppose that we have this data : http://www.info.univ-angers.fr/~gh/wstat/pg.dar

 ID TAILLE GROUPE
1 A01    130      0
2 A02    140      0
3 C01    162      0
4 C02    160      1
5 A03    136      0
6 C03    165      1

and we want to predic the group according to the height. The problematic can be the level of risk according to age, or the customer segment according the amounts of transaction, etc. 
Let's remind.
When we compute a linear model (let's assume just one predictor : simple linear model), we have : E(y) =Cste + a1x1. Linear regression like all regressions focuses on the conditional probability distribution of Y given X.

The first think generally do is to draw the groupe = f(taille), we got :

The idea of Generaliszed Model (logistic regression is a particular ) is to replace E(Y) by something else.
For our example, we are interested by the probability of a person to be in group 0 or 1.
So, Instead of E(y) =Cste + a1x1, we seek P(Groupe==1) = a0 +a1*Taille. But, to solve the roblem, which is exactly the same to the other hand, we have to transform left hand side using a bijection between the interval[0,1]. That means to seek a "link" function that can help us to work in R.

The most useful function in logistic regression is : logit(p) = log(p/1-p). But one can also use the inverse of normal distribution(probit), the log-log distribution, or poisson distribution.

The method used to perform logistic regression is the maximization of likelihod estimator (MLE)
Read this post.
We sum up :
- Suppose in a population from which we are sampling, each individual has the same probability p to be in groupe 1 or groupe 0
- The likelihood is the joint probability of the data L = Product(P ** {Gourpe = 1} *(1 - p)**{Groupe = 0})
** Means power
For instance, we use log-likelihood.



How to interpret the likelihood :?

When we try to assign the group for a new id, it's natural to assign the group which have the best probability according to height.

Apply the MLE and perform logistic regression is done by

> Test = fit.logis(y=don$GROUPE,x=don$TAILLE)
> Test
  coef.est std.err
a  -27.190   8.885
b    0.181   0.058

 We can get the same output using glm function with "binomial" option.

>viaglm

Call:  glm(formula = don$GROUPE ~ don$TAILLE, family = "binomial", data = don)

Coefficients:
(Intercept)   don$TAILLE  
   -27.2103       0.1812  

Degrees of Freedom: 29 Total (i.e. Null);  28 Residual
Null Deviance:     38.19 
Residual Deviance: 10.89  AIC: 14.89
 

So, we can see that our optimisation via optim function is quite equivalent to glm function.
Just have a look

May be this helps to understand how it works !

Mining the last French presidential debate

After reading this post (thanks to him), I think it could be interesting to replicate this with some specific up of french language and to see and we can perform rapid view of the debate between Sarkozy and Hollande of the last 2nd round of presidential election.

Key words : TextMining, Elections, France, Debate, 2nd Round

We use the packages qdap from (Tyler Rinker) and tm to perform textmining analysis and the classical package like ggplot or RColorBrewer make our  graphics look pretty.
For Hollande

Top words From hollande

For Sarkozy

Top words from Sarkozy

As we can see, I've a problem to manage french accent. If somebody have any idea... We can also perform a quick Gantt plot basing on Qdap package and get some information about who lead the debate No surprise about the winner.