Article de blog intéressant sur l'ajustement de spline 2D avec JAGS. La fonction jagam du package mgcv semble très intéressante. Décortiquer cette fonction, programmée par S.Wood lui-même (!), risque d'apporter pas mal d'infos super intéressantes!

J'ai passé pas mal de temps à résoudre la question, avant de me rendre compte qu'il y avait une meilleure solution au pb. Cela dit, comme ce n'est pas la première fois que je me pose la question, je stocke la solution ici.

J'ai un modèle de Poisson surdispersé simpliste décrivant la variable X_i, une variable correspondant au nombre d'individus dans une harde de cerfs i:

X_i ~ Poisson(lambda_i)
log(lambda_i) = a_0 + epsilon_i
epsilon_i ~ Normal(0,sigma)

Je parcours un circuit et je compte: (i) le nombre N de hardes sur le circuit, et (ii) le nombre total Z de cerfs sur toutes les hardes du circuit. Quelle est la distribution du nombre total de cerfs?

Mettons que je détecte 5 hardes. En fait j'ai ici 5 variables X_i i.i.d. suivant toutes la distribution de Poisson surdispersée décrite ci-dessus. La somme de N lois de Poisson de paramètre lambda_1, lambda_2,...,lambda_N est une loi de Poisson de paramètre (lambda_1 + lambda_2 + ... + lambda_N). Dans le cas présent, le paramètre lambda de la loi de Poisson sera:

(exp(a_0+epsilon_1)+exp(a_0+epsilon_2)+...+exp(a_0+epsilon_N)) = (exp(a_0)*(exp(epsilon_1)+exp(epsilon_2)+...+exp(epsilon_N))

Nous déplaçons le problème: quelle est la distribution suivie par (exp(epsilon_1)+exp(epsilon_2)+...+exp(epsilon_N))?

En fait, si epsilon_i suit une loi normale, exp(epsilon_i) suit une loi lognormale. Alors la question qui se pose est: quelle est la distribution suivie par une somme de distribution log-normale i.i.d.?

Dans l'article suivant:

@article{Fenton1960,
  title={The sum of log-normal probability distributions in scatter transmission systems},
  author={Fenton, Lawrence F},
  journal={Communications Systems, IRE Transactions on},
  volume={8},
  pages={57--67},
  year={1960}
}

L'auteur démontre que la somme de N lois log-normale peut être approchée, à la louche, par une log-normale, dont les paramètres mtilde et sigtilde (moyenne et écart-type de la log-variable) peuvent être obtenus facilement (en faisant en sorte que la moyenne et la variance de la somme des variables corresponde à la moyenne et la variance de la log-normale approchant cette somme).

J'ai essayé l'approche sous R. Je stocke le code ici (les paramètres de la fonction ont le même nom que dans la description plus haut):

## la fonction ci-dessous permet de calculer les paramètres m et
## sigma de la lognormale de synthèse à partir de l'écart-type
## des résidus de surdispersion epsilon_i
f <- function(N, sigma)
{
    sigtilde <- sqrt(log((1/N)*(exp(sigma^2)-1) + 1))
    mtilde <- log((N*exp((sigma^2)/2)/(exp((sigtilde^2)/2))))
    return(c(mtilde,sigtilde))
}

## Simulation sous R: pour un écart-type des résidus epsilon_i = 0.5
## et pour un groupe de 10 individus, on a une distribution
## quasi-identique

N <- 10
sigma <-0.5
aa <- sapply(1:1000, function(i) sum(exp(rnorm(N, mean=0, sd=sigma))))
bb <- exp(rnorm(1000, mean=f(N,sigma)[1], sd=f(N,sigma)[2]))
qqplot(aa,bb)
abline(0,1)

Ça marche pas mal, à condition de pas avoir un sigma trop grand. On peut essayer de jouer avec ce sigma et ce N pour tester l'approximation.

Remarque: Au final, j'ai préféré utiliser des résidus nu_i de surdispersion suivant une loi gamma, i.e. un modèle de la forme:

X_i ~ Poisson(lambda_i)
lambda_i = mu * nu_i
nu_i ~ Gamma(alpha,alpha)

Dans ce cas de figure, les résidus sont de moyenne égale à 1, et X_i suit une négative binomiale. L'avantage, c'est que la somme de variables suivant chacune une loi binomiale négative suit également une loi binomiale négative, et pour le coup, le résultat n'est pas une approximation. Mais bon, l'approximation de Fenton est quand même intéressante...

Toujours bon à rappeler: le DIC n'est pas idéal avec les occupancy models, du fait de la structure hiérarchique de tels modèles et du caractère pas clair du niveau "focus" auquel tirer les inférences.

A lire. M'a l'air intéressant

Des conseils pour se former aux statistiques pour ceux qui ne connaissent pas. Ya des refs intéressantes, et les commentaires sont plein de bons conseils.

La vache, ya plein d'infos, de webinaires, de formations sur ce site. Dès que j'ai un peu de temps, faut que je creuse ça

Un joli modèle de Weibull pour prédire l'espérance de vie du pape. Résultat: une chance sur deux de durer au moins 10 ans.
Bon après, le modèle est uniquement basé sur les durées de reignes des précédents papes, et pas du tout sur l'état de santé du bonhomme en particulier (qui n'est pas bon, si j'ai bien compris).
Pas grave, sur un plan statistique, le modèle est un exemple intéressant.

checklist et solutions toujours utiles

Intéressant, à creuser. via Guillaume.

Ben décidément, yen a des choses sur Nimble aujourd'hui