THEME : MODELE DE REGRESSION A INFLATION DE ZERO

Nous utilisons le jeu de données “bioChemists” du package R pscl comprenant des informations sur 915 étudiants diplômés en biochimie. Chaque étudiant est caractérisé par le nombre d’articles scientifiques produits au cours des trois dernières années de leur doctorat (“art”), leur genre (“fem” pour hommes ou femmes), leur statut marital (“mar” pour célibataire ou marié), le nombre d’enfants de moins de cinq ans à charge (“kid5”), le prestige du département de leur doctorat (“phd”), et le nombre d’articles produits par leur mentor académique (“ment”).

Dans cette étude, notre objectif est de déterminer le modèle optimal qui explique le nombre d’articles publiés en fonction des différentes variables explicatives disponibles. Nous allons explorer et comparer plusieurs modèles pour évaluer comment des facteurs tels que le genre, le statut marital, le nombre d’enfants, le prestige du département de doctorat, et l’influence du mentor académique influencent la production d’articles scientifiques chez les biochimistes diplômés. L’analyse vise à identifier les variables les plus significatives et à construire un modèle statistique robuste qui peut prédire efficacement la productivité académique dans ce contexte spécifique.

#installation et chargement de packages necessaires
#install.packages("lmtest")
#install.packages("car")
#install.packages("glmnet")
#install.packages("tseries")

library(MASS)
library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
library("lmtest")
## Le chargement a nécessité le package : zoo
## 
## Attachement du package : 'zoo'
## Les objets suivants sont masqués depuis 'package:base':
## 
##     as.Date, as.Date.numeric
library("car")
## Le chargement a nécessité le package : carData
library(glmnet)
## Le chargement a nécessité le package : Matrix
## Loaded glmnet 4.1-8
library(ResourceSelection)
## ResourceSelection 0.3-6   2023-06-27
data(bioChemists)
df <- bioChemists

head(df)
##   art   fem     mar kid5  phd ment
## 1   0   Men Married    0 2.52    7
## 2   0 Women  Single    0 2.05    6
## 3   0 Women  Single    0 3.75    6
## 4   0   Men Married    1 1.18    3
## 5   0 Women  Single    0 3.75   26
## 6   0 Women Married    2 3.59    2
summary(df)
##       art            fem           mar           kid5             phd       
##  Min.   : 0.000   Men  :494   Single :309   Min.   :0.0000   Min.   :0.755  
##  1st Qu.: 0.000   Women:421   Married:606   1st Qu.:0.0000   1st Qu.:2.260  
##  Median : 1.000                             Median :0.0000   Median :3.150  
##  Mean   : 1.693                             Mean   :0.4951   Mean   :3.103  
##  3rd Qu.: 2.000                             3rd Qu.:1.0000   3rd Qu.:3.920  
##  Max.   :19.000                             Max.   :3.0000   Max.   :4.620  
##       ment       
##  Min.   : 0.000  
##  1st Qu.: 3.000  
##  Median : 6.000  
##  Mean   : 8.767  
##  3rd Qu.:12.000  
##  Max.   :77.000

Matrice scatter plot

# Transformation des variables catégorielles en numériques pour pairs()
df$fem <- as.numeric(df$fem)
df$mar <- as.numeric(df$mar)

# Création de la matrice de scatter plots
pairs(df[, c("art", "fem", "mar", "kid5", "phd", "ment")], 
      main = "Matrice de Scatter Plots des Variables")

En examinant ce graphique de dispersion, nous constatons que les covariables ‘phd’ (prestige du département de doctorat) et ‘ment’ (nombre d’articles publiés par le mentor de thèse au cours des trois dernières années) semblent être corrélées.

