Freakonometrics

A la demande d'un étudiant en actuariat qui "cherche désespérément des codes tout faits sur internet, mais en vain" - que je ne vais pas satisfaire, faut pas rêver non plus* - je vais malgré tout mettre un petit billet sur l'utilisation des modèles GAM. Effectivement, dans quelques unes des formations R que j'ai pu organisé, on peut utiliser (rapidement) ces méthodes en provisionnement, au lieu d'introduire des facteurs de déroulé ou de date de survenance (de la même manière que sur les tables de mortalités prospectives ou des facteurs année et âge sont introduits).

  Les modèles GAM sont un outils fabuleux, mélant GLM (modèles linéaires généralisés), incluant les régressions binomiaux (probit ou logit) ou poissonniens, et régressions non-linéaires (à l'aide de fonctions splines). Ces modèles ont été introduits par Trevor Hastie et Rob Tibshirani en 1990 (sous la forme du livre chez Chapman). L'idée est de garder une forme linéaire, sauf que l'on s'autorise à avoir des transformations nonlinéaires de variables explicatives,

Contrairement aux approches classiques en régression nonlinéaire, on garde une forme additive, ce qui simplifie très fortement l'estimation. Au lieu d'avoir à estimer une fonction en dimension m, il est plus simple d'estimer m fonctions (univariées). Le code sous R est une simple extension de la fonction glm() qui est elle-même une extension de la fonction lm(). Mais autant ces fonctions étaient dans le package de base, autant il convient ici de charger la library(gam). La syntaxe est alors très simple,
> gam(formula, family = gaussian, data, weights, subset, na.action,  start, etastart, mustart, control = gam.control(...),...)
On retrouve la même syntaxe que la fonction 
> glm(formula, family = gaussian, data, weights, subset,    na.action, start = NULL, etastart, mustart,    offset, control = glm.control(...),  ...)
En fait, pour utiliser l'approche nonparamétrique, on utilise la syntaxe suivante (implémentée ici sur la base de données sur le prestige des métiers au Canada).
 > prestige<-read.table("http://perso.univ-rennes1.fr/arthur.charpentier/Prestige.txt",header=TRUE)
> head(prestige)
                    education income women prestige census type
GOV.ADMINISTRATORS      13.11  12351 11.16     68.8   1113 prof
GENERAL.MANAGERS        12.26  25879  4.02     69.1   1130 prof
ACCOUNTANTS             12.77   9271 15.70     63.4   1171 prof
PURCHASING.OFFICERS     11.42   8865  9.11     56.8   1175 prof
CHEMISTS                14.62   8403 11.68     73.5   2111 prof
PHYSICISTS              15.64  11030  5.13     77.6   2113 prof

> REG=gam(prestige~s(education)+s(income),data=prestige,
+ family=gaussian(link="identity"))
> plot(REG)

La dernière commande permet de visualiser la fonction des transformations non-linéaire utilisées. Malheureusement les deux graphiques sont imprimés à la suite, et comme la fenêtre graphique n'a pas d'historique, il est impossible de visualiser les deux. Une astuce possible est de générer un fichier pdf (qui contiendra ici deux pages),
 > pdf("graph-gam.pdf")
> plot(REG)
> dev.off()

Pour mémoire, pour retrouver l'endroit où est sauvé le fichier image, il suffit de demander où se trouve la working directory
> getwd()
[1] "H:/MES DOCUMENTS"

Sinon, on peut aussi tracer l'intervalle de confiance associée à l'estimation, tout simplement avec
 > pdf("graph-gam-se.pdf")
> plot(REG,se=TRUE)
> dev.off()
Voilà... ensuite pour faire du provisionnement avec des modèles GAM, on peut utiliser une régression Poissonnienne (au lieu du cas Gaussien utilisé ici), et mettre comme variable explicative la date de survenance et le delai (en tant que variables continues et non pas en tant que facteurs, comme pour les GLM).

Et sinon, on peut aussi utiliser la fonction predict() afin de représenter la surface de prédiciton.
 > Y=predict(REG,newdata=data.frame(income=rep(seq(0,25000,by=1000),26),education=rep(seq(6,16,by=10/25),each=26)))
> ZM=matrix(Y,26,26)
> XM=seq(0,25000,by=1000)
> YM=seq(6,16,by=10/25)
> persp(XM,YM,ZM,col="green",shade=TRUE,xlab="revenu",ylab="education",zlab="prestige")
>
> Y=predict(REG,newdata=data.frame(income=rep(seq(0,25000,by=100),251),education=rep(seq(6,16,by=10/250),each=251)))
> ZM=matrix(Y,251,251)
> XM=seq(0,25000,by=100)
> YM=seq(6,16,by=10/250)
> image(XM,YM,ZM,xlab="revenu",ylab="education",zlab="prestige",col=heat.colors(100))
> contour(XM,YM,ZM,add=TRUE) pdf("graph-gam-se.pdf")

* à l'attention des étudiants qui me contactent massivement: inutile de me demander de vous taper vos codes je ne le ferai pas... j'essaye toujours d'aider les étudiants qui me contactent (dans la limite de mon temps disponible), mais j'aide seulement, je ne tape pas à la place (ou alors ça s'appelle du conseil, et ça devient payant).