5.5 Applications en R

5.5.1 Algorithme d’estimation

Nous avons principalement couvert deux manières d’estimer les paramètres par maximum de vraisemblance:

  1. Utilisation d’un algorithme de type Newton-Raphson, nécessitant de trouver la condition de premier ordre, et une estimation de la matrice hessienne (le produit-croisé nous évite la double dérivée);

  2. Utilisation des fonctions d’un package R, comme glm(), gamlss() ou glm.nb(), qui ne nécessite que de définir le score du modèle. Les paramètres sont par la suite estimés automatiquement et une série de valeurs statistiques sont générées.

Entre l’approche i) qui demande pratiquement de tout faire, et l’approche ii) qui ne demande rien à l’utilisateur et qui fait presque tout, il existe une alternative à considérer.

Rappelons que l’estimateur par maximum de vraisemblance, noté MLE dans ces notes, est la valeur du paramètre qui maximise la fonction de vraisemblance, ou de logvraisemblance. Grâce à un optimisateur numérique, disponible dans la plupart des logiciels informatiques qui travaillent avec des données, il est possible de trouver les paramètres maximisant une fonction.

On doit tout d’abord utiliser le package optimx.

library(optimx)

Par la suite, nous devons simplement créer une fonction qui calcule la logvraisemblance d’un modèle en fonction de paramètres qu’on donne en entrée.

## Definition du score ##
score.glm <- as.formula(nb.sin ~ couleur + myopie + type_territoire + langue + alimentation)

# Generer la matrice des X, Y et expositions
Matrix.X <- model.matrix(score.glm, db.train)
Vect.Y <- as.matrix(db.train$nb.sin)
Vect.expo <- as.matrix(db.train$expo)

## logvraisemblance Poisson
max.loglike.Poisson <- function(parm){
  lambda <- Vect.expo*exp(Matrix.X %*% parm)
  return(-sum(Vect.Y*log(lambda) - lambda))
}

## initiation de beta ##
parm0 <- rep(0, ncol(Matrix.X)) 
parm0[1] <- log(sum(db.train$nb.sin)/sum(db.train$expo))

#optimisation 
optim.Poisson <- optimx(par = parm0, fn = max.loglike.Poisson, 
                        method = c("BFGS", "CG","Nelder-Mead"))

### Compare avec la fonction GLM ###
score.glm <- as.formula(nb.sin ~ couleur + myopie + type_territoire + langue + alimentation + offset(log(expo)))
Poisson1  <- glm(score.glm, family=poisson(link=log), data=db.train)

## Results
a <- cbind(logLik(Poisson1), t(-optim.Poisson$value - sum(lgamma(Vect.Y+1))))
b <- cbind(Poisson1$coefficients, t(coef(optim.Poisson)))
table <- rbind(a,b)
rownames(table)[1] <- 'Loglikelihood'

knitr::kable(table, digits = c(4, 4, 4, 4))
BFGS CG Nelder-Mead
Loglikelihood -142286.3267 -142286.3267 -142286.3268 -142308.9790
(Intercept) -1.2121 -1.2121 -1.2124 -1.2960
couleurRouge -0.0210 -0.0210 -0.0209 0.0577
myopieOui 0.0759 0.0759 0.0759 0.0805
type_territoireSemi-urbain -0.0955 -0.0955 -0.0955 -0.1047
type_territoireUrbain -0.1486 -0.1486 -0.1485 -0.1482
langueF 0.2503 0.2503 0.2503 0.2262
alimentationVégétalien -0.2431 -0.2431 -0.2429 -0.2450
alimentationVégétarien -0.4280 -0.4280 -0.4279 -0.4107

L’avantage d’utiliser optimx() est que nous n’avons plus à calculer la condition de premier ordre. Nous n’avons qu’à indiquer la fonction de logvraisemblance.

Par contre, il faut faire attention parce que le temps d’exécution de cette routine est beaucoup plus long que ce que nous avons vu depuis le début des notes de cours.

La fonction optimx() offre la possibilité d’utiliser diverses méthodes de maximisation. Dans l’exemple plus haut, nous avons utilisé “BFGS”, “CG” et “Nelder-Mead”. On voit que la dernière méthode génère de moins bons résultats puisque le logvraisemblance obtenu n’est pas aussi bon que les deux autres.

Je ne sais pas si “Nelder-Mead” est toujours moins bon que les deux autres approches, ou si tout cela dépend de la base de données utilisée, et des fonctions à maximiser. D’autres method sont aussi disponible. Il suffit de regarder la document sur optimx(). Choisir au moins trois méthodes, et prendre les paramètres qui produisent le meilleur logvraisemblance peut être une approche prudent, mais cela augmente le temps de calcul.

Le package optimx peut aussi être utilisé pour calculer la variance (ou l’écart-type) des estimateurs. Évidemment, le calcul des écart-types ne sera pas aussi précis que ce que nous avions fait plus tôt, car la fonction gHgen() calcul la matrice hessienne numériquement.

parm <- coef(optim.Poisson)[1,]
hessian <- gHgen(parm, max.loglike.Poisson)
std.err <- sqrt(diag(solve(hessian$Hn)))
table <- cbind(std.err, summary(Poisson1)$coefficients[, 2])
knitr::kable(table, digits = c(4, 4))
std.err
(Intercept) 0.0256 0.0256
couleurRouge 0.0189 0.0189
myopieOui 0.0135 0.0135
type_territoireSemi-urbain 0.0126 0.0126
type_territoireUrbain 0.0122 0.0122
langueF 0.0102 0.0102
alimentationVégétalien 0.0177 0.0177
alimentationVégétarien 0.0168 0.0168

On voit par contre que l’approximation numérique est très proche de ce que génère la fonction glm().

5.5.2 Autres modèles

En utilisant l’approche avec la fonction optimx(), et les fonctions calculant les logvraisemblances pour les modèles, il est ainsi possible d’estimer relativement facilement les paramètres et les écart-types de ces paramètres pour différents modèles.

En exercice, ou en devoir, vous pouvez ainsi estimer les modèles vus en début de chapitre.