Berechnung der logistischen Regressionsparameter mithilfe der verallgemeinerten Momentenmethode (GMM)

Ich möchte Koeffizienten zu einer Regression berechnen, die der logistischen Regression sehr ähnlich wenn könnte). Ich habe überlegt, GMM zur Berechnung der Koeffizienten zu verwenden, bin mir aber nicht sicher, welche Bedingungen ich im Moment verwenden soll.

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

Kann mir jemand dabei helfen?

Vielen Dank!

Unter der Annahme von hat dieses Modell die Bernoulli-Antwortvariable Y i mit$A\leq 1$$Y_i$

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

wobei (und möglicherweise A , abhängig davon, ob es als Konstante oder als Parameter behandelt wird) die angepassten Koeffizienten sind und X i die Daten für die Beobachtung i sind . Ich gehe davon aus, dass der Intercept-Term behandelt wird, indem der Datenmatrix eine Variable mit dem konstanten Wert 1 hinzugefügt wird.$b$$A$$X_i$$i$

Die momentanen Bedingungen sind:

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

Wir ersetzen dies durch das Beispielgegenstück der Bedingung unter der Annahme von Beobachtungen:$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$$

Dies wird praktisch gelöst, indem über alle möglichen Koeffizientenwerte b minimiert wird (im Folgenden wird der Nelder-Mead-Simplex verwendet , um diese Optimierung durchzuführen).$m'm$$b$

Aus einem hervorragenden Tutorial von R-Bloggern zu diesem Thema entlehnt, ist es ziemlich einfach, dies mit dem gmm-Paket in R zu implementieren. Als Beispiel wollen wir mit dem Iris-Dataset arbeiten und anhand der Länge und Breite der Kelchblätter und der Länge und Breite der Blütenblätter vorhersagen, ob eine Iris mehrfarbig ist. Ich gehe davon aus, dass konstant und in diesem Fall gleich 1 ist:$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

Hier sind die Koeffizienten, die mithilfe der logistischen Regression angepasst wurden:

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 *  

Das Hauptstück wir benutzen müssen gmm ist eine Funktion, kehrt die Momentbedingungen, nämlich Zeilen für jede Beobachtungi:$(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]
}

Wir können nun die Koeffizienten numerisch anpassen , indem wir die linearen Regressionskoeffizienten als geeigneten Anfangspunkt verwenden (wie im oben verlinkten Lernprogramm vorgeschlagen):$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 

Der Konvergenzcode 0 gibt die konvergierte Prozedur an und die Parameter sind mit denen identisch, die von der logistischen Regression zurückgegeben werden.

Ein kurzer Blick auf die gmm-Paketquelle (Funktionen momentEstim.baseGmm.iterativeund gmm:::.obj1die bereitgestellten Parameter) zeigt, dass das gmm-Paket wie oben angegeben minimiert . Der folgende äquivalente Code ruft die R- Funktion direkt auf und führt die gleiche Optimierung durch, die wir oben mit dem Aufruf von erreicht haben :$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