glmnet: Wie kann man die multinomiale Parametrisierung verstehen?

Folgendes Problem: Ich möchte eine kategoriale Antwortvariable mit einer (oder mehreren) kategorialen Variablen mithilfe von glmnet () vorhersagen.

Ich kann jedoch keinen Sinn für die Ausgabe machen, die glmnet mir gibt.

Ok, zuerst generieren wir zwei verwandte kategoriale Variablen:

Daten generieren

p <- 2 #number variables
mu <- rep(0,p)
sigma <- matrix(rep(0,p^2), ncol=p)
sigma[1,2] <- .8 #some relationship ..
diag(sigma) <- 1
sigma <- pmax(sigma, t(sigma))
n <- 100
set.seed(1)
library(MASS)
dat <- mvrnorm(n, mu, sigma)
#discretize
k <- 3 # number of categories
d <- apply(dat, 2, function(x) {
  q <- quantile(x, probs=seq(0,1, 1/k))[-c(1, k+1)]
  out <- numeric(length(x))
  for(i in 1:(k-1))
  {  out[x<q[k-i]] <- i } 
  return(out)
})
d <- data.frame(apply(d, 2, as.factor))
d[,2] <- relevel(d[,2], ref="0")
d[,1] <- relevel(d[,1], ref="0")
colnames(d) <- c("X1", "X2")

Wir bekommen:

> table(d)
   X2
X1   0  1  2
  0 22 11  1
  1  9 14 10
  2  3  8 22

Vorhersage: multinom ()

Dann lassen Sie uns X1 durch X2 mit dem multinom () aus dem nnet-Paket vorhersagen:

library(nnet)
mod1 <- multinom(X1~X2, data=d)
mod1

was uns gibt:

Call:
multinom(formula = X1 ~ X2, data = d)

Coefficients:
  (Intercept)      X21      X22
1  -0.8938246 1.134993 3.196476
2  -1.9924124 1.673949 5.083518

Manuelle Überprüfung

Lassen Sie uns nun prüfen, ob wir das manuell reproduzieren können:

tb <- table(d)
log(tb[2,1] / tb[1,1]) #intercept category1
[1] -0.8938179
log(tb[3,1] / tb[1,1]) #intercept category2
[1] -1.99243
log((tb[1,1]*tb[2,2]) / (tb[1,2]*tb[2,1])) #logodds-ratio cat X1 0vs1 in X2 0vs1
[1] 1.13498
#same for the three remaining log odds ratios

Wir produzieren die gleichen Zahlen, gut!

Vorhersage: glmnet ()

Jetzt machen wir dasselbe mit glmnet:

library(glmnet)
y <- d[,1]
X <- model.matrix(X1~X2, data=d)[,-1]
mod2 <- glmnet(X, y, family="multinomial", lambda=c(0))
coef(mod2, s=0) #meaning of coefficients unclear!
$`0`
3 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept)  0.9620216
X21         -1.1349130
X22         -3.1958293   

$`1`
3 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept) 0.06825755
X21         .         
X22         .         

$`2`
3 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept) -1.0302792
X21          0.5388814
X22          1.8870363

Beachten Sie, dass ich s = 0 setze, daher gibt es keine Regularisierung und die Parameter sollten genau die gleichen Informationen enthalten wie die Parameter der Funktion multinom ().

Trotzdem erhalten wir sehr unterschiedliche Parameter. Dies liegt an der unterschiedlichen Parametrisierung, die sie in glmnet verwenden, siehe z.

http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (Überschrift: Multinomial Models) oder das entsprechende Papier: http://www.jstatsoft.org/v33/i01/paper (Überschrift: 4. Regularized multinomiale Regression)

$P(Y=k | X)$

Bedingte Wahrscheinlichkeiten: multinom ()

Also berechne ich zuerst diese Wahrscheinlichkeiten aus multinom ():

p.fit <- predict(mod1, type="probs")
head(d)
head(p.fit)
ccp <- matrix(0,3,3)
ccp[,3] <- p.fit[1,]
ccp[,2] <- p.fit[2,]
ccp[,1] <- p.fit[4,]
ccp
           [,1]      [,2]       [,3]
[1,] 0.64705896 0.3333332 0.03030114
[2,] 0.26470416 0.4242450 0.30303140
[3,] 0.08823688 0.2424218 0.66666746
colSums(ccp) #sum to 1, ok; sorry for the awful code ...
[1] 1 1 1

