Как предсказание.cv.glmnet() вычисляет значение ссылки для биномиальных моделей?

Я думаю, что это ошибка кодирования с моей стороны, но я не могу понять, что я делаю неправильно. Я пытаюсь создать графики предельных эффектов для x, когда x + x ^ 2 находятся в биномиальной модели glmnet. Когда я делал прогнозы с помощью predict.cv.glmnet() для нового набора данных, казалось, что знаки коэффициентов модели поменялись местами. Итак, я сравнивал прогнозы при жестком кодировании формулы модели и функции прогнозирования. Я получаю разные прогнозы для новых данных, но не для данных обучения, и я не знаю, почему.

#############example############
library(tidyverse)
library(glmnet)
data(BinomialExample)
#head(x)
#head(y)


#training data
squareit<-function(x){x^2}
x%>%
  as.data.frame()%>%
  dplyr::select(V11, V23, V30)%>%#select three variables
  mutate_all(list(SQ = squareit))%>%#square them
  as.matrix()%>%{.->>x1}#back to matrix
# x1[1:5,]
#             V11         V23        V30    V11_SQ       V23_SQ    V30_SQ
# [1,]  0.5706039  0.01473546  0.9080937 0.3255888 0.0002171337 0.8246342
# [2,]  0.4138668 -0.81898245 -2.0492817 0.1712857 0.6707322586 4.1995554
# [3,] -0.4876853  0.03927982 -1.2089006 0.2378369 0.0015429039 1.4614407
# [4,]  0.4644753 -0.98728378  1.7796744 0.2157373 0.9747292699 3.1672408
# [5,] -0.5239595 -0.13288456 -1.1970731 0.2745336 0.0176583076 1.4329841

#example model + coefficients
set.seed(4484)
final.glmnet<-cv.glmnet(x1, y, type.measure="auc",alpha=1,family="binomial")
coef(final.glmnet,s="lambda.min")#print model coefficients
# (Intercept)  0.4097381
# V11         -0.4487949
# V23          0.3322128
# V30         -0.2924600
# V11_SQ      -0.3334376
# V23_SQ      -0.2867963
# V30_SQ       0.3864748

#get model formula - I'm sure there is a better way
mc<-data.frame(as.matrix(coef(final.glmnet,s="lambda.min")))%>%
  rownames_to_column(.,var="rowname")%>%mutate(rowname=recode(rowname,`(Intercept)`="1"))
paste("(",mc$rowname,"*",mc$X1,")",sep="",collapse = "+")
#model formula
# (1*0.409738094215867)+(V11*-0.44879489356345)+(V23*0.332212780157358)+
# (V30*-0.292459974168585)+(V11_SQ*-0.333437624504966)+
# (V23_SQ*-0.286796292157334)+(V30_SQ*0.386474813542322)

#####predict on training data -- same predictions!
#1 predict using predict.cv.glmnet()
pred_cv.glmnet<-predict(final.glmnet, s=final.glmnet$lambda.min, newx=x1[1:5,], type='link')
round(pred_cv.glmnet,3)
#2 predict using model formula
pred.model.formula<-data.frame(x1[1:5,])%>%
  mutate(link=(1*0.409738094215867)+(V11*-0.44879489356345)+
           (V23*0.332212780157358)+(V30*-0.292459974168585)+
           (V11_SQ*-0.333437624504966)+(V23_SQ*-0.286796292157334)+
           (V30_SQ*0.386474813542322))%>%
  dplyr::select(link)
round(pred.model.formula,3)
# 1         0.103
# 2         1.925
# 3         1.480
# 4         0.225
# 5         1.408

#new data - set up for maringal effects plot
colnames(x1)
test<-data.frame(V23=seq(min(x1[,2]),max(x1[,2]),0.1))%>%#let V23 vary from min-max by .1
  mutate(V23_SQ=V23^2, V11=0, V11_SQ=0, V30=0, V30_SQ=0)%>%#square V23 and set others to 0
as.matrix()
test[1:5,]
# > test[1:5,]
#           V23   V23_SQ V11 V11_SQ V30 V30_SQ
# [1,] -2.071573 4.291414   0      0   0      0
# [2,] -1.971573 3.887100   0      0   0      0
# [3,] -1.871573 3.502785   0      0   0      0
# [4,] -1.771573 3.138470   0      0   0      0
# [5,] -1.671573 2.794156   0      0   0      0

#####predict on new data -- different predictions!
#1 predict using predict.cv.glmnet()
pred_cv.glmnet<-predict(final.glmnet, s=final.glmnet$lambda.min, newx=test[1:5,], type='link')
round(pred_cv.glmnet,3)
# [1,] 2.765
# [2,] 2.586
# [3,] 2.413
# [4,] 2.247
# [5,] 2.088

#2 predict using model formula
pred.model.formula<-data.frame(test[1:5,])%>%
  mutate(link.longhand=(1*0.409738094215867)+(V11*-0.44879489356345)+
           (V23*0.332212780157358)+(V30*-0.292459974168585)+
           (V11_SQ*-0.333437624504966)+(V23_SQ*-0.286796292157334)+
           (V30_SQ*0.386474813542322))%>%
  dplyr::select(link.longhand)
round(pred.model.formula,3)
# 1        -1.509
# 2        -1.360
# 3        -1.217
# 4        -1.079
# 5        -0.947

person Kevin    schedule 25.01.2021    source источник


Ответы (1)


Ваше понимание того, как работает расчет, кажется правильным. Однако функция predict требует, чтобы матрица newx имела те же размеры, что и оригинал (т. е. столбцы должны быть в том же порядке). Когда столбцы перепутаны, результаты функции predict сойдут с ума, как вы видите в вашем примере.

Если вы добавите быструю корректировку, чтобы расположить матрицу test в том же порядке, что и матрица x1 (удобный ярлык — просто использовать dimnames из исходной матрицы в dplyr::select, как в строке 3 ниже)...

test <- data.frame(V23 = seq(min(x1[, 2]), max(x1[, 2]), 0.1)) %>% #let V23 vary from min-max by .1
    mutate(V23_SQ = V23^2, V11 = 0, V11_SQ = 0, V30 = 0, V30_SQ = 0) %>% #square V23 and set others to 0
    select(dimnames(x1)[[2]]) %>% #use the dimnames from x1 to select test in proper order!
    as.matrix()

... вы обнаружите, что получаете правильные результаты, которые совпадают с ожидаемыми.

#1 predict using predict.cv.glmnet()
pred_cv.glmnet <- predict(final.glmnet, s = 'lambda.min', newx = test[1:5,], type = 'link')
round(pred_cv.glmnet,3)
# [1,] -1.509
# [2,] -1.360
# [3,] -1.217
# [4,] -1.079
# [5,] -0.947

#2 predict using model formula
pred.model.formula <- data.frame(test[1:5,]) %>%
    mutate(link.longhand = ((1 * 0.409738094215867) +
               (V11 * -0.44879489356345) +
               (V23 * 0.332212780157358) +
               (V30 * -0.292459974168585) +
               (V11_SQ * -0.333437624504966) +
               (V23_SQ * -0.286796292157334) +
               (V30_SQ * 0.386474813542322))) %>%
    dplyr::select(link.longhand)
round(pred.model.formula, 3)
# link.longhand
# 1        -1.509
# 2        -1.360
# 3        -1.217
# 4        -1.079
# 5        -0.947
person Colin H    schedule 26.01.2021
comment
Большое спасибо за Вашу помощь! Я очень застрял. - person Kevin; 26.01.2021