4.4 Applications en R
4.4.1 Binomiale négative 2
Utilisons de nouveau la base de données du cours pour estimer un modèle de régression NB2.
Comme pour le modèle Poisson, un algorithme de type Newton-Rapshon peut être programmé pour estimer les paramètres de la régression NB2.
### Algorithme Newton-Raphson maison
## Initialisation
# Generer la matrice des X
score.glm <- as.formula(nb.sin ~ couleur + myopie + type_territoire + langue + alimentation)
Matrix.X <- model.matrix(score.glm, db.train)
# Generer le vecteur des Y
Vect.Y <- as.matrix(db.train$nb.sin)
Vect.expo <- as.matrix(db.train$expo)
# Valeur initiale des beta
Beta <- rep(0, ncol(Matrix.X))
Beta[1] <- log(mean(Vect.Y))
# Valeur initiale de alpha (attention doit être positif)
alpha <- 0.1
## Loops Newton-Raphson
max.U <- 1
arret <- 0
loop <- 0
while (arret < 1){
loop <- loop + 1
## lambda
lambda <- Vect.expo*exp(Matrix.X %*% Beta)
## Vecteur U (en deux parties: pour beta et alpha)
U1 <- t(Matrix.X) %*% ((Vect.Y - lambda)/(1 +lambda*alpha))
fct.1 <- 0
for(jj in 0:(max(Vect.Y)-1)){
fct.1 <- ifelse((Vect.Y-1) >= jj, fct.1 + 1/(jj + 1/alpha), fct.1)
}
U2a <- (log(1 +lambda*alpha) - fct.1)/alpha^2 + (Vect.Y - lambda)/(alpha*(1 +lambda*alpha))
U2 <- sum(U2a)
U <- rbind(U1, U2)
## Matrice H - au lieu de prendre les doubles dérivées, prenons OP
OP1 <- Matrix.X * as.vector((Vect.Y - lambda)/(1 + lambda*alpha))
OP2 <- cbind(OP1, U2a)
OP <- t(OP2) %*% OP2
## Update beta
parm <- as.vector(c(Beta, alpha))
parm <- parm + solve(OP)%*%U
Beta <- parm[1: ncol(Matrix.X)]
# attention alpha doit être positif
alpha <- max(0.01, parm[ncol(Matrix.X)+1])
## Criteres d'arret
arret <- ifelse(abs(max(U) - max.U) > 0.0001 & loop < 60, 0, 1)
max.U <- max(U)
}
## Print Results
rownames(parm)[ncol(Matrix.X)+1] <- 'alpha'
print(parm)## [,1]
## (Intercept) -1.19928847
## couleurRouge -0.02153588
## myopieOui 0.07534476
## type_territoireSemi-urbain -0.09520287
## type_territoireUrbain -0.14951787
## langueF 0.25203684
## alimentationVégétalien -0.24909662
## alimentationVégétarien -0.43673082
## alpha 0.63662213On peut remarquer que l’estimateur OP a été utilisé au lieu de l’hessienne, pour éviter le codage difficile de l’hessienne. On peut prendre la même astuce pour avoir un estimé de la variance des estimateurs.
Notez aussi que quelle manière le vecteur U a été construit: en collant (rbind) la dérivée du logvraisemblance par \(\beta\) à la dérivée du logvraisemblance par \(\alpha\). Cela fait en sorte que l’algorithme de Newton-Raphson cherche simultanément les valeurs des paramètres \(\beta\) et \(\alpha\) qui atteignent la racine de 0.
lambda <- Vect.expo*exp(Matrix.X %*% Beta)
fct.1 <- 0
for(jj in 0:(max(Vect.Y)-1)){
fct.1 <- ifelse((Vect.Y-1) >= jj, fct.1 + 1/(jj + 1/alpha), fct.1)
}
U2a <- (log(1 +lambda*alpha) - fct.1)/alpha^2 + (Vect.Y - lambda)/(alpha*(1 +lambda*alpha))
OP1 <- Matrix.X * as.vector((Vect.Y - lambda)/(1 + lambda*alpha))
OP2 <- cbind(OP1, U2a)
OP <- t(OP2) %*% OP2
varOP <- sqrt(diag(solve(OP)))
cbind(parm, varOP)## varOP
## (Intercept) -1.19928847 0.02731064
## couleurRouge -0.02153588 0.02016462
## myopieOui 0.07534476 0.01454167
## type_territoireSemi-urbain -0.09520287 0.01361203
## type_territoireUrbain -0.14951787 0.01303562
## langueF 0.25203684 0.01071314
## alimentationVégétalien -0.24909662 0.01927230
## alimentationVégétarien -0.43673082 0.01832446
## alpha 0.63662213 0.02000385On peut comparer avec certains packages R qui proposent la régression binomiale négative: la package MASS et le package gamlss.
library(MASS)
library(gamlss)On débute avec MASS. Il est à noter que le package donne la valeur de \(1/\alpha\) et non \(\alpha\), comme notre notation de la NB2.
score.nb <- as.formula(nb.sin ~ couleur + myopie + type_territoire + langue + alimentation + offset(log(expo)))
model.nb2.MASS <- glm.nb(score.nb, data=db.train)
print(summary(model.nb2.MASS),digits=3)##
## Call:
## glm.nb(formula = score.nb, data = db.train, init.theta = 1.570790499,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.861 -0.645 -0.573 -0.416 4.367
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1993 0.0274 -43.76 < 2e-16 ***
## couleurRouge -0.0215 0.0201 -1.07 0.28
## myopieOui 0.0753 0.0144 5.23 1.7e-07 ***
## type_territoireSemi-urbain -0.0952 0.0135 -7.06 1.6e-12 ***
## type_territoireUrbain -0.1495 0.0130 -11.47 < 2e-16 ***
## langueF 0.2520 0.0108 23.33 < 2e-16 ***
## alimentationVégétalien -0.2491 0.0191 -13.02 < 2e-16 ***
## alimentationVégétarien -0.4367 0.0182 -24.01 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.5708) family taken to be 1)
##
## Null deviance: 166034 on 270714 degrees of freedom
## Residual deviance: 163931 on 270707 degrees of freedom
## AIC: 283024
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.5708
## Std. Err.: 0.0492
##
## 2 x log-likelihood: -283006.2330On poursuit avec le package gamlss. Dans ce cas, il faut remarquer que:
- gamlss donne la valeur de \(\log(\alpha)\);
- gamlss note la distribution NB2 du cours comme une “NBI”;
- gamlss note la distribution NB1 du cours comme une “NBII”.
model.nb2.gamlss <- gamlss(score.nb,family=NBI, data=db.train)## GAMLSS-RS iteration 1: Global Deviance = 283006.2
## GAMLSS-RS iteration 2: Global Deviance = 283006.2
## GAMLSS-RS iteration 3: Global Deviance = 283006.2print(summary(model.nb2.gamlss),digits=3)## ******************************************************************
## Family: c("NBI", "Negative Binomial type I")
##
## Call: gamlss(formula = score.nb, family = NBI, data = db.train)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.19929 0.02741 -43.755 < 2e-16 ***
## couleurRouge -0.02154 0.02008 -1.072 0.284
## myopieOui 0.07534 0.01440 5.232 1.68e-07 ***
## type_territoireSemi-urbain -0.09520 0.01348 -7.062 1.65e-12 ***
## type_territoireUrbain -0.14952 0.01304 -11.468 < 2e-16 ***
## langueF 0.25204 0.01080 23.328 < 2e-16 ***
## alimentationVégétalien -0.24910 0.01914 -13.017 < 2e-16 ***
## alimentationVégétarien -0.43673 0.01819 -24.011 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.45158 0.03138 -14.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 270715
## Degrees of Freedom for the fit: 9
## Residual Deg. of Freedom: 270706
## at cycle: 3
##
## Global Deviance: 283006.2
## AIC: 283024.2
## SBC: 283118.8
## ******************************************************************
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1993 0.0274 -43.76 0.00e+00
## couleurRouge -0.0215 0.0201 -1.07 2.84e-01
## myopieOui 0.0753 0.0144 5.23 1.68e-07
## type_territoireSemi-urbain -0.0952 0.0135 -7.06 1.65e-12
## type_territoireUrbain -0.1495 0.0130 -11.47 1.94e-30
## langueF 0.2520 0.0108 23.33 3.00e-120
## alimentationVégétalien -0.2491 0.0191 -13.02 1.00e-38
## alimentationVégétarien -0.4367 0.0182 -24.01 2.92e-127
## (Intercept) -0.4516 0.0314 -14.39 6.19e-47On peut ainsi comparer les différents paramètres obtenus par les 3 méthodes.
MASS.parm <- c(model.nb2.MASS$coefficients, 1/model.nb2.MASS$theta)
gamlss.parm <- c(model.nb2.gamlss$mu.coefficients, exp(model.nb2.gamlss$sigma.coefficients))
cbind(parm, MASS.parm, gamlss.parm)## MASS.parm gamlss.parm
## (Intercept) -1.19928847 -1.19928847 -1.19928846
## couleurRouge -0.02153588 -0.02153588 -0.02153588
## myopieOui 0.07534476 0.07534476 0.07534476
## type_territoireSemi-urbain -0.09520287 -0.09520287 -0.09520287
## type_territoireUrbain -0.14951787 -0.14951787 -0.14951787
## langueF 0.25203684 0.25203684 0.25203684
## alimentationVégétalien -0.24909662 -0.24909662 -0.24909663
## alimentationVégétarien -0.43673082 -0.43673082 -0.43673083
## alpha 0.63662213 0.63662213 0.63662213On peut aussi comparer la variance des estimateurs.
MASS.stderr <- c(summary(model.nb2.MASS)$coefficients[, 2], model.nb2.MASS$SE.theta)
cbind(varOP, MASS.stderr)## varOP MASS.stderr
## (Intercept) 0.02731064 0.02740908
## couleurRouge 0.02016462 0.02008409
## myopieOui 0.01454167 0.01440132
## type_territoireSemi-urbain 0.01361203 0.01348186
## type_territoireUrbain 0.01303562 0.01303780
## langueF 0.01071314 0.01080385
## alimentationVégétalien 0.01927230 0.01913575
## alimentationVégétarien 0.01832446 0.01818895
## 0.02000385 0.04923304Pour le moment, je n’arrive pas à extraire directement les écart-types de la fonction gamlss. Mais on peut quand même connaitre ces valeurs en prenant le summary() de la fonction, comme on le voyait plus haut.
On peut voir que les écart-types des \(\beta\) sont similaires pour les trois méthodes. Pour ce qui en est du paramètre \(\alpha\), les deux autres méthodes ne génèrent pas la même forme alors les écart-types ne sont pas directement. comparables. La méthode delta serait bien utile ici.
4.4.2 Binomiale négative 1
L’application pour la binomiale négative 1 (NB1) est en devoir. Le même développement mathématique et informatique que celui développé pour la NB2 peut être fait. Pour s’assurer que les résultats obtenus sont bons, on pourra comparer avec le package gamlss avec la distribution NBII. Le package MASS ne propose pas la NB1.c