Quelle analyse

Bonjour tout le monde !

J'ai des données météo journalières avec le nombre d'individus de 14 espèces d'amphibiens (ex. 'tripal' et 'amphib') récupérés dans un crapodrome. Le crapodrome c'est une route avec des barrières sur les cotés et des seaux tout les 30 mètres pour éviter les collisions routières. Les seaux sont relevés quotidiennement.
Il y a eut trois saisons de 2016 à 2019 entre Décembre et Avril.

Code:

> head(crapo)
      date tmin tmoy tmax ventmoy ventmax  ori preci txhum pression  sol tripal amphib
1 20161215  8.4  9.5 10.4     7.9    20.9 S-SE   0.0    56   1017.5  0.2      0      0
2 20161216  6.7  9.1 12.8     3.2    11.2 S-SE   0.2    55   1022.8 52.3      0      0
3 20161217  2.6  6.9  9.6     1.4     6.5 N-NE   0.5    55   1037.5  0.0     10     12
4 20161218  5.7  6.9  8.1     3.2     7.9   NE   0.3    54   1036.2  0.0      6      8
5 20161219  4.9  6.4  7.7     4.7    11.2 N-NE   2.3    54   1025.9  0.0      3      3
6 20161220  3.8  4.7  6.4     7.9    24.1    S   1.8    54   1020.6  0.0      0      0

> str(crapo)
'data.frame': 321 obs. of  13 variables:
 $ date    : int  20161215 20161216 20161217 20161218 20161219 20161220 20161221 20161222 20161223 20161224 ...
 $ tmin    : num  8.4 6.7 2.6 5.7 4.9 3.8 6.3 5.7 3.3 5.4 ...
 $ tmoy    : num  9.5 9.1 6.9 6.9 6.4 4.7 9.6 8.8 8.2 8.2 ...
 $ tmax    : num  10.4 12.8 9.6 8.1 7.7 6.4 12.2 12.6 12 11 ...
 $ ventmoy : num  7.9 3.2 1.4 3.2 4.7 7.9 11.2 1.4 11.2 7.9 ...
 $ ventmax : num  20.9 11.2 6.5 7.9 11.2 24.1 28.8 9.7 27.4 16.2 ...
 $ ori     : Factor w/ 15 levels "E","E-NE","E-SE",..: 12 12 5 7 5 11 13 8 15 10 ...
 $ preci   : num  0 0.2 0.5 0.3 2.3 1.8 11.9 0.5 0.8 0.5 ...
 $ txhum   : int  56 55 55 54 54 54 54 56 55 55 ...
 $ pression: num  1018 1023 1038 1036 1026 ...
 $ sol     : num  0.2 52.3 0 0 0 0 0 41.4 2.4 20.8 ...
 $ tripal  : int  0 0 10 6 3 0 0 1 2 8 ...
 $ amphib  : int  0 0 12 8 3 0 2 19 2 14 ...

> summary(crapo)
      date               tmin            tmoy             tmax          ventmoy        ventmax           ori    
 Min.   :20161215   Min.   :-7.50   Min.   :-3.800   Min.   :-0.10   Min.   : 1.4   Min.   : 4.70   O-SO   :79  
 1st Qu.:20170305   1st Qu.: 2.20   1st Qu.: 5.400   1st Qu.: 8.80   1st Qu.: 7.9   1st Qu.:17.70   SO     :78  
 Median :20180206   Median : 4.20   Median : 7.700   Median :10.80   Median :12.2   Median :25.70   S-SE   :36  
 Mean   :20178787   Mean   : 4.13   Mean   : 7.271   Mean   :10.91   Mean   :12.8   Mean   :26.76   N-NE   :29  
 3rd Qu.:20190110   3rd Qu.: 6.50   3rd Qu.: 9.400   3rd Qu.:13.10   3rd Qu.:17.6   3rd Qu.:33.80   NE     :29  
 Max.   :20190331   Max.   :12.50   Max.   :15.300   Max.   :25.00   Max.   :29.8   Max.   :58.00   S-SO   :18  
                                                                                                    (Other):52  
     preci            txhum          pression         sol            tripal            amphib      
 Min.   : 0.000   Min.   : 7.00   Min.   : 983   Min.   : 0.00   Min.   :  0.000   Min.   :  0.00  
 1st Qu.: 0.000   1st Qu.:51.00   1st Qu.:1009   1st Qu.: 2.30   1st Qu.:  0.000   1st Qu.:  2.00  
 Median : 0.500   Median :78.00   Median :1018   Median :19.80   Median :  3.000   Median :  6.00  
 Mean   : 2.138   Mean   :69.02   Mean   :1017   Mean   :27.83   Mean   :  8.592   Mean   : 15.22  
 3rd Qu.: 2.800   3rd Qu.:88.00   3rd Qu.:1027   3rd Qu.:48.30   3rd Qu.: 10.000   3rd Qu.: 18.00  
 Max.   :25.900   Max.   :99.00   Max.   :1042   Max.   :85.70   Max.   :132.000   Max.   :198.00  

