Utilisation de la méthode des moments généralisés (GMM) pour calculer le paramètre de régression logistique

Je veux calculer les coefficients d'une régression très similaire à la régression logistique (en fait, régression logistique avec un autre coefficient: quandApourrait être donné). J'ai pensé à utiliser GMM pour calculer les coefficients, mais je ne sais pas quelles sont les conditions de moment que je devrais utiliser.

$$\frac{A}{1 + e^{- (b_0 + b_1 x_1 + b_2 x_2 + \ldots)}},$$ $A$

Est-ce que quelqu'un peut m'aider avec cela?

Merci!

En supposant , ce modèle a la variable de réponse de Bernoulli Y i avec$A\leq 1$$Y_i$

$$Pr(Y_i = 1) = \frac{A}{1+e^{-X_i'b}},$$

où (et éventuellement A , selon qu'il est traité comme une constante ou un paramètre) sont les coefficients ajustés et X i est les données pour l'observation i . Je suppose que le terme d'interception est géré en ajoutant une variable de valeur constante 1 à la matrice de données.$b$$A$$X_i$$i$

Les conditions de moment sont:

$$\begin{align*} \mathbb{E}\bigg[\bigg(Y_i-\frac{A}{1+e^{-X_i'b}}\bigg)X_i\bigg] &= 0. \end{align*}$$

Nous remplaçons cela par l'échantillon homologue de la condition, en supposant observations:$N$

$$m = \frac{1}{N}\sum_{i=1}^N \bigg[\bigg(Y_i-\frac{A}{1+e^{-X_i'b}}\bigg)X_i\bigg] = 0$$

Ceci est pratiquement résolu en minimisant sur toutes les valeurs de coefficient possibles b (ci-dessous, nous utiliserons le simplexe de Nelder-Mead pour effectuer cette optimisation).$m'm$$b$

Empruntant à un excellent tutoriel R-blogueurs sur le sujet , il est assez simple de l'implémenter dans R avec le package gmm. À titre d'exemple, travaillons avec le jeu de données iris, en prédisant si un iris est versicolore en fonction de sa longueur et de sa largeur sépale et de sa longueur et de sa largeur de pétale. Je suppose que est constant et égal à 1 dans ce cas:$A$

dat <- as.matrix(cbind(data.frame(IsVersicolor = as.numeric(iris$Species == "versicolor"), Intercept=1), iris[,1:4]))
head(dat)
#      IsVersicolor Intercept Sepal.Length Sepal.Width Petal.Length Petal.Width
# [1,]            0         1          5.1         3.5          1.4         0.2
# [2,]            0         1          4.9         3.0          1.4         0.2
# [3,]            0         1          4.7         3.2          1.3         0.2
# [4,]            0         1          4.6         3.1          1.5         0.2
# [5,]            0         1          5.0         3.6          1.4         0.2
# [6,]            0         1          5.4         3.9          1.7         0.4

Voici les coefficients ajustés à l'aide de la régression logistique:

summary(glm(IsVersicolor~., data=as.data.frame(dat[,-2]), family="binomial"))
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)    7.3785     2.4993   2.952 0.003155 ** 
# Sepal.Length  -0.2454     0.6496  -0.378 0.705634    
# Sepal.Width   -2.7966     0.7835  -3.569 0.000358 ***
# Petal.Length   1.3136     0.6838   1.921 0.054713 .  
# Petal.Width   -2.7783     1.1731  -2.368 0.017868 *  

La pièce principale dont nous avons besoin pour utiliser gmm est une fonction qui renvoie les conditions de moment, à savoir les lignes pour chaque observationi:$(Y_i-\frac{A}{1+e^{-X_i'b}})X_i$$i$

moments <- function(b, X) {
  A <- 1
  as.vector(X[,1] - A / (1 + exp(-(X[,-1] %*% cbind(b))))) * X[,-1]
}

$b$

init.coef <- lm(IsVersicolor~., data=as.data.frame(dat[,-2]))$coefficients
library(gmm)
fitted <- gmm(moments, x = dat, t0 = init.coef, type = "iterative", crit = 1e-19,
              wmatrix = "optimal", method = "Nelder-Mead",
              control = list(reltol = 1e-19, maxit = 20000))
fitted
#  (Intercept)  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width  
#      7.37849      -0.24536      -2.79657       1.31364      -2.77834  
# 
# Convergence code =  0 

Le code de convergence de 0 indique la procédure convergée et les paramètres sont identiques à ceux retournés par régression logistique.

momentEstim.baseGmm.iterativegmm:::.obj1$m'm$optimgmm

gmm.objective <- function(theta, x, momentFun) {
  avg.moment <- colMeans(momentFun(theta, x))
  sum(avg.moment^2)
}
optim(init.coef, gmm.objective, x=dat, momentFun=moments,
      control = list(reltol = 1e-19, maxit = 20000))$par
#  (Intercept) Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
#    7.3784866   -0.2453567   -2.7965681    1.3136433   -2.7783439