Как заявил Томас Грейф в: подгонка GLM (логистическая регрессия) к SQL

Оригинальный вопрос:

Мы часто оцениваем данные в базе данных напрямую для простых моделей, таких как линейная или логистическая регрессия. Всегда немного сложно правильно перевести все коэффициенты из R в SQL. Я подумал, что могу сделать некоторый перевод R в SQL для результата glm. Для числовых переменных это довольно просто:

library(rpart)

fit <- glm(Kyphosis ~ ., data = kyphosis, family = binomial())

coefs <- fit$coef[2:length(fit$coef)]
expr <- paste0('1/(1 + exp(-(',fit$coef[1], '+', paste0('(', 
           coefs, '*', names(coefs), ')', collapse = '+'),')))')

print(expr)

a <- with(kyphosis, eval(parse(text = expr)))
b <- predict(fit, kyphosis, type = 'response')
names(b) <- NULL
all.equal(a, b)

Сгенерированное выражение:

1/(1 + exp(-(-2.03693352129613+      (0.0109304821420485*Age)+   (0.410601186932733*Number)+(-0.206510049753697*Start)))).

Есть ли способ заставить это работать для факторных переменных? Я хотел бы поставить коэффициенты в случае ... когда ... тогда ... в конце предложения. Предположим, у нас есть следующая модель:

kyphosis$factor_variable <- rep(LETTERS[1:5],20)[1:81]
fit <- glm(Kyphosis ~ ., data = kyphosis, family = binomial())

Я просматриваю структуру фита, но не вижу ничего полезного. Единственный вариант для анализа имен (fit $ coef)?

Вот ссылка на лучший ответ на данный момент... https://stackoverflow.com/a/33659431/6497137

Возможное решение

glm_to_sql <- function(glmmodel) {
  xlev <- data.frame(unlist(glmmodel$xlevels))
  xlev$xlevrowname <- rownames(xlev)
  rownames(xlev) <- NULL
  colnames(xlev)[1] <- "xlevel"
  if (nrow(xlev)==0){xlev <- data.frame(xlevrowname=character(0), xlevel=character(0), stringsAsFactors=F)}

  modcoeffs <- data.frame(unlist(glmmodel$coefficients))
  modcoeffs$coeffname <- rownames(modcoeffs)
  rownames(modcoeffs) <- NULL
  colnames(modcoeffs)[1] <- "coeffvalue"

  coeffmatrix <- sqldf("select a.*,b.*,'' as sqlstr, 
                       substr(coeffname,1,instr(coeffname, xlevel)-1) as varname 
                       from modcoeffs a left join xlev b on coeffname like '%' || xlevel and xlevrowname like substr(coeffname,1,instr(coeffname, xlevel)-1) || '%'")

  for (i in 1:nrow(coeffmatrix)) {
    if(coeffmatrix$coeffname[i] == "(Intercept)") 
    {
      coeffmatrix$sqlstr[i] <- coeffmatrix$coeffvalue[i]
    } else if (is.na(coeffmatrix$xlevel[i]) ) {    
      coeffmatrix$sqlstr[i] <- paste("(",coeffmatrix$coeffvalue[i],"*",coeffmatrix$coeffname[i],")")
    } else {
      coeffmatrix$sqlstr[i] <- paste("(case when ",coeffmatrix$varname[i],"='",coeffmatrix$xlevel[i], "' THEN ",coeffmatrix$coeffvalue[i]," ELSE 0 END)",sep="")
    }

    if (i==1){x.sql0 <- coeffmatrix$sqlstr[i]} else {x.sql0 <- paste(x.sql0,"+",coeffmatrix$sqlstr[i])}
  }

  if (glmmodel$family$link == "logit") {
    x.sql <- paste("1/(1 + exp(-(",x.sql0,")))")  
  } else if (glmmodel$family$link == "identity") {
    x.sql <- x.sql0
  }

  return(x.sql)
}

Проблема

Соединение sqldf не идеально:

where varname is null or length(varname) >0 ## additional filter  

Это не избавит вас от всех углов. Если переменная (например, человек) оканчивается на «n», а другая переменная (например, выживший) равна y/n, тогда она вычтет «n» из человека и соединит его со всеми другими переменными y/n.

У кого-нибудь есть потенциальное обходное решение?

РЕДАКТИРОВАТЬ: пример

library(sqldf)
ID <- seq(1,50,  1)

cabin <- as.numeric(as.character((seq(1,25.5,  .5))))

str(cabin)

Defect <-     c(1,0,1,0,0,1,0,1,0,1,0,1,0,0,0,0,1,0,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,0,1,0,1,1,0,0,0,1,0,1,0,0,0,0,0)

Pre_register <- c("Y", "N", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y",
             "N", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y",
             "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", "Y", 
             "Y", "N", "N", "Y", "N", "N", "Y", "N")

length(Pre_register)
length(cabin)
length(ID)

x <- data.frame(cbind(ID, cabin, Pre_register, Defect))

x$cabin <- as.numeric(as.character(x$cabin))

str(x)

glm_ex <- glm(Defect ~ cabin + Pre_register ,
           data=x,
           family=binomial(link="logit"))

summary(glm_ex)

И вот результат:

> glm_to_sql(glm_ex)

[1] "1/(1 + exp(-( 0,97216 + (случай, когда FLT_REV_Jan_Sep_2015='Y' THEN Round(-1.95327, 3) ELSE 0 END) + (случай, когда ='N' THEN Round(- 1.93112, 3) ИНАЧЕ 0 КОНЕЦ) )))"

Обратите внимание, что в операторе case есть пробел, равный «N». Эта часть неверна и связана с логикой glm_to_sql.

Это соединение, в котором кабина оканчивается на "n", перепутанную с Y/N. Это гораздо меньший пример.

РЕДАКТИРОВАТЬ2:

Проходим через glm_to_sql:

xlev <- data.frame(unlist(glm_ex$xlevels))

xlev$xlevrowname <- rownames(xlev)

rownames(xlev) <- NULL

colnames(xlev)[1] <- "xlevel"

if (nrow(xlev)==0){xlev <- data.frame(xlevrowname=character(0), xlevel=character(0), stringsAsFactors=F)}

xlev

modcoeffs <- data.frame(unlist(glm_ex$coefficients))

modcoeffs$coeffname <- rownames(modcoeffs)

rownames(modcoeffs) <- NULL

colnames(modcoeffs)[1] <- "coeffvalue"

modcoeffs

Вот где существует проблема:

coeffmatrix <- sqldf("select a.*,b.*,'' as sqlstr, 
                   substr(coeffname,1,instr(coeffname, xlevel)-1) as varname 
                 from modcoeffs a left join xlev b on coeffname like '%' || xlevel and xlevrowname like substr(coeffname,1,instr(coeffname, xlevel)-1) || '%'")

Выход:

   coeffvalue     coeffname xlevel   xlevrowname sqlstr      varname
1 -0.51243845   (Intercept)   <NA>          <NA>                <NA>
2 -0.04240967         cabin      N Pre_register1                    
3  1.17625756 Pre_registerY      Y Pre_register2        Pre_register

Проблема существует в строке 2 вывода - кабина ассоциируется с Y/N из уровней Pre_register Y/N, а кабина, оканчивающаяся на букву n, превращается в уровень.