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

13

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.

A1+e(b0+b1x1+b2x2+),
A

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

Merci!

user5497
la source
Lorsque vous dites " pourrait être donné", voulez-vous dire qu'il est spécifié par l'utilisateur ou estimé par le modèle? A
Macro
d'une manière ou d'une autre. Je peux le mettre en entrée (egA = 0,25) ou être l'un des coefficients à trouver
user5497
Cela varie-t-il d'un sujet à l'autre (c.-à-d. S'agit-il de données) ou s'agit-il d'une constante fixe dans toutes les observations?
Macro
corrigé sur toutes les observations (comme b0, b1, ...)
user5497
2
Pourquoi ne pas utiliser le maximum de vraisemblance au lieu de GMM?
Macro

Réponses:

6

En supposant , ce modèle a la variable de réponse de Bernoulli Y i avecA1Yi

Pr(Yi=1)=A1+eXib,

(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.bUNEXii

Les conditions de moment sont:

E[(YiA1+eXib)Xi]=0.

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

m=1Ni=1N[(YiA1+eXib)Xi]=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).mmb

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:(YiA1+eXib)Xii

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:::.obj1mmoptimgmm

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 
josliber
la source