Berechnung der logistischen Regressionsparameter mithilfe der verallgemeinerten Momentenmethode (GMM)

13

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.

EIN1+e-(b0+b1x1+b2x2+),
EIN

Kann mir jemand dabei helfen?

Vielen Dank!

user5497
quelle
Wenn Sie " könnte gegeben werden" sagen, meinen Sie damit, dass es vom Benutzer angegeben oder vom Modell geschätzt wird? EIN
Makro
in jedem Fall. Ich kann es als Eingabe eingeben (zB A = 0,25) oder es kann einer der zu findenden Koeffizienten sein
user5497
Unterscheidet es sich von Subjekt zu Subjekt (dh sind es Daten) oder ist es eine feste Konstante für alle Beobachtungen?
Makro
Fix für alle Beobachtungen (wie b0, b1, ...)
user5497
2
Warum nicht Maximum Likelihood anstelle von GMM verwenden?
Makro

Antworten:

6

Unter der Annahme von hat dieses Modell die Bernoulli-Antwortvariable Y i mitEIN1Y.ich

Pr(Y.ich=1)=EIN1+e-Xichb,

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.bEINXichich

Die momentanen Bedingungen sind:

E[(Y.ich-EIN1+e-Xichb)Xich]=0.

Wir ersetzen dies durch das Beispielgegenstück der Bedingung unter der Annahme von Beobachtungen:N

m=1Nich=1N[(Y.ich-EIN1+e-Xichb)Xich]=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).mmb

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:EIN

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.ich-EIN1+e-Xichb)Xichich

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 :mmoptimgmm

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
quelle