Problème suivant: Je veux prédire une variable de réponse catégorielle avec une (ou plusieurs) variables catégorielles en utilisant glmnet ().
Cependant, je ne peux pas comprendre la sortie que me donne glmnet.
Ok, commençons par générer deux variables catégorielles liées:
Générer des données
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")
On a:
> table(d)
X2
X1 0 1 2
0 22 11 1
1 9 14 10
2 3 8 22
Prédiction: multinom ()
Prédisons alors X1 par X2 en utilisant le multinom () du package nnet:
library(nnet)
mod1 <- multinom(X1~X2, data=d)
mod1
ce qui nous donne:
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
Vérification manuelle
Voyons maintenant si nous pouvons reproduire cela manuellement:
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
Nous produisons les mêmes numéros, bon!
Prédiction: glmnet ()
Maintenant, faisons de même avec 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
Notez que je mets s = 0, donc il n'y a pas de régularisation et les paramètres doivent contenir exactement les mêmes informations que les paramètres de la fonction multinom ().
Pourtant, nous obtenons des paramètres très différents. Cela est dû aux différents paramétrages qu'ils utilisent dans glmnet, voir par exemple:
http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (en-tête: Modèles multinomiaux) ou le document correspondant: http://www.jstatsoft.org/v33/i01/paper (en-tête: 4. Régularisé) régression multinomiale)
Probabilités conditionnelles: multinom ()
J'ai donc d'abord calculé ces probabilités à partir de 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
Comme nous avons un modèle saturé ici, cela devrait être le même que ce que nous pouvons calculer à partir des données:
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
ce qui est effectivement le cas.
Probabilités conditionnelles: glmnet ()
Maintenant la même chose avec 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
Les probabilités conditionnelles des cellules sont quelque peu liées mais différentes. De plus, ils ne résument pas à un.
Nous avons donc deux problèmes ici:
a) les probabilités conditionnelles ne totalisent pas 1 et
b) les paramètres ne décrivent pas ce que nous voyons dans les données: par exemple, dans la ligne 2, il y a des différences entre les colonnes, mais glmnet estime les deux coefficients (et non l'ordonnée à l'origine) comme nuls.
J'ai utilisé un problème de régression linéaire et j'ai comparé glm et glmnet à s = 0 pour m'assurer que s = 0 signifie une régularisation nulle (les solutions étaient presque identiques).
Toute aide et idée serait très appréciée!
Pour s'assurer que la somme des probabilités de choix est 1, tous les paramètres de l'alternative de référence doivent être nuls. Donc, je pense que le résultat de glmnet () est étrange.
Q: Pourquoi glmnet peut-il être calculé des paramètres pour toutes les catégories?
la source