# Calcul de la matrice de corrélation pour les variables numériques
numeric_vars <- df[, c("kid5", "phd", "ment")]
cor(numeric_vars)
##             kid5         phd       ment
## kid5  1.00000000 -0.02650604 0.03702032
## phd  -0.02650604  1.00000000 0.26041413
## ment  0.03702032  0.26041413 1.00000000
# Calcul de la corrélation de Pearson entre les variables 'phd' et 'ment'
cor.test(bioChemists$phd, bioChemists$ment, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  bioChemists$phd and bioChemists$ment
## t = 8.1498, df = 913, p-value = 1.195e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1989622 0.3198262
## sample estimates:
##       cor 
## 0.2604141

cor(phd,ment) = 0.2604

Le test de corrélation de Pearson montre une corrélation positive (0.260) et hautement significative (p-value ≈ 1.195e-15) entre le prestige du département de doctorat et le nombre d’articles publiés par le mentor de thèse au cours des trois dernières années, ce qui est tout a fait normal. Dans la suite, nous pourrions envisager de ne prendre en compte qu’une seule des deux variables pour l’ajustement des modèles.

Distribution du nombre d’article pubilé

hist(df$art, breaks = 20, main = "Distribution du nombre d'article publié", xlab = "Nombre d'articles", ylab = "Frequency")

Nous notons d’apres cet histogramme une abondance de zéros dans les comptes.

Vérification de la dispertion

mean(df$art)
## [1] 1.692896
var(df$art)
## [1] 3.709742
  • La variance dépasse largement la moyenne( trois fois environ la moyenne),nous en déduisons une surdispersion dans les données. Un modèle Poisson pourrait ne pas être approprié.

Vérification de l’excès de zéro

  • Comparons le nombre de zéros observés avec le nombre attendu selon une distribution de Poisson.

    # Compte de zéros observé
zero_count <- sum(df$art == 0)
total_count <- nrow(df)
proportion_zeros_observed <- zero_count / total_count

# le nombre de zéros observé
table(df$art == 0)
    ## 
## FALSE  TRUE 
##   640   275
    # nombre d'articles publiés
mean_art <- mean(df$art)

# Probabilité de zéro dans une distribution Poisson
prob_zero_poisson <- dpois(0, lambda = mean_art)

proportion_zeros_observed
    ## [1] 0.3005464
    prob_zero_poisson
    ## [1] 0.1839859

Interpretation des résultats

30% de zéros observés contre environ 18% attendus selon la distribution de Poisson. La distribution de posson n’est de ce fait pas adapté pour modéliser ses données. la proportion de zéros observée est significativement plus élevée que la probabilité de zéro attendue dans une distribution Poisson, cela suggère un excès de zéros dans les données.

Dans ce cas, un modèle ZIP (Zero-Inflated Poisson) ou ZINB (Zero-Inflated Negative Binomial) serait plus approprié, car ces modèles prennent en compte la présence excessive de zéros en modélisant deux processus : l’un générant les zéros en excès et l’autre générant les comptes positifs selon une distribution Poisson ou Negative Binomial.

Diagnostic

Nous ajustons en premier lieux un modèle binomial négatif afin tenir compte de la surdispertion

nb_model <- glm.nb(art ~ fem + mar + kid5 + phd + ment, data = df)

# Résumé du modèle ( Test de Wald )
summary(nb_model)
## 
## Call:
## glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = df, 
##     init.theta = 2.264387695, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.322073   0.222090   1.450 0.147004    
## fem         -0.216418   0.072636  -2.979 0.002887 ** 
## mar          0.150489   0.082097   1.833 0.066791 .  
## kid5        -0.176415   0.052813  -3.340 0.000837 ***
## phd          0.015271   0.035873   0.426 0.670326    
## ment         0.029082   0.003214   9.048  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2644) family taken to be 1)
## 
##     Null deviance: 1109.0  on 914  degrees of freedom
## Residual deviance: 1004.3  on 909  degrees of freedom
## AIC: 3135.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.264 
##           Std. Err.:  0.271 
## 
##  2 x log-likelihood:  -3121.917

Interprétation de quelques coefficients

nb_model <- glm.nb(art ~ fem + mar + kid5 + ment, data = df)

# Résumé du modèle ( Test de Wald )
summary(nb_model)
## 
## Call:
## glm.nb(formula = art ~ fem + mar + kid5 + ment, data = df, init.theta = 2.26411693, 
##     link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.373057   0.186978   1.995 0.046022 *  
## fem         -0.216673   0.072624  -2.983 0.002850 ** 
## mar          0.146944   0.081765   1.797 0.072312 .  
## kid5        -0.176797   0.052826  -3.347 0.000818 ***
## ment         0.029430   0.003108   9.470  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2641) family taken to be 1)
## 
##     Null deviance: 1108.9  on 914  degrees of freedom
## Residual deviance: 1004.4  on 910  degrees of freedom
## AIC: 3134.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.264 
##           Std. Err.:  0.271 
## 
##  2 x log-likelihood:  -3122.096

Le AIC de ce nouveau modèle (3134.1) est inférieur à celui du modèle précédent, indiquant que le nouveau modèle s’ajuste mieux aux données (AIC précédent : 3135.9).

Nous comparerons ce modèle aux modèles ZIP et ZINB, que nous utiliserons pour diagnostiquer nos données présentant une surdispersion et une inflation de zéros.

Ajustement d’un modèle ZIP

zip_model <- zeroinfl(art ~ fem + mar + kid5 + phd + ment, data = df, dist = "poisson")
# Sélection de variables par élimination pas à pas basée sur l'AIC
step_model <- stepAIC(zip_model, direction = "both")
## Start:  AIC=3233.55
## art ~ fem + mar + kid5 + phd + ment
## 
##        Df    AIC
## - phd   2 3229.6
## <none>    3233.5
## - mar   2 3234.8
## - fem   2 3243.6
## - kid5  2 3245.3
## - ment  2 3340.4
## 
## Step:  AIC=3229.59
## art ~ fem + mar + kid5 + ment
## 
##        Df    AIC
## <none>    3229.6
## - mar   2 3231.0
## + phd   2 3233.5
## - fem   2 3239.6
## - kid5  2 3241.3
## - ment  2 3344.2

Le modèle ZIP ajusté avec l’élimination pas à pas basée sur le critère AIC montre que la variable phd (prestige du département de doctorat) peut être retirée sans augmenter l’AIC. Cela indique qu’elle n’améliore pas significativement le modèle. Ainsi, le modèle final retenu inclut les variables fem, mar, kid5et ment, qui sont jugées les plus pertinentes pour expliquer le nombre d’articles publiés par les biochimistes doctorants.

Modèle ZIP sélectionné

# Nous laissons tomber la variable phd 
zip_model <- zeroinfl(art ~ fem + mar + kid5 + ment, data = df, dist = "poisson")
#summary(zip_model)
# Nous intégrons uniquement la variable ment dans la composante logit vu quelle est #la seule étant significative de la composante logit du modèle avec toutes les variables
print("AIC avec toutes les var de la composante logit:")
## [1] "AIC avec toutes les var de la composante logit:"
AIC(zip_model)
## [1] 3229.594
zip_model <- zeroinfl(art ~ fem + mar + kid5 + ment|ment, data = df, dist = "poisson")
print("AIC sans les var non significatifs de la composante logit:")
## [1] "AIC sans les var non significatifs de la composante logit:"
AIC(zip_model)
## [1] 3225.517
# Résumé du nouveau modèle
summary(zip_model)
## 
## Call:
## zeroinfl(formula = art ~ fem + mar + kid5 + ment | ment, data = df, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3222 -0.8700 -0.2588  0.5475  7.1875 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.69291    0.15170   4.568 4.93e-06 ***
## fem         -0.21826    0.05878  -3.713 0.000205 ***
## mar          0.13483    0.06587   2.047 0.040670 *  
## kid5        -0.16277    0.04337  -3.753 0.000175 ***
## ment         0.01819    0.00221   8.227  < 2e-16 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.68569    0.20548  -3.337 0.000847 ***
## ment        -0.13007    0.04023  -3.233 0.001224 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 26 
## Log-likelihood: -1606 on 7 Df

Interpretation

Modèle de Comptage (Poisson avec lien log)

  • Intercept (Intercept) : En moyenne, lorsque toutes les autres variables sont nulles, environ 1.999 articles sont publiés.
  • fem (Genre - Women) : Les femmes publient en moyenne environ 80.3% moins d’articles que les hommes, toutes choses étant égales par ailleurs.
  • mar (Statut matrimonial - Married) : Les biochimistes mariés publient en moyenne environ 14.5% plus d’articles que ceux célibataires.
  • kid5 (Nombre d’enfants de moins de cinq ans) : Chaque enfant de moins de cinq ans supplémentaire est associé à une diminution moyenne de 16.5% dans le nombre d’articles publiés.
  • ment (Nombre d’articles publiés par le mentor) : Chaque augmentation d’une unité dans le nombre d’articles publiés par le mentor est associée à une augmentation moyenne de 1.8% dans le nombre d’articles publiés par le biochimiste.

Modèle d’Inflation de Zéro (Binomial avec lien logit)

  • ment (Nombre d’articles publiés par le mentor) : Le coefficient pour ment est -0.13007 (p = 0.00122), ce qui signifie que chaque augmentation d’une unité dans le nombre d’articles publiés par le mentor est associée à une diminution d’environ 13.01% de la probabilité que ses étudiants ne publient aucun article, soulignant son impact positif sur leur productivité scientifique.

Ajustement d’un modèle ZINB

Un modèle ZINB tiendra compte non seulement de l’inflation de zéro mais aussi de la surdispersion dans les données.

zinb_model <- zeroinfl(art ~ fem + mar + kid5 + phd + ment, data = df, dist = "negbin")
# Sélection de variables par élimination pas à pas basée sur l'AIC
step_model <- stepAIC(zinb_model, direction = "both")
## Start:  AIC=3125.98
## art ~ fem + mar + kid5 + phd + ment
## Warning in value[[3L]](cond): le système est numériquement singulier :
## conditionnement de la réciproque = 5.3869e-34FALSE
##        Df    AIC
## - phd   2 3122.0
## <none>    3126.0
## - mar   2 3128.0
## - fem   2 3131.2
## - kid5  2 3134.3
## - ment  2 3211.6
## 
## Step:  AIC=3122
## art ~ fem + mar + kid5 + ment
## 
##        Df    AIC
## <none>    3122.0
## - mar   2 3124.1
## + phd   2 3126.0
## - fem   2 3127.4
## - kid5  2 3130.3
## - ment  2 3216.3
zinb_model <- zeroinfl(art ~ fem + mar + kid5 + ment, data = df, dist = "negbin")
summary(zinb_model)
## 
## Call:
## zeroinfl(formula = art ~ fem + mar + kid5 + ment, data = df, dist = "negbin")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2939 -0.7625 -0.2912  0.4424  6.4196 
## 
## Count model coefficients (negbin with log link):
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.511184   0.196642   2.600  0.00933 ** 
## fem         -0.194715   0.075316  -2.585  0.00973 ** 
## mar          0.098105   0.083871   1.170  0.24211    
## kid5        -0.151707   0.054234  -2.797  0.00515 ** 
## ment         0.024749   0.003405   7.269 3.61e-13 ***
## Log(theta)   0.977494   0.135333   7.223 5.09e-13 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.4931     2.1900   0.225  0.82185   
## fem           0.6585     0.8235   0.800  0.42395   
## mar          -1.4691     0.8916  -1.648  0.09939 . 
## kid5          0.6217     0.4381   1.419  0.15586   
## ment         -0.8772     0.3127  -2.806  0.00502 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 2.6578 
## Number of iterations in BFGS optimization: 39 
## Log-likelihood: -1550 on 11 Df

Interprétation

Modèle de Comptage (Binomial négatif avec lien log)

  • Intercept (Intercept) : Lorsque toutes les autres variables sont nulles, le nombre attendu d’articles publiés est exp(0.511184) ≈ 1.67 articles.

  • fem (Genre - Women) : Les femmes publient en moyenne environ 17.7% (exp(-0.194715)) moins d’articles que les hommes, toutes choses étant égales par ailleurs.

  • kid5 : Chaque enfant de moins de cinq ans supplémentaire est associé à une diminution moyenne de 14.1% (exp(-0.151707)) dans le nombre d’articles publiés.

  • ment : Chaque augmentation d’une unité dans le nombre d’articles publiés par le mentor est associée à une augmentation moyenne de 2.5% (exp(0.024749)) dans le nombre d’articles publiés par le biochimiste.

    Log(theta) : Le paramètre de dispersion (theta) est 2.6578, suggérant une surdispersion dans les données.

Modèle d’Inflation de Zéro (Binomial avec lien logit)

  • ment (Nombre d’articles publiés par le mentor) : Chaque augmentation d’une unité dans le nombre d’articles publiés par le mentor est associée à une diminution de la probabilité de faire partie du groupe à zéros gonflés, avec un coefficient de -0.8772. Cet effet est statistiquement significatif (p = 0.00502), ce qui signifie que les mentors plus productifs réduisent la probabilité que leurs étudiants ne publient aucun article.

Nous remarquons que seule la variable ment est significative pour la composante zero, nous sommes donc tenté d’ajuster un modèle n’incluant qu’elle seule pour cette composante afin d’avoir un modèle plus performant.

print("AIC avec toutes les var de la composante logit:")
## [1] "AIC avec toutes les var de la composante logit:"
AIC(zinb_model)
## [1] 3121.997
zinb_model_ment <- zeroinfl(art ~ fem + mar + kid5 + ment| ment, data = df, dist = "negbin")

print("AIC avec avec ment uniquement dans la composante logit:")
## [1] "AIC avec avec ment uniquement dans la composante logit:"
AIC(zinb_model_ment)
## [1] 3122.545

Le modèle incluant uniquement la variable ment pour la composante logit n’améliore pas le critère AIC.

Test de rapport de vraisemblance (modèle vs intercept only pour le count )

# Modèle Zero-Inflated Negative Binomial (ZINB) avec uniquement intercept
zinb_model_intercept <- zeroinfl(art ~ fem + mar + kid5 + ment | 1, data = df, dist = "negbin")

# Test de rapport de vraisemblance entre les deux modèles ZINB
lrtest(zinb_model, zinb_model_intercept)
## Likelihood ratio test
## 
## Model 1: art ~ fem + mar + kid5 + ment
## Model 2: art ~ fem + mar + kid5 + ment | 1
##   #Df LogLik Df  Chisq Pr(>Chisq)    
## 1  11  -1550                         
## 2   7  -1561 -4 22.102  0.0001913 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modèle Zero-Inflated Poisson (ZIP) avec uniquement intercept
zip_model_intercept <- zeroinfl(art ~ fem + mar + kid5 + ment | 1, data = df, dist = "poisson")

# Test de rapport de vraisemblance entre les deux modèles ZIP
lrtest(zip_model, zip_model_intercept)
## Likelihood ratio test
## 
## Model 1: art ~ fem + mar + kid5 + ment | ment
## Model 2: art ~ fem + mar + kid5 + ment | 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   7 -1605.8                         
## 2   6 -1620.8 -1 30.058  4.192e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La statistique du test de rapport de vraisemblance (Chisq) est très significative pour les deux paires de modèles, avec des valeurs de p (Pr(>Chisq)) bien en dessous de 0.001. Cela indique que les modèles avec les memes régresseurs pour les deux composantes sont significativement meilleurs que ceux avec uniquement l’intercept pour la coposante zéro. Les résultats des tests de rapport de vraisemblance montrent que l’inclusion de toutes les variables (sauf phd) dans la composante de zéros améliore significativement l’ajustement du modèle pour les modèles ZINB et ZIP.

Comparaison des Modèles

# Ajustement du modèle Poisson
poisson_model <- glm(art ~ fem + mar + kid5 + ment, data = df, family = poisson)
# Calculer les AIC pour chaque modèle
aic_poisson <- AIC(poisson_model)
aic_nb_model <- AIC(nb_model)
aic_zip <- AIC(zip_model)
aic_zinb <- AIC(zinb_model)

# Afficher les AIC
aic_values <- data.frame(
  Model = c("Poisson", "Negative Binomial", "Zero-Inflated Poisson", "Zero-Inflated Negative Binomial"),
  AIC = c(aic_poisson, aic_nb_model, aic_zip, aic_zinb)
)
print(aic_values)
##                             Model      AIC
## 1                         Poisson 3312.349
## 2               Negative Binomial 3134.096
## 3           Zero-Inflated Poisson 3225.517
## 4 Zero-Inflated Negative Binomial 3121.997

Test de Vuong pour comparer les modèles avec inflation de zéros à leurs homologues

Le test de Vuong est utilisé pour comparer les modèles à inflation de zéros à leurs contreparties sans inflation de zéros.

vuong_test = vuong(zip_model, poisson_model)
## Vuong Non-Nested Hypothesis Test-Statistic: 
## (test-statistic is asymptotically distributed N(0,1) under the
##  null that the models are indistinguishible)
## -------------------------------------------------------------
##               Vuong z-statistic             H_A    p-value
## Raw                    4.107180 model1 > model2 2.0026e-05
## AIC-corrected          3.926310 model1 > model2 4.3130e-05
## BIC-corrected          3.490509 model1 > model2 0.00024105
print(vuong_test)
## NULL
vuong_test <- vuong(zinb_model, nb_model)
## Vuong Non-Nested Hypothesis Test-Statistic: 
## (test-statistic is asymptotically distributed N(0,1) under the
##  null that the models are indistinguishible)
## -------------------------------------------------------------
##               Vuong z-statistic             H_A p-value
## Raw                    2.240175 model1 > model2 0.01254
## AIC-corrected          1.226485 model1 > model2 0.11001
## BIC-corrected         -1.215962 model2 > model1 0.11200
print(vuong_test)
## NULL
# Créer un tableau des résultats du test de Vuong (Raw uniquement)
results_raw <- data.frame(
  Comparison = c("ZIP vs Poisson", "ZINB vs NB"),
  Vuong_z_statistic = c(4.107180, 2.240175),
  Hypothesis_Alternative = c("model1 > model2", "model1 > model2"),
  p_value = c(2.0026e-05, 0.01254)
)

# Imprimer le tableau des résultats
print(results_raw)
##       Comparison Vuong_z_statistic Hypothesis_Alternative    p_value
## 1 ZIP vs Poisson          4.107180        model1 > model2 2.0026e-05
## 2     ZINB vs NB          2.240175        model1 > model2 1.2540e-02

  1. ZIP vs Poisson: Le test de Vuong avec une statistique z de 4.107180 et une p-value de 2.0026e-05 indique que le modèle ZIP est significativement meilleur que le modèle Poisson.

  2. ZINB vs NB: Le test de Vuong avec une statistique z de 2.240175 et une p-value de 0.01254 montre que le modèle ZINB est significativement meilleur que le modèle NB.

Ce tableau et cette conclusion montrent clairement la supériorité des modèles ZIP et ZINB par rapport à leurs équivalents simples (Poisson et NB) selon le test de Vuong, basé sur les résultats Raw.

Évaluation de la Précision des Modèles pour la Prédiction des Zéros

# Calculer le pourcentage de zéros observés dans les données
observed_zeros_percentage <- mean(df$art == 0) * 100

# Prédire le pourcentage de zéros pour chaque modèle
zinb_zero_probs <- predict(zinb_model, type = "prob")[, 1]
zip_zero_probs <- predict(zip_model, type = "prob")[, 1]

# Calculer les probabilités de zéros pour les modèles Poisson et NB
poisson_means <- predict(poisson_model, type = "response")
nb_means <- predict(nb_model, type = "response")
nb_theta <- nb_model$theta

poisson_zero_probs <- dpois(0, lambda = poisson_means)
nb_zero_probs <- dnbinom(0, mu = nb_means, size = nb_theta)

# Moyenne des probabilités de zéros prédites pour chaque modèle
predicted_zeros_percentage <- c(
  mean(zinb_zero_probs) * 100,
  mean(nb_zero_probs) * 100,
  mean(poisson_zero_probs) * 100,
  mean(zip_zero_probs) * 100
)

# Créer un tableau récapitulatif des résultats
results <- data.frame(
  Model = c("ZINB", "NB", "Poisson", "ZIP"),
  Observed_Counts = rep(observed_zeros_percentage, 4),  # Remplacer "*" par les valeurs observées si disponibles
  Expected_Counts = round(predicted_zeros_percentage, 2)  # Arrondir les prédictions à deux décimales
)

# Imprimer les résultats
print(results)
##     Model Observed_Counts Expected_Counts
## 1    ZINB        30.05464           31.20
## 2      NB        30.05464           30.35
## 3 Poisson        30.05464           20.90
## 4     ZIP        30.05464           29.79

Les modèles ZINB et NB prédisent des pourcentages de zéros proches de l’observé (31.20% et 30.35% respectivement). Cela indique leur efficacité pour capturer l’overdispersion et l’inflation des zéros pour le ZINB. Le modèle Poisson sous-estime les zéros à 20.90%, tandis que le modèle ZIP prédit 29.79%, se rapprochant de l’observé (30.05%).

Nous sommes quand même tenté d’évaluer la prédicion de zéro du modèle ZINB n’incluant que l’intercepte pour la composante zéro, puisque ’après nos analyses passées, toutes les variables de cette composante sont non significatives sauf la varible ment, mais le modèle n’incluant que la variable ment pour cette composante n’améliore pas le AIC.

Prédiction du modèle intercepté pour la composante zéro

zinb_zero_probs <- predict(zinb_model_intercept, type = "prob")[, 1]
# Imprimer les résultats
print("Le pourcentage de zéros prédit par le modèle ZINB intercepté est:")
## [1] "Le pourcentage de zéros prédit par le modèle ZINB intercepté est:"
# Moyenne des probabilités de zéros prédites pour chaque modèle
predicted_zeros_percentage <-mean(zinb_zero_probs) * 100
predicted_zeros_percentage
## [1] 30.35516

En conclusion

Le modèle incluant uniquement la variable ment pour la composante logit n’améliore pas significativement le critère AIC. Cela indique que la productivité d’un mentor de thèse n’a pas un impact aussi significatif sur la probabilité qu’un étudiant ne publie aucun article (de façon global). De plus le modèle intercepté prédit mieux les zéros que tous les autres modèles de l’analyse: 30.35% de zéros prédits pour 30.05% observés.

Par conséquent notre choix du meilleur modèle se porte sur le modèle ZINB intercepté.

Pertinence

# Test de Hosmer Lemeshow 
hoslem.test(df$art, fitted(zinb_model_intercept), g = 10)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  df$art, fitted(zinb_model_intercept)
## X-squared = -120.63, df = 8, p-value = 1

Vérification des résidus pour détecter les schémas systématiques.

# Obtenir les valeurs ajustées et les résidus pour chaque modèle
fitted_values_zip <- fitted(zip_model)
residuals_zip <- residuals(zip_model, type = "pearson")

fitted_values_zinb <- fitted(zinb_model_intercept)
residuals_zinb <- residuals(zinb_model_intercept, type = "pearson")

fitted_values_poisson <- fitted(poisson_model)
residuals_poisson <- residuals(poisson_model, type = "pearson")

fitted_values_nb <- fitted(nb_model)
residuals_nb <- residuals(nb_model, type = "pearson")

# Créer un data frame pour les valeurs ajustées et les résidus
data <- data.frame(
  Model = rep(c("ZIP", "ZINB", "Poisson", "Negative Binomial"), each = length(fitted_values_zip)),
  Fitted_Values = c(fitted_values_zip, fitted_values_zinb, fitted_values_poisson, fitted_values_nb),
  Residuals = c(residuals_zip, residuals_zinb, residuals_poisson, residuals_nb)
)

# Installer ggplot2 si nécessaire
# install.packages("ggplot2")
library(ggplot2)

# Tracer les graphiques de valeurs ajustées par rapport aux résidus avec ggplot2
ggplot(data, aes(x = Fitted_Values, y = Residuals, color = Model)) +
  geom_point() +
  labs(x = "Fitted Values", y = "Residuals", title = "Fitted Values vs Residuals by Model") +
  theme_minimal() +
  facet_wrap(~ Model, scales = "free") +
  theme(legend.position = "top")

En examinant les graphiques de résidus, il est clair qu’ajuster un modèle prenant en compte la surdispersion et l’inflation des zéros est une bonne idée. De plus, l’analyse du graphique de résidus pour le modèle binomial négatif à inflation de zéros montre des résidus plus faibles. Cela suggère qu’il est préférable d’ajuster un modèle binomial négatif avec inflation de zéros art ~ fem + mar + kid5 + ment | 1, c’est-à-dire un modèle sans régresseurs pour la partie zéro du modèle (puisqu’aucun d’entre eux ne semble significatif). Ce que nous avons fait plutôt.

Conclusion

L’analyse des données sur le nombre d’articles publiés par les doctorants en biochimie révèle la présence de surdispersion et d’inflation de zéros. Les modèles de Poisson et Binomial Négatif ne sont pas adaptés à ces données, car ils sous-estiment considérablement le nombre de zéros observés.

Les modèles ZINB et ZIP, qui tiennent compte de l’inflation des zéros, se révèlent être significativement meilleurs que les modèles Poisson et Binomial Négatif selon le test de Vuong. Ils prédisent un pourcentage de zéros proche de l’observé, indiquant leur capacité à capturer les caractéristiques particulières des données.

L’analyse des variables explicatives suggère que les variables fem, mar, kid5 et ment sont les plus pertinentes pour expliquer le nombre d’articles publiés. La variable phd n’améliore pas significativement l’ajustement du modèle et peut être supprimée vu qu’elle présente déja une forte correlation avec la variable ment.

Concernant la composante zéro du modèle, seul le nombre d’articles publiés par le mentor (“ment”) semble avoir un impact significatif. Cependant, l’inclusion de cette variable dans la composante zéro n’améliore pas le critère AIC. Le modèle ZINB intercepté (sans variable dans la composante zéro) prédit avec précision le pourcentage de zéros observé, ce qui en fait le meilleur choix pour modéliser les données.

En conclusion, le modèle ZINB , incluant les variables fem, mar, kid5 et ment dans la composante de comptageet uniquement l’intercept pour la composante zéro représente le modèle le plus performant pour ces données. Ce modèle capture efficacement la surdispersion et l’inflation de zéros, prédit avec précision l’inflation de zéro dans les données.

Des recherches supplémentaires pourraient être menées pour explorer l’impact potentiel d’autres variables non incluses dans l’analyse actuelle.

L’analyse pourrait être étendue à d’autres disciplines scientifiques pour évaluer la généralisabilité des conclusions.