Masalah berikut: Saya ingin memprediksi variabel respons kategoris dengan satu (atau lebih) variabel kategorikal menggunakan glmnet ().
Namun, saya tidak dapat memahami output yang diberikan glmnet kepada saya.
Ok, pertama mari kita buat dua variabel kategori terkait:
Hasilkan Data
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")
Kita mendapatkan:
> table(d)
X2
X1 0 1 2
0 22 11 1
1 9 14 10
2 3 8 22
Prediksi: multinom ()
Kemudian mari kita prediksi X1 oleh X2 menggunakan multinom () dari paket nnet:
library(nnet)
mod1 <- multinom(X1~X2, data=d)
mod1
yang memberi kita:
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
Pemeriksaan manual
Sekarang mari kita periksa, apakah kita dapat mereproduksi itu secara manual:
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
Kami menghasilkan angka yang sama, bagus!
Prediksi: glmnet ()
Sekarang mari kita lakukan hal yang sama dengan 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
Perhatikan bahwa saya menetapkan s = 0, sehingga tidak ada regularisasi dan parameter harus berisi informasi yang sama persis dengan parameter fungsi multinom ().
Namun, kami mendapatkan parameter yang sangat berbeda. Ini disebabkan oleh berbagai parameterisasi yang mereka gunakan di glmnet, lihat misalnya:
http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (judul: model multinomial) atau kertas yang sesuai: http://www.jstatsoft.org/v33/i01/paper (heading: 4. Diatur secara teratur regresi multinomial)
Peluang bersyarat: multinom ()
Jadi saya pertama-tama saya menghitung probabilitas ini dari 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
Karena kita memiliki model jenuh di sini, ini harus sama dengan apa yang dapat kita hitung dari data:
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
yang memang demikian.
Peluang bersyarat: glmnet ()
Sekarang sama dari 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
Probabilitas kondisional sel agak terkait tetapi berbeda. Juga mereka tidak meringkas satu.
Jadi kita punya dua masalah di sini:
a) probabilitas bersyarat tidak meringkaskan hingga 1 dan
b) parameter tidak menggambarkan apa yang kita lihat dalam data: misalnya di baris 2 ada perbedaan di seluruh kolom, tetapi glmnet memperkirakan kedua koefisien (bukan intersep) sebagai nol.
Saya menggunakan masalah regresi linier dan membandingkan glm dan glmnet dengan s = 0 untuk memastikan bahwa s = 0 berarti nol regularisasi (solusinya hampir identik).
Bantuan & ide apa pun akan sangat dihargai!
Untuk memastikan bahwa jumlah probabilitas pilihan adalah 1, semua parameter dari alternatif referensi harus nol. Jadi, saya pikir hasil glmnet () aneh.
T terkait: Mengapa glmnet dapat dihitung parameter untuk semua kategori?
sumber