Siguiente problema: quiero predecir una variable de respuesta categórica con uno (o más) variables categóricas mediante glmnet().
Sin embargo, no puedo hacer sentido de la salida de glmnet me da.
Ok, primero que vamos a generar dos variables categóricas:
Generar Datos
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")
Obtenemos:
> table(d)
X2
X1 0 1 2
0 22 11 1
1 9 14 10
2 3 8 22
Predicción: multinom()
A continuación, vamos a predecir X1 por X2 utilizando el multinom() de la nnet paquete:
library(nnet)
mod1 <- multinom(X1~X2, data=d)
mod1
lo que nos da:
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
Manual de verificación
Ahora vamos a ver, si podemos reproducir manualmente:
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
Producimos los mismos números, ¡bueno!
Predicción: glmnet()
Ahora vamos a hacer lo mismo con 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
Tenga en cuenta que he conjunto s=0, por lo tanto no es de regularización y los parámetros que debe contener exactamente la misma información que los parámetros de la multinom() función.
Aún así, tenemos muy diferentes parámetros. Esto es debido a las diferentes parametrización utilizan en glmnet, véase por ejemplo:
http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (título: modelos Multinomiales) o el papel correspondiente: http://www.jstatsoft.org/v33/i01/paper (título: 4. Regularización de la regresión multinomial)
Pero no importa cómo, exactamente, una parametriza, se debe obtener el mismo $P(Y=k | X)$, la probabilidad de la categoría k condicional en X.
Probabilidades condicionales: multinom()
Así que primero tengo que calcular estas probabilidades 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
Como hemos saturado el modelo de aquí, este debe ser el mismo que lo podemos calcular a partir de los datos:
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
que es de hecho el caso.
Probabilidades condicionales: glmnet()
Ahora, el mismo de 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
La celda de probabilidades condicionales son algo diferentes pero relacionadas. También que no suma de a uno.
Así que tenemos dos problemas aquí:
a) las probabilidades condicionales no suma 1 y
b) los parámetros que no se describir lo que vemos en los datos: por ejemplo, en la fila 2, existen diferencias a través de las columnas, pero glmnet las estimaciones de ambos coeficientes (no la intercepción) como cero.
He utilizado una regresión lineal problema y en comparación con glm y glmnet con s=0 para asegurarse de que s=0, significa cero de regularización (las soluciones eran casi idénticos).
Cualquier ayuda y de ideas, sería muy apreciada!
