Interlude de programmation — Emilie Champagne

Avis aux lecteurs n’utilisant jamais les statistiques : ce texte n’est pas pour vous. J’y décris un problème auquel je me suis heurtée et la solution à laquelle je suis arrivée. Ça peut avoir l’air plate, mais peut-être que ça aidera quelqu’un, un jour.

J’analyse présentement des données qui ont une distribution binomiale négative. La distribution binomiale négative et la distribution de poisson, contrairement à la distribution normale, ne sont pas symétriques. En fait, elles s’étirent comme une queue de poisson.

La distribution négative binomiale, sous forme de GIF (Merci Wikipedia https://en.wikipedia.org/wiki/Negative_binomial_distribution)

Le problème : l’analyse statistique que j’utilise (modèle linéaire généralisé mixte) a comme supposition que nos données se présentent sous la forme d’une distribution normale.

La solution : certaines fonctions des logiciels statistiques ont des fonctions particulières pour tenir compte de ces distributions non normales. Dans le logiciel SAS, c’est la fonction GLIMMIX. Dans R, on a l’option glmer.nb.

Le nouveau problème : mes modèles ne convergent pas (1), autant avec SAS qu’avec le logiciel R.

Nouvelle solution : la fonction glmmadmb (package glmmADMB, R) permet aussi de faire des modèles avec distribution binomiale négative. Et ça marche!

Le nouveau nouveau problème : Les autres fonctions que j’utilise, comme une me permettant de faire de la sélection de modèle, ne sont pas compatibles avec glmmadmb. Aussi, c’est dur à dire glmmadmb.

Nouvelle nouvelle solution (comme un modèle, les problèmes, c'est itératif) : j’ai écrit mon propre code pour faire un tableau de sélection de modèles. Je me suis fortement inspirée du package AICcmodavg et du cours que j'ai suivi avec son auteur, Marc Mazerolle (2). Probablement pas la fonction la plus élégante, mais la voici :

#Tout d'abord, créer une liste des modèles (Modlist) pour la sélection ainsi qu'une liste de noms de modèles (Modnames)

table.sel <- function(Modlist, Modnames){
    Out <- data.frame(cbind(Modnames, unlist(lapply(Modlist, "[[", "npar")), unlist(lapply(Modlist, FUN= AIC)),
            unlist(lapply(Modlist, "[[", "loglik"))), stringsAsFactors = FALSE)
    colnames(Out) [1:4] <- c("Names", "K", "AIC", "LogL")
    order <- order(Out$AIC)
    Out <- Out[order, ]
    Out$AIC <- as.numeric(Out$AIC)
    Out$LogL <- as.numeric(Out$LogL)
    Out$K <- as.numeric(Out$K)
    Out$Delta_AIC <- abs(ave(Out$AIC,FUN=function(x)x[1]-x))
    Out$ModelLik<-exp(-0.5*Out$Delta_AIC)
    Out$AICWt<-Out$ModelLik/sum(Out$ModelLik)
    print(Out)
}

tab <- table.sel(Modlist= Modlist, Modnames= Modnames)

Il est possible assez aisément de transformer le code pour sélectionner des AICc

Ajout (29-06-2015): Finalement, j'ai du changer de types d'analyses que je faisais et j'ai donc tout fait ça...pour rien ! Mais peut-être que mon expérience vous inspirera/vous sera utile !

(1) Si vous ne connaissez pas les statistiques et que vous lisez toujours, les modèles linéaires généralisés sont une équation. Le logiciel statistique essaie, avec plusieurs essais, de résoudre l’équation. Si après plusieurs itérations il n’y réussit pas, on dit que le modèle ne converge pas.

(2) Plutôt que d'utiliser mon code, utilisez ses excellentes fonctions. À moins d'être pris avec une distribution binomiale négative.