Da wir hier ein gesättigtes Modell haben, sollte dies das gleiche sein, was wir aus den Daten berechnen können:

emp <- table(d)/100
cemp <- apply(emp, 2, function(x) {
  x / sum(x)
})
cemp 
   X2
             0         1          2
  0 0.64705882 0.3333333 0.03030303
  1 0.26470588 0.4242424 0.30303030
  2 0.08823529 0.2424242 0.66666667

was in der Tat der Fall ist.

Bedingte Wahrscheinlichkeiten: glmnet ()

Nun das gleiche von glmnet:

c1 <- coef(mod2, s=0)
c <-matrix(rapply(c1, function(x) { as.matrix(x)}, how="unlist"), 3,3, byrow=T)

ccp2 <- matrix(0,3,3)
config <- rbind(c(0,0), c(1,0), c(0,1))

for(l in 1:3) #loop through categories
{
  denom <- numeric(3)
  for(i in 1:3) # loop through possible predictor combinations
  { 
    x1 <- config[i, 1]
    x2 <- config[i, 2]
    denom[i] <- exp(c[l,1] + x1 * c[l,2]  + x2 * c[l,3])
  }
  ccp2[l,1] <- denom[1] / sum(denom)
  ccp2[l,2] <- denom[2] / sum(denom)
  ccp2[l,3] <- denom[3] / sum(denom)
}
ccp2
          [,1]      [,2]       [,3]
[1,] 0.7340082 0.2359470 0.03004484
[2,] 0.3333333 0.3333333 0.33333333
[3,] 0.1073668 0.1840361 0.70859708
colSums(ccp2)
[1] 1.1747083 0.7533165 1.0719753

Die zellbedingten Wahrscheinlichkeiten sind etwas verwandt, aber unterschiedlich. Auch summieren sie sich nicht zu einem.

Wir haben hier also zwei Probleme:

a) Die bedingten Wahrscheinlichkeiten summieren sich nicht zu 1 und

b) Die Parameter beschreiben nicht, was wir in den Daten sehen: z. B. gibt es in Zeile 2 Unterschiede zwischen den Spalten, aber glmnet schätzt beide Koeffizienten (nicht den Achsenabschnitt) als Null.

Ich habe ein lineares Regressionsproblem verwendet und glm und glmnet mit s = 0 verglichen, um sicherzustellen, dass s = 0 eine Regularisierung von Null bedeutet (die Lösungen waren fast identisch).

Jede Hilfe und Ideen wäre sehr dankbar!

In Bezug auf die Parameter von Multinom und glmnet fand ich diese Antwort hilfreich. Kann ich glm-Algorithmen verwenden, um eine multinomiale logistische Regression durchzuführen?

insbesondere: "Ja, mit einem Poisson GLM (logarithmisches lineares Modell) können Sie multinomiale Modelle anpassen. Daher sind multinomiale logistische oder logarithmische lineare Poisson-Modelle gleichwertig."

Also werde ich die Reparametrisierung der glmnet-Koeffizienten zu Multinom-Koeffizienten zeigen.

n.subj=1000
x1 <- rnorm(n.subj)
x2 <- rnorm(n.subj)
prob <- matrix(c(rep(1,n.subj), exp(3+2*x1+x2), exp(-1+x1-3*x2)), , ncol=3)
prob <- sweep(prob, 1, apply(prob, 1, sum), "/")

y = c()
for (i in 1:n.subj)
  y[i] <- sample(3, 1, replace = T, prob = prob[i,])

multinom(y~x1+x2)

x <- cbind(x1,x2); y2 <- factor(y)
fit <- glmnet(x, y2, family="multinomial", lambda=0, type.multinomial =     "grouped")
cf <- coef(fit)

cf[[2]]@x - cf[[1]]@x   # for the category 2
cf[[3]]@x - cf[[1]]@x   # for the category 3

Hoffe das hilft. Aber ich glaube nicht, dass ich die Äquivalenz von Generalized Linear Model (Poisson) und multinomialem Logistikmodell rein und raus verstehe.

Sagen Sie mir, ob es eine gute und lesbare und "leicht" verständliche Quelle gibt.

Um sicherzustellen, dass die Summe der Auswahlwahrscheinlichkeiten 1 ist, müssen alle Parameter der Referenzalternative Null sein. Daher finde ich das Ergebnis von glmnet () seltsam.

Verwandte F: Warum kann glmnet Parameter für alle Kategorien berechnet werden?