La question : peut-on prédire le déplacements des amphibiens avec les données météo ?
L'analyse : je pense à un GLM
Mon problème : les comptages ne suivent pas une loi de poisson


Distribution théorique (rose) et observée (bleue)


Y a-t-il quelque chose à faire ? d'autres analyses ?



Merci d'avance

Avant tout chose, d'où provient la distribution observée que vous nous montrez ?

Il s'agit de regarder la distribution des données dans chaque modalité des variables explicatives (non poolé). Est-ce que vous faites ?

Si tel est le cas, vous avez apparemment un wagon de zéros, et c'est probablement pour ceci que vos données ne sont pas conformes à des lois de Poisson. Vous pourriez jeter un oeil sur les méthodes dites "zero inflated Poisson" qui combinent intelligemment un processus de Poisson (pour les comptages) avec un loi binomiale (pour la proportion de zéros). Je pense que ça vaudrait le cout d'essayer.

HTH, Eric.

Salut Eric,

Merci pour ta réponse !

Eric Wajnberg a écrit:Avant tout chose, d'où provient la distribution observée que vous nous montrez ?

Il s'agit de la distribution observée des tritons palmés (tripal)

Eric Wajnberg a écrit:Il s'agit de regarder la distribution des données dans chaque modalité des variables explicatives (non poolé). Est-ce que vous faites ?

Mais mes variables explicatives sont quantitatives ...

Eric Wajnberg a écrit: Si tel est le cas, vous avez apparemment un wagon de zéros, et c'est probablement pour ceci que vos données ne sont pas conformes à des lois de Poisson.

Code:

> table(crapo$tripal)

  0  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  24  25  26  27  28  30  31  35  36  37  39  44
 85  40  28  22  17  10  11  12  11  3  11  4  3  7  5  3  1  4  4  3  1  1  2  2  2  1  2  1  2  1  1  1  2  2  2
 45  46  48  49  52  57  62  67  72  74 110 132
  1  1  2  2  1  1  1  1  1  1  1  1

Vous pensez ?! Parce que si j'enlève les zéros -> la moyenne de tripal augmente -> le paramètre lambda augmente -> le décalage augmente.
Je pensais que c'était surtout lié aux valeurs extrêmes mais en supprimant les deux plus grandes cela ne change pas beaucoup lambda. Et puis ca devient un peu de la cuisine ...

Je n'ai pas dit qu'il fallait enlever les zéros. Relisez ma réponse.

Pour le problème de la comparaison de la distribution des loi observées à une loi de Poisson, si vous laissez l'influence de vos variables explicatives, vous ne pouvez évidement pas retomber sur une loi de Poisson, simplement parce que vos variables explicatives modifient les valeurs observées. Le graphique de distribution que vous fournissez ici n'a guère de sens, ni d'utilité, donc.

Il n'y pas pas trop de solutions, et - d'une manière générale - il est toujours difficile de savoir quelles distributions suivent vos données. Bref, si vous avez un comptage, alors vous pouvez penser que vous êtes dans le cas poissonien, et une régression log-linéaire (GLM) doit faire l'affaire. Vous avez toujours la ressource de prendre l'option "family=quasipoisson" et de jeter un oeil sur le paramètre de dispersion que l'ajustement va estimer, et qui vous indiquera - a posteriori - si vous être bien dans un cas poissonien ou pas.

HTH, Eric.

Bonjour Eric,

Eric Wajnberg a écrit:Pour le problème de la comparaison de la distribution des loi observées à une loi de Poisson, si vous laissez l'influence de vos variables explicatives, vous ne pouvez évidement pas retomber sur une loi de Poisson, simplement parce que vos variables explicatives modifient les valeurs observées.

Il ne me semble pas avoir laisser l'influence des variables explicatives.

Code:

set.seed(1234) # permet de simuler toujours les mêmes comptages.
theoretic_count <-rbin(321,3)

# on incorpore ces comptages théoriques dans un data frame
tc_df <-data.frame(theoretic_count)
# on plot simultanémaent les comptages observés et les comptages théoriques
library(ggplot2)
ggplot(crapo,aes([color=#ff0000]tripal[/color]))+
  geom_bar(fill="#1E90FF")+
  geom_bar(data=tc_df, aes(theoretic_count,fill="#1E90FF", alpha=0.5))+
  theme_classic()+
  theme(legend.position="none")

Pour la ZIP j'ai un souci ...

Code:

crapo = structure(list(tmoy = c(9.5, 9.1, 6.9, 6.9, 6.4, 4.7, 9.6, 8.8,
8.2, 8.2, 10.1, 9.4, 4.1, 2.5, 3.1, 1.8, -1.5, 4.1, 4.2, 2.8,
4.8, 4.8, 1.4, 4.2, 8.9, 8.4, 8.6, 9.9, 6.8, 5.6), tmax = c(10.4,
12.8, 9.6, 8.1, 7.7, 6.4, 12.2, 12.6, 12, 11, 12.1, 12.4, 7.7,
5.6, 9, 5.1, 4.9, 6.8, 6.4, 8.8, 7.8, 9.2, 6.9, 8.2, 10.7, 10.1,
10.7, 11.5, 10, 8.3), ventmoy = c(7.9, 3.2, 1.4, 3.2, 4.7, 7.9,
11.2, 1.4, 11.2, 7.9, 9.7, 6.5, 4.7, 3.2, 3.2, 1.4, 1.4, 9.7,
9.7, 1.4, 4.7, 3.2, 3.2, 4.7, 3.2, 14.4, 6.5, 7.9, 27.4, 14.4
), ventmax = c(20.9, 11.2, 6.5, 7.9, 11.2, 24.1, 28.8, 9.7, 27.4,
16.2, 20.9, 17.6, 17.6, 11.2, 11.2, 4.7, 6.5, 19.4, 22.7, 7.9,
16.2, 7.9, 9.7, 11.2, 11.2, 28.8, 22.7, 28.8, 54.7, 38.5), ori = structure(c(12L,
12L, 5L, 7L, 5L, 11L, 13L, 8L, 15L, 10L, 15L, 4L, 1L, 12L, 12L,
11L, 13L, 15L, 7L, 13L, 10L, 5L, 11L, 15L, 10L, 10L, 6L, 10L,
15L, 6L), .Label = c("E", "E-NE", "E-SE", "N", "N-NE", "N-NO",
"NE", "NO", "O", "O-SO", "S", "S-SE", "S-SO", "SE", "SO"), class = "factor"),
    preci = c(0, 0.2, 0.5, 0.3, 2.3, 1.8, 11.9, 0.5, 0.8, 0.5,
    0, 1.3, 0, 0.2, 0.3, 0.2, 0.3, 0.3, 0, 0.5, 0, 0, 0.2, 0.3,
    1.5, 2.3, 1.3, 2.2, 7.1, 2.1), txhum = c(56L, 55L, 55L, 54L,
    54L, 54L, 54L, 56L, 55L, 55L, 55L, 55L, 53L, 51L, 50L, 50L,
    49L, 48L, 49L, 47L, 47L, 46L, 45L, 46L, 49L, 50L, 51L, 52L,
    52L, 48L), pression = c(1017.5, 1022.8, 1037.5, 1036.2, 1025.9,
    1020.6, 1022.8, 1029.4, 1031.5, 1032.9, 1030.6, 1034.8, 1041.9,
    1038.7, 1036.1, 1036.6, 1032.6, 1020.6, 1025, 1030.8, 1028.7,
    1032.8, 1037.2, 1035.8, 1033.8, 1024.8, 1015.9, 1018, 1011,
    1016.2), sol = c(0.2, 52.3, 0, 0, 0, 0, 0, 41.4, 2.4, 20.8,
    0, 43.9, 40.5, 16.4, 65.4, 0, 59.7, 12.9, 0.6, 64.5, 0, 22.7,
    57.5, 0, 3.5, 10.5, 15.9, 3.7, 1.2, 36.5), tripal = c(0L,
    0L, 10L, 6L, 3L, 0L, 0L, 1L, 2L, 8L, 15L, 18L, 0L, 0L, 0L,
    2L, 3L, 0L, 5L, 2L, 0L, 10L, 2L, 1L, 3L, 17L, 5L, 18L, 10L,
    2L)), class = "data.frame", row.names = c(NA, 30L))

> mzip <- zeroinfl(tripal ~ tmin + tmoy + tmax + ventmoy + ventmax + ori + preci + txhum + pression + sol | tmin + tmoy + tmax + ventmoy + ventmax + ori + preci + txhum + pression + sol ,data=crapo)
[color=#ff0000]Warning message:
In value[[3L]](cond) :
  routine Lapack dgesv : le système est exactement singulier : U[32,32] = 0FALSE[/color]

> print(summary(mzip))

Call:
zeroinfl(formula = tripal ~ tmin + tmoy + tmax + ventmoy + ventmax + ori + preci + txhum + pression + sol | tmin + tmoy +
    tmax + ventmoy + ventmax + ori + preci + txhum + pression + sol, data = crapo)

Pearson residuals:
    Min      1Q  Median      3Q     Max
-3.1325 -1.1457 -0.5894  0.5694 13.2141

Count model coefficients (poisson with log link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.2514268         NA      NA       NA
tmin         0.0004383         NA      NA       NA
tmoy         0.3944597         NA      NA       NA
tmax        -0.2172592         NA      NA       NA
ventmoy     -0.0078335         NA      NA       NA
ventmax      0.0081198         NA      NA       NA
oriE-NE      1.0852193         NA      NA       NA
oriE-SE      2.1518525         NA      NA       NA
oriN         1.7972461         NA      NA       NA
oriN-NE      1.5591143         NA      NA       NA
oriN-NO      0.7622988         NA      NA       NA
oriNE        1.6045508         NA      NA       NA
oriNO        1.5663493         NA      NA       NA
oriO         0.0828207         NA      NA       NA
oriO-SO      1.6066215         NA      NA       NA
oriS         2.0638507         NA      NA       NA
oriS-SE      1.4392252         NA      NA       NA
oriS-SO      2.5209954         NA      NA       NA
oriSE        2.8224036         NA      NA       NA
oriSO        1.9302062         NA      NA       NA
preci        0.0476095         NA      NA       NA
txhum        0.0068328         NA      NA       NA
pression     0.0041875         NA      NA       NA
sol          0.0029493         NA      NA       NA

Zero-inflation model coefficients (binomial with logit link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  20.086859         NA      NA       NA
tmin         -0.245618         NA      NA       NA
tmoy         -0.061699         NA      NA       NA
tmax          0.087050         NA      NA       NA
ventmoy      -0.112970         NA      NA       NA
ventmax       0.062839         NA      NA       NA
oriE-NE      14.605326         NA      NA       NA
oriE-SE     -17.076171         NA      NA       NA
oriN         14.956201         NA      NA       NA
oriN-NE      15.267773         NA      NA       NA
oriN-NO      12.556253         NA      NA       NA
oriNE        14.644322         NA      NA       NA
oriNO        14.507729         NA      NA       NA
oriO        -13.385202         NA      NA       NA
oriO-SO      14.258431         NA      NA       NA
oriS         14.177546         NA      NA       NA
oriS-SE      15.283597         NA      NA       NA
oriS-SO      14.680558         NA      NA       NA
oriSE       -15.700968         NA      NA       NA
oriSO        14.324057         NA      NA       NA
preci        -0.037658         NA      NA       NA
txhum        -0.010720         NA      NA       NA
pression     -0.034038         NA      NA       NA
sol          -0.008373         NA      NA       NA

Number of iterations in BFGS optimization: 77
Log-likelihood: -1622 on 48 Df

letolah a écrit:Il ne me semble pas avoir laisser l'influence des variables explicatives.

Lorsque vous faites un histogramme de la distribution des données de "tripal", vous prenez toutes les données globalement, mais ces données sont sous l'influence des autres variables (que vous voulez tester). Prenons un exemple plus simple. Imaginons que l'on veuille comparer la taille des individus mâles et femelles dans une population, et que l'on utilise un test-t sur données supposément suivre une loi normale. On doit alors regarder la distribution des tailles séparément chez les mâles et chez les femelles. Faire une histogramme globalement des tailles mélangerait deux distributions, et n'apporterait aucune information pertinente.

Pour votre régression zero-inflated, vous avez beaucoup de régresseurs. Je commencerais par faire un essai avec moins de régresseurs (un ou deux) pour voir si le problème ne vient pas de là..

HTH, Eric.